# HG changeset patch
# User dbateman
# Date 1165436356 0
# Node ID 3c92b8d892dd62c20b33bdeeaa1cb0f15d2c6d42
# Parent  cb8c62c78b42f7176bf9b96a339ed98efc07af71
[project @ 2006-12-06 20:19:14 by dbateman]

diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -41,6 +41,7 @@
 #include "CmplxDET.h"
 #include "CmplxSCHUR.h"
 #include "CmplxSVD.h"
+#include "CmplxCHOL.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
 #include "lo-ieee.h"
@@ -142,6 +143,13 @@
 			     F77_CHAR_ARG_LEN_DECL);
 
   F77_RET_T
+  F77_FUNC (ztrtri, ZTRTRI) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, 
+			     const octave_idx_type&, const Complex*, 
+			     const octave_idx_type&, octave_idx_type& 
+			     F77_CHAR_ARG_LEN_DECL
+			     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
   F77_FUNC (ztrcon, ZTRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, 
 			     F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, 
 			     const Complex*, const octave_idx_type&, double&,
@@ -959,19 +967,94 @@
 {
   octave_idx_type info;
   double rcond;
-  return inverse (info, rcond, 0, 0);
+  MatrixType mattype (*this);
+  return inverse (mattype, info, rcond, 0, 0);
+}
+
+ComplexMatrix
+ComplexMatrix::inverse (MatrixType &mattype) const
+{
+  octave_idx_type info;
+  double rcond;
+  return inverse (mattype, info, rcond, 0, 0);
+}
+
+ComplexMatrix
+ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const
+{
+  double rcond;
+  return inverse (mattype, info, rcond, 0, 0);
 }
 
 ComplexMatrix
-ComplexMatrix::inverse (octave_idx_type& info) const
+ComplexMatrix::tinverse (MatrixType &mattype, octave_idx_type& info,
+			 double& rcond, int force, int calc_cond) const
 {
-  double rcond;
-  return inverse (info, rcond, 0, 0);
+  ComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+
+  if (nr != nc || nr == 0 || nc == 0)
+    (*current_liboctave_error_handler) ("inverse requires square matrix");
+  else
+    {
+      int typ = mattype.type ();
+      char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
+      char udiag = 'N';
+      retval = *this;
+      Complex *tmp_data = retval.fortran_vec ();
+
+      F77_XFCN (ztrtri, ZTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),
+				 F77_CONST_CHAR_ARG2 (&udiag, 1),
+				 nr, tmp_data, nr, info 
+				 F77_CHAR_ARG_LEN (1)
+				 F77_CHAR_ARG_LEN (1)));
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in ztrtri");
+      else
+	{
+	  // Throw-away extra info LAPACK gives so as to not change output.
+	  rcond = 0.0;
+	  if (info != 0) 
+	    info = -1;
+	  else if (calc_cond) 
+	    {
+	      octave_idx_type ztrcon_info = 0;
+	      char job = '1';
+
+	      OCTAVE_LOCAL_BUFFER (Complex, cwork, 2 * nr);
+	      OCTAVE_LOCAL_BUFFER (double, rwork, nr);
+
+	      F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
+					 F77_CONST_CHAR_ARG2 (&uplo, 1),
+					 F77_CONST_CHAR_ARG2 (&udiag, 1),
+					 nr, tmp_data, nr, rcond, 
+					 cwork, rwork, ztrcon_info 
+					 F77_CHAR_ARG_LEN (1)
+					 F77_CHAR_ARG_LEN (1)
+					 F77_CHAR_ARG_LEN (1)));
+
+	      if (f77_exception_encountered)
+		(*current_liboctave_error_handler) 
+		  ("unrecoverable error in ztrcon");
+
+	      if (ztrcon_info != 0) 
+		info = -1;
+	    }
+	}
+
+      if (info == -1 && ! force)
+	retval = *this; // Restore matrix contents.
+    }
+
+  return retval;
 }
 
 ComplexMatrix
-ComplexMatrix::inverse (octave_idx_type& info, double& rcond, int force, 
-			int calc_cond) const
+ComplexMatrix::finverse (MatrixType &mattype, octave_idx_type& info,
+			 double& rcond, int force, int calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -1062,12 +1145,45 @@
 		info = -1;
 	    }
 	}
+
+      if (info != 0)
+	mattype.mark_as_rectangular();
     }
   
   return retval;
 }
 
 ComplexMatrix
+ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info,
+			double& rcond, int force, int calc_cond) const
+{
+  int typ = mattype.type (false);
+  ComplexMatrix ret;
+
+  if (typ == MatrixType::Unknown)
+    typ = mattype.type (*this);
+
+  if (typ == MatrixType::Upper || typ == MatrixType::Lower)
+    ret = tinverse (mattype, info, rcond, force, calc_cond);
+  else if (typ != MatrixType::Rectangular)
+    {
+      if (mattype.is_hermitian ())
+	{
+	  ComplexCHOL chol (*this, info);
+	  if (info == 0)
+	    ret = chol.inverse ();
+	  else
+	    mattype.mark_as_unsymmetric ();
+	}
+
+      if (!mattype.is_hermitian ())
+	ret = finverse(mattype, info, rcond, force, calc_cond);
+    }
+
+  return ret;
+}
+
+ComplexMatrix
 ComplexMatrix::pseudo_inverse (double tol) const
 {
   ComplexMatrix retval;
diff --git a/liboctave/CMatrix.h b/liboctave/CMatrix.h
--- a/liboctave/CMatrix.h
+++ b/liboctave/CMatrix.h
@@ -135,9 +135,19 @@
 
   ComplexColumnVector column (octave_idx_type i) const;
 
+private:
+  ComplexMatrix tinverse (MatrixType &mattype, octave_idx_type& info,
+			  double& rcond, int force, int calc_cond) const;
+
+  ComplexMatrix finverse (MatrixType &mattype, octave_idx_type& info,
+			  double& rcond, int force, int calc_cond) const;
+
+public:
   ComplexMatrix inverse (void) const;
-  ComplexMatrix inverse (octave_idx_type& info) const;
-  ComplexMatrix inverse (octave_idx_type& info, double& rcond, int force = 0,
+  ComplexMatrix inverse (MatrixType &mattype) const;
+  ComplexMatrix inverse (MatrixType &mattype, octave_idx_type& info) const;
+  ComplexMatrix inverse (MatrixType &mattype, octave_idx_type& info,
+			 double& rcond, int force = 0, 
 			 int calc_cond = 1) const;
 
   ComplexMatrix pseudo_inverse (double tol = 0.0) const;
diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc
--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -180,12 +180,33 @@
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
 
-  if (is_square () && nr > 0)
+  if (nr == nc && nr > 0)
     {
-      for (octave_idx_type i = 0; i < nr; i++)
-	for (octave_idx_type j = i; j < nc; j++)
-	  if (elem (i, j) != conj (elem (j, i)))
-	    return false;
+      for (octave_idx_type j = 0; j < nc; j++)
+	{
+	  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+	    {
+	      octave_idx_type ri = ridx(i);
+
+	      if (ri != j)
+		{
+		  bool found = false;
+
+		  for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++)
+		    {
+		      if (ridx(k) == j)
+			{
+			  if (data(i) == conj(data(k)))
+			    found = true;
+			  break;
+			}
+		    }
+
+		  if (! found)
+		    return false;
+		}
+	    }
+	}
 
       return true;
     }
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,17 @@
+2006-12-06  David Bateman  <dbateman@free.fr>
+
+	* dSparse.cc (SparseMatrix::is_symmetric(void) const): Faster
+	implementation.
+	* CSparse.cc (SparseComplexMatrix::is_symmetric(void) const): ditto.
+
+	* dMatrrix.cc (finverse): Old inverse method renamed inverse.
+        (tinverse): New method for triangular matrices.
+        (inverse): New function with matrix type probing.
+        * dMatrix.h (finverse, tinverse, inverse): New and modified
+        declarations.
+        * CMatrix.cc: ditto.
+        * CMatrix.h: ditto.
+
 2006-12-06  John W. Eaton  <jwe@octave.org>
 
 	* strptime.c (day_of_the_week): Use code from current glibc sources.
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -37,6 +37,7 @@
 #include "dbleDET.h"
 #include "dbleSCHUR.h"
 #include "dbleSVD.h"
+#include "dbleCHOL.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
 #include "lo-ieee.h"
@@ -138,6 +139,12 @@
 			     F77_CHAR_ARG_LEN_DECL);
 
   F77_RET_T
