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