Mercurial > octave-nkf

changeset 7798:42c42c640108

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improved matrix_type check
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 21 May 2008 17:14:08 +0200
parents f42c6f8d6d8e
children 199181592240
files liboctave/ChangeLog liboctave/MatrixType.cc src/ChangeLog src/DLD-FUNCTIONS/matrix_type.cc
diffstat 4 files changed, 187 insertions(+), 269 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog	Wed May 14 23:28:41 2008 +0200
+++ b/liboctave/ChangeLog	Wed May 21 17:14:08 2008 +0200
@@ -66,6 +66,20 @@
 	FloatHESS, FloatComplexHESS, FloatQR, FloatComplexQR, QRP,
 	ComplexQRP, FloatQRP, FloatComplexQRP):	Declare classes.
 	
+2008-05-21  Jaroslav Hajek <highegg@gmail.com>
+
+	* MatrixType.cc (matrix_real_probe, matrix_complex_probe):
+	New template functions.
+	(MatrixType::MatrixType (const Matrix&),
+	MatrixType::MatrixType (const FloatMatrix&)):
+	just call matrix_real_probe.
+	(MatrixType::MatrixType (const ComplexMatrix&),
+	MatrixType::MatrixType (const FloatComplexMatrix&)):
+	just call matrix_complex_probe.
+
+	* MatrixType.cc (MatrixType::MatrixType (matrix_type, bool)):
+	add missing test for Unknown.
+
 2008-05-20  David Bateman  <dbateman@free.fr>
 
 	* Array.cc (Array<T> Array<T>::transpose () const): Modify for tiled
--- a/liboctave/MatrixType.cc	Wed May 14 23:28:41 2008 +0200
+++ b/liboctave/MatrixType.cc	Wed May 21 17:14:08 2008 +0200
@@ -55,13 +55,15 @@
     }
 }
 
-MatrixType::MatrixType (const Matrix &a)
-  : typ (MatrixType::Unknown),
-    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
-    dense (false), full (true), nperm (0), perm (0)
+template<class T> 
+MatrixType::matrix_type 
+matrix_real_probe (const MArray2<T>& a)
 {
+  MatrixType::matrix_type typ;
   octave_idx_type nrows = a.rows ();
   octave_idx_type ncols = a.cols ();
+
+  const T zero = 0;
  
   if (ncols == nrows)
     {
@@ -69,35 +71,30 @@
       bool lower = true;
       bool hermitian = true;
 
-      for (octave_idx_type j = 0; j < ncols; j++)
+      // do the checks for lower/upper/hermitian all in one pass.
+      ColumnVector diag(ncols);
+
+      for (octave_idx_type j = 0; 
+           j < ncols && upper; j++)
 	{
-	  if (j < nrows)
-	    {
-	      if (a.elem (j,j) == 0.)
-		{
-		  upper = false;
-		  lower = false;
-		  hermitian = false;
-		  break;
-		}
-	      if (a.elem (j,j) < 0.)
-		hermitian = false;
-	    }      
-	  for (octave_idx_type i = 0; i < j; i++)
-	    if (lower && a.elem (i,j) != 0.)
-	      {
-		lower = false;
-		break;
-	      }
-	  for (octave_idx_type i = j+1; i < nrows; i++)
-	    {
-	      if (hermitian && a.elem (i, j) != a.elem (j, i))
-		hermitian = false;
-	      if (upper && a.elem (i,j) != 0)
-		upper = false;
-	    }
-	  if (!upper && !lower && !hermitian)
-	    break;
+          T d = a.elem (j,j);
+          upper = upper && (d != zero);
+          lower = lower && (d != zero);
+          hermitian = hermitian && (d > zero);
+          diag(j) = d;
+        }
+
+      for (octave_idx_type j = 0; 
+           j < ncols && (upper || lower || hermitian); j++)
+	{
+          for (octave_idx_type i = 0; i < j; i++)
+            {
+              double aij = a.elem (i,j), aji = a.elem (j,i);
+              lower = lower && (aij == zero);
+              upper = upper && (aji == zero);
+              hermitian = hermitian && (aij == aji 
+                                        && aij*aij < diag(i)*diag(j));
+            }
 	}
 
       if (upper)
@@ -106,62 +103,58 @@
 	typ = MatrixType::Lower;
       else if (hermitian)
 	typ = MatrixType::Hermitian;
-      else if (ncols == nrows)
+      else 
 	typ = MatrixType::Full;
     }
   else
     typ = MatrixType::Rectangular;
+
+  return typ;
 }
 
