function [L,esterr] = logm(A)
%LOGM   Matrix logarithm.
%   L = LOGM(A) is the matrix logarithm of A, the inverse of EXPM(A).
%   Complex results are produced if A has negative eigenvalues.
%   A warning message is printed if the computed expm(L) is not close to A.
%
%   [L,esterr] = logm(A) does not print any warning message, but 
%   returns an estimate of the relative residual, norm(expm(L)-A)/norm(A).
%
%   If A is real symmetric or complex Hermitian, then so is LOGM(A).
%
%   Some matrices, like A = [0 1; 0 0], do not have any logarithms,
%   real or complex, and LOGM cannot be expected to produce one.
%
%   See also EXPM, SQRTM, FUNM.

%   CBM, 7-12-92, 2-1-96.
%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 5.8 $  $Date: 2000/06/01 02:04:17 $

% First try Parlett's method directly.
[L,esterr] = funm(A,'log');
tol = 1000*eps;
% If funm's error estimate is small, accept the result.
if esterr >= tol
   % Use norm of residual instead of funm's crude error estimate.
   esterr = norm(expm(L)-A,1)/norm(A,1);
   if esterr >= tol | ~isfinite(esterr)
      [n,n] = size(A);
      I = eye(n,n);
      % Try again with a not-quite-random rotation.
      [J,K] = meshgrid(1:n);
      R = orth(I + J + K);
      [L,ignore] = funm(R*A*R','log');
      L = R'*L*R;
      if norm(imag(L),1) <= 1000*tol*norm(L,1), L = real(L); end
      % One step of improvement.
      E = expm(L);
      L = L - (E - A)/E;
      esterr = norm(expm(L)-A,1)/norm(A,1);
   end
end
if (esterr >= tol | ~isfinite(esterr)) & (nargout < 2)
   warning(['LOGM appears inaccurate.  esterr = ',num2str(esterr)])
end
