## Copyright (C) 1993, 1995, 1996, 1997, 1999, 2000, 2005, 2006, 2007
##               John W. Eaton
##
## This file is part of Octave.
##
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{S}, @var{iters}] =} logm (@var{A}, @var{opt_iters})
## Compute the matrix logarithm of the square matrix @var{a}.  Utilizes a Pade 
## approximant and the identity
##
## @code{ logm(@var{A}) = 2^k * logm(@var{A}^(1 / 2^k)) }.
##
## Optional argument @var{opt_iters} is the number of square roots computed
## and defaults to 100.  Optional output @var{iters} is the number of square
## roots actually computed.
##
## The Pade approximation was adapted from the mftoobox by D. Higham.
##
## Reference: D. Higham, Functions of Matrices: Theory and Computation 
## 	 (SIAM, 2008.)
##
## @end deftypefn
 
## Author: Richard T. Guy with significant code modified from mftoolbox (D. Higham)
## Created: April 29, 2010 

function [S,iters] = logm(A,opt_iters)
 
if (nargin == 0)
  print_usage();
  return;
elseif (nargin < 2)
  opt_iters = 100;
elseif (nargin > 2)
  print_usage();
  return;
endif

if (! issquare(A))
  error("logm: argument must be a square matrix.\n");
  return;
endif

tic;
[u,s] = schur(A);
toc

tic;
[U,S] = rsf2csf (u, s);
toc

if (any(diag(S) < 0))
  warning("Matrix has nonegative eigenvalues.  Principal matrix logarithm is not defined.  Computing non-principal logarithm.");
endif

k = 0;
## Algorithm 11.9 in "Function of matrices", by N. Higham
tic;
theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1];
p = 0;
while (k < opt_iters)
  tau = norm (S - eye (size (S)),1);
  if (tau <= theta (7))
    p = p + 1;
    j (1) = min (find (tau <=  theta));
    j (2) = min (find (tau / 2 <=  theta));
    if (j (1) - j (2) <= 1 || p == 2)
      m = j (1);
      break
    endif
  endif
  k = k + 1;
  S = sqrtm (S);  # TODO: There are faster ways to do this.
endwhile
toc

if (k >= opt_iters)
  warning("Maximum number of square roots exceeded.  Results may still be accurate.");
endif

tic;
S = logm_pade_pf(S-eye(size(S)),m);
toc

tic;
S = 2^k * U * S * U';
toc

if (nargout == 2)
  iters = k;
endif

endfunction

################## ANCILLARY FUNCTIONS ################################
######	Taken from the mfttoolbox (GPL 3) by D. Higham.
######  Reference: 
######      D. Higham, Functions of Matrices: Theory and Computation 
######		(SIAM, 2008.).
#######################################################################

function S = logm_pade_pf(A,m)

##LOGM_PADE_PF   Evaluate Pade approximant to matrix log by partial fractions.
##   Y = LOGM_PADE_PF(A,M) evaluates the [M/M] Pade approximation to
##   LOG(EYE(SIZE(A))+A) using a partial fraction expansion.

[nodes,wts] = gauss_legendre(m);
## Convert from [-1,1] to [0,1].
nodes = (nodes + 1)/2;
wts = wts/2;

n = length(A);
S = zeros(n);
for j=1:m
  S = S + wts(j)*(A/(eye(n) + nodes(j)*A));
endfor

endfunction

######################################################################

function [x,w] = gauss_legendre(n)

## GAUSS_LEGENDRE  Nodes and weights for Gauss-Legendre quadrature.
##   [X,W] = GAUSS_LEGENDRE(N) computes the nodes X and weights W
##   for N-point Gauss-Legendre quadrature.

## Reference:
## G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature
## rules, Math. Comp., 23(106):221-230, 1969.
i = 1:n-1;
v = i./sqrt((2*i).^2-1);
[V,D] = eig( diag(v,-1)+diag(v,1) );
x = diag(D);
w = 2*(V(1,:)'.^2);
endfunction


%!assert(norm(logm([1 -1;0 1]) - [0 -1; 0 0]) < 1e-5);
%!assert(norm(expm(logm([-1 2 ; 4 -1])) - [-1 2 ; 4 -1]) < 1e-5);
%!
%!error logm ();
%!error logm (1, 2, 3);
%!error <logm: argument must be a square matrix.> logm([1 0;0 1; 2 2]);s