-MatrixType::MatrixType (const ComplexMatrix &a)
-  : typ (MatrixType::Unknown),
-    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
-    dense (false), full (true), nperm (0), perm (0)
+template<class T> 
+MatrixType::matrix_type 
+matrix_complex_probe (const MArray2<T>& a)
 {
+  MatrixType::matrix_type typ;
   octave_idx_type nrows = a.rows ();
   octave_idx_type ncols = a.cols ();
 
+  const typename T::value_type zero = 0;
+
   if (ncols == nrows)
     {
       bool upper = true;
       bool lower = true;
       bool hermitian = true;
 
-      for (octave_idx_type j = 0; j < ncols; j++)
+      // do the checks for lower/upper/hermitian all in one pass.
+      ColumnVector diag(ncols);
+
+      for (octave_idx_type j = 0; 
+           j < ncols && upper; j++)
 	{
-	  if (j < ncols)
-	    {
-	      if (imag(a.elem (j,j)) == 0. && 
-		  real(a.elem (j,j)) == 0.)
-		{
-		  upper = false;
-		  lower = false;
-		  hermitian = false;
-		  break;
-		}
+          T d = a.elem (j,j);
+          upper = upper && (d != zero);
+          lower = lower && (d != zero);
+          hermitian = hermitian && (d.real() > zero && d.imag() == zero);
+          diag (j) = d.real();
+        }
 
-	      if (imag(a.elem (j,j)) != 0. || 
-		  real(a.elem (j,j)) < 0.)
-		    hermitian = false;
-	    }
-	  for (octave_idx_type i = 0; i < j; i++)
-	    if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0))
-	      {
-		lower = false;
-		break;
-	      }
-	  for (octave_idx_type i = j+1; i < nrows; i++)
-	    {
-	      if (hermitian && a.elem (i, j) != conj(a.elem (j, i)))
-		hermitian = false;
-	      if (upper && (real(a.elem (i,j)) != 0 || 
-			    imag(a.elem (i,j)) != 0))
-		upper = false;
-	    }
-	  if (!upper && !lower && !hermitian)
-	    break;
+      for (octave_idx_type j = 0; 
+           j < ncols && (upper || lower || hermitian); j++)
+	{
+          for (octave_idx_type i = 0; i < j; i++)
+            {
+              T aij = a.elem (i,j), aji = a.elem (j,i);
+              lower = lower && (aij == zero);
+              upper = upper && (aji == zero);
+              hermitian = hermitian && (aij == std::conj (aji)
+                                        && std::norm (aij) < diag(i)*diag(j));
+            }
 	}
 
+
       if (upper)
 	typ = MatrixType::Upper;
       else if (lower)
@@ -173,6 +166,24 @@
     }
   else
     typ = MatrixType::Rectangular;
+
+  return typ;
+}
+
+MatrixType::MatrixType (const Matrix &a)
+  : typ (MatrixType::Unknown),
+    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
+    dense (false), full (true), nperm (0), perm (0)
+{
+  typ = matrix_real_probe (a);
+}
+
+MatrixType::MatrixType (const ComplexMatrix &a)
+  : typ (MatrixType::Unknown),
+    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
+    dense (false), full (true), nperm (0), perm (0)
+{
+  typ = matrix_complex_probe (a);
 }
 
 