+  F77_FUNC (dtrtri, DTRTRI) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, 
+			     const octave_idx_type&, const double*, 
+			     const octave_idx_type&, octave_idx_type&
+			     F77_CHAR_ARG_LEN_DECL
+			     F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
   F77_FUNC (dtrcon, DTRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, 
 			     F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, 
 			     const double*, const octave_idx_type&, double&,
@@ -628,18 +635,95 @@
 {
   octave_idx_type info;
   double rcond;
-  return inverse (info, rcond, 0, 0);
+  MatrixType mattype (*this);
+  return inverse (mattype, info, rcond, 0, 0);
+}
+
+Matrix
+Matrix::inverse (MatrixType& mattype) const
+{
+  octave_idx_type info;
+  double rcond;
+  return inverse (mattype, info, rcond, 0, 0);
+}
+
+Matrix
+Matrix::inverse (MatrixType &mattype, octave_idx_type& info) const
+{
+  double rcond;
+  return inverse (mattype, info, rcond, 0, 0);
 }
 
 Matrix
-Matrix::inverse (octave_idx_type& info) const
+Matrix::tinverse (MatrixType &mattype, octave_idx_type& info, double& rcond, 
+		  int force, int calc_cond) const
 {
-  double rcond;
-  return inverse (info, rcond, 0, 0);
+  Matrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+
+  if (nr != nc || nr == 0 || nc == 0)
+    (*current_liboctave_error_handler) ("inverse requires square matrix");
+  else
+    {
+      int typ = mattype.type ();
+      char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
+      char udiag = 'N';
+      retval = *this;
+      double *tmp_data = retval.fortran_vec ();
+
+      F77_XFCN (dtrtri, DTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),
+				 F77_CONST_CHAR_ARG2 (&udiag, 1),
+				 nr, tmp_data, nr, info 
+				 F77_CHAR_ARG_LEN (1)
+				 F77_CHAR_ARG_LEN (1)));
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in dtrtri");
+      else
+	{
+	  // Throw-away extra info LAPACK gives so as to not change output.
+	  rcond = 0.0;
+	  if (info != 0) 
+	    info = -1;
+	  else if (calc_cond) 
+	    {
+	      octave_idx_type dtrcon_info = 0;
+	      char job = '1';
+
+	      OCTAVE_LOCAL_BUFFER (double, work, 3 * nr);
+	      OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr);
+
+	      F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
+					 F77_CONST_CHAR_ARG2 (&uplo, 1),
+					 F77_CONST_CHAR_ARG2 (&udiag, 1),
+					 nr, tmp_data, nr, rcond, 
+					 work, iwork, dtrcon_info 
+					 F77_CHAR_ARG_LEN (1)
+					 F77_CHAR_ARG_LEN (1)
+					 F77_CHAR_ARG_LEN (1)));
+
+	      if (f77_exception_encountered)
+		(*current_liboctave_error_handler) 
+		  ("unrecoverable error in dtrcon");
+
+	      if (dtrcon_info != 0) 
+		info = -1;
+	    }
+	}
+
+      if (info == -1 && ! force)
+	retval = *this; // Restore matrix contents.
+    }
+
+  return retval;
 }
 
+
 Matrix
-Matrix::inverse (octave_idx_type& info, double& rcond, int force, int calc_cond) const
+Matrix::finverse (MatrixType &mattype, octave_idx_type& info, double& rcond, 
+		  int force, int calc_cond) const
 {
   Matrix retval;
 
@@ -730,12 +814,45 @@
 		info = -1;
 	    }
 	}
+
+      if (info != 0)
+	mattype.mark_as_rectangular();
     }
 
   return retval;
 }
 
 Matrix
