# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1227104799 -3600
# Node ID 9813c07ca946de9bd60f97968005ea1c407482b8
# Parent  64cf956a109cc310f236fdd0d4f275e06d8b67ff
make det take advantage of matrix type

diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -56,6 +56,7 @@
 #include "mx-cm-s.h"
 #include "mx-inlines.cc"
 #include "oct-cmplx.h"
+#include "oct-norm.h"
 
 #if defined (HAVE_FFTW3)
 #include "oct-fftw.h"
@@ -1574,74 +1575,132 @@
 ComplexDET
 ComplexMatrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const
 {
-  ComplexDET retval;
+  MatrixType mattype (*this);
+  return determinant (mattype, info, rcon, calc_cond);
+}
+
+ComplexDET
+ComplexMatrix::determinant (MatrixType& mattype,
+                            octave_idx_type& info, double& rcon, int calc_cond) const
+{
+  ComplexDET retval (1.0);
 
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
 
-  if (nr == 0 || nc == 0)
-    {
-      retval = ComplexDET (1.0, 0);
-    }
+  if (nr != nc)
+    (*current_liboctave_error_handler) ("matrix must be square");
   else
     {
-      Array<octave_idx_type> ipvt (nr);
-      octave_idx_type *pipvt = ipvt.fortran_vec ();
-
-      ComplexMatrix atmp = *this;
-      Complex *tmp_data = atmp.fortran_vec ();
-
-      info = 0;
-
-      // Calculate the norm of the matrix, for later use.
-      double anorm = 0;
-      if (calc_cond) 
-	anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
-
-      F77_XFCN (zgetrf, ZGETRF, (nr, nc, tmp_data, nr, pipvt, info));
-
-      // Throw-away extra info LAPACK gives so as to not change output.
-      rcon = 0.0;
-      if (info != 0) 
-	{
-	  info = -1;
-	  retval = ComplexDET ();
-	} 
-      else 
-	{
-	  if (calc_cond) 
-	    {
-	      // Now calc the condition number for non-singular matrix.
-	      char job = '1';
-	      Array<Complex> z (2*nr);
-	      Complex *pz = z.fortran_vec ();
-	      Array<double> rz (2*nr);
-	      double *prz = rz.fortran_vec ();
-
-	      F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
-					 nc, tmp_data, nr, anorm, 
-					 rcon, pz, prz, info
-					 F77_CHAR_ARG_LEN (1)));
-	    }
-
-	  if (info != 0) 
-	    {
-	      info = -1;
-	      retval = ComplexDET ();
-	    } 
-	  else 
-	    {
-              retval = ComplexDET (1.0);
-              
-	      for (octave_idx_type i = 0; i < nc; i++) 
-		{
-                  Complex c = atmp(i,i);
-                  retval *= (ipvt(i) != (i+1)) ? -c : c;
+      int typ = mattype.type ();
+
+      if (typ == MatrixType::Lower || typ == MatrixType::Upper)
+        {
+          for (octave_idx_type i = 0; i < nc; i++) 
+            retval *= elem (i,i);
+        }
+      else if (typ == MatrixType::Hermitian)
+        {
+          ComplexMatrix atmp = *this;
+          Complex *tmp_data = atmp.fortran_vec ();
+
+          info = 0;
+          double anorm = 0;
+          if (calc_cond) anorm = xnorm (*this, 1);
+
+
+          char job = 'L';
+          F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 
+                                     tmp_data, nr, info
+                                     F77_CHAR_ARG_LEN (1)));
+
+          if (info != 0) 
+            {
+              rcon = 0.0;
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+          else 
+            {
+              Array<Complex> z (2 * nc);
+              Complex *pz = z.fortran_vec ();
+              Array<double> rz (nc);
+              double *prz = rz.fortran_vec ();
+
+              F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                         nr, tmp_data, nr, anorm,
+                                         rcon, pz, prz, info
+                                         F77_CHAR_ARG_LEN (1)));
+
+              if (info != 0) 
+                rcon = 0.0;
+
+              for (octave_idx_type i = 0; i < nc; i++) 
+                retval *= elem (i,i);
+
+              retval = retval.square ();
+            }
+        }
+      else if (typ != MatrixType::Full)
+        (*current_liboctave_error_handler) ("det: invalid dense matrix type");
+
+      if (typ == MatrixType::Full)
+        {
+          Array<octave_idx_type> ipvt (nr);
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          ComplexMatrix atmp = *this;
+          Complex *tmp_data = atmp.fortran_vec ();
+
+          info = 0;
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm = 0;
+          if (calc_cond) anorm = xnorm (*this, 1);
+
+          F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
+
+          // Throw-away extra info LAPACK gives so as to not change output.
+          rcon = 0.0;
+          if (info != 0) 
+            {
+              info = -1;
+              retval = ComplexDET ();
+            } 
+          else 
+            {
+              if (calc_cond) 
+                {
+                  // Now calc the condition number for non-singular matrix.
+                  char job = '1';
+                  Array<Complex> z (2 * nc);
+                  Complex *pz = z.fortran_vec ();
+                  Array<double> rz (2 * nc);
+                  double *prz = rz.fortran_vec ();
+
+                  F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                             nc, tmp_data, nr, anorm, 
+                                             rcon, pz, prz, info
+                                             F77_CHAR_ARG_LEN (1)));
                 }
