function [D, fX] = dsqrtm(X, H)
% Fréchet derivative of the matrix square root.
%
% function [D, fX] = dsqrtm(X, H)
%
% Computes the directional derivative (the Fréchet derivative) of sqrtm at
% X along H (square matrices).
%
% Thus, D = lim_(t -> 0) (sqrtm(X + tH) - sqrtm(X)) / t.
%
% The second output is fX = sqrtm(X), though it may be less accurate.
%
% Note: the adjoint of dsqrtm(X, .) is dsqrtm(X', .), which is a fact often
% useful to derive gradients of matrix functions involving sqrtm(X).
% (This is wrt the inner product inner = @(A, B) real(trace(A'*B))).
% 
% See also: dfunm dlogm dexpm

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 3, 2015.
% Contributors:
% Change log:
%
%   June 14, 2019 (NB): now also outputs sqrtm(X) as a by-product.
    
    % Note: following Higham, 'Functions of Matrices', 2008, page 58: this
    % could also be computed as fX = sqrtm(X), then solving the Sylvester
    % equation fX*D + D*fX = H, e.g. via , D = sylvester(fX, fX, H).
    % If X has special structure (e.g., if it is symmetric or Hermitian),
    % then this may be faster and more accurate. This should be tested
    % before considering a replacement.

    [D, fX] = dfunm(@sqrtm, X, H);
    
end
