# HG changeset patch
# User jwe
# Date 1175699871 0
# Node ID e978a9233cf64c0252d72be17aaa3eb72c46ca33
# Parent  0f233b5b96a1d661c1bf6277d5f7668ffbc35f98
[project @ 2007-04-04 15:16:46 by jwe]

diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -1186,9 +1186,15 @@
     {
       if (mattype.is_hermitian ())
 	{
-	  ComplexCHOL chol (*this, info);
+	  ComplexCHOL chol (*this, info, calc_cond);
 	  if (info == 0)
-	    ret = chol.inverse ();
+	    {
+	      if (calc_cond)
+		rcond = chol.rcond();
+	      else
+		rcond = 1.0;
+	      ret = chol.inverse ();
+	    }
 	  else
 	    mattype.mark_as_unsymmetric ();
 	}
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,26 @@
+2007-04-04  David Bateman  <dbateman@free.fr>
+
+	* dMatrix.cc (Matrix::inverse): If calc_cond is true, calculate
+	the condition number for positive definite matrices.
+	* CMatrix.cc (ComplexMatrix::inverse): Ditto.
+	* dbleChol.h (CHOL(const Matrix&, bool)): New arg, calc_cond.
+	(CHOL(const Matrix&, octave_idx_type&, bool): Ditto.
+	(octave_idx_type init (const Matrix&, bool)): Ditto.
+	(CHOL(const CHOL&)): Copy xrcond.
+	(CHOL& operator = (const CHOL&)): Copy xrcond.
+	(xrcond): New private data member.
+	* CmplxCHOL.h (ComplexCHOL(const ComplexMatrix&, bool)): New arg,
+	calc_cond.
+	(ComplexCHOL(const ComplexMatrix&, octave_idx_type&, bool): Ditto
+	(octave_idx_type init (const ComplexMatrix&, bool)): Ditto.
+	(ComplexCHOL(const ComplexCHOL&)): Copy xrcond.
+	(ComplexCHOL& operator = (const ComplexCHOL&)): Copy xrcond.
+	(xrcond): New private data member.
+	* dbleCHOL.cc (CHOL::init(const Matrix&, bool)): If calc_cond is
+	true, calculate the condition number with dpocon.
+	* CmplxCHOL.cc (ComplexCHOL::init(const ComplexMatrix&, bool)): If
+	calc_cond is true, calculate the condition number with zpocon.
+
 2007-04-03  John W. Eaton  <jwe@octave.org>
 
 	* intNDArray.cc (intNDArray): Delete spurious semicolon.
diff --git a/liboctave/CmplxCHOL.cc b/liboctave/CmplxCHOL.cc
--- a/liboctave/CmplxCHOL.cc
+++ b/liboctave/CmplxCHOL.cc
@@ -25,6 +25,8 @@
 #include <config.h>
 #endif
 
+#include "dMatrix.h"
+#include "dRowVector.h"
 #include "CmplxCHOL.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
@@ -39,10 +41,16 @@
   F77_FUNC (zpotri, ZPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
 			     Complex*, const octave_idx_type&, octave_idx_type&
 			     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+			     Complex*, const octave_idx_type&, const double&,
+			     double&, Complex*, double*, 
+			     octave_idx_type& F77_CHAR_ARG_LEN_DECL);
 }
 
 octave_idx_type
-ComplexCHOL::init (const ComplexMatrix& a)
+ComplexCHOL::init (const ComplexMatrix& a, bool calc_cond)
 {
   octave_idx_type a_nr = a.rows ();
   octave_idx_type a_nc = a.cols ();
@@ -60,6 +68,11 @@
   chol_mat = a;
   Complex *h = chol_mat.fortran_vec ();
 
+  // Calculate the norm of the matrix, for later use.
+  double anorm = 0;
+  if (calc_cond) 
+    anorm = chol_mat.abs().sum().row(static_cast<octave_idx_type>(0)).max();
+
   F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
 			     F77_CHAR_ARG_LEN (1)));
 
@@ -67,13 +80,39 @@
     (*current_liboctave_error_handler) ("unrecoverable error in zpotrf");
   else
     {
-      // If someone thinks of a more graceful way of doing this (or
-      // faster for that matter :-)), please let me know!
+      xrcond = 0.0;
+      if (info != 0)
+	info = -1;
+      else if (calc_cond) 
+	{
+	  octave_idx_type zpocon_info = 0;
+
+	  // Now calculate the condition number for non-singular matrix.
+	  Array<Complex> z (2*n);
+	  Complex *pz = z.fortran_vec ();
+	  Array<double> rz (n);
+	  double *prz = rz.fortran_vec ();
+	  F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
+				     n, anorm, xrcond, pz, prz, zpocon_info
+				     F77_CHAR_ARG_LEN (1)));
 
-      if (n > 1)
-	for (octave_idx_type j = 0; j < a_nc; j++)
-	  for (octave_idx_type i = j+1; i < a_nr; i++)
-	    chol_mat.xelem (i, j) = 0.0;
+	  if (f77_exception_encountered)
+	    (*current_liboctave_error_handler) 
+	      ("unrecoverable error in zpocon");
+
+	  if (zpocon_info != 0) 
+	    info = -1;
+	}
+      else
+	{
+	  // If someone thinks of a more graceful way of doing this (or
+	  // faster for that matter :-)), please let me know!
+
+	  if (n > 1)
+	    for (octave_idx_type j = 0; j < a_nc; j++)
+	      for (octave_idx_type i = j+1; i < a_nr; i++)
+		chol_mat.xelem (i, j) = 0.0;
+	}
     }
 
   return info;