-	    }
-	}
+
+              if (info != 0) 
+                {
+                  info = -1;
+                  retval = ComplexDET ();
+                } 
+              else 
+                {
+                  for (octave_idx_type i = 0; i < nc; i++) 
+                    {
+                      Complex c = atmp(i,i);
+                      retval *= (ipvt(i) != (i+1)) ? -c : c;
+                    }
+                }
+            }
+        }
     }
-  
+
   return retval;
 }
 
diff --git a/liboctave/CMatrix.h b/liboctave/CMatrix.h
--- a/liboctave/CMatrix.h
+++ b/liboctave/CMatrix.h
@@ -176,6 +176,8 @@
   ComplexDET determinant (void) const;
   ComplexDET determinant (octave_idx_type& info) const;
   ComplexDET determinant (octave_idx_type& info, double& rcon, int calc_cond = 1) const;
+  ComplexDET determinant (MatrixType &mattype, octave_idx_type& info, 
+                          double& rcon, int calc_cond = 1) const;
 
   double rcond (void) const;
   double rcond (MatrixType &mattype) const;
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,13 @@
+2008-11-19  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DET.h (base_det<T>::square): New member function.
+	* dMatrix.cc (Matrix::determinant),
+	fMatrix.cc (FloatMatrix::determinant),
+	CMatrix.cc (ComplexMatrix::determinant),
+	fCMatrix.cc (FloatComplexMatrix::determinant):
+	Allow taking MatrixType argument.
+	* dMatrix.h, fMatrix.h, CMatrix.h, fCMatrix.h: Update decls.
+
 2008-11-19  Jaroslav Hajek  <highegg@gmail.com>
 
 	* DET.h: New source.