@@ -181,57 +192,7 @@
     sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
     dense (false), full (true), nperm (0), perm (0)
 {
-  octave_idx_type nrows = a.rows ();
-  octave_idx_type ncols = a.cols ();
- 
-  if (ncols == nrows)
-    {
-      bool upper = true;
-      bool lower = true;
-      bool hermitian = true;
-
-      for (octave_idx_type j = 0; j < ncols; j++)
-	{
-	  if (j < nrows)
-	    {
-	      if (a.elem (j,j) == 0.)
-		{
-		  upper = false;
-		  lower = false;
-		  hermitian = false;
-		  break;
-		}
-	      if (a.elem (j,j) < 0.)
-		hermitian = false;
-	    }      
-	  for (octave_idx_type i = 0; i < j; i++)
-	    if (lower && a.elem (i,j) != 0.)
-	      {
-		lower = false;
-		break;
-	      }
-	  for (octave_idx_type i = j+1; i < nrows; i++)
-	    {
-	      if (hermitian && a.elem (i, j) != a.elem (j, i))
-		hermitian = false;
-	      if (upper && a.elem (i,j) != 0)
-		upper = false;
-	    }
-	  if (!upper && !lower && !hermitian)
-	    break;
-	}
-
-      if (upper)
-	typ = MatrixType::Upper;
-      else if (lower)
-	typ = MatrixType::Lower;
-      else if (hermitian)
-	typ = MatrixType::Hermitian;
-      else if (ncols == nrows)
-	typ = MatrixType::Full;
-    }
-  else
-    typ = MatrixType::Rectangular;
+  typ = matrix_real_probe (a);
 }
 
 MatrixType::MatrixType (const FloatComplexMatrix &a)
@@ -239,61 +200,7 @@
     sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
     dense (false), full (true), nperm (0), perm (0)
 {
-  octave_idx_type nrows = a.rows ();
-  octave_idx_type ncols = a.cols ();
-
-  if (ncols == nrows)
-    {
-      bool upper = true;
-      bool lower = true;
-      bool hermitian = true;
-
-      for (octave_idx_type j = 0; j < ncols; j++)
-	{
-	  if (j < ncols)
-	    {
-	      if (imag(a.elem (j,j)) == 0. && 
-		  real(a.elem (j,j)) == 0.)
-		{
-		  upper = false;
-		  lower = false;
-		  hermitian = false;
-		  break;
-		}
-
-	      if (imag(a.elem (j,j)) != 0. || 
-		  real(a.elem (j,j)) < 0.)
-		    hermitian = false;
-	    }
-	  for (octave_idx_type i = 0; i < j; i++)
-	    if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0))
-	      {
-		lower = false;
-		break;
-	      }
-	  for (octave_idx_type i = j+1; i < nrows; i++)
-	    {
-	      if (hermitian && a.elem (i, j) != conj(a.elem (j, i)))
-		hermitian = false;
-	      if (upper && (real(a.elem (i,j)) != 0 || 
-			    imag(a.elem (i,j)) != 0))
-		upper = false;
-	    }
-	  if (!upper && !lower && !hermitian)
-	    break;
-	}
-
-      if (upper)
-	typ = MatrixType::Upper;
-      else if (lower)
-	typ = MatrixType::Lower;
-      else if (hermitian)
-	typ = MatrixType::Hermitian;
-      else if (ncols == nrows)
-	typ = MatrixType::Full;
-    }
-  else
-    typ = MatrixType::Rectangular;
+  typ = matrix_complex_probe (a);
 }
 
 MatrixType::MatrixType (const SparseMatrix &a)
