Mercurial > hg > octave-lojdl
changeset 6207:3c92b8d892dd
[project @ 2006-12-06 20:19:14 by dbateman]
author | dbateman |
---|---|
date | Wed, 06 Dec 2006 20:19:16 +0000 |
parents | cb8c62c78b42 |
children | 323be5eeed1f |
files | liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/CSparse.cc liboctave/ChangeLog liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/dSparse.cc src/ChangeLog src/DLD-FUNCTIONS/inv.cc src/xpow.cc |
diffstat | 10 files changed, 358 insertions(+), 30 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -41,6 +41,7 @@ #include "CmplxDET.h" #include "CmplxSCHUR.h" #include "CmplxSVD.h" +#include "CmplxCHOL.h" #include "f77-fcn.h" #include "lo-error.h" #include "lo-ieee.h" @@ -142,6 +143,13 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T + F77_FUNC (ztrtri, ZTRTRI) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const Complex*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T F77_FUNC (ztrcon, ZTRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, const Complex*, const octave_idx_type&, double&, @@ -959,19 +967,94 @@ { octave_idx_type info; double rcond; - return inverse (info, rcond, 0, 0); + MatrixType mattype (*this); + return inverse (mattype, info, rcond, 0, 0); +} + +ComplexMatrix +ComplexMatrix::inverse (MatrixType &mattype) const +{ + octave_idx_type info; + double rcond; + return inverse (mattype, info, rcond, 0, 0); +} + +ComplexMatrix +ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const +{ + double rcond; + return inverse (mattype, info, rcond, 0, 0); } ComplexMatrix -ComplexMatrix::inverse (octave_idx_type& info) const +ComplexMatrix::tinverse (MatrixType &mattype, octave_idx_type& info, + double& rcond, int force, int calc_cond) const { - double rcond; - return inverse (info, rcond, 0, 0); + ComplexMatrix retval; + + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + + if (nr != nc || nr == 0 || nc == 0) + (*current_liboctave_error_handler) ("inverse requires square matrix"); + else + { + int typ = mattype.type (); + char uplo = (typ == MatrixType::Lower ? 'L' : 'U'); + char udiag = 'N'; + retval = *this; + Complex *tmp_data = retval.fortran_vec (); + + F77_XFCN (ztrtri, ZTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in ztrtri"); + else + { + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + info = -1; + else if (calc_cond) + { + octave_idx_type ztrcon_info = 0; + char job = '1'; + + OCTAVE_LOCAL_BUFFER (Complex, cwork, 2 * nr); + OCTAVE_LOCAL_BUFFER (double, rwork, nr); + + F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, rcond, + cwork, rwork, ztrcon_info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in ztrcon"); + + if (ztrcon_info != 0) + info = -1; + } + } + + if (info == -1 && ! force) + retval = *this; // Restore matrix contents. + } + + return retval; } ComplexMatrix -ComplexMatrix::inverse (octave_idx_type& info, double& rcond, int force, - int calc_cond) const +ComplexMatrix::finverse (MatrixType &mattype, octave_idx_type& info, + double& rcond, int force, int calc_cond) const { ComplexMatrix retval; @@ -1062,12 +1145,45 @@ info = -1; } } + + if (info != 0) + mattype.mark_as_rectangular(); } return retval; } ComplexMatrix +ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info, + double& rcond, int force, int calc_cond) const +{ + int typ = mattype.type (false); + ComplexMatrix ret; + + if (typ == MatrixType::Unknown) + typ = mattype.type (*this); + + if (typ == MatrixType::Upper || typ == MatrixType::Lower) + ret = tinverse (mattype, info, rcond, force, calc_cond); + else if (typ != MatrixType::Rectangular) + { + if (mattype.is_hermitian ()) + { + ComplexCHOL chol (*this, info); + if (info == 0) + ret = chol.inverse (); + else + mattype.mark_as_unsymmetric (); + } + + if (!mattype.is_hermitian ()) + ret = finverse(mattype, info, rcond, force, calc_cond); + } + + return ret; +} + +ComplexMatrix ComplexMatrix::pseudo_inverse (double tol) const { ComplexMatrix retval;
--- a/liboctave/CMatrix.h +++ b/liboctave/CMatrix.h @@ -135,9 +135,19 @@ ComplexColumnVector column (octave_idx_type i) const; +private: + ComplexMatrix tinverse (MatrixType &mattype, octave_idx_type& info, + double& rcond, int force, int calc_cond) const; + + ComplexMatrix finverse (MatrixType &mattype, octave_idx_type& info, + double& rcond, int force, int calc_cond) const; + +public: ComplexMatrix inverse (void) const; - ComplexMatrix inverse (octave_idx_type& info) const; - ComplexMatrix inverse (octave_idx_type& info, double& rcond, int force = 0, + ComplexMatrix inverse (MatrixType &mattype) const; + ComplexMatrix inverse (MatrixType &mattype, octave_idx_type& info) const; + ComplexMatrix inverse (MatrixType &mattype, octave_idx_type& info, + double& rcond, int force = 0, int calc_cond = 1) const; ComplexMatrix pseudo_inverse (double tol = 0.0) const;
--- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -180,12 +180,33 @@ octave_idx_type nr = rows (); octave_idx_type nc = cols (); - if (is_square () && nr > 0) + if (nr == nc && nr > 0) { - for (octave_idx_type i = 0; i < nr; i++) - for (octave_idx_type j = i; j < nc; j++) - if (elem (i, j) != conj (elem (j, i))) - return false; + for (octave_idx_type j = 0; j < nc; j++) + { + for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) + { + octave_idx_type ri = ridx(i); + + if (ri != j) + { + bool found = false; + + for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++) + { + if (ridx(k) == j) + { + if (data(i) == conj(data(k))) + found = true; + break; + } + } + + if (! found) + return false; + } + } + } return true; }
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,17 @@ +2006-12-06 David Bateman <dbateman@free.fr> + + * dSparse.cc (SparseMatrix::is_symmetric(void) const): Faster + implementation. + * CSparse.cc (SparseComplexMatrix::is_symmetric(void) const): ditto. + + * dMatrrix.cc (finverse): Old inverse method renamed inverse. + (tinverse): New method for triangular matrices. + (inverse): New function with matrix type probing. + * dMatrix.h (finverse, tinverse, inverse): New and modified + declarations. + * CMatrix.cc: ditto. + * CMatrix.h: ditto. + 2006-12-06 John W. Eaton <jwe@octave.org> * strptime.c (day_of_the_week): Use code from current glibc sources.
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -37,6 +37,7 @@ #include "dbleDET.h" #include "dbleSCHUR.h" #include "dbleSVD.h" +#include "dbleCHOL.h" #include "f77-fcn.h" #include "lo-error.h" #include "lo-ieee.h" @@ -138,6 +139,12 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T + F77_FUNC (dtrtri, DTRTRI) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const double*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + F77_RET_T F77_FUNC (dtrcon, DTRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, const double*, const octave_idx_type&, double&, @@ -628,18 +635,95 @@ { octave_idx_type info; double rcond; - return inverse (info, rcond, 0, 0); + MatrixType mattype (*this); + return inverse (mattype, info, rcond, 0, 0); +} + +Matrix +Matrix::inverse (MatrixType& mattype) const +{ + octave_idx_type info; + double rcond; + return inverse (mattype, info, rcond, 0, 0); +} + +Matrix +Matrix::inverse (MatrixType &mattype, octave_idx_type& info) const +{ + double rcond; + return inverse (mattype, info, rcond, 0, 0); } Matrix -Matrix::inverse (octave_idx_type& info) const +Matrix::tinverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + int force, int calc_cond) const { - double rcond; - return inverse (info, rcond, 0, 0); + Matrix retval; + + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + + if (nr != nc || nr == 0 || nc == 0) + (*current_liboctave_error_handler) ("inverse requires square matrix"); + else + { + int typ = mattype.type (); + char uplo = (typ == MatrixType::Lower ? 'L' : 'U'); + char udiag = 'N'; + retval = *this; + double *tmp_data = retval.fortran_vec (); + + F77_XFCN (dtrtri, DTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dtrtri"); + else + { + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + info = -1; + else if (calc_cond) + { + octave_idx_type dtrcon_info = 0; + char job = '1'; + + OCTAVE_LOCAL_BUFFER (double, work, 3 * nr); + OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr); + + F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, rcond, + work, iwork, dtrcon_info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dtrcon"); + + if (dtrcon_info != 0) + info = -1; + } + } + + if (info == -1 && ! force) + retval = *this; // Restore matrix contents. + } + + return retval; } + Matrix -Matrix::inverse (octave_idx_type& info, double& rcond, int force, int calc_cond) const +Matrix::finverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + int force, int calc_cond) const { Matrix retval; @@ -730,12 +814,45 @@ info = -1; } } + + if (info != 0) + mattype.mark_as_rectangular(); } return retval; } Matrix +Matrix::inverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + int force, int calc_cond) const +{ + int typ = mattype.type (false); + Matrix ret; + + if (typ == MatrixType::Unknown) + typ = mattype.type (*this); + + if (typ == MatrixType::Upper || typ == MatrixType::Lower) + ret = tinverse (mattype, info, rcond, force, calc_cond); + else if (typ != MatrixType::Rectangular) + { + if (mattype.is_hermitian ()) + { + CHOL chol (*this, info); + if (info == 0) + ret = chol.inverse (); + else + mattype.mark_as_unsymmetric (); + } + + if (!mattype.is_hermitian ()) + ret = finverse(mattype, info, rcond, force, calc_cond); + } + + return ret; +} + +Matrix Matrix::pseudo_inverse (double tol) const { SVD result (*this, SVD::economy);
--- a/liboctave/dMatrix.h +++ b/liboctave/dMatrix.h @@ -107,10 +107,19 @@ ColumnVector column (octave_idx_type i) const; +private: + Matrix tinverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + int force, int calc_cond) const; + + Matrix finverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + int force, int calc_cond) const; + +public: Matrix inverse (void) const; - Matrix inverse (octave_idx_type& info) const; - Matrix inverse (octave_idx_type& info, double& rcond, int force = 0, - int calc_cond = 1) const; + Matrix inverse (MatrixType &mattype) const; + Matrix inverse (MatrixType &mattype, octave_idx_type& info) const; + Matrix inverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + int force = 0, int calc_cond = 1) const; Matrix pseudo_inverse (double tol = 0.0) const;
--- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -170,12 +170,36 @@ bool SparseMatrix::is_symmetric (void) const { - if (is_square () && rows () > 0) + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + + if (nr == nc && nr > 0) { - for (octave_idx_type i = 0; i < rows (); i++) - for (octave_idx_type j = i+1; j < cols (); j++) - if (elem (i, j) != elem (j, i)) - return false; + for (octave_idx_type j = 0; j < nc; j++) + { + for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) + { + octave_idx_type ri = ridx(i); + + if (ri != j) + { + bool found = false; + + for (octave_idx_type k = cidx(ri); k < cidx(ri+1); k++) + { + if (ridx(k) == j) + { + if (data(i) == data(k)) + found = true; + break; + } + } + + if (! found) + return false; + } + } + } return true; }
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,10 @@ +2006-12-04 David Bateman <dbateman@free.fr> + + * xpow.cc (xpow (const Matrix&, double)): Add matrix type probing + to matrix inverse. + (xpow (const ComplexMatrix&, double)): ditto. + * DLD-FUNCTIONS/inv.cc (Finv): Add matrix type probing. + 2006-12-06 John W. Eaton <jwe@octave.org> 2006-12-06 Michael Goffioul <michael.goffioul@swing.be>
--- a/src/DLD-FUNCTIONS/inv.cc +++ b/src/DLD-FUNCTIONS/inv.cc @@ -74,10 +74,14 @@ if (! error_state) { + MatrixType mattyp = args(0).matrix_type (); + octave_idx_type info; double rcond = 0.0; - Matrix result = m.inverse (info, rcond, 1); + Matrix result = m.inverse (mattyp, info, rcond, 1); + + args(0).matrix_type (mattyp); if (nargout > 1) retval(1) = rcond; @@ -97,10 +101,14 @@ if (! error_state) { + MatrixType mattyp = args(0).matrix_type (); + octave_idx_type info; double rcond = 0.0; - ComplexMatrix result = m.inverse (info, rcond, 1); + ComplexMatrix result = m.inverse (mattyp, info, rcond, 1); + + args(0).matrix_type (mattyp); if (nargout > 1) retval(1) = rcond;
--- a/src/xpow.cc +++ b/src/xpow.cc @@ -192,8 +192,9 @@ octave_idx_type info; double rcond = 0.0; + MatrixType mattype (a); - atmp = a.inverse (info, rcond, 1); + atmp = a.inverse (mattype, info, rcond, 1); if (info == -1) warning ("inverse: matrix singular to machine\ @@ -388,8 +389,9 @@ octave_idx_type info; double rcond = 0.0; + MatrixType mattype (a); - atmp = a.inverse (info, rcond, 1); + atmp = a.inverse (mattype, info, rcond, 1); if (info == -1) warning ("inverse: matrix singular to machine\