diff --git a/liboctave/DET.h b/liboctave/DET.h
--- a/liboctave/DET.h
+++ b/liboctave/DET.h
@@ -34,7 +34,7 @@
 {
 public:
 
-  base_det (T c = 0, int e = 0) 
+  base_det (T c = 1, int e = 0) 
     { 
       c2 = xlog2 (c, e2); 
       e2 += e; 
@@ -65,6 +65,8 @@
   T value () const { return c2 * static_cast<T> (std::ldexp (1.0, e2)); }
   operator T () const { return value (); }
 
+  base_det square () const { return base_det (c2*c2, e2+e2); }
+
   void operator *= (T t)
     {
       int e;
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -51,6 +51,7 @@
 #include "mx-dm-m.h"
 #include "mx-inlines.cc"
 #include "oct-cmplx.h"
+#include "oct-norm.h"
 #include "quit.h"
 
 #if defined (HAVE_FFTW3)
@@ -1240,72 +1241,130 @@
 DET
 Matrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const
 {
-  DET retval;
+  MatrixType mattype (*this);
+  return determinant (mattype, info, rcon, calc_cond);
+}
+
+DET
+Matrix::determinant (MatrixType& mattype,
+                     octave_idx_type& info, double& rcon, int calc_cond) const
+{
+  DET retval (1.0);
 
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
 
-  if (nr == 0 || nc == 0)
-    {
-      retval = DET (1.0, 0);
-    }
+  if (nr != nc)
+    (*current_liboctave_error_handler) ("matrix must be square");
   else
     {
-      Array<octave_idx_type> ipvt (nr);
-      octave_idx_type *pipvt = ipvt.fortran_vec ();
-
-      Matrix atmp = *this;
-      double *tmp_data = atmp.fortran_vec ();
-
-      info = 0;
-
-      // Calculate the norm of the matrix, for later use.
-      double anorm = 0;
-      if (calc_cond) 
-	anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
-
-      F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
-
-      // Throw-away extra info LAPACK gives so as to not change output.
-      rcon = 0.0;
-      if (info != 0) 
-	{
-	  info = -1;
-	  retval = DET ();
-	} 
-      else 
-	{
-	  if (calc_cond) 
-	    {
-	      // Now calc the condition number for non-singular matrix.
-	      char job = '1';
-	      Array<double> z (4 * nc);
-	      double *pz = z.fortran_vec ();
-	      Array<octave_idx_type> iz (nc);
-	      octave_idx_type *piz = iz.fortran_vec ();
-
-	      F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
-					 nc, tmp_data, nr, anorm, 
-					 rcon, pz, piz, info
-					 F77_CHAR_ARG_LEN (1)));
-	    }
-
-	  if (info != 0) 
-	    {
-	      info = -1;
-	      retval = DET ();
-	    } 
-	  else 
-	    {
-              retval = DET (1.0);
-              
-	      for (octave_idx_type i = 0; i < nc; i++) 
-		{
-                  double c = atmp(i,i);
-                  retval *= (ipvt(i) != (i+1)) ? -c : c;
+      int typ = mattype.type ();
+
+      if (typ == MatrixType::Lower || typ == MatrixType::Upper)
+        {
+          for (octave_idx_type i = 0; i < nc; i++) 
+            retval *= elem (i,i);
+        }
+      else if (typ == MatrixType::Hermitian)
+        {
+          Matrix atmp = *this;
+          double *tmp_data = atmp.fortran_vec ();
+
+          info = 0;
+          double anorm = 0;
+          if (calc_cond) anorm = xnorm (*this, 1);
+
+
+          char job = 'L';
+          F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 
+                                     tmp_data, nr, info
+                                     F77_CHAR_ARG_LEN (1)));
+
+          if (info != 0) 
+            {
+              rcon = 0.0;
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+          else 
+            {
+              Array<double> z (3 * nc);
+              double *pz = z.fortran_vec ();
+              Array<octave_idx_type> iz (nc);
+              octave_idx_type *piz = iz.fortran_vec ();
+
+              F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                         nr, tmp_data, nr, anorm,
+                                         rcon, pz, piz, info
+                                         F77_CHAR_ARG_LEN (1)));
+
+              if (info != 0) 
+                rcon = 0.0;
+
+              for (octave_idx_type i = 0; i < nc; i++) 
+                retval *= elem (i,i);
+
+              retval = retval.square ();
+            }
+        }
+      else if (typ != MatrixType::Full)
+        (*current_liboctave_error_handler) ("det: invalid dense matrix type");
+
+      if (typ == MatrixType::Full)
+        {
+          Array<octave_idx_type> ipvt (nr);
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          Matrix atmp = *this;
+          double *tmp_data = atmp.fortran_vec ();
+
+          info = 0;
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm = 0;
+          if (calc_cond) anorm = xnorm (*this, 1);
+
+          F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
+
+          // Throw-away extra info LAPACK gives so as to not change output.
+          rcon = 0.0;
+          if (info != 0) 
+            {
+              info = -1;
+              retval = DET ();
+            } 
+          else 
+            {
+              if (calc_cond) 
+                {
+                  // Now calc the condition number for non-singular matrix.
+                  char job = '1';
+                  Array<double> z (4 * nc);
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (nc);
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                             nc, tmp_data, nr, anorm, 
+                                             rcon, pz, piz, info
+                                             F77_CHAR_ARG_LEN (1)));
                 }
-	    }
-	}
+
+              if (info != 0) 
+                {
+                  info = -1;
+                  retval = DET ();
+                } 
+              else 
+                {
+                  for (octave_idx_type i = 0; i < nc; i++) 
+                    {
+                      double c = atmp(i,i);
+                      retval *= (ipvt(i) != (i+1)) ? -c : c;
+                    }
+                }
+            }
+        }
     }
 
   return retval;
diff --git a/liboctave/dMatrix.h b/liboctave/dMatrix.h
--- a/liboctave/dMatrix.h
+++ b/liboctave/dMatrix.h
@@ -145,6 +145,8 @@
   DET determinant (void) const;
   DET determinant (octave_idx_type& info) const;
   DET determinant (octave_idx_type& info, double& rcon, int calc_cond = 1) const;
+  DET determinant (MatrixType &mattype, octave_idx_type& info, 
+                   double& rcon, int calc_cond = 1) const;
 
   double rcond (void) const;
   double rcond (MatrixType &mattype) const;
