changeset 13098:8b7e6f359cee stable

logm.m: Return real matrix when all eigenvalues are real (Bug #32121). * logm.m: Remove complex numbers of order eps() which may have entered return value through numeric roundoff.
author Rik <octave@nomad.inbox5.com>
date Sun, 04 Sep 2011 16:46:16 -0700
parents 70d32160c90b
children c7512d0d52e8
files scripts/linear-algebra/logm.m
diffstat 1 files changed, 9 insertions(+), 1 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/linear-algebra/logm.m	Fri Sep 02 21:43:10 2011 +0100
+++ b/scripts/linear-algebra/logm.m	Sun Sep 04 16:46:16 2011 -0700
@@ -63,11 +63,14 @@
     [u, s] = rsf2csf (u, s);
   endif
 
-  if (any (diag (s) < 0))
+  eigv = diag (s);
+  if (any (eigv < 0))
     warning ("Octave:logm:non-principal",
              "logm: principal matrix logarithm is not defined for matrices with negative eigenvalues; computing non-principal logarithm");
   endif
 
+  real_eig = all (eigv >= 0);
+
   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];
@@ -100,6 +103,11 @@
 
   s = 2^k * u * s * u';
 
+  ## Remove small complex values (O(eps)) which may have entered calculation
+  if (real_eig)
+    s = real (s);
+  endif
+
   if (nargout == 2)
     iters = k;
   endif