Mercurial > octave-nkf
changeset 8336:9813c07ca946
make det take advantage of matrix type
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Wed, 19 Nov 2008 15:26:39 +0100 |
| parents | 64cf956a109c |
| children | e02242c54c49 |
| files | liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/ChangeLog liboctave/DET.h liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/fCMatrix.cc liboctave/fCMatrix.h liboctave/fMatrix.cc liboctave/fMatrix.h |
| diffstat | 10 files changed, 496 insertions(+), 239 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/CMatrix.cc Wed Nov 19 15:26:39 2008 +0100 @@ -56,6 +56,7 @@ #include "mx-cm-s.h" #include "mx-inlines.cc" #include "oct-cmplx.h" +#include "oct-norm.h" #if defined (HAVE_FFTW3) #include "oct-fftw.h" @@ -1574,74 +1575,132 @@ ComplexDET ComplexMatrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const { - ComplexDET retval; + MatrixType mattype (*this); + return determinant (mattype, info, rcon, calc_cond); +} + +ComplexDET +ComplexMatrix::determinant (MatrixType& mattype, + octave_idx_type& info, double& rcon, int calc_cond) const +{ + ComplexDET retval (1.0); octave_idx_type nr = rows (); octave_idx_type nc = cols (); - if (nr == 0 || nc == 0) - { - retval = ComplexDET (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 (); - - ComplexMatrix atmp = *this; - Complex *tmp_data = atmp.fortran_vec (); - - info = 0; - - // Calculate the norm of the matrix, for later use. - double anorm = 0; - if (calc_cond) - anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); - - F77_XFCN (zgetrf, ZGETRF, (nr, nc, 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 = ComplexDET (); - } - else - { - if (calc_cond) - { - // Now calc the condition number for non-singular matrix. - char job = '1'; - Array<Complex> z (2*nr); - Complex *pz = z.fortran_vec (); - Array<double> rz (2*nr); - double *prz = rz.fortran_vec (); - - F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcon, pz, prz, info - F77_CHAR_ARG_LEN (1))); - } - - if (info != 0) - { - info = -1; - retval = ComplexDET (); - } - else - { - retval = ComplexDET (1.0); - - for (octave_idx_type i = 0; i < nc; i++) - { - Complex 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) + { + ComplexMatrix atmp = *this; + Complex *tmp_data = atmp.fortran_vec (); + + info = 0; + double anorm = 0; + if (calc_cond) anorm = xnorm (*this, 1); + + + char job = 'L'; + F77_XFCN (zpotrf, ZPOTRF, (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<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (nc); + double *prz = rz.fortran_vec (); + + F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcon, pz, prz, 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 (); + + ComplexMatrix atmp = *this; + Complex *tmp_data = atmp.fortran_vec (); + + info = 0; + + // Calculate the norm of the matrix, for later use. + double anorm = 0; + if (calc_cond) anorm = xnorm (*this, 1); + + F77_XFCN (zgetrf, ZGETRF, (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 = ComplexDET (); + } + else + { + if (calc_cond) + { + // Now calc the condition number for non-singular matrix. + char job = '1'; + Array<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (2 * nc); + double *prz = rz.fortran_vec (); + + F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, prz, info + F77_CHAR_ARG_LEN (1))); } - } - } + + if (info != 0) + { + info = -1; + retval = ComplexDET (); + } + else + { + for (octave_idx_type i = 0; i < nc; i++) + { + Complex c = atmp(i,i); + retval *= (ipvt(i) != (i+1)) ? -c : c; + } + } + } + } } - + return retval; }
--- a/liboctave/CMatrix.h Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/CMatrix.h Wed Nov 19 15:26:39 2008 +0100 @@ -176,6 +176,8 @@ ComplexDET determinant (void) const; ComplexDET determinant (octave_idx_type& info) const; ComplexDET determinant (octave_idx_type& info, double& rcon, int calc_cond = 1) const; + ComplexDET determinant (MatrixType &mattype, octave_idx_type& info, + double& rcon, int calc_cond = 1) const; double rcond (void) const; double rcond (MatrixType &mattype) const;
--- a/liboctave/ChangeLog Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/ChangeLog Wed Nov 19 15:26:39 2008 +0100 @@ -1,3 +1,13 @@ +2008-11-19 Jaroslav Hajek <highegg@gmail.com> + + * DET.h (base_det<T>::square): New member function. + * dMatrix.cc (Matrix::determinant), + fMatrix.cc (FloatMatrix::determinant), + CMatrix.cc (ComplexMatrix::determinant), + fCMatrix.cc (FloatComplexMatrix::determinant): + Allow taking MatrixType argument. + * dMatrix.h, fMatrix.h, CMatrix.h, fCMatrix.h: Update decls. + 2008-11-19 Jaroslav Hajek <highegg@gmail.com> * DET.h: New source.
--- a/liboctave/DET.h Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/DET.h Wed Nov 19 15:26:39 2008 +0100 @@ -34,7 +34,7 @@ { public: - base_det (T c = 0, int e = 0) + base_det (T c = 1, int e = 0) { c2 = xlog2 (c, e2); e2 += e; @@ -65,6 +65,8 @@ T value () const { return c2 * static_cast<T> (std::ldexp (1.0, e2)); } operator T () const { return value (); } + base_det square () const { return base_det (c2*c2, e2+e2); } + void operator *= (T t) { int e;
--- a/liboctave/dMatrix.cc Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/dMatrix.cc Wed Nov 19 15:26:39 2008 +0100 @@ -51,6 +51,7 @@ #include "mx-dm-m.h" #include "mx-inlines.cc" #include "oct-cmplx.h" +#include "oct-norm.h" #include "quit.h" #if defined (HAVE_FFTW3) @@ -1240,72 +1241,130 @@ DET Matrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const { - DET retval; + MatrixType mattype (*this); + return determinant (mattype, info, rcon, calc_cond); +} + +DET +Matrix::determinant (MatrixType& mattype, + octave_idx_type& info, double& rcon, int calc_cond) const +{ + DET retval (1.0); octave_idx_type nr = rows (); octave_idx_type nc = cols (); - if (nr == 0 || nc == 0) - { - retval = DET (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 (); - - Matrix atmp = *this; - double *tmp_data = atmp.fortran_vec (); - - info = 0; - - // Calculate the norm of the matrix, for later use. - double anorm = 0; - if (calc_cond) - anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); - - F77_XFCN (dgetrf, DGETRF, (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 = DET (); - } - else - { - if (calc_cond) - { - // Now calc the condition number for non-singular matrix. - char job = '1'; - Array<double> z (4 * nc); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dgecon, DGECON, (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 = DET (); - } - else - { - retval = DET (1.0); - - for (octave_idx_type i = 0; i < nc; i++) - { - double 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) + { + Matrix atmp = *this; + double *tmp_data = atmp.fortran_vec (); + + info = 0; + double anorm = 0; + if (calc_cond) anorm = xnorm (*this, 1); + + + char job = 'L'; + F77_XFCN (dpotrf, DPOTRF, (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<double> z (3 * nc); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dpocon, DPOCON, (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 (); + + Matrix atmp = *this; + double *tmp_data = atmp.fortran_vec (); + + info = 0; + + // Calculate the norm of the matrix, for later use. + double anorm = 0; + if (calc_cond) anorm = xnorm (*this, 1); + + F77_XFCN (dgetrf, DGETRF, (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 = DET (); + } + else + { + if (calc_cond) + { + // Now calc the condition number for non-singular matrix. + char job = '1'; + Array<double> z (4 * nc); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dgecon, DGECON, (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 = DET (); + } + else + { + for (octave_idx_type i = 0; i < nc; i++) + { + double c = atmp(i,i); + retval *= (ipvt(i) != (i+1)) ? -c : c; + } + } + } + } } return retval;
--- a/liboctave/dMatrix.h Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/dMatrix.h Wed Nov 19 15:26:39 2008 +0100 @@ -145,6 +145,8 @@ DET determinant (void) const; DET determinant (octave_idx_type& info) const; DET determinant (octave_idx_type& info, double& rcon, int calc_cond = 1) const; + DET determinant (MatrixType &mattype, octave_idx_type& info, + double& rcon, int calc_cond = 1) const; double rcond (void) const; double rcond (MatrixType &mattype) const;
--- a/liboctave/fCMatrix.cc Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/fCMatrix.cc Wed Nov 19 15:26:39 2008 +0100 @@ -3,6 +3,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -54,6 +55,7 @@ #include "mx-fcm-fs.h" #include "mx-inlines.cc" #include "oct-cmplx.h" +#include "oct-norm.h" #if defined (HAVE_FFTW3) #include "oct-fftw.h" @@ -1567,74 +1569,132 @@ FloatComplexDET FloatComplexMatrix::determinant (octave_idx_type& info, float& rcon, int calc_cond) const { - FloatComplexDET retval; + MatrixType mattype (*this); + return determinant (mattype, info, rcon, calc_cond); +} + +FloatComplexDET +FloatComplexMatrix::determinant (MatrixType& mattype, + octave_idx_type& info, float& rcon, int calc_cond) const +{ + FloatComplexDET retval (1.0); octave_idx_type nr = rows (); octave_idx_type nc = cols (); - if (nr == 0 || nc == 0) - { - retval = FloatComplexDET (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 (); - - FloatComplexMatrix atmp = *this; - FloatComplex *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 (cgetrf, CGETRF, (nr, nc, 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 = FloatComplexDET (); - } - else - { - if (calc_cond) - { - // Now calc the condition number for non-singular matrix. - char job = '1'; - Array<FloatComplex> z (2*nr); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (2*nr); - float *prz = rz.fortran_vec (); - - F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcon, pz, prz, info - F77_CHAR_ARG_LEN (1))); - } - - if (info != 0) - { - info = -1; - retval = FloatComplexDET (); - } - else - { - retval = FloatComplexDET (1.0); - - for (octave_idx_type i = 0; i < nc; i++) - { - FloatComplex 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) + { + FloatComplexMatrix atmp = *this; + FloatComplex *tmp_data = atmp.fortran_vec (); + + info = 0; + float anorm = 0; + if (calc_cond) anorm = xnorm (*this, 1); + + + char job = 'L'; + F77_XFCN (cpotrf, CPOTRF, (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<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcon, pz, prz, 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 (); + + FloatComplexMatrix atmp = *this; + FloatComplex *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 (cgetrf, CGETRF, (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 = FloatComplexDET (); + } + else + { + if (calc_cond) + { + // Now calc the condition number for non-singular matrix. + char job = '1'; + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (2 * nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, prz, info + F77_CHAR_ARG_LEN (1))); } - } - } + + if (info != 0) + { + info = -1; + retval = FloatComplexDET (); + } + else + { + for (octave_idx_type i = 0; i < nc; i++) + { + FloatComplex c = atmp(i,i); + retval *= (ipvt(i) != (i+1)) ? -c : c; + } + } + } + } } - + return retval; }
--- a/liboctave/fCMatrix.h Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/fCMatrix.h Wed Nov 19 15:26:39 2008 +0100 @@ -176,6 +176,8 @@ FloatComplexDET determinant (void) const; FloatComplexDET determinant (octave_idx_type& info) const; FloatComplexDET determinant (octave_idx_type& info, float& rcon, int calc_cond = 1) const; + FloatComplexDET determinant (MatrixType &mattype, octave_idx_type& info, + float& rcon, int calc_cond = 1) const; float rcond (void) const; float rcond (MatrixType &mattype) const;
--- 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;
--- a/liboctave/fMatrix.h Wed Nov 19 11:23:07 2008 +0100 +++ b/liboctave/fMatrix.h Wed Nov 19 15:26:39 2008 +0100 @@ -145,6 +145,8 @@ FloatDET determinant (void) const; FloatDET determinant (octave_idx_type& info) const; FloatDET determinant (octave_idx_type& info, float& rcon, int calc_cond = 1) const; + FloatDET determinant (MatrixType &mattype, octave_idx_type& info, + float& rcon, int calc_cond = 1) const; float rcond (void) const; float rcond (MatrixType &mattype) const;