@@ -560,54 +467,49 @@
 			      typ == MatrixType::Tridiagonal || 
 			      typ == MatrixType::Banded))
 	{
-	  // Check for symmetry, with positive real diagonal, which
-	  // has a very good chance of being symmetric positive
-	  // definite..
 	  bool is_herm = true;
 
-	  for (octave_idx_type j = 0; j < ncols; j++)
-	    {
-	      bool diag_positive = false;
-
-	      for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
-		{
-		  octave_idx_type ri = a.ridx(i);
+          // first, check whether the diagonal is positive & extract it
+          ColumnVector diag (ncols);
 
-		  if (ri == j)
-		    {
-		      if (a.data(i) == std::abs(a.data(i)))
-			diag_positive = true;
-		      else
-			break;
-		    }
-		  else
-		    {
-		      bool found = false;
+	  for (octave_idx_type j = 0; is_herm && j < ncols; j++)
+            {
+              is_herm = false;
+              for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
+                {
+                  if (a.ridx(i) == j)
+                    {
+                      double d = a.data(i);
+                      is_herm = d > 0.;
+                      diag(j) = d;
+                      break;
+                    }
+                }
+            }
+
+
+          // next, check symmetry and 2x2 positiveness
 
-		      for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++)
-			{
-			  if (a.ridx(k) == j)
-			    {
-			      if (a.data(i) == a.data(k))
-				found = true;
-			      break;
-			    }
-			}
-
-		      if (! found)
-			{
-			  is_herm = false;
-			  break;
-			}
-		    }
-		}
-
-	      if (! diag_positive || ! is_herm)
-		{
-		  is_herm = false;
-		  break;
-		} 
-	    }
+          for (octave_idx_type j = 0; is_herm && j < ncols; j++)
+            for (octave_idx_type i = a.cidx(j); is_herm && i < a.cidx(j+1); i++)
+              {
+                octave_idx_type k = a.ridx(i);
+                is_herm = k == j;
+                if (is_herm) 
+                  continue;
+                double d = a.data(i);
+                if (d*d < diag(j)*diag(k))
+                  {
+                    for (octave_idx_type l = a.cidx(k); l < a.cidx(k+1); l++)
+                      {
+                        if (a.ridx(l) == j)
+                          {
+                            is_herm = a.data(l) == d;
+                            break;
+                          }
+                      }
+                  }
+              }
 
 	  if (is_herm)
 	    {
@@ -886,54 +788,50 @@
 			      typ == MatrixType::Tridiagonal || 
 			      typ == MatrixType::Banded))
 	{
-	  // Check for symmetry, with positive real diagonal, which
-	  // has a very good chance of being symmetric positive
-	  // definite..
 	  bool is_herm = true;
 
-	  for (octave_idx_type j = 0; j < ncols; j++)
-	    {
-	      bool diag_positive = false;
-
-	      for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
-		{
-		  octave_idx_type ri = a.ridx(i);
+          // first, check whether the diagonal is positive & extract it
+          ColumnVector diag (ncols);
 
-		  if (ri == j)
-		    {
-		      if (a.data(i) == std::abs(a.data(i)))
-			diag_positive = true;
-		      else
-			break;
-		    }
-		  else
-		    {
-		      bool found = false;
+	  for (octave_idx_type j = 0; is_herm && j < ncols; j++)
+            {
+              is_herm = false;
+              for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
+                {
+                  if (a.ridx(i) == j)
+                    {
+                      Complex d = a.data(i);
+                      is_herm = d.real() > 0. && d.imag() == 0.;
+                      diag(j) = d.real();
+                      break;
+                    }
+                }
+            }
+
+          // next, check symmetry and 2x2 positiveness
 
-		      for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++)
-			{
-			  if (a.ridx(k) == j)
-			    {
-			      if (a.data(i) == conj(a.data(k)))
-				found = true;
-			      break;
-			    }
-			}
+          for (octave_idx_type j = 0; is_herm && j < ncols; j++)
+            for (octave_idx_type i = a.cidx(j); is_herm && i < a.cidx(j+1); i++)
+              {
+                octave_idx_type k = a.ridx(i);
+                is_herm = k == j;
+                if (is_herm) 
+                  continue;
+                Complex d = a.data(i);
+                if (std::norm (d) < diag(j)*diag(k))
+                  {
+                    d = std::conj (d);
+                    for (octave_idx_type l = a.cidx(k); l < a.cidx(k+1); l++)
+                      {
+                        if (a.ridx(l) == j)
+                          {
+                            is_herm = a.data(l) == d;
+                            break;
+                          }
+                      }
+                  }
+              }
 
-		      if (! found)
-			{
-			  is_herm = false;
-			  break;
-			}
-		    }
-		}
-
-	      if (! diag_positive || ! is_herm)
-		{
-		  is_herm = false;
-		  break;
-		} 
-	    }
 
 	  if (is_herm)
 	    {
@@ -953,10 +851,11 @@
     bandden (0), upper_band (0), lower_band (0),
     dense (false), full (_full), nperm (0), perm (0)
 {
-  if (t == MatrixType::Full || t == MatrixType::Diagonal ||
-      t == MatrixType::Permuted_Diagonal || t == MatrixType::Upper ||
-      t == MatrixType::Lower || t == MatrixType::Tridiagonal ||
-      t == MatrixType::Tridiagonal_Hermitian || t == MatrixType::Rectangular)
+  if (t == MatrixType::Unknown || t == MatrixType::Full 
+      || t == MatrixType::Diagonal || t == MatrixType::Permuted_Diagonal 
+      || t == MatrixType::Upper || t == MatrixType::Lower 
+      || t == MatrixType::Tridiagonal || t == MatrixType::Tridiagonal_Hermitian
+      || t == MatrixType::Rectangular)
     typ = t;
   else
     (*current_liboctave_warning_handler) ("Invalid matrix type");
--- a/src/ChangeLog	Wed May 14 23:28:41 2008 +0200
+++ b/src/ChangeLog	Wed May 21 17:14:08 2008 +0200
@@ -12,6 +12,11 @@
 	(Feps): Modify behavior for a single numerical argument to give
 	difference to next largest value in the class of the type passed.
 
+2008-05-21  Jaroslav Hajek <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/matrix_type.cc: Fix tests relying on the
+	older more optimistic hermitian check.
+
 2008-05-21  John W. Eaton  <jwe@octave.org>
 
 	* pt-idx.h (tree_index_expression::tree_index_expression (int, int)): 
--- a/src/DLD-FUNCTIONS/matrix_type.cc	Wed May 14 23:28:41 2008 +0200
+++ b/src/DLD-FUNCTIONS/matrix_type.cc	Wed May 21 17:14:08 2008 +0200
@@ -512,9 +512,9 @@
 %!test
 %! bnd=spparms("bandden");
 %! spparms("bandden",0.5);
-%! a = spdiags(randn(10,3),[-1,0,1],10,10);
+%! a = spdiags(rand(10,3)-0.5,[-1,0,1],10,10);
 %! assert(matrix_type(a),"Tridiagonal");
-%! assert(matrix_type(abs(a')+abs(a)),"Tridiagonal Positive Definite");
+%! assert(matrix_type(a'+a+2*speye(10)),"Tridiagonal Positive Definite");
 %! spparms("bandden",bnd);
 %!test
 %! bnd=spparms("bandden");
@@ -551,14 +551,14 @@
 %! bnd=spparms("bandden");
 %! spparms("bandden",0.5);
 %! assert(matrix_type(spdiags(1i*randn(10,3),[-1,0,1],10,10)),"Tridiagonal");
-%! a = 1i*randn(9,1);a=[[a;0],ones(10,1),[0;-a]];
+%! a = 1i*(rand(9,1)-0.5);a=[[a;0],ones(10,1),[0;-a]];
 %! assert(matrix_type(spdiags(a,[-1,0,1],10,10)),"Tridiagonal Positive Definite");
 %! spparms("bandden",bnd);
 %!test
 %! bnd=spparms("bandden");
 %! spparms("bandden",0.5);
 %! assert(matrix_type(spdiags(1i*randn(10,4),[-2:1],10,10)),"Banded");
-%! a = 1i*randn(9,2);a=[[a;[0,0]],ones(10,1),[[0;-a(:,2)],[0;0;-a(1:8,1)]]];
+%! a = 1i*(rand(9,2)-0.5);a=[[a;[0,0]],ones(10,1),[[0;-a(:,2)],[0;0;-a(1:8,1)]]];
 %! assert(matrix_type(spdiags(a,[-2:2],10,10)),"Banded Positive Definite");
 %! spparms("bandden",bnd);
 %!test