changeset 8402:2176f2b4599e

Fix sparse cholesky inversion
author David Bateman <dbateman@free.fr>
date Fri, 12 Dec 2008 23:18:20 +0100
parents 712cfdc2e417
children 87cca636a6c6
files liboctave/ChangeLog liboctave/sparse-base-chol.cc src/ChangeLog src/DLD-FUNCTIONS/chol.cc
diffstat 4 files changed, 30 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog	Fri Dec 12 13:35:10 2008 +0100
+++ b/liboctave/ChangeLog	Fri Dec 12 23:18:20 2008 +0100
@@ -24,6 +24,11 @@
 	* CmplxAEPBAL.h, CmplxAEPBAL.cc: Rebase ComplexAEPBAL on base_aepbal.
 	* fCmplxAEPBAL.h, fCmplxAEPBAL.cc: Rebase FloatComplexAEPBAL on base_aepbal.
 
+2008-12-12  David Bateman  <dbateman@free.fr>
+
+	* sparse-base-chol.cc (inverse): Fix inversion based on cholesky 
+	factorization.
+	
 2008-12-08  Jaroslav Hajek  <highegg@gmail.com>
 
 	* idx-vector.cc (idx_vector::idx_vector_rep::idx_vector_rep (const 
--- a/liboctave/sparse-base-chol.cc	Fri Dec 12 13:35:10 2008 +0100
+++ b/liboctave/sparse-base-chol.cc	Fri Dec 12 23:18:20 2008 +0100
@@ -277,15 +277,15 @@
   double rcond2;
   octave_idx_type info;
   MatrixType mattype (MatrixType::Upper);
-  chol_type linv = L().transpose().inverse(mattype, info, rcond2, 1, 0);
+  chol_type linv = L().hermitian().inverse(mattype, info, rcond2, 1, 0);
 
   if (perms.length() == n)
     {
       p_type Qc = Q();
-      retval = Qc * linv.transpose() * linv * Qc.transpose();
+      retval = Qc * linv * linv.hermitian() * Qc.transpose();
     }
   else
-    retval = linv.transpose() * linv;
+    retval = linv * linv.hermitian ();
 #endif
   return retval;
 }
--- a/src/ChangeLog	Fri Dec 12 13:35:10 2008 +0100
+++ b/src/ChangeLog	Fri Dec 12 23:18:20 2008 +0100
@@ -46,6 +46,10 @@
 
 	* DLD-FUNCTIONS/balance.cc (Fbalance): Exploit the new AEPBAL functionality.
 
+2008-12-12  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/chol.cc (Fcholinv): Add test.
+
 2008-12-08  Jaroslav Hajek  <highegg@gmail.com>
 	
 	* xpow.cc ( xpow (const DiagMatrix& a, double b), 
--- a/src/DLD-FUNCTIONS/chol.cc	Fri Dec 12 13:35:10 2008 +0100
+++ b/src/DLD-FUNCTIONS/chol.cc	Fri Dec 12 23:18:20 2008 +0100
@@ -469,6 +469,24 @@
   return retval;
 }
 
+/*
+
+%!test
+%! A = [2,0.2;0.2,1];
+%! issymmetric(A)
+%! min(eig(A))
+%! Ainv = inv(A);
+%! Ainv1 = cholinv(A);
+%! Ainv2 = inv(sparse(A));
+%! Ainv3 = cholinv(sparse(A));
+%! Ainv4 = spcholinv(sparse(A));
+%! assert (norm(Ainv-Ainv1),1e-10)
+%! assert (norm(Ainv-Ainv2),1e-10)
+%! assert (norm(Ainv-Ainv3),1e-10)
+%! assert (norm(Ainv-Ainv4),1e-10)
+
+*/
+
 DEFUN_DLD (chol2inv, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {} chol2inv (@var{u})\n\