--- a/scripts/ChangeLog	Tue Jul 27 10:46:11 2010 -0700
+++ b/scripts/ChangeLog	Wed Jul 28 09:22:41 2010 +0200
@@ -1,3 +1,7 @@
+2010-07-28  Jaroslav Hajek  <highegg@gmail.com>
+
+	* linear-algebra/logm.m: Rewrite. Thanks to M. Caliari and R. T. Guy.
+
 2010-07-26  Rik <octave@nomad.inbox5.com>
 	* deprecated/complement.m, deprecated/intwarning.m, general/arrayfun.m,
 	general/circshift.m, general/colon.m, general/common_size.m,

--- a/scripts/linear-algebra/logm.m	Tue Jul 27 10:46:11 2010 -0700
+++ b/scripts/linear-algebra/logm.m	Wed Jul 28 09:22:41 2010 +0200
@@ -1,35 +1,133 @@
-## Copyright (C) 2003, 2005, 2006, 2007 John W. Eaton
+## Copyright (C) 2010 Richard T. Guy <guyrt7@wfu.edu>
+## Copyright (C) 2010 Marco Caliari <marco.caliari@univr.it>
+## Copyright (C) 2008 N.J. Higham
 ##
 ## This file is part of Octave.
 ##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {} logm (@var{a})
-## Compute the matrix logarithm of the square matrix @var{a}.  Note that
-## this is currently implemented in terms of an eigenvalue expansion and
-## needs to be improved to be more robust.
+## @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.
+##
 ## @end deftypefn
 
-function B = logm (A)
+## Reference: N. J. Higham, Functions of Matrices: Theory and Computation 
+## 	 (SIAM, 2008.)
+##
 
-  if (nargin != 1)
+function [s, iters] = logm (a, opt_iters)
+ 
+  if (nargin == 0)
+    print_usage ();
+  elseif (nargin < 2)
+    opt_iters = 100;
+  elseif (nargin > 2)
     print_usage ();
   endif
 
-  [V, D] = eig (A);
-  B = V * diag (log (diag (D))) * inv (V);
+  if (! issquare (a))
+    error ("logm: argument must be a square matrix.");
+  endif
+
+  [u, s] = schur (a);
+
+  if (isreal (a))
+    [u, s] = rsf2csf (u, s);
+  endif
+
+  if (any (diag (s) < 0))
+    warning ("Octave:logm:non-principal",
+    ["logm: 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
+  theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1];
+  p = 0;
+  m = 7;
+  while (k < opt_iters)
+    tau = norm (s - eye (size (s)),1);
+    if (tau <= theta (7))
+      p = p + 1;
+      j(1) = find (tau <= theta, 1);
+      j(2) = find (tau / 2 <= theta, 1);
+      if (j(1) - j(2) <= 1 || p == 2)
+        m = j(1);
+        break
+      endif
+    endif
+    k = k + 1;
+    s = sqrtm (s);
+  endwhile
+
+  if (k >= opt_iters)
+    warning ("logm: Maximum number of square roots exceeded.  Results may still be accurate.");
+  endif
+
+  s = logm_pade_pf (s - eye (size (s)), m);
+
+  s = 2^k * u * s * u';
+
+  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.).
+#######################################################################
+
+##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.
+
+function s = logm_pade_pf(a,m)
+  [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
+
+######################################################################
+## 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.
+
+function [x,w] = gauss_legendre(n)
+  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]);