# HG changeset patch
# User Jaroslav Hajek
# Date 1211382848 -7200
# Node ID e0fdf5deced75945bf75268540bafa6529357630
# Parent 39c1026191e93605b94bad19be3bae85c80e10ff
improved matrix_type check
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,17 @@
+2008-05-21 Jaroslav Hajek
+
+ * 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
* Array.cc (Array Array::transpose () const): Modify for tiled
diff --git a/liboctave/MatrixType.cc b/liboctave/MatrixType.cc
--- a/liboctave/MatrixType.cc
+++ b/liboctave/MatrixType.cc
@@ -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
+MatrixType::matrix_type
+matrix_real_probe (const MArray2& 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,36 +71,89 @@
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)
+ typ = MatrixType::Upper;
+ else if (lower)
+ typ = MatrixType::Lower;
+ else if (hermitian)
+ typ = MatrixType::Hermitian;
+ else
+ typ = MatrixType::Full;
+ }
+ else
+ typ = MatrixType::Rectangular;
+
+ return typ;
+}
+
+template
+MatrixType::matrix_type
+matrix_complex_probe (const MArray2& 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;
+
+ // do the checks for lower/upper/hermitian all in one pass.
+ ColumnVector diag(ncols);
+
+ for (octave_idx_type j = 0;
+ j < ncols && upper; j++)
+ {
+ 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();
+ }
+
+ 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;
@@ -111,6 +166,16 @@
}
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)
@@ -118,61 +183,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);
}
@@ -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;
+ // first, check whether the diagonal is positive & extract it
+ ColumnVector diag (ncols);
- for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
- {
- octave_idx_type ri = a.ridx(i);
+ 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;
+ }
+ }
+ }
- 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 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;
- }
- }
+ // next, check symmetry and 2x2 positiveness
- 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;
+ // first, check whether the diagonal is positive & extract it
+ ColumnVector diag (ncols);
- for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
- {
- octave_idx_type ri = a.ridx(i);
+ 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;
+ }
+ }
+ }
- if (ri == j)
- {
- if (a.data(i) == std::abs(a.data(i)))
- diag_positive = true;
- else
- break;
- }
- else
- {
- bool found = false;
+ // 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");
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@
+2008-05-21 Jaroslav Hajek
+
+ * DLD-FUNCTIONS/matrix_type.cc: Fix tests relying on the
+ older more optimistic hermitian check.
+
2008-05-21 John W. Eaton
* pt-idx.h (tree_index_expression::tree_index_expression (int, int)):
diff --git a/src/DLD-FUNCTIONS/matrix_type.cc b/src/DLD-FUNCTIONS/matrix_type.cc
--- a/src/DLD-FUNCTIONS/matrix_type.cc
+++ b/src/DLD-FUNCTIONS/matrix_type.cc
@@ -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