diff --git a/liboctave/CmplxCHOL.h b/liboctave/CmplxCHOL.h
--- a/liboctave/CmplxCHOL.h
+++ b/liboctave/CmplxCHOL.h
@@ -36,26 +36,31 @@
 
   ComplexCHOL (void) : chol_mat () { }
 
-  ComplexCHOL (const ComplexMatrix& a) { init (a); }
+  ComplexCHOL (const ComplexMatrix& a, bool calc_cond = false) { init (a, calc_cond); }
 
-  ComplexCHOL (const ComplexMatrix& a, octave_idx_type& info)
+  ComplexCHOL (const ComplexMatrix& a, octave_idx_type& info, bool calc_cond = false)
     {
-      info = init (a);
+      info = init (a, calc_cond);
     }
 
   ComplexCHOL (const ComplexCHOL& a)
-    : chol_mat (a.chol_mat) { }
+    : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
 
   ComplexCHOL& operator = (const ComplexCHOL& a)
     {
       if (this != &a)
-	chol_mat = a.chol_mat;
+	{
+	  chol_mat = a.chol_mat;
+	  xrcond = a.xrcond;
+	}
 
       return *this;
     }
 
   ComplexMatrix chol_matrix (void) const { return chol_mat; }
 
+  double rcond (void) const { return xrcond; }
+
   ComplexMatrix inverse (void) const;
 
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ComplexCHOL& a);
@@ -64,7 +69,9 @@
 
   ComplexMatrix chol_mat;
 
-  octave_idx_type init (const ComplexMatrix& a);
+  double xrcond;
+
+  octave_idx_type init (const ComplexMatrix& a, bool calc_cond);
 };
 
 ComplexMatrix OCTAVE_API chol2inv (const ComplexMatrix& r);
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -855,9 +855,15 @@
     {
       if (mattype.is_hermitian ())
 	{
-	  CHOL chol (*this, info);
+	  CHOL chol (*this, info, calc_cond);
 	  if (info == 0)
-	    ret = chol.inverse ();
+	    {
+	      if (calc_cond)
+		rcond = chol.rcond ();
+	      else
+		rcond = 1.0;
+	      ret = chol.inverse ();
+	    }
 	  else
 	    mattype.mark_as_unsymmetric ();
 	}
diff --git a/liboctave/dbleCHOL.cc b/liboctave/dbleCHOL.cc
--- a/liboctave/dbleCHOL.cc
+++ b/liboctave/dbleCHOL.cc
@@ -25,6 +25,7 @@
 #include <config.h>
 #endif
 
+#include "dRowVector.h"
 #include "dbleCHOL.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
@@ -40,10 +41,16 @@
   F77_FUNC (dpotri, DPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
 			     double*, const octave_idx_type&, octave_idx_type&
 			     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+			     double*, const octave_idx_type&, const double&,
+			     double&, double*, octave_idx_type*, 
+			     octave_idx_type& F77_CHAR_ARG_LEN_DECL);
 }
 
 octave_idx_type
-CHOL::init (const Matrix& a)
+CHOL::init (const Matrix& a, bool calc_cond)
 {
   octave_idx_type a_nr = a.rows ();
   octave_idx_type a_nc = a.cols ();
@@ -60,6 +67,11 @@
   chol_mat = a;
   double *h = chol_mat.fortran_vec ();
 
+  // Calculate the norm of the matrix, for later use.
+  double anorm = 0;
+  if (calc_cond) 
+    anorm = chol_mat.abs().sum().row(static_cast<octave_idx_type>(0)).max();
+
   F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
 			     n, h, n, info
 			     F77_CHAR_ARG_LEN (1)));
