Mercurial > hg > octave-lojdl
changeset 8335:64cf956a109c
templatize & fix DET
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Wed, 19 Nov 2008 11:23:07 +0100 |
parents | 70dd33450061 |
children | 9813c07ca946 |
files | liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/CSparse.cc liboctave/CSparse.h liboctave/ChangeLog liboctave/CmplxDET.cc liboctave/CmplxDET.h liboctave/DET.h liboctave/Makefile.in liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/dSparse.cc liboctave/dSparse.h liboctave/dbleDET.cc liboctave/dbleDET.h liboctave/fCMatrix.cc liboctave/fCMatrix.h liboctave/fCmplxDET.cc liboctave/fCmplxDET.h liboctave/fMatrix.cc liboctave/fMatrix.h liboctave/floatDET.cc liboctave/floatDET.h liboctave/mx-defs.h liboctave/mx-ext.h src/ChangeLog src/DLD-FUNCTIONS/det.cc |
diffstat | 27 files changed, 158 insertions(+), 964 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -40,7 +40,7 @@ #include "Array-util.h" #include "CMatrix.h" #include "CmplxAEPBAL.h" -#include "CmplxDET.h" +#include "DET.h" #include "CmplxSCHUR.h" #include "CmplxSVD.h" #include "CmplxCHOL.h" @@ -1631,33 +1631,13 @@ } else { - Complex c = 1.0; - int e = 0; - + retval = ComplexDET (1.0); + for (octave_idx_type i = 0; i < nc; i++) { - if (ipvt(i) != (i+1)) - c = -c; - - c *= atmp(i,i); - - if (c == 0.0) - break; - - while (std::abs(c) < 0.5) - { - c *= 2.0; - e--; - } - - while (std::abs(c) >= 2.0) - { - c /= 2.0; - e++; - } - } - - retval = ComplexDET (c, e); + Complex c = atmp(i,i); + retval *= (ipvt(i) != (i+1)) ? -c : c; + } } } }
--- a/liboctave/CMatrix.h +++ b/liboctave/CMatrix.h @@ -31,6 +31,7 @@ #include "mx-defs.h" #include "mx-op-defs.h" #include "oct-cmplx.h" +#include "DET.h" class OCTAVE_API
--- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -1114,10 +1114,7 @@ if (nr == 0 || nc == 0 || nr != nc) { - Complex d[2]; - d[0] = 1.0; - d[1] = 0.0; - retval = ComplexDET (d); + retval = ComplexDET (1.0); } else { @@ -1201,13 +1198,10 @@ { UMFPACK_ZNAME (report_numeric) (Numeric, control); - Complex d[2]; - double d_exponent; - - status = UMFPACK_ZNAME (get_determinant) - (reinterpret_cast<double *> (&d[0]), 0, &d_exponent, - Numeric, info); - d[1] = d_exponent; + double c10[2], e10; + + status = UMFPACK_ZNAME (get_determinant) (c10, 0, &e10, + Numeric, info); if (status < 0) { @@ -1218,7 +1212,7 @@ UMFPACK_ZNAME (report_info) (control, info); } else - retval = ComplexDET (d); + retval = ComplexDET (Complex (c10[0], c10[1]), e10, 10); UMFPACK_ZNAME (free_numeric) (&Numeric); }
--- a/liboctave/CSparse.h +++ b/liboctave/CSparse.h @@ -31,7 +31,7 @@ #include "CColVector.h" #include "oct-cmplx.h" -#include "CmplxDET.h" +#include "DET.h" #include "MSparse.h" #include "MSparse-defs.h" #include "Sparse-op-defs.h"
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,19 @@ +2008-11-19 Jaroslav Hajek <highegg@gmail.com> + + * DET.h: New source. + * CmplxDET.cc, CmplxDET.h, dbleDET.cc, dbleDET.h, fCmplxDET.cc, + fCmplxDET.h, floatDET.cc, floatDET.h: Remove. + * Makefile.in: Reflect changes. + * mx-defs.h: Remove DET decls. + * mx-ext.h, dMatrix.h, fMatrix.h, CMatrix.h, fCMatrix.h, + dSparse.h, CSparse.h: Include only DET.h. + * dMatrix.cc (Matrix::determinant), + fMatrix.cc (FloatMatrix::determinant), + CMatrix.cc (ComplexMatrix::determinant), + fCMatrix.cc (FloatComplexMatrix::determinant), + dSparse.cc (SparseMatrix::determinant), + CSparse.cc (SparseComplexMatrix::determinant): Use new class. + 2008-11-18 David Bateman <dbateman@free.fr> * file-ops.cc (std::string file_ops::tilde_expand (const
deleted file mode 100644 --- a/liboctave/CmplxDET.cc +++ /dev/null @@ -1,86 +0,0 @@ -/* - -Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2006, - 2007 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <cassert> -#include <cfloat> - -#include "CmplxDET.h" -#include "lo-mappers.h" -#include "lo-math.h" -#include "oct-cmplx.h" - -bool -ComplexDET::value_will_overflow (void) const -{ - return base2 - ? (e2 + 1 > xlog2 (DBL_MAX) ? 1 : 0) - : (e10 + 1 > log10 (DBL_MAX) ? 1 : 0); -} - -bool -ComplexDET::value_will_underflow (void) const -{ - return base2 - ? (e2 - 1 < xlog2 (DBL_MIN) ? 1 : 0) - : (e10 - 1 < log10 (DBL_MIN) ? 1 : 0); -} - -void -ComplexDET::initialize10 (void) -{ - if (c2 != 0.0) - { - double etmp = e2 / xlog2 (static_cast<double>(10)); - e10 = static_cast<int> (xround (etmp)); - etmp -= e10; - c10 = c2 * pow (10.0, etmp); - } -} - -void -ComplexDET::initialize2 (void) -{ - if (c10 != 0.0) - { - double etmp = e10 / log10 (2.0); - e2 = static_cast<int> (xround (etmp)); - etmp -= e2; - c2 = c10 * xexp2 (etmp); - } -} - -Complex -ComplexDET::value (void) const -{ - return base2 ? c2 * xexp2 (static_cast<double>(e2)) : c10 * pow (10.0, e10); -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
deleted file mode 100644 --- a/liboctave/CmplxDET.h +++ /dev/null @@ -1,121 +0,0 @@ -/* - -Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, - 2007 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#if !defined (octave_ComplexDET_h) -#define octave_ComplexDET_h 1 - -#include <iostream> - -#include "oct-cmplx.h" - -// FIXME -- we could use templates here; compare with dbleDET.h - -class -OCTAVE_API -ComplexDET -{ -friend class ComplexMatrix; -friend class SparseComplexMatrix; - -public: - - ComplexDET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { } - - ComplexDET (const ComplexDET& a) - : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2) - { } - - ComplexDET& operator = (const ComplexDET& a) - { - if (this != &a) - { - c2 = a.c2; - e2 = a.e2; - - c10 = a.c10; - e10 = a.e10; - - base2 = a.base2; - } - return *this; - } - - bool value_will_overflow (void) const; - bool value_will_underflow (void) const; - - // These two functions were originally defined in base 10, so we are - // preserving that interface here. - - Complex coefficient (void) const { return coefficient10 (); } - int exponent (void) const { return exponent10 (); } - - Complex coefficient10 (void) const { return c10; } - int exponent10 (void) const { return e10; } - - Complex coefficient2 (void) const { return c2; } - int exponent2 (void) const { return e2; } - - Complex value (void) const; - - friend std::ostream& operator << (std::ostream& os, const ComplexDET& a); - -private: - - // Constructed this way, we assume base 2. - - ComplexDET (const Complex& c, int e) - : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true) - { - initialize10 (); - } - - // Original interface had only this constructor and it was assumed - // to be base 10, so we are preserving that interface here. - - ComplexDET (const Complex *d) - : c2 (0), c10 (d[0]), e2 (0), e10 (static_cast<int> (d[1].real ())), - base2 (false) - { - initialize2 (); - } - - void initialize2 (void); - void initialize10 (void); - - Complex c2; - Complex c10; - - int e2; - int e10; - - // TRUE means the original values were provided in base 2. - bool base2; -}; - -#endif - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
new file mode 100644 --- /dev/null +++ b/liboctave/DET.h @@ -0,0 +1,87 @@ +/* + +Copyright (C) 2008 Jaroslav Hajek + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#if !defined (octave_DET_h) +#define octave_DET_h 1 + +#include <cmath> +#include "oct-cmplx.h" +#include "lo-mappers.h" + +template <class T> +class +OCTAVE_API +base_det +{ +public: + + base_det (T c = 0, int e = 0) + { + c2 = xlog2 (c, e2); + e2 += e; + } + + base_det (T c, double e, double b) + { + e *= xlog2 (b); + e2 = e; + c *= xexp2 (e - e2); + int f; + c2 = xlog2 (c, f); + e2 += f; + } + + base_det (const base_det& a) : c2(a.c2), e2(a.e2) { } + + base_det& operator = (const base_det& a) + { + c2 = a.c2; + e2 = a.e2; + return *this; + } + + T coef (void) const { return c2; } + int exp (void) const { return e2; } + + T value () const { return c2 * static_cast<T> (std::ldexp (1.0, e2)); } + operator T () const { return value (); } + + void operator *= (T t) + { + int e; + c2 *= xlog2 (t, e); + e2 += e; + } + +private: + + T c2; + int e2; +}; + +// Provide the old types by typedefs. +typedef base_det<double> DET; +typedef base_det<float> FloatDET; +typedef base_det<Complex> ComplexDET; +typedef base_det<FloatComplex> FloatComplexDET; + +#endif
--- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -46,10 +46,10 @@ base-lu.h dim-vector.h mx-base.h mx-op-defs.h \ mx-defs.h mx-ext.h CColVector.h CDiagMatrix.h CMatrix.h \ CNDArray.h CRowVector.h CmplxAEPBAL.h CmplxCHOL.h \ - CmplxDET.h CmplxGEPBAL.h CmplxHESS.h CmplxLU.h CmplxQR.h CmplxQRP.h \ + CmplxGEPBAL.h CmplxHESS.h CmplxLU.h CmplxQR.h CmplxQRP.h \ CmplxSCHUR.h CmplxSVD.h EIG.h fEIG.h boolMatrix.h boolNDArray.h \ chMatrix.h chNDArray.h dColVector.h dDiagMatrix.h dMatrix.h \ - dNDArray.h dRowVector.h dbleAEPBAL.h dbleCHOL.h dbleDET.h \ + dNDArray.h dRowVector.h dbleAEPBAL.h dbleCHOL.h DET.h \ dbleGEPBAL.h dbleHESS.h dbleLU.h dbleQR.h dbleQRP.h dbleSCHUR.h \ dbleSVD.h boolSparse.h CSparse.h dSparse.h MSparse-defs.h MSparse.h \ Sparse.h sparse-base-lu.h SparseCmplxLU.h SparsedbleLU.h \ @@ -62,9 +62,9 @@ fCColVector.h fCRowVector.h fCDiagMatrix.h fCMatrix.h fCNDArray.h \ fColVector.h fRowVector.h fDiagMatrix.h fMatrix.h fNDArray.h \ fCmplxAEPBAL.h fCmplxGEPBAL.h fCmplxHESS.h fCmplxCHOL.h \ - fCmplxDET.h fCmplxLU.h fCmplxSCHUR.h fCmplxSVD.h fCmplxQR.h \ + fCmplxLU.h fCmplxSCHUR.h fCmplxSVD.h fCmplxQR.h \ fCmplxQRP.h floatAEPBAL.h \ - floatCHOL.h floatDET.h floatGEPBAL.h floatHESS.h floatLU.h \ + floatCHOL.h floatGEPBAL.h floatHESS.h floatLU.h \ floatSCHUR.h floatSVD.h floatQR.h floatQRP.h MX_OP_INC := $(shell $(AWK) -f $(srcdir)/mk-ops.awk prefix=mx list_h_files=1 $(srcdir)/mx-ops) @@ -116,12 +116,12 @@ MATRIX_SRC := Array-util.cc CColVector.cc \ CDiagMatrix.cc CMatrix.cc CNDArray.cc CRowVector.cc \ - CmplxAEPBAL.cc CmplxCHOL.cc CmplxDET.cc CmplxGEPBAL.cc CmplxHESS.cc \ + CmplxAEPBAL.cc CmplxCHOL.cc CmplxGEPBAL.cc CmplxHESS.cc \ CmplxLU.cc CmplxQR.cc CmplxQRP.cc CmplxSCHUR.cc CmplxSVD.cc \ EIG.cc fEIG.cc boolMatrix.cc boolNDArray.cc chMatrix.cc \ chNDArray.cc dColVector.cc dDiagMatrix.cc dMatrix.cc \ dNDArray.cc dRowVector.cc dbleAEPBAL.cc dbleCHOL.cc \ - dbleDET.cc dbleGEPBAL.cc dbleHESS.cc dbleLU.cc dbleQR.cc dbleQRP.cc \ + dbleGEPBAL.cc dbleHESS.cc dbleLU.cc dbleQR.cc dbleQRP.cc \ dbleSCHUR.cc dbleSVD.cc boolSparse.cc CSparse.cc dSparse.cc \ MSparse.cc Sparse.cc SparseCmplxLU.cc SparsedbleLU.cc \ SparseCmplxCHOL.cc SparsedbleCHOL.cc \ @@ -130,9 +130,9 @@ int32NDArray.cc uint32NDArray.cc int64NDArray.cc uint64NDArray.cc \ fCColVector.cc fCRowVector.cc fCDiagMatrix.cc fCMatrix.cc fCNDArray.cc \ fColVector.cc fRowVector.cc fDiagMatrix.cc fMatrix.cc fNDArray.cc \ - fCmplxAEPBAL.cc fCmplxCHOL.cc fCmplxDET.cc fCmplxGEPBAL.cc \ + fCmplxAEPBAL.cc fCmplxCHOL.cc fCmplxGEPBAL.cc \ fCmplxHESS.cc fCmplxLU.cc fCmplxSCHUR.cc fCmplxSVD.cc fCmplxQR.cc \ - fCmplxQRP.cc floatAEPBAL.cc floatCHOL.cc floatDET.cc \ + fCmplxQRP.cc floatAEPBAL.cc floatCHOL.cc \ floatGEPBAL.cc floatHESS.cc floatLU.cc \ floatSCHUR.cc floatSVD.cc floatQR.cc floatQRP.cc
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -36,7 +36,7 @@ #include "byte-swap.h" #include "dMatrix.h" #include "dbleAEPBAL.h" -#include "dbleDET.h" +#include "DET.h" #include "dbleSCHUR.h" #include "dbleSVD.h" #include "dbleCHOL.h" @@ -1297,33 +1297,13 @@ } else { - double c = 1.0; - int e = 0; - + retval = DET (1.0); + for (octave_idx_type i = 0; i < nc; i++) { - if (ipvt(i) != (i+1)) - c = -c; - - c *= atmp(i,i); - - if (c == 0.0) - break; - - while (fabs (c) < 0.5) - { - c *= 2.0; - e--; - } - - while (fabs (c) >= 2.0) - { - c /= 2.0; - e++; - } - } - - retval = DET (c, e); + double c = atmp(i,i); + retval *= (ipvt(i) != (i+1)) ? -c : c; + } } } }
--- a/liboctave/dMatrix.h +++ b/liboctave/dMatrix.h @@ -30,6 +30,7 @@ #include "mx-defs.h" #include "mx-op-defs.h" +#include "DET.h" class OCTAVE_API
--- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -1188,10 +1188,7 @@ if (nr == 0 || nc == 0 || nr != nc) { - double d[2]; - d[0] = 1.0; - d[1] = 0.0; - retval = DET (d); + retval = DET (1.0); } else { @@ -1270,10 +1267,9 @@ { UMFPACK_DNAME (report_numeric) (Numeric, control); - double d[2]; - - status = UMFPACK_DNAME (get_determinant) (&d[0], - &d[1], Numeric, info); + double c10, e10; + + status = UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric, info); if (status < 0) { @@ -1284,7 +1280,7 @@ UMFPACK_DNAME (report_info) (control, info); } else - retval = DET (d); + retval = DET (c10, e10, 10); UMFPACK_DNAME (free_numeric) (&Numeric); }
--- a/liboctave/dSparse.h +++ b/liboctave/dSparse.h @@ -30,7 +30,7 @@ #include "dColVector.h" #include "CColVector.h" -#include "dbleDET.h" +#include "DET.h" #include "MSparse.h" #include "MSparse-defs.h" #include "Sparse-op-defs.h"
deleted file mode 100644 --- a/liboctave/dbleDET.cc +++ /dev/null @@ -1,84 +0,0 @@ -/* - -Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2006, - 2007 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <cfloat> - -#include "dbleDET.h" -#include "lo-mappers.h" -#include "lo-math.h" - -bool -DET::value_will_overflow (void) const -{ - return base2 - ? (e2 + 1 > xlog2 (DBL_MAX) ? 1 : 0) - : (e10 + 1 > log10 (DBL_MAX) ? 1 : 0); -} - -bool -DET::value_will_underflow (void) const -{ - return base2 - ? (e2 - 1 < xlog2 (DBL_MIN) ? 1 : 0) - : (e10 - 1 < log10 (DBL_MIN) ? 1 : 0); -} - -void -DET::initialize10 (void) -{ - if (c2 != 0.0) - { - double etmp = e2 / xlog2 (static_cast<double>(10)); - e10 = static_cast<int> (xround (etmp)); - etmp -= e10; - c10 = c2 * pow (10.0, etmp); - } -} - -void -DET::initialize2 (void) -{ - if (c10 != 0.0) - { - double etmp = e10 / log10 (2.0); - e2 = static_cast<int> (xround (etmp)); - etmp -= e2; - c2 = c10 * xexp2 (etmp); - } -} - -double -DET::value (void) const -{ - return base2 ? c2 * xexp2 (static_cast<double>(e2)) : c10 * pow (10.0, e10); -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
deleted file mode 100644 --- a/liboctave/dbleDET.h +++ /dev/null @@ -1,118 +0,0 @@ -/* - -Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, - 2007 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#if !defined (octave_DET_h) -#define octave_DET_h 1 - -#include <iostream> - -// FIXME -- we could use templates here; compare with CmplxDET.h - -class -OCTAVE_API -DET -{ -friend class Matrix; -friend class SparseMatrix; - -public: - - DET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { } - - DET (const DET& a) - : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2) - { } - - DET& operator = (const DET& a) - { - if (this != &a) - { - c2 = a.c2; - e2 = a.e2; - - c10 = a.c10; - e10 = a.e10; - - base2 = a.base2; - } - return *this; - } - - bool value_will_overflow (void) const; - bool value_will_underflow (void) const; - - // These two functions were originally defined in base 10, so we are - // preserving that interface here. - - double coefficient (void) const { return coefficient10 (); } - int exponent (void) const { return exponent10 (); } - - double coefficient10 (void) const { return c10; } - int exponent10 (void) const { return e10; } - - double coefficient2 (void) const { return c2; } - int exponent2 (void) const { return e2; } - - double value (void) const; - - friend std::ostream& operator << (std::ostream& os, const DET& a); - -private: - - // Constructed this way, we assume base 2. - - DET (double c, int e) - : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true) - { - initialize10 (); - } - - // Original interface had only this constructor and it was assumed - // to be base 10, so we are preserving that interface here. - - DET (const double *d) - : c2 (0), c10 (d[0]), e2 (0), e10 (static_cast<int> (d[1])), base2 (false) - { - initialize2 (); - } - - void initialize2 (void); - void initialize10 (void); - - double c2; - double c10; - - int e2; - int e10; - - // TRUE means the original values were provided in base 2. - bool base2; -}; - -#endif - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/liboctave/fCMatrix.cc +++ b/liboctave/fCMatrix.cc @@ -38,7 +38,7 @@ #include "Array-util.h" #include "fCMatrix.h" -#include "fCmplxDET.h" +#include "DET.h" #include "fCmplxSCHUR.h" #include "fCmplxSVD.h" #include "fCmplxCHOL.h" @@ -1624,33 +1624,13 @@ } else { - FloatComplex c = 1.0; - int e = 0; - + retval = FloatComplexDET (1.0); + for (octave_idx_type i = 0; i < nc; i++) { - if (ipvt(i) != (i+1)) - c = -c; - - c *= atmp(i,i); - - if (c == static_cast<float> (0.0)) - break; - - while (std::abs(c) < 0.5) - { - c *= 2.0; - e--; - } - - while (std::abs(c) >= 2.0) - { - c /= 2.0; - e++; - } - } - - retval = FloatComplexDET (c, e); + FloatComplex c = atmp(i,i); + retval *= (ipvt(i) != (i+1)) ? -c : c; + } } } }
--- a/liboctave/fCMatrix.h +++ b/liboctave/fCMatrix.h @@ -31,6 +31,7 @@ #include "mx-defs.h" #include "mx-op-defs.h" #include "oct-cmplx.h" +#include "DET.h" class OCTAVE_API
deleted file mode 100644 --- a/liboctave/fCmplxDET.cc +++ /dev/null @@ -1,86 +0,0 @@ -/* - -Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2006, - 2007 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <cassert> -#include <cfloat> - -#include "fCmplxDET.h" -#include "lo-mappers.h" -#include "lo-math.h" -#include "oct-cmplx.h" - -bool -FloatComplexDET::value_will_overflow (void) const -{ - return base2 - ? (e2 + 1 > xlog2 (FLT_MAX) ? 1 : 0) - : (e10 + 1 > log10 (FLT_MAX) ? 1 : 0); -} - -bool -FloatComplexDET::value_will_underflow (void) const -{ - return base2 - ? (e2 - 1 < xlog2 (FLT_MIN) ? 1 : 0) - : (e10 - 1 < log10 (FLT_MIN) ? 1 : 0); -} - -void -FloatComplexDET::initialize10 (void) -{ - if (c2 != static_cast<float> (0.0)) - { - float etmp = e2 / xlog2 (static_cast<float>(10)); - e10 = static_cast<int> (xround (etmp)); - etmp -= e10; - c10 = c2 * static_cast<float> (pow (10.0f, etmp)); - } -} - -void -FloatComplexDET::initialize2 (void) -{ - if (c10 != static_cast<float> (0.0)) - { - float etmp = e10 / log10 (2.0); - e2 = static_cast<int> (xround (etmp)); - etmp -= e2; - c2 = c10 * xexp2 (etmp); - } -} - -FloatComplex -FloatComplexDET::value (void) const -{ - return base2 ? c2 * xexp2 (static_cast<float>(e2)) : c10 * static_cast<float> (pow (10.0, e10)); -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
deleted file mode 100644 --- a/liboctave/fCmplxDET.h +++ /dev/null @@ -1,120 +0,0 @@ -/* - -Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, - 2007 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#if !defined (octave_FloatComplexDET_h) -#define octave_FloatComplexDET_h 1 - -#include <iostream> - -#include "oct-cmplx.h" - -// FIXME -- we could use templates here; compare with dbleDET.h - -class -OCTAVE_API -FloatComplexDET -{ -friend class FloatComplexMatrix; - -public: - - FloatComplexDET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { } - - FloatComplexDET (const FloatComplexDET& a) - : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2) - { } - - FloatComplexDET& operator = (const FloatComplexDET& a) - { - if (this != &a) - { - c2 = a.c2; - e2 = a.e2; - - c10 = a.c10; - e10 = a.e10; - - base2 = a.base2; - } - return *this; - } - - bool value_will_overflow (void) const; - bool value_will_underflow (void) const; - - // These two functions were originally defined in base 10, so we are - // preserving that interface here. - - FloatComplex coefficient (void) const { return coefficient10 (); } - int exponent (void) const { return exponent10 (); } - - FloatComplex coefficient10 (void) const { return c10; } - int exponent10 (void) const { return e10; } - - FloatComplex coefficient2 (void) const { return c2; } - int exponent2 (void) const { return e2; } - - FloatComplex value (void) const; - - friend std::ostream& operator << (std::ostream& os, const FloatComplexDET& a); - -private: - - // Constructed this way, we assume base 2. - - FloatComplexDET (const FloatComplex& c, int e) - : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true) - { - initialize10 (); - } - - // Original interface had only this constructor and it was assumed - // to be base 10, so we are preserving that interface here. - - FloatComplexDET (const FloatComplex *d) - : c2 (0), c10 (d[0]), e2 (0), e10 (static_cast<int> (d[1].real ())), - base2 (false) - { - initialize2 (); - } - - void initialize2 (void); - void initialize10 (void); - - FloatComplex c2; - FloatComplex c10; - - int e2; - int e10; - - // TRUE means the original values were provided in base 2. - bool base2; -}; - -#endif - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/liboctave/fMatrix.cc +++ b/liboctave/fMatrix.cc @@ -35,7 +35,7 @@ #include "Array-util.h" #include "byte-swap.h" #include "fMatrix.h" -#include "floatDET.h" +#include "DET.h" #include "floatSCHUR.h" #include "floatSVD.h" #include "floatCHOL.h" @@ -1296,33 +1296,13 @@ } else { - float c = 1.0; - int e = 0; - + retval = FloatDET (1.0); + for (octave_idx_type i = 0; i < nc; i++) { - if (ipvt(i) != (i+1)) - c = -c; - - c *= atmp(i,i); - - if (c == 0.0) - break; - - while (fabs (c) < 0.5) - { - c *= 2.0; - e--; - } - - while (fabs (c) >= 2.0) - { - c /= 2.0; - e++; - } - } - - retval = FloatDET (c, e); + float c = atmp(i,i); + retval *= (ipvt(i) != (i+1)) ? -c : c; + } } } }
--- a/liboctave/fMatrix.h +++ b/liboctave/fMatrix.h @@ -30,6 +30,7 @@ #include "mx-defs.h" #include "mx-op-defs.h" +#include "DET.h" class OCTAVE_API
deleted file mode 100644 --- a/liboctave/floatDET.cc +++ /dev/null @@ -1,84 +0,0 @@ -/* - -Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2006, - 2007 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <cfloat> - -#include "floatDET.h" -#include "lo-mappers.h" -#include "lo-math.h" - -bool -FloatDET::value_will_overflow (void) const -{ - return base2 - ? (e2 + 1 > xlog2 (FLT_MAX) ? 1 : 0) - : (e10 + 1 > log10 (FLT_MAX) ? 1 : 0); -} - -bool -FloatDET::value_will_underflow (void) const -{ - return base2 - ? (e2 - 1 < xlog2 (FLT_MIN) ? 1 : 0) - : (e10 - 1 < log10 (FLT_MIN) ? 1 : 0); -} - -void -FloatDET::initialize10 (void) -{ - if (c2 != 0.0) - { - float etmp = e2 / xlog2 (static_cast<float>(10)); - e10 = static_cast<int> (xround (etmp)); - etmp -= e10; - c10 = c2 * powf (10.0, etmp); - } -} - -void -FloatDET::initialize2 (void) -{ - if (c10 != 0.0) - { - float etmp = e10 / log10 (2.0); - e2 = static_cast<int> (xround (etmp)); - etmp -= e2; - c2 = c10 * xexp2 (etmp); - } -} - -float -FloatDET::value (void) const -{ - return base2 ? c2 * xexp2 (static_cast<float>(e2)) : c10 * powf (10.0, e10); -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
deleted file mode 100644 --- a/liboctave/floatDET.h +++ /dev/null @@ -1,117 +0,0 @@ -/* - -Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, - 2007 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#if !defined (octave_FloatDET_h) -#define octave_FloatDET_h 1 - -#include <iostream> - -// FIXME -- we could use templates here; compare with CmplxFloatDET.h - -class -OCTAVE_API -FloatDET -{ -friend class FloatMatrix; - -public: - - FloatDET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { } - - FloatDET (const FloatDET& a) - : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2) - { } - - FloatDET& operator = (const FloatDET& a) - { - if (this != &a) - { - c2 = a.c2; - e2 = a.e2; - - c10 = a.c10; - e10 = a.e10; - - base2 = a.base2; - } - return *this; - } - - bool value_will_overflow (void) const; - bool value_will_underflow (void) const; - - // These two functions were originally defined in base 10, so we are - // preserving that interface here. - - float coefficient (void) const { return coefficient10 (); } - int exponent (void) const { return exponent10 (); } - - float coefficient10 (void) const { return c10; } - int exponent10 (void) const { return e10; } - - float coefficient2 (void) const { return c2; } - int exponent2 (void) const { return e2; } - - float value (void) const; - - friend std::ostream& operator << (std::ostream& os, const FloatDET& a); - -private: - - // Constructed this way, we assume base 2. - - FloatDET (float c, int e) - : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true) - { - initialize10 (); - } - - // Original interface had only this constructor and it was assumed - // to be base 10, so we are preserving that interface here. - - FloatDET (const float *d) - : c2 (0), c10 (d[0]), e2 (0), e10 (static_cast<int> (d[1])), base2 (false) - { - initialize2 (); - } - - void initialize2 (void); - void initialize10 (void); - - float c2; - float c10; - - int e2; - int e10; - - // TRUE means the original values were provided in base 2. - bool base2; -}; - -#endif - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/liboctave/mx-defs.h +++ b/liboctave/mx-defs.h @@ -70,11 +70,6 @@ class FloatCHOL; class FloatComplexCHOL; -class DET; -class ComplexDET; -class FloatDET; -class FloatComplexDET; - class EIG; class HESS;
--- a/liboctave/mx-ext.h +++ b/liboctave/mx-ext.h @@ -31,10 +31,7 @@ // Result of a Determinant calculation. -#include "dbleDET.h" -#include "CmplxDET.h" -#include "floatDET.h" -#include "fCmplxDET.h" +#include "DET.h" // Result of a Cholesky Factorization