Mercurial > octave-nkf
changeset 7792:39c1026191e9
add missing files from single-precision merge
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Wed, 21 May 2008 09:36:46 -0400 |
| parents | 975e9540be2c |
| children | 96ba591be50f |
| files | liboctave/CmplxGEPBAL.cc liboctave/CmplxGEPBAL.h liboctave/dbleGEPBAL.cc liboctave/dbleGEPBAL.h liboctave/fCmplxAEPBAL.cc liboctave/fCmplxAEPBAL.h liboctave/fCmplxGEPBAL.cc liboctave/fCmplxGEPBAL.h liboctave/fCmplxHESS.cc liboctave/fCmplxHESS.h liboctave/fCmplxQR.cc liboctave/fCmplxQR.h liboctave/fCmplxQRP.cc liboctave/fCmplxQRP.h liboctave/floatAEPBAL.cc liboctave/floatAEPBAL.h liboctave/floatGEPBAL.cc liboctave/floatGEPBAL.h liboctave/floatHESS.cc liboctave/floatHESS.h liboctave/floatQR.cc liboctave/floatQR.h liboctave/floatQRP.cc liboctave/floatQRP.h |
| diffstat | 24 files changed, 2870 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxGEPBAL.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,130 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, + 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 <string> +#include <vector> + +#include "CmplxGEPBAL.h" +#include "Array-util.h" +#include "f77-fcn.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type& N, + Complex* A, const octave_idx_type& LDA, Complex* B, + const octave_idx_type& LDB, octave_idx_type& ILO, octave_idx_type& IHI, + double* LSCALE, double* RSCALE, + double* WORK, octave_idx_type& INFO + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type& N, const octave_idx_type& ILO, + const octave_idx_type& IHI, const double* LSCALE, + const double* RSCALE, octave_idx_type& M, double* V, + const octave_idx_type& LDV, octave_idx_type& INFO + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + +} + +octave_idx_type +ComplexGEPBALANCE::init (const ComplexMatrix& a, const ComplexMatrix& b, + const std::string& balance_job) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + { + (*current_liboctave_error_handler) ("ComplexGEPBALANCE requires square matrix"); + return -1; + } + + if (a.dims() != b.dims ()) + { + gripe_nonconformant ("ComplexGEPBALANCE", n, n, b.rows(), b.cols()); + return -1; + } + + octave_idx_type info; + octave_idx_type ilo; + octave_idx_type ihi; + + OCTAVE_LOCAL_BUFFER (double, plscale, n); + OCTAVE_LOCAL_BUFFER (double, prscale, n); + OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n); + + balanced_mat = a; + Complex *p_balanced_mat = balanced_mat.fortran_vec (); + balanced_mat2 = b; + Complex *p_balanced_mat2 = balanced_mat2.fortran_vec (); + + char job = balance_job[0]; + + F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1), + n, p_balanced_mat, n, p_balanced_mat2, + n, ilo, ihi, plscale, prscale, pwork, info + F77_CHAR_ARG_LEN (1))); + + balancing_mat = Matrix (n, n, 0.0); + balancing_mat2 = Matrix (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + { + OCTAVE_QUIT; + balancing_mat.elem (i ,i) = 1.0; + balancing_mat2.elem (i ,i) = 1.0; + } + + double *p_balancing_mat = balancing_mat.fortran_vec (); + double *p_balancing_mat2 = balancing_mat2.fortran_vec (); + + // first left + F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 ("L", 1), + n, ilo, ihi, plscale, prscale, + n, p_balancing_mat, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + // then right + F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 ("R", 1), + n, ilo, ihi, plscale, prscale, + n, p_balancing_mat2, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/CmplxGEPBAL.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,91 @@ +/* + +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_ComplexGEPBALANCE_h) +#define octave_ComplexGEPBALANCE_h 1 + +#include <iostream> +#include <string> + +#include "CMatrix.h" +#include "dMatrix.h" + +class +OCTAVE_API +ComplexGEPBALANCE +{ +public: + + ComplexGEPBALANCE (void) : balanced_mat (), balancing_mat () { } + + ComplexGEPBALANCE (const ComplexMatrix& a, const ComplexMatrix& b, const std::string& balance_job) + { + init (a, b, balance_job); + } + + ComplexGEPBALANCE (const ComplexGEPBALANCE& a) + : balanced_mat (a.balanced_mat), balanced_mat2 (a.balanced_mat2), + balancing_mat (a.balancing_mat), balancing_mat2 (a.balancing_mat2) { } + + ComplexGEPBALANCE& operator = (const ComplexGEPBALANCE& a) + { + if (this != &a) + { + balanced_mat = a.balanced_mat; + balanced_mat2 = a.balanced_mat2; + balancing_mat = a.balancing_mat; + balancing_mat2 = a.balancing_mat2; + } + return *this; + } + + ~ComplexGEPBALANCE (void) { } + + ComplexMatrix balanced_matrix (void) const { return balanced_mat; } + + ComplexMatrix balanced_matrix2 (void) const { return balanced_mat2; } + + Matrix balancing_matrix (void) const { return balancing_mat; } + + Matrix balancing_matrix2 (void) const { return balancing_mat2; } + + friend std::ostream& operator << (std::ostream& os, const ComplexGEPBALANCE& a); + +private: + + ComplexMatrix balanced_mat; + ComplexMatrix balanced_mat2; + Matrix balancing_mat; + Matrix balancing_mat2; + + octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, + const std::string& balance_job); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleGEPBAL.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,130 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, + 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 <string> +#include <vector> + +#include "dbleGEPBAL.h" +#include "Array-util.h" +#include "f77-fcn.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type& N, + double* A, const octave_idx_type& LDA, double* B, + const octave_idx_type& LDB, octave_idx_type& ILO, octave_idx_type& IHI, + double* LSCALE, double* RSCALE, + double* WORK, octave_idx_type& INFO + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type& N, const octave_idx_type& ILO, + const octave_idx_type& IHI, const double* LSCALE, + const double* RSCALE, octave_idx_type& M, double* V, + const octave_idx_type& LDV, octave_idx_type& INFO + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + +} + +octave_idx_type +GEPBALANCE::init (const Matrix& a, const Matrix& b, + const std::string& balance_job) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + { + (*current_liboctave_error_handler) ("GEPBALANCE requires square matrix"); + return -1; + } + + if (a.dims() != b.dims ()) + { + gripe_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols()); + return -1; + } + + octave_idx_type info; + octave_idx_type ilo; + octave_idx_type ihi; + + OCTAVE_LOCAL_BUFFER (double, plscale, n); + OCTAVE_LOCAL_BUFFER (double, prscale, n); + OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n); + + balanced_mat = a; + double *p_balanced_mat = balanced_mat.fortran_vec (); + balanced_mat2 = b; + double *p_balanced_mat2 = balanced_mat2.fortran_vec (); + + char job = balance_job[0]; + + F77_XFCN (dggbal, DGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1), + n, p_balanced_mat, n, p_balanced_mat2, + n, ilo, ihi, plscale, prscale, pwork, info + F77_CHAR_ARG_LEN (1))); + + balancing_mat = Matrix (n, n, 0.0); + balancing_mat2 = Matrix (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + { + OCTAVE_QUIT; + balancing_mat.elem (i ,i) = 1.0; + balancing_mat2.elem (i ,i) = 1.0; + } + + double *p_balancing_mat = balancing_mat.fortran_vec (); + double *p_balancing_mat2 = balancing_mat2.fortran_vec (); + + // first left + F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 ("L", 1), + n, ilo, ihi, plscale, prscale, + n, p_balancing_mat, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + // then right + F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 ("R", 1), + n, ilo, ihi, plscale, prscale, + n, p_balancing_mat2, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/dbleGEPBAL.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,90 @@ +/* + +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_GEPBALANCE_h) +#define octave_GEPBALANCE_h 1 + +#include <iostream> +#include <string> + +#include "dMatrix.h" + +class +OCTAVE_API +GEPBALANCE +{ +public: + + GEPBALANCE (void) : balanced_mat (), balancing_mat () { } + + GEPBALANCE (const Matrix& a, const Matrix& b, const std::string& balance_job) + { + init (a, b, balance_job); + } + + GEPBALANCE (const GEPBALANCE& a) + : balanced_mat (a.balanced_mat), balanced_mat2 (a.balanced_mat2), + balancing_mat (a.balancing_mat), balancing_mat2 (a.balancing_mat2) { } + + GEPBALANCE& operator = (const GEPBALANCE& a) + { + if (this != &a) + { + balanced_mat = a.balanced_mat; + balanced_mat2 = a.balanced_mat2; + balancing_mat = a.balancing_mat; + balancing_mat2 = a.balancing_mat2; + } + return *this; + } + + ~GEPBALANCE (void) { } + + Matrix balanced_matrix (void) const { return balanced_mat; } + + Matrix balanced_matrix2 (void) const { return balanced_mat2; } + + Matrix balancing_matrix (void) const { return balancing_mat; } + + Matrix balancing_matrix2 (void) const { return balancing_mat2; } + + friend std::ostream& operator << (std::ostream& os, const GEPBALANCE& a); + +private: + + Matrix balanced_mat; + Matrix balanced_mat2; + Matrix balancing_mat; + Matrix balancing_mat2; + + octave_idx_type init (const Matrix& a, const Matrix& b, + const std::string& balance_job); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxAEPBAL.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,102 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, + 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 <string> + +#include "fCmplxAEPBAL.h" +#include "fMatrix.h" +#include "f77-fcn.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (cgebal, CGEBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + octave_idx_type&, octave_idx_type&, float*, + octave_idx_type& F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (cgebak, CGEBAK) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + const octave_idx_type&, float*, + const octave_idx_type&, FloatComplex*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); +} + +octave_idx_type +FloatComplexAEPBALANCE::init (const FloatComplexMatrix& a, + const std::string& balance_job) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + { + (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); + return -1; + } + + octave_idx_type info; + octave_idx_type ilo; + octave_idx_type ihi; + + Array<float> scale (n); + float *pscale = scale.fortran_vec (); + + balanced_mat = a; + FloatComplex *p_balanced_mat = balanced_mat.fortran_vec (); + + char job = balance_job[0]; + + F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), + n, p_balanced_mat, n, ilo, ihi, + pscale, info + F77_CHAR_ARG_LEN (1))); + + balancing_mat = FloatComplexMatrix (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + balancing_mat.elem (i, i) = 1.0; + + FloatComplex *p_balancing_mat = balancing_mat.fortran_vec (); + + char side = 'R'; + + F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, pscale, n, + p_balancing_mat, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxAEPBAL.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,80 @@ +/* + +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_FloatComplexAEPBALANCE_h) +#define octave_FloatComplexAEPBALANCE_h 1 + +#include <iostream> +#include <string> + +#include "fCMatrix.h" + +class +OCTAVE_API +FloatComplexAEPBALANCE +{ +public: + + FloatComplexAEPBALANCE (void) : balanced_mat (), balancing_mat () { } + + FloatComplexAEPBALANCE (const FloatComplexMatrix& a, const std::string& balance_job) + { + init (a, balance_job); + } + + FloatComplexAEPBALANCE (const FloatComplexAEPBALANCE& a) + : balanced_mat (a.balanced_mat), balancing_mat (a.balancing_mat) { } + + FloatComplexAEPBALANCE& operator = (const FloatComplexAEPBALANCE& a) + { + if (this != &a) + { + balanced_mat = a.balanced_mat; + balancing_mat = a.balancing_mat; + } + return *this; + } + + ~FloatComplexAEPBALANCE (void) { } + + FloatComplexMatrix balanced_matrix (void) const { return balanced_mat; } + + FloatComplexMatrix balancing_matrix (void) const { return balancing_mat; } + + friend std::ostream& operator << (std::ostream& os, const FloatComplexAEPBALANCE& a); + +private: + + FloatComplexMatrix balanced_mat; + FloatComplexMatrix balancing_mat; + + octave_idx_type init (const FloatComplexMatrix& a, const std::string& balance_job); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxGEPBAL.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,130 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, + 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 <string> +#include <vector> + +#include "fCmplxGEPBAL.h" +#include "Array-util.h" +#include "f77-fcn.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (cggbal, CGGBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type& N, + FloatComplex* A, const octave_idx_type& LDA, FloatComplex* B, + const octave_idx_type& LDB, octave_idx_type& ILO, octave_idx_type& IHI, + float* LSCALE, float* RSCALE, + float* WORK, octave_idx_type& INFO + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (sggbak, SGGBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type& N, const octave_idx_type& ILO, + const octave_idx_type& IHI, const float* LSCALE, + const float* RSCALE, octave_idx_type& M, float* V, + const octave_idx_type& LDV, octave_idx_type& INFO + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + +} + +octave_idx_type +FloatComplexGEPBALANCE::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, + const std::string& balance_job) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + { + (*current_liboctave_error_handler) ("FloatComplexGEPBALANCE requires square matrix"); + return -1; + } + + if (a.dims() != b.dims ()) + { + gripe_nonconformant ("FloatComplexGEPBALANCE", n, n, b.rows(), b.cols()); + return -1; + } + + octave_idx_type info; + octave_idx_type ilo; + octave_idx_type ihi; + + OCTAVE_LOCAL_BUFFER (float, plscale, n); + OCTAVE_LOCAL_BUFFER (float, prscale, n); + OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n); + + balanced_mat = a; + FloatComplex *p_balanced_mat = balanced_mat.fortran_vec (); + balanced_mat2 = b; + FloatComplex *p_balanced_mat2 = balanced_mat2.fortran_vec (); + + char job = balance_job[0]; + + F77_XFCN (cggbal, CGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1), + n, p_balanced_mat, n, p_balanced_mat2, + n, ilo, ihi, plscale,prscale, pwork, info + F77_CHAR_ARG_LEN (1))); + + balancing_mat = FloatMatrix (n, n, 0.0); + balancing_mat2 = FloatMatrix (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + { + OCTAVE_QUIT; + balancing_mat.elem (i ,i) = 1.0; + balancing_mat2.elem (i ,i) = 1.0; + } + + float *p_balancing_mat = balancing_mat.fortran_vec (); + float *p_balancing_mat2 = balancing_mat2.fortran_vec (); + + // first left + F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 ("L", 1), + n, ilo, ihi, plscale, prscale, + n, p_balancing_mat, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + // then right + F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 ("R", 1), + n, ilo, ihi, plscale, prscale, + n, p_balancing_mat2, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxGEPBAL.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,91 @@ +/* + +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_FloatComplexGEPBALANCE_h) +#define octave_FloatComplexGEPBALANCE_h 1 + +#include <iostream> +#include <string> + +#include "fMatrix.h" +#include "fCMatrix.h" + +class +OCTAVE_API +FloatComplexGEPBALANCE +{ +public: + + FloatComplexGEPBALANCE (void) : balanced_mat (), balancing_mat () { } + + FloatComplexGEPBALANCE (const FloatComplexMatrix& a, const FloatComplexMatrix& b, const std::string& balance_job) + { + init (a, b, balance_job); + } + + FloatComplexGEPBALANCE (const FloatComplexGEPBALANCE& a) + : balanced_mat (a.balanced_mat), balanced_mat2 (a.balanced_mat2), + balancing_mat (a.balancing_mat), balancing_mat2 (a.balancing_mat2) { } + + FloatComplexGEPBALANCE& operator = (const FloatComplexGEPBALANCE& a) + { + if (this != &a) + { + balanced_mat = a.balanced_mat; + balanced_mat2 = a.balanced_mat2; + balancing_mat = a.balancing_mat; + balancing_mat2 = a.balancing_mat2; + } + return *this; + } + + ~FloatComplexGEPBALANCE (void) { } + + FloatComplexMatrix balanced_matrix (void) const { return balanced_mat; } + + FloatComplexMatrix balanced_matrix2 (void) const { return balanced_mat2; } + + FloatMatrix balancing_matrix (void) const { return balancing_mat; } + + FloatMatrix balancing_matrix2 (void) const { return balancing_mat2; } + + friend std::ostream& operator << (std::ostream& os, const FloatComplexGEPBALANCE& a); + +private: + + FloatComplexMatrix balanced_mat; + FloatComplexMatrix balanced_mat2; + FloatMatrix balancing_mat; + FloatMatrix balancing_mat2; + + octave_idx_type init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, + const std::string& balance_job); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxHESS.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,127 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 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 "fCmplxHESS.h" +#include "f77-fcn.h" +#include "lo-error.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (cgebal, CGEBAL) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, FloatComplex*, const octave_idx_type&, + octave_idx_type&, octave_idx_type&, float*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (cgehrd, CGEHRD) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, FloatComplex*, + FloatComplex*, const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (cunghr, CUNGHR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, FloatComplex*, + FloatComplex*, const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (cgebak, CGEBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, float*, + const octave_idx_type&, FloatComplex*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); +} + +octave_idx_type +FloatComplexHESS::init (const FloatComplexMatrix& a) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + if (a_nr != a_nc) + { + (*current_liboctave_error_handler) + ("FloatComplexHESS requires square matrix"); + return -1; + } + + char job = 'N'; + char side = 'R'; + + octave_idx_type n = a_nc; + octave_idx_type lwork = 32 * n; + octave_idx_type info; + octave_idx_type ilo; + octave_idx_type ihi; + + hess_mat = a; + FloatComplex *h = hess_mat.fortran_vec (); + + Array<float> scale (n); + float *pscale = scale.fortran_vec (); + + F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), + n, h, n, ilo, ihi, pscale, info + F77_CHAR_ARG_LEN (1))); + + Array<FloatComplex> tau (n-1); + FloatComplex *ptau = tau.fortran_vec (); + + Array<FloatComplex> work (lwork); + FloatComplex *pwork = work.fortran_vec (); + + F77_XFCN (cgehrd, CGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info)); + + unitary_hess_mat = hess_mat; + FloatComplex *z = unitary_hess_mat.fortran_vec (); + + F77_XFCN (cunghr, CUNGHR, (n, ilo, ihi, z, n, ptau, pwork, + lwork, info)); + + F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, pscale, n, z, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + // If someone thinks of a more graceful way of + // doing this (or faster for that matter :-)), + // please let me know! + + if (n > 2) + for (octave_idx_type j = 0; j < a_nc; j++) + for (octave_idx_type i = j+2; i < a_nr; i++) + hess_mat.elem (i, j) = 0; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxHESS.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,81 @@ +/* + +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_FloatComplexHESS_h) +#define octave_FloatComplexHESS_h 1 + +#include <iostream> + +#include "fCMatrix.h" + +class +OCTAVE_API +FloatComplexHESS +{ +public: + + FloatComplexHESS (void) : hess_mat (), unitary_hess_mat () { } + + FloatComplexHESS (const FloatComplexMatrix& a) { init (a); } + + FloatComplexHESS (const FloatComplexMatrix& a, octave_idx_type& info) { info = init (a); } + + FloatComplexHESS (const FloatComplexHESS& a) + : hess_mat (a.hess_mat), unitary_hess_mat (a.unitary_hess_mat) { } + + FloatComplexHESS& operator = (const FloatComplexHESS& a) + { + if (this != &a) + { + hess_mat = a.hess_mat; + unitary_hess_mat = a.unitary_hess_mat; + } + return *this; + } + + ~FloatComplexHESS (void) { } + + FloatComplexMatrix hess_matrix (void) const { return hess_mat; } + + FloatComplexMatrix unitary_hess_matrix (void) const + { + return unitary_hess_mat; + } + + friend std::ostream& operator << (std::ostream& os, const FloatComplexHESS& a); + +private: + + FloatComplexMatrix hess_mat; + FloatComplexMatrix unitary_hess_mat; + + octave_idx_type init (const FloatComplexMatrix& a); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxQR.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,307 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 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/>. + +*/ + +// updating/downdating by Jaroslav Hajek 2008 + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "fCmplxQR.h" +#include "f77-fcn.h" +#include "lo-error.h" +#include "Range.h" +#include "idx-vector.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (cgeqrf, CGEQRF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*, + const octave_idx_type&, FloatComplex*, FloatComplex*, + const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (cungqr, CUNGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, FloatComplex*, + FloatComplex*, const octave_idx_type&, octave_idx_type&); + + // these come from qrupdate + + F77_RET_T + F77_FUNC (cqr1up, CQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + FloatComplex*, FloatComplex*, const FloatComplex*, const FloatComplex*); + + F77_RET_T + F77_FUNC (cqrinc, CQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const FloatComplex*, FloatComplex*, const octave_idx_type&, const FloatComplex*); + + F77_RET_T + F77_FUNC (cqrdec, CQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const FloatComplex*, FloatComplex*, const octave_idx_type&); + + F77_RET_T + F77_FUNC (cqrinr, CQRINR) (const octave_idx_type&, const octave_idx_type&, + const FloatComplex*, FloatComplex*, const FloatComplex*, FloatComplex*, + const octave_idx_type&, const FloatComplex*); + + F77_RET_T + F77_FUNC (cqrder, CQRDER) (const octave_idx_type&, const octave_idx_type&, + const FloatComplex*, FloatComplex*, const FloatComplex*, FloatComplex *, + const octave_idx_type&); + + F77_RET_T + F77_FUNC (cqrshc, CQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + FloatComplex*, FloatComplex*, const octave_idx_type&, const octave_idx_type&); +} + +FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& a, QR::type qr_type) + : q (), r () +{ + init (a, qr_type); +} + +void +FloatComplexQR::init (const FloatComplexMatrix& a, QR::type qr_type) +{ + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + if (m == 0 || n == 0) + { + (*current_liboctave_error_handler) + ("FloatComplexQR must have non-empty matrix"); + return; + } + + octave_idx_type min_mn = m < n ? m : n; + + Array<FloatComplex> tau (min_mn); + FloatComplex *ptau = tau.fortran_vec (); + + octave_idx_type lwork = 32*n; + Array<FloatComplex> work (lwork); + FloatComplex *pwork = work.fortran_vec (); + + octave_idx_type info = 0; + + FloatComplexMatrix A_fact; + if (m > n && qr_type != QR::economy) + { + A_fact.resize (m, m); + A_fact.insert (a, 0, 0); + } + else + A_fact = a; + + FloatComplex *tmp_data = A_fact.fortran_vec (); + + F77_XFCN (cgeqrf, CGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info)); + + if (qr_type == QR::raw) + { + for (octave_idx_type j = 0; j < min_mn; j++) + { + octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; + for (octave_idx_type i = limit + 1; i < m; i++) + A_fact.elem (i, j) *= tau.elem (j); + } + + r = A_fact; + + if (m > n) + r.resize (m, n); + } + else + { + octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; + + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); + else + r.resize (m, n, 0.0); + + for (octave_idx_type j = 0; j < n; j++) + { + octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; + for (octave_idx_type i = 0; i <= limit; i++) + r.elem (i, j) = A_fact.elem (i, j); + } + + lwork = 32 * n2; + work.resize (lwork); + FloatComplex *pwork2 = work.fortran_vec (); + + F77_XFCN (cungqr, CUNGQR, (m, n2, min_mn, tmp_data, m, ptau, + pwork2, lwork, info)); + + q = A_fact; + q.resize (m, n2); + } +} + +FloatComplexQR::FloatComplexQR (const FloatComplexMatrix &q, const FloatComplexMatrix& r) +{ + if (q.columns () != r.rows ()) + { + (*current_liboctave_error_handler) ("QR dimensions mismatch"); + return; + } + + this->q = q; + this->r = r; +} + +void +FloatComplexQR::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () == m && v.length () == n) + F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), + u.data (), v.data ())); + else + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); +} + +void +FloatComplexQR::insert_col (const FloatComplexMatrix& u, octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () != m) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("QR insert index out of range"); + else + { + FloatComplexMatrix r1 (m,n+1); + + F77_XFCN (cqrinc, CQRINC, (m, n, k, q.fortran_vec (), r.data (), + r1.fortran_vec (), j+1, u.data ())); + + r = r1; + } +} + +void +FloatComplexQR::delete_col (octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type k = r.rows (); + octave_idx_type n = r.columns (); + + if (k < m && k < n) + (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + else if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("QR delete index out of range"); + else + { + FloatComplexMatrix r1 (k, n-1); + + F77_XFCN (cqrdec, CQRDEC, (m, n, k, q.fortran_vec (), r.data (), + r1.fortran_vec (), j+1)); + + r = r1; + } +} + +void +FloatComplexQR::insert_row (const FloatComplexMatrix& u, octave_idx_type j) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square () || u.length () != n) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 0 || j > m) + (*current_liboctave_error_handler) ("QR insert index out of range"); + else + { + FloatComplexMatrix q1 (m+1, m+1); + FloatComplexMatrix r1 (m+1, n); + + F77_XFCN (cqrinr, CQRINR, (m, n, q.data (), q1.fortran_vec (), + r.data (), r1.fortran_vec (), j+1, u.data ())); + + q = q1; + r = r1; + } +} + +void +FloatComplexQR::delete_row (octave_idx_type j) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square ()) + (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + else if (j < 0 || j > m-1) + (*current_liboctave_error_handler) ("QR delete index out of range"); + else + { + FloatComplexMatrix q1 (m-1, m-1); + FloatComplexMatrix r1 (m-1, n); + + F77_XFCN (cqrder, CQRDER, (m, n, q.data (), q1.fortran_vec (), + r.data (), r1.fortran_vec (), j+1 )); + + q = q1; + r = r1; + } +} + +void +FloatComplexQR::shift_cols (octave_idx_type i, octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type k = r.rows (); + octave_idx_type n = r.columns (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("QR shift index out of range"); + else + F77_XFCN (cqrshc, CQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1)); +} + +void +FloatComplexQR::economize (void) +{ + octave_idx_type r_nc = r.columns (); + + if (r.rows () > r_nc) + { + q.resize (q.rows (), r_nc); + r.resize (r_nc, r_nc); + } +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxQR.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,96 @@ +/* + +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/>. + +*/ + +// updating/downdating by Jaroslav Hajek 2008 + +#if !defined (octave_FloatComplexQR_h) +#define octave_FloatComplexQR_h 1 + +#include <iostream> + +#include "fCMatrix.h" +#include "fCColVector.h" +#include "fCRowVector.h" +#include "dbleQR.h" + +class +OCTAVE_API +FloatComplexQR +{ +public: + + FloatComplexQR (void) : q (), r () { } + + FloatComplexQR (const FloatComplexMatrix&, QR::type = QR::std); + + FloatComplexQR (const FloatComplexMatrix& q, const FloatComplexMatrix& r); + + FloatComplexQR (const FloatComplexQR& a) : q (a.q), r (a.r) { } + + FloatComplexQR& operator = (const FloatComplexQR& a) + { + if (this != &a) + { + q = a.q; + r = a.r; + } + return *this; + } + + ~FloatComplexQR (void) { } + + void init (const FloatComplexMatrix&, QR::type = QR::std); + + FloatComplexMatrix Q (void) const { return q; } + + FloatComplexMatrix R (void) const { return r; } + + void update (const FloatComplexMatrix& u, const FloatComplexMatrix& v); + + void insert_col (const FloatComplexMatrix& u, octave_idx_type j); + + void delete_col (octave_idx_type j); + + void insert_row (const FloatComplexMatrix& u, octave_idx_type j); + + void delete_row (octave_idx_type j); + + void shift_cols (octave_idx_type i, octave_idx_type j); + + void economize(); + + friend std::ostream& operator << (std::ostream&, const FloatComplexQR&); + +protected: + + FloatComplexMatrix q; + FloatComplexMatrix r; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxQRP.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,138 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 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 "fCmplxQRP.h" +#include "f77-fcn.h" +#include "lo-error.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (cgeqpf, CGEQPF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*, + const octave_idx_type&, octave_idx_type*, FloatComplex*, FloatComplex*, + float*, octave_idx_type&); + + F77_RET_T + F77_FUNC (cungqr, CUNGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, FloatComplex*, + FloatComplex*, const octave_idx_type&, octave_idx_type&); +} + +// It would be best to share some of this code with FloatComplexQR class... + +FloatComplexQRP::FloatComplexQRP (const FloatComplexMatrix& a, QR::type qr_type) + : FloatComplexQR (), p () +{ + init (a, qr_type); +} + +void +FloatComplexQRP::init (const FloatComplexMatrix& a, QR::type qr_type) +{ + assert (qr_type != QR::raw); + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + if (m == 0 || n == 0) + { + (*current_liboctave_error_handler) + ("FloatComplexQR must have non-empty matrix"); + return; + } + + octave_idx_type min_mn = m < n ? m : n; + Array<FloatComplex> tau (min_mn); + FloatComplex *ptau = tau.fortran_vec (); + + octave_idx_type lwork = 3*n > 32*m ? 3*n : 32*m; + Array<FloatComplex> work (lwork); + FloatComplex *pwork = work.fortran_vec (); + + octave_idx_type info = 0; + + FloatComplexMatrix A_fact = a; + if (m > n && qr_type != QR::economy) + A_fact.resize (m, m, 0.0); + + FloatComplex *tmp_data = A_fact.fortran_vec (); + + Array<float> rwork (2*n); + float *prwork = rwork.fortran_vec (); + + Array<octave_idx_type> jpvt (n, 0); + octave_idx_type *pjpvt = jpvt.fortran_vec (); + + // Code to enforce a certain permutation could go here... + + F77_XFCN (cgeqpf, CGEQPF, (m, n, tmp_data, m, pjpvt, ptau, pwork, + prwork, info)); + + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + if (qr_type == QR::economy) + { + p.resize (1, n, 0.0); + for (octave_idx_type j = 0; j < n; j++) + p.elem (0, j) = jpvt.elem (j); + } + else + { + p.resize (n, n, 0.0); + for (octave_idx_type j = 0; j < n; j++) + p.elem (jpvt.elem (j) - 1, j) = 1.0; + } + + octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; + + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); + else + r.resize (m, n, 0.0); + + for (octave_idx_type j = 0; j < n; j++) + { + octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; + for (octave_idx_type i = 0; i <= limit; i++) + r.elem (i, j) = A_fact.elem (i, j); + } + + F77_XFCN (cungqr, CUNGQR, (m, n2, min_mn, tmp_data, m, ptau, + pwork, lwork, info)); + + q = A_fact; + q.resize (m, n2); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/fCmplxQRP.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,73 @@ +/* + +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_FloatComplexQRP_h) +#define octave_FloatComplexQRP_h 1 + +#include <iostream> + +#include "fCmplxQR.h" +#include "fMatrix.h" + +class +OCTAVE_API +FloatComplexQRP : public FloatComplexQR +{ +public: + + FloatComplexQRP (void) : FloatComplexQR (), p () { } + + FloatComplexQRP (const FloatComplexMatrix&, QR::type = QR::std); + + FloatComplexQRP (const FloatComplexQRP& a) : FloatComplexQR (a), p (a.p) { } + + FloatComplexQRP& operator = (const FloatComplexQRP& a) + { + if (this != &a) + { + FloatComplexQR::operator = (a); + p = a.p; + } + return *this; + } + + ~FloatComplexQRP (void) { } + + void init (const FloatComplexMatrix&, QR::type = QR::std); + + FloatMatrix P (void) const { return p; } + + friend std::ostream& operator << (std::ostream&, const FloatComplexQRP&); + +private: + + FloatMatrix p; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatAEPBAL.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,99 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, + 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 <string> + +#include "floatAEPBAL.h" +#include "f77-fcn.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (sgebal, SGEBAL) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, float*, const octave_idx_type&, octave_idx_type&, + octave_idx_type&, float*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (sgebak, SGEBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, float*, + const octave_idx_type&, float*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); +} + +octave_idx_type +FloatAEPBALANCE::init (const FloatMatrix& a, const std::string& balance_job) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + { + (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); + return -1; + } + + octave_idx_type info; + octave_idx_type ilo; + octave_idx_type ihi; + + Array<float> scale (n); + float *pscale = scale.fortran_vec (); + + balanced_mat = a; + float *p_balanced_mat = balanced_mat.fortran_vec (); + + char job = balance_job[0]; + + F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), + n, p_balanced_mat, n, ilo, ihi, pscale, info + F77_CHAR_ARG_LEN (1))); + + balancing_mat = Matrix (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + balancing_mat.elem (i ,i) = 1.0; + + float *p_balancing_mat = balancing_mat.fortran_vec (); + + char side = 'R'; + + F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, pscale, n, + p_balancing_mat, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatAEPBAL.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,80 @@ +/* + +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_FloatAEPBALANCE_h) +#define octave_FloatAEPBALANCE_h 1 + +#include <iostream> +#include <string> + +#include "fMatrix.h" + +class +OCTAVE_API +FloatAEPBALANCE +{ +public: + + FloatAEPBALANCE (void) : balanced_mat (), balancing_mat () { } + + FloatAEPBALANCE (const FloatMatrix& a,const std::string& balance_job) + { + init (a, balance_job); + } + + FloatAEPBALANCE (const FloatAEPBALANCE& a) + : balanced_mat (a.balanced_mat), balancing_mat (a.balancing_mat) { } + + FloatAEPBALANCE& operator = (const FloatAEPBALANCE& a) + { + if (this != &a) + { + balanced_mat = a.balanced_mat; + balancing_mat = a.balancing_mat; + } + return *this; + } + + ~FloatAEPBALANCE (void) { } + + FloatMatrix balanced_matrix (void) const { return balanced_mat; } + + FloatMatrix balancing_matrix (void) const { return balancing_mat; } + + friend std::ostream& operator << (std::ostream& os, const FloatAEPBALANCE& a); + +private: + + FloatMatrix balanced_mat; + FloatMatrix balancing_mat; + + octave_idx_type init (const FloatMatrix& a, const std::string& balance_job); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatGEPBAL.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,130 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, + 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 <string> +#include <vector> + +#include "floatGEPBAL.h" +#include "Array-util.h" +#include "f77-fcn.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (sggbal, SGGBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type& N, + float* A, const octave_idx_type& LDA, float* B, + const octave_idx_type& LDB, octave_idx_type& ILO, octave_idx_type& IHI, + float* LSCALE, float* RSCALE, + float* WORK, octave_idx_type& INFO + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (sggbak, SGGBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type& N, const octave_idx_type& ILO, + const octave_idx_type& IHI, const float* LSCALE, + const float* RSCALE, octave_idx_type& M, float* V, + const octave_idx_type& LDV, octave_idx_type& INFO + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + +} + +octave_idx_type +FloatGEPBALANCE::init (const FloatMatrix& a, const FloatMatrix& b, + const std::string& balance_job) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + { + (*current_liboctave_error_handler) ("FloatGEPBALANCE requires square matrix"); + return -1; + } + + if (a.dims() != b.dims ()) + { + gripe_nonconformant ("FloatGEPBALANCE", n, n, b.rows(), b.cols()); + return -1; + } + + octave_idx_type info; + octave_idx_type ilo; + octave_idx_type ihi; + + OCTAVE_LOCAL_BUFFER (float, plscale, n); + OCTAVE_LOCAL_BUFFER (float, prscale, n); + OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n); + + balanced_mat = a; + float *p_balanced_mat = balanced_mat.fortran_vec (); + balanced_mat2 = b; + float *p_balanced_mat2 = balanced_mat2.fortran_vec (); + + char job = balance_job[0]; + + F77_XFCN (sggbal, SGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1), + n, p_balanced_mat, n, p_balanced_mat2, + n, ilo, ihi, plscale, prscale, pwork, info + F77_CHAR_ARG_LEN (1))); + + balancing_mat = FloatMatrix (n, n, 0.0); + balancing_mat2 = FloatMatrix (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + { + OCTAVE_QUIT; + balancing_mat.elem (i ,i) = 1.0; + balancing_mat2.elem (i ,i) = 1.0; + } + + float *p_balancing_mat = balancing_mat.fortran_vec (); + float *p_balancing_mat2 = balancing_mat2.fortran_vec (); + + // first left + F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 ("L", 1), + n, ilo, ihi, plscale, prscale, + n, p_balancing_mat, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + // then right + F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 ("R", 1), + n, ilo, ihi, plscale, prscale, + n, p_balancing_mat2, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatGEPBAL.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,90 @@ +/* + +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_FloatGEPBALANCE_h) +#define octave_FloatGEPBALANCE_h 1 + +#include <iostream> +#include <string> + +#include "fMatrix.h" + +class +OCTAVE_API +FloatGEPBALANCE +{ +public: + + FloatGEPBALANCE (void) : balanced_mat (), balancing_mat () { } + + FloatGEPBALANCE (const FloatMatrix& a, const FloatMatrix& b, const std::string& balance_job) + { + init (a, b, balance_job); + } + + FloatGEPBALANCE (const FloatGEPBALANCE& a) + : balanced_mat (a.balanced_mat), balanced_mat2 (a.balanced_mat2), + balancing_mat (a.balancing_mat), balancing_mat2 (a.balancing_mat2) { } + + FloatGEPBALANCE& operator = (const FloatGEPBALANCE& a) + { + if (this != &a) + { + balanced_mat = a.balanced_mat; + balanced_mat2 = a.balanced_mat2; + balancing_mat = a.balancing_mat; + balancing_mat2 = a.balancing_mat2; + } + return *this; + } + + ~FloatGEPBALANCE (void) { } + + FloatMatrix balanced_matrix (void) const { return balanced_mat; } + + FloatMatrix balanced_matrix2 (void) const { return balanced_mat2; } + + FloatMatrix balancing_matrix (void) const { return balancing_mat; } + + FloatMatrix balancing_matrix2 (void) const { return balancing_mat2; } + + friend std::ostream& operator << (std::ostream& os, const FloatGEPBALANCE& a); + +private: + + FloatMatrix balanced_mat; + FloatMatrix balanced_mat2; + FloatMatrix balancing_mat; + FloatMatrix balancing_mat2; + + octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b, + const std::string& balance_job); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatHESS.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,128 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 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 "floatHESS.h" +#include "f77-fcn.h" +#include "lo-error.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (sgebal, SGEBAL) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, float*, const octave_idx_type&, octave_idx_type&, + octave_idx_type&, float*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (sgehrd, SGEHRD) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + float*, const octave_idx_type&, float*, float*, + const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (sorghr, SORGHR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + float*, const octave_idx_type&, float*, float*, + const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (sgebak, SGEBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, float*, + const octave_idx_type&, float*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); +} + +octave_idx_type +FloatHESS::init (const FloatMatrix& a) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + if (a_nr != a_nc) + { + (*current_liboctave_error_handler) ("FloatHESS requires square matrix"); + return -1; + } + + char job = 'N'; + char side = 'R'; + + octave_idx_type n = a_nc; + octave_idx_type lwork = 32 * n; + octave_idx_type info; + octave_idx_type ilo; + octave_idx_type ihi; + + hess_mat = a; + float *h = hess_mat.fortran_vec (); + + Array<float> scale (n); + float *pscale = scale.fortran_vec (); + + F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), + n, h, n, ilo, ihi, pscale, info + F77_CHAR_ARG_LEN (1))); + + Array<float> tau (n-1); + float *ptau = tau.fortran_vec (); + + Array<float> work (lwork); + float *pwork = work.fortran_vec (); + + F77_XFCN (sgehrd, SGEHRD, (n, ilo, ihi, h, n, ptau, pwork, + lwork, info)); + + unitary_hess_mat = hess_mat; + float *z = unitary_hess_mat.fortran_vec (); + + F77_XFCN (sorghr, SORGHR, (n, ilo, ihi, z, n, ptau, pwork, + lwork, info)); + + F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, pscale, n, z, + n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + // If someone thinks of a more graceful way of doing + // this (or faster for that matter :-)), please let + // me know! + + if (n > 2) + for (octave_idx_type j = 0; j < a_nc; j++) + for (octave_idx_type i = j+2; i < a_nr; i++) + hess_mat.elem (i, j) = 0; + + return info; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatHESS.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,78 @@ +/* + +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_FloatHESS_h) +#define octave_FloatHESS_h 1 + +#include <iostream> + +#include "fMatrix.h" + +class +OCTAVE_API +FloatHESS +{ +public: + + FloatHESS (void) : hess_mat (), unitary_hess_mat () { } + + FloatHESS (const FloatMatrix& a) { init (a); } + + FloatHESS (const FloatMatrix& a, octave_idx_type& info) { info = init (a); } + + FloatHESS (const FloatHESS& a) + : hess_mat (a.hess_mat), unitary_hess_mat (a.unitary_hess_mat) { } + + FloatHESS& operator = (const FloatHESS& a) + { + if (this != &a) + { + hess_mat = a.hess_mat; + unitary_hess_mat = a.unitary_hess_mat; + } + return *this; + } + + ~FloatHESS (void) { } + + FloatMatrix hess_matrix (void) const { return hess_mat; } + + FloatMatrix unitary_hess_matrix (void) const { return unitary_hess_mat; } + + friend std::ostream& operator << (std::ostream& os, const FloatHESS& a); + +private: + + FloatMatrix hess_mat; + FloatMatrix unitary_hess_mat; + + octave_idx_type init (const FloatMatrix& a); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatQR.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,297 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 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 "floatQR.h" +#include "f77-fcn.h" +#include "lo-error.h" +#include "Range.h" +#include "idx-vector.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (sgeqrf, SGEQRF) (const octave_idx_type&, const octave_idx_type&, float*, const octave_idx_type&, + float*, float*, const octave_idx_type&, octave_idx_type&); + + F77_RET_T + F77_FUNC (sorgqr, SORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, float*, + const octave_idx_type&, float*, float*, const octave_idx_type&, octave_idx_type&); + + // these come from qrupdate + + F77_RET_T + F77_FUNC (sqr1up, SQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + float*, float*, const float*, const float*); + + F77_RET_T + F77_FUNC (sqrinc, SQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + float*, const float*, float*, const octave_idx_type&, const float*); + + F77_RET_T + F77_FUNC (sqrdec, SQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + float*, const float*, float*, const octave_idx_type&); + + F77_RET_T + F77_FUNC (sqrinr, SQRINR) (const octave_idx_type&, const octave_idx_type&, + const float*, float*, const float*, float*, + const octave_idx_type&, const float*); + + F77_RET_T + F77_FUNC (sqrder, SQRDER) (const octave_idx_type&, const octave_idx_type&, + const float*, float*, const float*, float *, + const octave_idx_type&); + + F77_RET_T + F77_FUNC (sqrshc, SQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + float*, float*, const octave_idx_type&, const octave_idx_type&); +} + +FloatQR::FloatQR (const FloatMatrix& a, QR::type qr_type) + : q (), r () +{ + init (a, qr_type); +} + +void +FloatQR::init (const FloatMatrix& a, QR::type qr_type) +{ + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + if (m == 0 || n == 0) + { + (*current_liboctave_error_handler) ("QR must have non-empty matrix"); + return; + } + + octave_idx_type min_mn = m < n ? m : n; + Array<float> tau (min_mn); + float *ptau = tau.fortran_vec (); + + octave_idx_type lwork = 32*n; + Array<float> work (lwork); + float *pwork = work.fortran_vec (); + + octave_idx_type info = 0; + + FloatMatrix A_fact = a; + if (m > n && qr_type != QR::economy) + A_fact.resize (m, m, 0.0); + + float *tmp_data = A_fact.fortran_vec (); + + F77_XFCN (sgeqrf, SGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info)); + + if (qr_type == QR::raw) + { + for (octave_idx_type j = 0; j < min_mn; j++) + { + octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; + for (octave_idx_type i = limit + 1; i < m; i++) + A_fact.elem (i, j) *= tau.elem (j); + } + + r = A_fact; + + if (m > n) + r.resize (m, n); + } + else + { + octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; + + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); + else + r.resize (m, n, 0.0); + + for (octave_idx_type j = 0; j < n; j++) + { + octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; + for (octave_idx_type i = 0; i <= limit; i++) + r.elem (i, j) = tmp_data[m*j+i]; + } + + lwork = 32 * n2; + work.resize (lwork); + float *pwork2 = work.fortran_vec (); + + F77_XFCN (sorgqr, SORGQR, (m, n2, min_mn, tmp_data, m, ptau, + pwork2, lwork, info)); + + q = A_fact; + q.resize (m, n2); + } +} + +FloatQR::FloatQR (const FloatMatrix& q, const FloatMatrix& r) +{ + if (q.columns () != r.rows ()) + { + (*current_liboctave_error_handler) ("QR dimensions mismatch"); + return; + } + + this->q = q; + this->r = r; +} + +void +FloatQR::update (const FloatMatrix& u, const FloatMatrix& v) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () == m && v.length () == n) + F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), + u.data (), v.data ())); + else + (*current_liboctave_error_handler) ("QR update dimensions mismatch"); +} + +void +FloatQR::insert_col (const FloatMatrix& u, octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type n = r.columns (); + octave_idx_type k = q.columns (); + + if (u.length () != m) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 0 || j > n) + (*current_liboctave_error_handler) ("QR insert index out of range"); + else + { + FloatMatrix r1 (m, n+1); + + F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), r.data (), + r1.fortran_vec (), j+1, u.data ())); + + r = r1; + } +} + +void +FloatQR::delete_col (octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type k = r.rows (); + octave_idx_type n = r.columns (); + + if (k < m && k < n) + (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + else if (j < 0 || j > n-1) + (*current_liboctave_error_handler) ("QR delete index out of range"); + else + { + FloatMatrix r1 (k, n-1); + + F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), r.data (), + r1.fortran_vec (), j+1)); + + r = r1; + } +} + +void +FloatQR::insert_row (const FloatMatrix& u, octave_idx_type j) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square () || u.length () != n) + (*current_liboctave_error_handler) ("QR insert dimensions mismatch"); + else if (j < 0 || j > m) + (*current_liboctave_error_handler) ("QR insert index out of range"); + else + { + FloatMatrix q1 (m+1, m+1); + FloatMatrix r1 (m+1, n); + + F77_XFCN (sqrinr, SQRINR, (m, n, q.data (), q1.fortran_vec (), + r.data (), r1.fortran_vec (), j+1, u.data ())); + + q = q1; + r = r1; + } +} + +void +FloatQR::delete_row (octave_idx_type j) +{ + octave_idx_type m = r.rows (); + octave_idx_type n = r.columns (); + + if (! q.is_square ()) + (*current_liboctave_error_handler) ("QR delete dimensions mismatch"); + else if (j < 0 || j > m-1) + (*current_liboctave_error_handler) ("QR delete index out of range"); + else + { + FloatMatrix q1 (m-1, m-1); + FloatMatrix r1 (m-1, n); + + F77_XFCN (sqrder, SQRDER, (m, n, q.data (), q1.fortran_vec (), + r.data (), r1.fortran_vec (), j+1 )); + + q = q1; + r = r1; + } +} + +void +FloatQR::shift_cols (octave_idx_type i, octave_idx_type j) +{ + octave_idx_type m = q.rows (); + octave_idx_type k = r.rows (); + octave_idx_type n = r.columns (); + + if (i < 0 || i > n-1 || j < 0 || j > n-1) + (*current_liboctave_error_handler) ("QR shift index out of range"); + else + F77_XFCN (sqrshc, SQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1)); +} + +void +FloatQR::economize (void) +{ + octave_idx_type r_nc = r.columns (); + + if (r.rows () > r_nc) + { + q.resize (q.rows (), r_nc); + r.resize (r_nc, r_nc); + } +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatQR.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,96 @@ +/* + +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/>. + +*/ + +// updating/downdating by Jaroslav Hajek 2008 + +#if !defined (octave_FloatQR_h) +#define octave_FloatQR_h 1 + +#include <iostream> + +#include "fMatrix.h" +#include "fColVector.h" +#include "fRowVector.h" +#include "dbleQR.h" + +class +OCTAVE_API +FloatQR +{ +public: + + FloatQR (void) : q (), r () { } + + FloatQR (const FloatMatrix&, QR::type = QR::std); + + FloatQR (const FloatMatrix& q, const FloatMatrix& r); + + FloatQR (const FloatQR& a) : q (a.q), r (a.r) { } + + FloatQR& operator = (const FloatQR& a) + { + if (this != &a) + { + q = a.q; + r = a.r; + } + return *this; + } + + ~FloatQR (void) { } + + void init (const FloatMatrix&, QR::type); + + FloatMatrix Q (void) const { return q; } + + FloatMatrix R (void) const { return r; } + + void update (const FloatMatrix& u, const FloatMatrix& v); + + void insert_col (const FloatMatrix& u, octave_idx_type j); + + void delete_col (octave_idx_type j); + + void insert_row (const FloatMatrix& u, octave_idx_type j); + + void delete_row (octave_idx_type j); + + void shift_cols (octave_idx_type i, octave_idx_type j); + + void economize (void); + + friend std::ostream& operator << (std::ostream&, const FloatQR&); + +protected: + + FloatMatrix q; + FloatMatrix r; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatQRP.cc Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,132 @@ +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 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 "floatQRP.h" +#include "f77-fcn.h" +#include "lo-error.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (sgeqpf, SGEQPF) (const octave_idx_type&, const octave_idx_type&, float*, + const octave_idx_type&, octave_idx_type*, float*, float*, octave_idx_type&); + + F77_RET_T + F77_FUNC (sorgqr, SORGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + float*, const octave_idx_type&, float*, float*, + const octave_idx_type&, octave_idx_type&); +} + +// It would be best to share some of this code with QR class... + +FloatQRP::FloatQRP (const FloatMatrix& a, QR::type qr_type) + : FloatQR (), p () +{ + init (a, qr_type); +} + +void +FloatQRP::init (const FloatMatrix& a, QR::type qr_type) +{ + assert (qr_type != QR::raw); + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + if (m == 0 || n == 0) + { + (*current_liboctave_error_handler) ("QR must have non-empty matrix"); + return; + } + + octave_idx_type min_mn = m < n ? m : n; + Array<float> tau (min_mn); + float *ptau = tau.fortran_vec (); + + octave_idx_type lwork = 3*n > 32*m ? 3*n : 32*m; + Array<float> work (lwork); + float *pwork = work.fortran_vec (); + + octave_idx_type info = 0; + + FloatMatrix A_fact = a; + if (m > n && qr_type != QR::economy) + A_fact.resize (m, m, 0.0); + + float *tmp_data = A_fact.fortran_vec (); + + Array<octave_idx_type> jpvt (n, 0); + octave_idx_type *pjpvt = jpvt.fortran_vec (); + + // Code to enforce a certain permutation could go here... + + F77_XFCN (sgeqpf, SGEQPF, (m, n, tmp_data, m, pjpvt, ptau, pwork, info)); + + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + if (qr_type == QR::economy) + { + p.resize (1, n, 0.0); + for (octave_idx_type j = 0; j < n; j++) + p.elem (0, j) = jpvt.elem (j); + } + else + { + p.resize (n, n, 0.0); + for (octave_idx_type j = 0; j < n; j++) + p.elem (jpvt.elem (j) - 1, j) = 1.0; + } + + octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; + + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); + else + r.resize (m, n, 0.0); + + for (octave_idx_type j = 0; j < n; j++) + { + octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; + for (octave_idx_type i = 0; i <= limit; i++) + r.elem (i, j) = A_fact.elem (i, j); + } + + F77_XFCN (sorgqr, SORGQR, (m, n2, min_mn, tmp_data, m, ptau, + pwork, lwork, info)); + + q = A_fact; + q.resize (m, n2); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/floatQRP.h Wed May 21 09:36:46 2008 -0400 @@ -0,0 +1,74 @@ +/* + +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_FloatQRP_h) +#define octave_FloatQRP_h 1 + +#include <iostream> + +#include "floatQR.h" +#include "fMatrix.h" + +class +OCTAVE_API +FloatQRP : public FloatQR +{ +public: + + FloatQRP (void) : FloatQR (), p () { } + + FloatQRP (const FloatMatrix&, QR::type = QR::std); + + FloatQRP (const FloatQRP& a) : FloatQR (a), p (a.p) { } + + FloatQRP& operator = (const FloatQRP& a) + { + if (this != &a) + { + FloatQR::operator = (a); + p = a.p; + } + + return *this; + } + + ~FloatQRP (void) { } + + void init (const FloatMatrix&, QR::type = QR::std); + + FloatMatrix P (void) const { return p; } + + friend std::ostream& operator << (std::ostream&, const FloatQRP&); + +protected: + + FloatMatrix p; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/