+Matrix::inverse (MatrixType &mattype, octave_idx_type& info, double& rcond, 
+		 int force, int calc_cond) const
+{
+  int typ = mattype.type (false);
+  Matrix ret;
+
+  if (typ == MatrixType::Unknown)
+    typ = mattype.type (*this);
+
+  if (typ == MatrixType::Upper || typ == MatrixType::Lower)
+    ret = tinverse (mattype, info, rcond, force, calc_cond);
+  else if (typ != MatrixType::Rectangular)
+    {
+      if (mattype.is_hermitian ())
+	{
+	  CHOL chol (*this, info);
+	  if (info == 0)
+	    ret = chol.inverse ();
+	  else
+	    mattype.mark_as_unsymmetric ();
+	}
+
+      if (!mattype.is_hermitian ())
+	ret = finverse(mattype, info, rcond, force, calc_cond);
+    }
+
+  return ret;
+}
+
+Matrix
 Matrix::pseudo_inverse (double tol) const
 {
   SVD result (*this, SVD::economy);
diff --git a/liboctave/dMatrix.h b/liboctave/dMatrix.h
--- a/liboctave/dMatrix.h
+++ b/liboctave/dMatrix.h
@@ -107,10 +107,19 @@
 
   ColumnVector column (octave_idx_type i) const;
 
+private:
+  Matrix tinverse (MatrixType &mattype, octave_idx_type& info, double& rcond, 
+		   int force, int calc_cond) const;
+
+  Matrix finverse (MatrixType &mattype, octave_idx_type& info, double& rcond, 
+		   int force, int calc_cond) const;
+
+public:
   Matrix inverse (void) const;
-  Matrix inverse (octave_idx_type& info) const;
-  Matrix inverse (octave_idx_type& info, double& rcond, int force = 0, 
-		  int calc_cond = 1) const;
+  Matrix inverse (MatrixType &mattype) const;
+  Matrix inverse (MatrixType &mattype, octave_idx_type& info) const;
+  Matrix inverse (MatrixType &mattype, octave_idx_type& info, double& rcond,
+		  int force = 0, int calc_cond = 1) const;
 
   Matrix pseudo_inverse (double tol = 0.0) const;
 
diff --git a/liboctave/dSparse.cc b/liboctave/dSparse.cc
--- a/liboctave/dSparse.cc
+++ b/liboctave/dSparse.cc
@@ -170,12 +170,36 @@
 bool
 SparseMatrix::is_symmetric (void) const
 {
-  if (is_square () && rows () > 0)
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+
+  if (nr == nc && nr > 0)
     {
-      for (octave_idx_type i = 0; i < rows (); i++)
-	for (octave_idx_type j = i+1; j < cols (); j++)
-	  if (elem (i, j) != elem (j, i))
-	    return false;
+      for (octave_idx_type j = 0; j < nc; j++)
+	{
+	  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+	    {
+	      octave_idx_type ri = ridx(i);
+
+	      if (ri != j)
+		{
+		  bool found = false;
+
+		  for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++)
+		    {
+		      if (ridx(k) == j)
+			{
+			  if (data(i) == data(k))
+			    found = true;
+			  break;
+			}
+		    }
+
+		  if (! found)
+		    return false;
+		}
+	    }
+	}
 
       return true;
     }
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,10 @@
+2006-12-04  David Bateman  <dbateman@free.fr>
+
+	* xpow.cc (xpow (const Matrix&, double)): Add matrix type probing
+        to matrix inverse.
+        (xpow (const ComplexMatrix&, double)): ditto.
+        * DLD-FUNCTIONS/inv.cc (Finv): Add matrix type probing.
+
 2006-12-06  John W. Eaton  <jwe@octave.org>
 
 2006-12-06  Michael Goffioul  <michael.goffioul@swing.be>
diff --git a/src/DLD-FUNCTIONS/inv.cc b/src/DLD-FUNCTIONS/inv.cc
--- a/src/DLD-FUNCTIONS/inv.cc
+++ b/src/DLD-FUNCTIONS/inv.cc
@@ -74,10 +74,14 @@
 
       if (! error_state)
 	{
+	  MatrixType mattyp = args(0).matrix_type ();
+
 	  octave_idx_type info;
 	  double rcond = 0.0;
 
-	  Matrix result = m.inverse (info, rcond, 1);
+	  Matrix result = m.inverse (mattyp, info, rcond, 1);
+
+	  args(0).matrix_type (mattyp);
 
 	  if (nargout > 1)
 	    retval(1) = rcond;
@@ -97,10 +101,14 @@
 
       if (! error_state)
 	{
+	  MatrixType mattyp = args(0).matrix_type ();
+
 	  octave_idx_type info;
 	  double rcond = 0.0;
 
-	  ComplexMatrix result = m.inverse (info, rcond, 1);
+	  ComplexMatrix result = m.inverse (mattyp, info, rcond, 1);
+
+	  args(0).matrix_type (mattyp);
 
 	  if (nargout > 1)
 	    retval(1) = rcond;
diff --git a/src/xpow.cc b/src/xpow.cc
--- a/src/xpow.cc
+++ b/src/xpow.cc
@@ -192,8 +192,9 @@
 
 		  octave_idx_type info;
 		  double rcond = 0.0;
+		  MatrixType mattype (a);
 
-		  atmp = a.inverse (info, rcond, 1);
+		  atmp = a.inverse (mattype, info, rcond, 1);
 
 		  if (info == -1)
 		    warning ("inverse: matrix singular to machine\
@@ -388,8 +389,9 @@
 
 		  octave_idx_type info;
 		  double rcond = 0.0;
+		  MatrixType mattype (a);
 
-		  atmp = a.inverse (info, rcond, 1);
+		  atmp = a.inverse (mattype, info, rcond, 1);
 
 		  if (info == -1)
 		    warning ("inverse: matrix singular to machine\