Mercurial > octave

changeset 22242:5417dad26a25

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
chol2inv: "inv (Y'*Y)" instead of "(inv(Y))' * inv(Y)" (bug #36437) * liboctave/numeric/sparse-chol.cc (chol2inv): compute "inv (Y'*Y)" instead of "(inv(Y))' * inv(Y)".
author Barbara Locsi <locsi.barbara@gmail.com>
date Tue, 09 Aug 2016 23:25:09 +0100
parents 71d86e88589f
children 654de580bdb3
files liboctave/numeric/sparse-chol.cc
diffstat 1 files changed, 7 insertions(+), 5 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/numeric/sparse-chol.cc	Sat Aug 06 19:38:32 2016 +0200
+++ b/liboctave/numeric/sparse-chol.cc	Tue Aug 09 23:25:09 2016 +0100
@@ -546,21 +546,23 @@
   int typ = mattype.type (false);
   double rcond;
   octave_idx_type info;
-  chol_type rinv;
+  chol_type rtra, multip;
 
   if (typ == MatrixType::Upper)
     {
-      rinv = r.inverse (mattype, info, rcond, true, false);
-      retval = rinv.transpose () * rinv;
+      rtra = r.transpose ();
+      multip = (rtra*r);
     }
   else if (typ == MatrixType::Lower)
     {
-      rinv = r.transpose ().inverse (mattype, info, rcond, true, false);
-      retval = rinv.transpose () * rinv;
+      rtra = r.transpose ();
+      multip = (r*rtra);
     }
   else
     (*current_liboctave_error_handler) ("U must be a triangular matrix");
 
+  MatrixType mattypenew (multip);
+  retval = multip.inverse (mattypenew, info, rcond, true, false);
   return retval;
 }