home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / MATFUN.DI$ / SQRTM.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.4 KB  |  41 lines
  1. function [S,esterr] = sqrtm(A)
  2. %SQRTM    Matrix square root.
  3. %    S = SQRTM(A) is the matrix square root of A.
  4. %    Complex results are produced if A has negative eigenvalues.
  5. %    A warning message is printed if the computed S*S is not close to A.
  6. %
  7. %    [S,esterr] = sqrtm(A) does not print any warning message, but
  8. %    returns an estimate of the relative residual, norm(S*S-A)/norm(A).
  9. %
  10. %    If A is real, symmetric and positive definite, or complex, Hermitian
  11. %    and positive definite, then so is the computed matrix square root.
  12. %
  13. %    Some matrices, like A = [0 1; 0 0], do not have any square roots,
  14. %    real or complex, and SQRTM cannot be expected to produce one.
  15. %
  16. %    See also EXPM, LOGM, FUNM.
  17.  
  18. %    CBM, 7-12-92.
  19. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  20.  
  21. % First try Parlett's method directly.
  22. [S,esterr] = funm(A,'sqrt');
  23. tol = 1000*eps;
  24. % If funm's error estimate is small, accept the result.
  25. if esterr >= tol
  26.    % Use norm of residual instead of funm's crude error estimate.
  27.    esterr = norm(S*S-A,1)/norm(A,1);
  28.    if esterr >= tol
  29.       [n,n] = size(A);
  30.       % Try again with a not-quite-random rotation.
  31.       R = orth(eye(n,n) + magic(n)/(n*(n*n+1)/2));
  32.       [S,ignore] = funm(R*A*R','sqrt');
  33.       S = R'*S*R;
  34.       esterr = norm(S*S-A,1)/norm(A,1);
  35.    end
  36. end
  37. if (esterr >= tol) & (nargout < 2)
  38.    disp(' ')
  39.    disp(['Warning: SQRTM appears inaccurate.  esterr = ',num2str(esterr)])
  40. end
  41.