Mercurial > hg > octave-lojdl
changeset 7798:42c42c640108
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 +++ b/liboctave/ChangeLog @@ -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 +++ 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<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 +++ b/src/ChangeLog @@ -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 +++ 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