diff --git a/liboctave/fCMatrix.cc b/liboctave/fCMatrix.cc
--- a/liboctave/fCMatrix.cc
+++ b/liboctave/fCMatrix.cc
@@ -3,6 +3,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
               2003, 2004, 2005, 2006, 2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -54,6 +55,7 @@
 #include "mx-fcm-fs.h"
 #include "mx-inlines.cc"
 #include "oct-cmplx.h"
+#include "oct-norm.h"
 
 #if defined (HAVE_FFTW3)
 #include "oct-fftw.h"
@@ -1567,74 +1569,132 @@
 FloatComplexDET
 FloatComplexMatrix::determinant (octave_idx_type& info, float& rcon, int calc_cond) const
 {
-  FloatComplexDET retval;
+  MatrixType mattype (*this);
+  return determinant (mattype, info, rcon, calc_cond);
+}
+
+FloatComplexDET
+FloatComplexMatrix::determinant (MatrixType& mattype,
+                                 octave_idx_type& info, float& rcon, int calc_cond) const
+{
+  FloatComplexDET retval (1.0);
 
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
 
-  if (nr == 0 || nc == 0)
-    {
-      retval = FloatComplexDET (1.0, 0);
-    }
+  if (nr != nc)
+    (*current_liboctave_error_handler) ("matrix must be square");
   else
     {
-      Array<octave_idx_type> ipvt (nr);
-      octave_idx_type *pipvt = ipvt.fortran_vec ();
-
-      FloatComplexMatrix atmp = *this;
-      FloatComplex *tmp_data = atmp.fortran_vec ();
-
-      info = 0;
-
-      // Calculate the norm of the matrix, for later use.
-      float anorm = 0;
-      if (calc_cond) 
-	anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
-
-      F77_XFCN (cgetrf, CGETRF, (nr, nc, tmp_data, nr, pipvt, info));
-
-      // Throw-away extra info LAPACK gives so as to not change output.
-      rcon = 0.0;
-      if (info != 0) 
-	{
-	  info = -1;
-	  retval = FloatComplexDET ();
-	} 
-      else 
-	{
-	  if (calc_cond) 
-	    {
-	      // Now calc the condition number for non-singular matrix.
-	      char job = '1';
-	      Array<FloatComplex> z (2*nr);
-	      FloatComplex *pz = z.fortran_vec ();
-	      Array<float> rz (2*nr);
-	      float *prz = rz.fortran_vec ();
-
-	      F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
-					 nc, tmp_data, nr, anorm, 
-					 rcon, pz, prz, info
-					 F77_CHAR_ARG_LEN (1)));
-	    }
-
-	  if (info != 0) 
-	    {
-	      info = -1;
-	      retval = FloatComplexDET ();
-	    } 
-	  else 
-	    {
-              retval = FloatComplexDET (1.0);
-              
-	      for (octave_idx_type i = 0; i < nc; i++) 
-		{
-                  FloatComplex c = atmp(i,i);
-                  retval *= (ipvt(i) != (i+1)) ? -c : c;
+      int typ = mattype.type ();
+
+      if (typ == MatrixType::Lower || typ == MatrixType::Upper)
+        {
+          for (octave_idx_type i = 0; i < nc; i++) 
+            retval *= elem (i,i);
+        }
+      else if (typ == MatrixType::Hermitian)
+        {
+          FloatComplexMatrix atmp = *this;
+          FloatComplex *tmp_data = atmp.fortran_vec ();
+
+          info = 0;
+          float anorm = 0;
+          if (calc_cond) anorm = xnorm (*this, 1);
+
+
+          char job = 'L';
+          F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 
+                                     tmp_data, nr, info
+                                     F77_CHAR_ARG_LEN (1)));
+
+          if (info != 0) 
+            {
+              rcon = 0.0;
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+          else 
+            {
+              Array<FloatComplex> z (2 * nc);
+              FloatComplex *pz = z.fortran_vec ();
+              Array<float> rz (nc);
+              float *prz = rz.fortran_vec ();
+
+              F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                         nr, tmp_data, nr, anorm,
+                                         rcon, pz, prz, info
+                                         F77_CHAR_ARG_LEN (1)));
+
+              if (info != 0) 
+                rcon = 0.0;
+
+              for (octave_idx_type i = 0; i < nc; i++) 
+                retval *= elem (i,i);
+
+              retval = retval.square ();
+            }
+        }
+      else if (typ != MatrixType::Full)
+        (*current_liboctave_error_handler) ("det: invalid dense matrix type");
+
+      if (typ == MatrixType::Full)
+        {
+          Array<octave_idx_type> ipvt (nr);
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          FloatComplexMatrix atmp = *this;
+          FloatComplex *tmp_data = atmp.fortran_vec ();
+
+          info = 0;
+
+          // Calculate the norm of the matrix, for later use.
+          float anorm = 0;
+          if (calc_cond) anorm = xnorm (*this, 1);
+
+          F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
+
+          // Throw-away extra info LAPACK gives so as to not change output.
+          rcon = 0.0;
+          if (info != 0) 
+            {
+              info = -1;
+              retval = FloatComplexDET ();
+            } 
+          else 
+            {
+              if (calc_cond) 
+                {
+                  // Now calc the condition number for non-singular matrix.
+                  char job = '1';
+                  Array<FloatComplex> z (2 * nc);
+                  FloatComplex *pz = z.fortran_vec ();
+                  Array<float> rz (2 * nc);
+                  float *prz = rz.fortran_vec ();
+
+                  F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                             nc, tmp_data, nr, anorm, 
+                                             rcon, pz, prz, info
+                                             F77_CHAR_ARG_LEN (1)));
                 }
