function [X, arg2, condest] = sqrtm(A)
%SQRTM     Matrix square root.
%   X = SQRTM(A) is the principal square root of the matrix A, i.e. X*X = A.
%          
%   X is the unique square root for which every eigenvalue has nonnegative
%   real part.  If A has any eigenvalues with negative real parts then a
%   complex result is produced.  If A is singular then A may not have a
%   square root.  A warning is printed if exactly singularity is detected.
%          
%   With two output arguments, [X, RESNORM] = SQRTM(A) does not print any
%   warning, and returns the residual, norm(A-X^2,'fro')/norm(A,'fro').
%
%   With three output arguments, [X, ALPHA, CONDEST] = SQRTM(A) returns a
%   stability factor ALPHA and an estimate CONDEST of the matrix square root
%   condition number of X.  The residual NORM(A-X^2,'fro')/NORM(A,'fro') is
%   bounded by (n+1)*ALPHA*EPS and the Frobenius norm relative error in X is
%   bounded approximately by N*ALPHA*CONDEST*EPS, where N = MAX(SIZE(A)).
%
%   See also EXPM, SQRTM, FUNM.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.6 $  $Date: 1998/10/15 15:15:57 $
%
%   References:
%   N. J. Higham, Computing real square roots of a real
%       matrix, Linear Algebra and Appl., 88/89 (1987), pp. 405-430.
%   A. Bjorck and S. Hammarling, A Schur method for the square root of a
%       matrix, Linear Algebra and Appl., 52/53 (1983), pp. 127-140.

n = max(size(A));
[Q, T] = schur(A);        % T is real/complex according to A.
[Q, T] = rsf2csf(Q, T);   % T is now complex Schur form.

% Compute upper triangular square root R of T, a column at a time.

warns = warning;
warning('off');

R = zeros(n);
for j=1:n
    R(j,j) = sqrt(T(j,j));
    for i=j-1:-1:1
        s = 0;
        for k=i+1:j-1
            s = s + R(i,k)*R(k,j);
        end
        R(i,j) = (T(i,j) - s)/(R(i,i) + R(j,j));
    end
end

warning(warns);

if nargout > 2
   arg2 = norm(R,'fro')^2 / norm(T,'fro');
end

X = Q*R*Q';

nzeig = any(diag(T)==0);

if nzeig & (nargout ~= 2)
   warning('Matrix is singular and may not have a square root.')
end

if nargout == 2
   arg2 = norm(X*X-A,'fro')/norm(A,'fro');
end

if nargout > 2

   if nzeig
      condest = inf;
   else

      % Power method to get condition number estimate.

      tol = 1e-2;
      x = ones(n^2,1);    % Starting vector.
      cnt = 1;
      e = 1;
      e0 = 0;
      while abs(e-e0) > tol*e & cnt <= 6
         x = x/norm(x);
         x0 = x;
         e0 = e;
         Sx = tksolve(R, x);
         x = tksolve(R, Sx, 'T');
         e = sqrt(real(x0'*x));  % sqrt of Rayleigh quotient.
         % fprintf('cnt = %2.0f, e = %9.4e\n', cnt, e)
         cnt = cnt+1;
      end

      condest = e*norm(A,'fro')/norm(X,'fro');

   end

end

%====================================

function x = tksolve(R, b, tran)
%TKSOLVE     Solves triangular Kronecker system.
%            x = TKSOLVE(R, b, TRAN) solves
%                  A*x = b  if TRAN = '',
%                 A'*x = b  if TRAN = 'T',
%            where A = KRON(EYE,R) + KRON(TRANSPOSE(R),EYE).
%            Default: TRAN = ''.

if nargin < 3, tran = ''; end

n = max(size(R));
x = zeros(n^2,1);

I = eye(n);

if isempty(tran)

   % Forward substitution.
   for i = 1:n
       temp = b(n*(i-1)+1:n*i);
       for j = 1:i-1
           temp = temp - R(j,i)*x(n*(j-1)+1:n*j);
       end
       x(n*(i-1)+1:n*i) = (R + R(i,i)*I) \ temp;
   end

elseif strcmp(tran,'T')

   % Back substitution.
   for i = n:-1:1
       temp = b(n*(i-1)+1:n*i);
       for j = i+1:n
           temp = temp - conj(R(i,j))*x(n*(j-1)+1:n*j);
       end
       x(n*(i-1)+1:n*i) = (R' + conj(R(i,i))*I) \ temp;
   end

end

return
