--- a/liboctave/fMatrix.cc	Wed Nov 19 11:23:07 2008 +0100
+++ b/liboctave/fMatrix.cc	Wed Nov 19 15:26:39 2008 +0100
@@ -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;