@@ -68,13 +80,39 @@
     (*current_liboctave_error_handler) ("unrecoverable error in dpotrf");
   else
     {
-      // If someone thinks of a more graceful way of doing this (or
-      // faster for that matter :-)), please let me know!
+      xrcond = 0.0;
+      if (info != 0)
+	info = -1;
+      else if (calc_cond) 
+	{
+	  octave_idx_type dpocon_info = 0;
+
+	  // Now calculate the condition number for non-singular matrix.
+	  Array<double> z (3*n);
+	  double *pz = z.fortran_vec ();
+	  Array<octave_idx_type> iz (n);
+	  octave_idx_type *piz = iz.fortran_vec ();
+	  F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
+				     n, anorm, xrcond, pz, piz, dpocon_info
+				     F77_CHAR_ARG_LEN (1)));
 
-      if (n > 1)
-	for (octave_idx_type j = 0; j < a_nc; j++)
-	  for (octave_idx_type i = j+1; i < a_nr; i++)
-	    chol_mat.xelem (i, j) = 0.0;
+	  if (f77_exception_encountered)
+	    (*current_liboctave_error_handler) 
+	      ("unrecoverable error in dpocon");
+
+	  if (dpocon_info != 0) 
+	    info = -1;
+	}
+      else
+	{
+	  // If someone thinks of a more graceful way of doing this (or
+	  // faster for that matter :-)), please let me know!
+
+	  if (n > 1)
+	    for (octave_idx_type j = 0; j < a_nc; j++)
+	      for (octave_idx_type i = j+1; i < a_nr; i++)
+		chol_mat.xelem (i, j) = 0.0;
+	}
     }
 
   return info;
@@ -91,27 +129,32 @@
   if (r_nr == r_nc)
     {
       octave_idx_type n = r_nc;
-      octave_idx_type info;
+      octave_idx_type info = 0;
 
       Matrix tmp = r;
+      double *v = tmp.fortran_vec();
 
-      F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
-				 tmp.fortran_vec (), n, info
-				 F77_CHAR_ARG_LEN (1)));
-
-      if (f77_exception_encountered)
-	(*current_liboctave_error_handler) ("unrecoverable error in dpotri");
-      else
+      if (info == 0)
 	{
-	  // If someone thinks of a more graceful way of doing this (or
-	  // faster for that matter :-)), please let me know!
+	  F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
+				     v, n, info
+				     F77_CHAR_ARG_LEN (1)));
 
-	  if (n > 1)
-	    for (octave_idx_type j = 0; j < r_nc; j++)
-	      for (octave_idx_type i = j+1; i < r_nr; i++)
-		tmp.xelem (i, j) = tmp.xelem (j, i);
+	  if (f77_exception_encountered)
+	    (*current_liboctave_error_handler) 
+	      ("unrecoverable error in dpotri");
+	  else
+	    {
+	      // If someone thinks of a more graceful way of doing this (or
+	      // faster for that matter :-)), please let me know!
 
-	  retval = tmp;
+	      if (n > 1)
+		for (octave_idx_type j = 0; j < r_nc; j++)
+		  for (octave_idx_type i = j+1; i < r_nr; i++)
+		    tmp.xelem (i, j) = tmp.xelem (j, i);
+
+	      retval = tmp;
+	    }
 	}
     }
   else
diff --git a/liboctave/dbleCHOL.h b/liboctave/dbleCHOL.h
--- a/liboctave/dbleCHOL.h
+++ b/liboctave/dbleCHOL.h
@@ -36,22 +36,27 @@
 
   CHOL (void) : chol_mat () { }
 
-  CHOL (const Matrix& a) { init (a); }
+  CHOL (const Matrix& a, bool calc_cond = false) { init (a, calc_cond); }
 
-  CHOL (const Matrix& a, octave_idx_type& info) { info = init (a); }
+  CHOL (const Matrix& a, octave_idx_type& info, bool calc_cond = false) 
+    { info = init (a, calc_cond); }
 
-  CHOL (const CHOL& a) : chol_mat (a.chol_mat) { }
+  CHOL (const CHOL& a) : chol_mat (a.chol_mat), xrcond (a.xrcond) { }
 
   CHOL& operator = (const CHOL& a)
     {
       if (this != &a)
-	chol_mat = a.chol_mat;
-
+	{
+	  chol_mat = a.chol_mat;
+	  xrcond = a.xrcond;
+	}
       return *this;
     }
 
   Matrix chol_matrix (void) const { return chol_mat; }
 
+  double rcond (void) const { return xrcond; }
+
   // Compute the inverse of a matrix using the Cholesky factorization.
   Matrix inverse (void) const;
 
@@ -61,7 +66,9 @@
 
   Matrix chol_mat;
 
-  octave_idx_type init (const Matrix& a);
+  double xrcond;
+
+  octave_idx_type init (const Matrix& a, bool calc_cond);
 };
 
 Matrix OCTAVE_API chol2inv (const Matrix& r);