Mercurial > hg > octave-lyh

changeset 13097:52e4aa30d5b2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 5553412c6614
children bf61bc523805
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
+++ b/scripts/linear-algebra/logm.m
@@ -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