-	    }
-	}
+
+              if (info != 0) 
+                {
+                  info = -1;
+                  retval = FloatComplexDET ();
+                } 
+              else 
+                {
+                  for (octave_idx_type i = 0; i < nc; i++) 
+                    {
+                      FloatComplex c = atmp(i,i);
+                      retval *= (ipvt(i) != (i+1)) ? -c : c;
+                    }
+                }
+            }
+        }
     }
-  
+
   return retval;
 }
 
diff --git a/liboctave/fCMatrix.h b/liboctave/fCMatrix.h
--- a/liboctave/fCMatrix.h
+++ b/liboctave/fCMatrix.h
@@ -176,6 +176,8 @@
   FloatComplexDET determinant (void) const;
   FloatComplexDET determinant (octave_idx_type& info) const;
   FloatComplexDET determinant (octave_idx_type& info, float& rcon, int calc_cond = 1) const;
+  FloatComplexDET determinant (MatrixType &mattype, octave_idx_type& info, 
+                               float& rcon, int calc_cond = 1) const;
 
   float rcond (void) const;
   float rcond (MatrixType &mattype) const;
diff --git a/liboctave/fMatrix.cc b/liboctave/fMatrix.cc
--- a/liboctave/fMatrix.cc
+++ b/liboctave/fMatrix.cc
@@ -50,6 +50,7 @@
 #include "mx-fdm-fm.h"
 #include "mx-inlines.cc"
 #include "oct-cmplx.h"
+#include "oct-norm.h"
 #include "quit.h"
 
 #if defined (HAVE_FFTW3)
