function [X, alpha, condest] = sqrtm2(A)
%SQRTM2    Matrix square root.
%          X = SQRTM2(A) is a square root of the matrix A (A = X*X).
%          X is the unique square root for which every eigenvalue has
%          nonnegative real part (the principal square root).
%          If A is real with a negative real eigenvalue then a complex
%          result will be produced.
%          If A is singular then A may not have a square root.
%          [X, ALPHA, CONDEST] = SQRTM2(A) returns a stability factor ALPHA and
%          an estimate CONDEST for 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)).

%          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.

nzeig = length(find(diag(T)==0));
if nzeig, warning('Matrix is singular and may not have a square root.'), end

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
        if R(i,i) + R(j,j) ~= 0
           R(i,j) = (T(i,j) - s)/(R(i,i) + R(j,j));
        end
    end
end

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

X = Q*R*Q';

if nargout > 2

   if nzeig
      condest = inf;
   else

      % Power method to get condition number estimate.
      warns = warning;
      warning('off');

      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 = solve(R, x);
         x = solve(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');
      warning(warns);

   end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = solve(R, b, tran)
%SOLVE       Solves Kronecker system.
%            x = SOLVE(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