@@ -1239,72 +1240,130 @@
 FloatDET
 FloatMatrix::determinant (octave_idx_type& info, float& rcon, int calc_cond) const
 {
-  FloatDET retval;
+  MatrixType mattype (*this);
+  return determinant (mattype, info, rcon, calc_cond);
+}
+
+FloatDET
+FloatMatrix::determinant (MatrixType& mattype,
+                          octave_idx_type& info, float& rcon, int calc_cond) const
+{
+  FloatDET retval (1.0);
 
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
 
-  if (nr == 0 || nc == 0)
-    {
-      retval = FloatDET (1.0, 0);
-    }
+  if (nr != nc)
+    (*current_liboctave_error_handler) ("matrix must be square");
   else
     {
-      Array<octave_idx_type> ipvt (nr);
-      octave_idx_type *pipvt = ipvt.fortran_vec ();
-
-      FloatMatrix atmp = *this;
-      float *tmp_data = atmp.fortran_vec ();
-
-      info = 0;
-
-      // Calculate the norm of the matrix, for later use.
-      float anorm = 0;
-      if (calc_cond) 
-	anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
-
-      F77_XFCN (sgetrf, SGETRF, (nr, nr, tmp_data, nr, pipvt, info));
-
-      // Throw-away extra info LAPACK gives so as to not change output.
-      rcon = 0.0;
-      if (info != 0) 
-	{
-	  info = -1;
-	  retval = FloatDET ();
-	} 
-      else 
-	{
-	  if (calc_cond) 
-	    {
-	      // Now calc the condition number for non-singular matrix.
-	      char job = '1';
-	      Array<float> z (4 * nc);
-	      float *pz = z.fortran_vec ();
-	      Array<octave_idx_type> iz (nc);
-	      octave_idx_type *piz = iz.fortran_vec ();
-
-	      F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
-					 nc, tmp_data, nr, anorm, 
-					 rcon, pz, piz, info
-					 F77_CHAR_ARG_LEN (1)));
-	    }
-
-	  if (info != 0) 
-	    {
-	      info = -1;
-	      retval = FloatDET ();
-	    } 
-	  else 
-	    {
-              retval = FloatDET (1.0);
-              
-	      for (octave_idx_type i = 0; i < nc; i++) 
-		{
-                  float c = atmp(i,i);
-                  retval *= (ipvt(i) != (i+1)) ? -c : c;
+      int typ = mattype.type ();
+
+      if (typ == MatrixType::Lower || typ == MatrixType::Upper)
+        {
+          for (octave_idx_type i = 0; i < nc; i++) 
+            retval *= elem (i,i);
+        }
+      else if (typ == MatrixType::Hermitian)
+        {
+          FloatMatrix atmp = *this;
+          float *tmp_data = atmp.fortran_vec ();
+
+          info = 0;
+          float anorm = 0;
+          if (calc_cond) anorm = xnorm (*this, 1);
+
+
+          char job = 'L';
+          F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 
+                                     tmp_data, nr, info
+                                     F77_CHAR_ARG_LEN (1)));
+
+          if (info != 0) 
+            {
+              rcon = 0.0;
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+          else 
+            {
+              Array<float> z (3 * nc);
+              float *pz = z.fortran_vec ();
+              Array<octave_idx_type> iz (nc);
+              octave_idx_type *piz = iz.fortran_vec ();
+
+              F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                         nr, tmp_data, nr, anorm,
+                                         rcon, pz, piz, info
+                                         F77_CHAR_ARG_LEN (1)));
+
+              if (info != 0) 
+                rcon = 0.0;
+
+              for (octave_idx_type i = 0; i < nc; i++) 
+                retval *= elem (i,i);
+
+              retval = retval.square ();
+            }
+        }
+      else if (typ != MatrixType::Full)
+        (*current_liboctave_error_handler) ("det: invalid dense matrix type");
+
+      if (typ == MatrixType::Full)
+        {
+          Array<octave_idx_type> ipvt (nr);
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          FloatMatrix atmp = *this;
+          float *tmp_data = atmp.fortran_vec ();
+
+          info = 0;
+
+          // Calculate the norm of the matrix, for later use.
+          float anorm = 0;
+          if (calc_cond) anorm = xnorm (*this, 1);
+
+          F77_XFCN (sgetrf, SGETRF, (nr, nr, tmp_data, nr, pipvt, info));
+
+          // Throw-away extra info LAPACK gives so as to not change output.
+          rcon = 0.0;
+          if (info != 0) 
+            {
+              info = -1;
+              retval = FloatDET ();
+            } 
+          else 
+            {
+              if (calc_cond) 
+                {
+                  // Now calc the condition number for non-singular matrix.
+                  char job = '1';
+                  Array<float> z (4 * nc);
+                  float *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (nc);
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                             nc, tmp_data, nr, anorm, 
+                                             rcon, pz, piz, info
+                                             F77_CHAR_ARG_LEN (1)));
                 }
-	    }
-	}
+
+              if (info != 0) 
+                {
+                  info = -1;
+                  retval = FloatDET ();
+                } 
+              else 
+                {
+                  for (octave_idx_type i = 0; i < nc; i++) 
+                    {
+                      float c = atmp(i,i);
+                      retval *= (ipvt(i) != (i+1)) ? -c : c;
+                    }
+                }
+            }
+        }
     }
 
   return retval;
diff --git a/liboctave/fMatrix.h b/liboctave/fMatrix.h
--- a/liboctave/fMatrix.h
+++ b/liboctave/fMatrix.h
@@ -145,6 +145,8 @@
   FloatDET determinant (void) const;
   FloatDET determinant (octave_idx_type& info) const;
   FloatDET determinant (octave_idx_type& info, float& rcon, int calc_cond = 1) const;
+  FloatDET determinant (MatrixType &mattype, octave_idx_type& info, 
+                        float& rcon, int calc_cond = 1) const;
 
   float rcond (void) const;
   float rcond (MatrixType &mattype) const;