Mercurial > octave
changeset 21268:f08ae27289e4
better use of templates for balance classes
* aepbalance.h, aepbalance.cc: New files generated from base-aepbal.h,
CmplxAEPBAL.cc, CmplxAEPBAL.h, dbleAEPBAL.cc, dbleAEPBAL.h,
fCmplxAEPBAL.cc, fCmplxAEPBAL.h, floatAEPBAL.cc, and floatAEPBAL.h and
making them templates.
* gepbalance.h, gepbalance.cc: New files generate from CmplxGEPBAL.cc,
CmplxGEPBAL.h, dbleGEPBAL.cc, dbleGEPBAL.h, fCmplxGEPBAL.cc,
fCmplxGEPBAL.h, floatGEPBAL.cc, and floatGEPBAL.h and making them
templates.
* liboctave/numeric/module.mk: Update.
* balance.cc, CMatrix.h, dMatrix.h, fCMatrix.h, fMatrix.h, mx-defs.h,
mx-ext.h: Use new classes.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 16 Feb 2016 00:32:29 -0500 |
| parents | f5b8c3aca5f8 |
| children | 3c8a3d35661a |
| files | libinterp/corefcn/balance.cc liboctave/array/CMatrix.h liboctave/array/dMatrix.h liboctave/array/fCMatrix.h liboctave/array/fMatrix.h liboctave/numeric/CmplxAEPBAL.cc liboctave/numeric/CmplxAEPBAL.h liboctave/numeric/CmplxGEPBAL.cc liboctave/numeric/CmplxGEPBAL.h liboctave/numeric/aepbalance.cc liboctave/numeric/aepbalance.h liboctave/numeric/base-aepbal.h liboctave/numeric/dbleAEPBAL.cc liboctave/numeric/dbleAEPBAL.h liboctave/numeric/dbleGEPBAL.cc liboctave/numeric/dbleGEPBAL.h liboctave/numeric/fCmplxAEPBAL.cc liboctave/numeric/fCmplxAEPBAL.h liboctave/numeric/fCmplxGEPBAL.cc liboctave/numeric/fCmplxGEPBAL.h liboctave/numeric/floatAEPBAL.cc liboctave/numeric/floatAEPBAL.h liboctave/numeric/floatGEPBAL.cc liboctave/numeric/floatGEPBAL.h liboctave/numeric/gepbalance.cc liboctave/numeric/gepbalance.h liboctave/numeric/module.mk liboctave/operators/mx-defs.h liboctave/operators/mx-ext.h |
| diffstat | 29 files changed, 930 insertions(+), 1612 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/balance.cc Mon Feb 15 22:33:45 2016 -0500 +++ b/libinterp/corefcn/balance.cc Tue Feb 16 00:32:29 2016 -0500 @@ -29,14 +29,12 @@ #include <string> -#include "CmplxAEPBAL.h" -#include "fCmplxAEPBAL.h" -#include "dbleAEPBAL.h" -#include "floatAEPBAL.h" -#include "CmplxGEPBAL.h" -#include "fCmplxGEPBAL.h" -#include "dbleGEPBAL.h" -#include "floatGEPBAL.h" +#include "CMatrix.h" +#include "aepbalance.h" +#include "dMatrix.h" +#include "fCMatrix.h" +#include "fMatrix.h" +#include "gepbalance.h" #include "quit.h" #include "defun.h" @@ -152,7 +150,7 @@ { if (complex_case) { - FloatComplexAEPBALANCE result (fcaa, noperm, noscal); + aepbalance<FloatComplexMatrix> result (fcaa, noperm, noscal); if (nargout == 0 || nargout == 1) retval = ovl (result.balanced_matrix ()); @@ -166,7 +164,7 @@ } else { - FloatAEPBALANCE result (faa, noperm, noscal); + aepbalance<FloatMatrix> result (faa, noperm, noscal); if (nargout == 0 || nargout == 1) retval = ovl (result.balanced_matrix ()); @@ -183,7 +181,7 @@ { if (complex_case) { - ComplexAEPBALANCE result (caa, noperm, noscal); + aepbalance<ComplexMatrix> result (caa, noperm, noscal); if (nargout == 0 || nargout == 1) retval = ovl (result.balanced_matrix ()); @@ -197,7 +195,7 @@ } else { - AEPBALANCE result (aa, noperm, noscal); + aepbalance<Matrix> result (aa, noperm, noscal); if (nargout == 0 || nargout == 1) retval = ovl (result.balanced_matrix ()); @@ -251,7 +249,7 @@ { if (complex_case) { - FloatComplexGEPBALANCE result (fcaa, fcbb, bal_job); + gepbalance<FloatComplexMatrix> result (fcaa, fcbb, bal_job); switch (nargout) { @@ -276,7 +274,7 @@ } else { - FloatGEPBALANCE result (faa, fbb, bal_job); + gepbalance<FloatMatrix> result (faa, fbb, bal_job); switch (nargout) { @@ -304,7 +302,7 @@ { if (complex_case) { - ComplexGEPBALANCE result (caa, cbb, bal_job); + gepbalance<ComplexMatrix> result (caa, cbb, bal_job); switch (nargout) { @@ -329,7 +327,7 @@ } else { - GEPBALANCE result (aa, bb, bal_job); + gepbalance<Matrix> result (aa, bb, bal_job); switch (nargout) {
--- a/liboctave/array/CMatrix.h Mon Feb 15 22:33:45 2016 -0500 +++ b/liboctave/array/CMatrix.h Tue Feb 16 00:32:29 2016 -0500 @@ -44,6 +44,12 @@ typedef ComplexColumnVector column_vector_type; typedef ComplexRowVector row_vector_type; + typedef ColumnVector real_column_vector_type; + typedef RowVector real_row_vector_type; + + typedef Matrix real_matrix_type; + typedef ComplexMatrix complex_matrix_type; + typedef void (*solve_singularity_handler) (double rcon); ComplexMatrix (void) : ComplexNDArray () { }
--- a/liboctave/array/dMatrix.h Mon Feb 15 22:33:45 2016 -0500 +++ b/liboctave/array/dMatrix.h Tue Feb 16 00:32:29 2016 -0500 @@ -43,6 +43,12 @@ typedef ColumnVector column_vector_type; typedef RowVector row_vector_type; + typedef ColumnVector real_column_vector_type; + typedef RowVector real_row_vector_type; + + typedef Matrix real_matrix_type; + typedef ComplexMatrix complex_matrix_type; + typedef void (*solve_singularity_handler) (double rcon); Matrix (void) : NDArray () { }
--- a/liboctave/array/fCMatrix.h Mon Feb 15 22:33:45 2016 -0500 +++ b/liboctave/array/fCMatrix.h Tue Feb 16 00:32:29 2016 -0500 @@ -44,6 +44,12 @@ typedef FloatComplexColumnVector column_vector_type; typedef FloatComplexRowVector row_vector_type; + typedef FloatColumnVector real_column_vector_type; + typedef FloatRowVector real_row_vector_type; + + typedef FloatMatrix real_matrix_type; + typedef FloatComplexMatrix complex_matrix_type; + typedef void (*solve_singularity_handler) (float rcon); FloatComplexMatrix (void) : FloatComplexNDArray () { }
--- a/liboctave/array/fMatrix.h Mon Feb 15 22:33:45 2016 -0500 +++ b/liboctave/array/fMatrix.h Tue Feb 16 00:32:29 2016 -0500 @@ -43,6 +43,12 @@ typedef FloatColumnVector column_vector_type; typedef FloatRowVector row_vector_type; + typedef FloatColumnVector real_column_vector_type; + typedef FloatRowVector real_row_vector_type; + + typedef FloatMatrix real_matrix_type; + typedef FloatComplexMatrix complex_matrix_type; + typedef void (*solve_singularity_handler) (float rcon); FloatMatrix (void) : FloatNDArray () { }
--- a/liboctave/numeric/CmplxAEPBAL.cc Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,102 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -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/>. - -*/ - -#ifdef HAVE_CONFIG_H -# include <config.h> -#endif - -#include <string> - -#include "CmplxAEPBAL.h" -#include "dMatrix.h" -#include "f77-fcn.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, Complex*, - const octave_idx_type&, octave_idx_type&, - octave_idx_type&, double*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (zgebak, ZGEBAK) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - const octave_idx_type&, const double*, - const octave_idx_type&, Complex*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); -} - -ComplexAEPBALANCE::ComplexAEPBALANCE (const ComplexMatrix& a, - bool noperm, bool noscal) - : base_aepbal<ComplexMatrix, ColumnVector> () -{ - octave_idx_type n = a.cols (); - - if (a.rows () != n) - (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); - - octave_idx_type info; - - scale = ColumnVector (n); - double *pscale = scale.fortran_vec (); - - balanced_mat = a; - Complex *p_balanced_mat = balanced_mat.fortran_vec (); - - job = noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B'); - - F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), - n, p_balanced_mat, n, ilo, ihi, - pscale, info - F77_CHAR_ARG_LEN (1))); -} - -ComplexMatrix -ComplexAEPBALANCE::balancing_matrix (void) const -{ - octave_idx_type n = balanced_mat.rows (); - ComplexMatrix balancing_mat (n, n, 0.0); - for (octave_idx_type i = 0; i < n; i++) - balancing_mat.elem (i, i) = 1.0; - - Complex *p_balancing_mat = balancing_mat.fortran_vec (); - const double *pscale = scale.fortran_vec (); - - octave_idx_type info; - - char side = 'R'; - - F77_XFCN (zgebak, ZGEBAK, (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 balancing_mat; -}
--- a/liboctave/numeric/CmplxAEPBAL.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2008-2009 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_CmplxAEPBAL_h) -#define octave_CmplxAEPBAL_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "base-aepbal.h" -#include "CMatrix.h" -#include "dColVector.h" - -class -OCTAVE_API -ComplexAEPBALANCE : public base_aepbal<ComplexMatrix, ColumnVector> -{ -public: - - ComplexAEPBALANCE (void) : base_aepbal<ComplexMatrix, ColumnVector> () { } - - ComplexAEPBALANCE (const ComplexMatrix& a, bool noperm = false, - bool noscal = false); - - ComplexAEPBALANCE (const ComplexAEPBALANCE& a) - : base_aepbal<ComplexMatrix, ColumnVector> (a) { } - - ComplexMatrix balancing_matrix (void) const; -}; - -#endif
--- a/liboctave/numeric/CmplxGEPBAL.cc Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,123 +0,0 @@ -/* - -Copyright (C) 1994-2015 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" -#include "oct-locbuf.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"); - - if (a.dims () != b.dims ()) - err_nonconformant ("ComplexGEPBALANCE", n, n, b.rows(), b.cols()); - - 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; -}
--- a/liboctave/numeric/CmplxGEPBAL.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,91 +0,0 @@ -/* - -Copyright (C) 1994-2015 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_CmplxGEPBAL_h) -#define octave_CmplxGEPBAL_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "CMatrix.h" -#include "dMatrix.h" - -class -OCTAVE_API -ComplexGEPBALANCE -{ -public: - - ComplexGEPBALANCE (void) - : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () - { } - - ComplexGEPBALANCE (const ComplexMatrix& a, const ComplexMatrix& b, - const std::string& balance_job) - : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () - { - 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/aepbalance.cc Tue Feb 16 00:32:29 2016 -0500 @@ -0,0 +1,297 @@ +/* + +Copyright (C) 1994-2015 John W. Eaton +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/>. + +*/ + +#ifdef HAVE_CONFIG_H +# include <config.h> +#endif + +#include <string> + +#include "CColVector.h" +#include "CMatrix.h" +#include "aepbalance.h" +#include "dColVector.h" +#include "dMatrix.h" +#include "f77-fcn.h" +#include "fCColVector.h" +#include "fCMatrix.h" +#include "fColVector.h" +#include "fMatrix.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type&, + octave_idx_type&, double*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + const octave_idx_type&, const double*, + const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + 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&, const float*, + const octave_idx_type&, float*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, Complex*, + const octave_idx_type&, octave_idx_type&, + octave_idx_type&, double*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zgebak, ZGEBAK) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + const octave_idx_type&, const double*, + const octave_idx_type&, Complex*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + 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&, const float*, + const octave_idx_type&, FloatComplex*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); +} + +static inline char +get_job (bool noperm, bool noscal) +{ + return noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B'); +} + +template <> +aepbalance<Matrix>::aepbalance (const Matrix& a, bool noperm, bool noscal) + : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal)) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + (*current_liboctave_error_handler) ("aepbalance: requires square matrix"); + + scale = ColumnVector (n); + + octave_idx_type info; + + F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n, + balanced_mat.fortran_vec (), n, ilo, ihi, + scale.fortran_vec (), info + F77_CHAR_ARG_LEN (1))); +} + +template <> +Matrix +aepbalance<Matrix>::balancing_matrix (void) const +{ + octave_idx_type n = balanced_mat.rows (); + Matrix balancing_mat (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + balancing_mat.elem (i ,i) = 1.0; + + octave_idx_type info; + + char side = 'R'; + + F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, scale.data (), n, + balancing_mat.fortran_vec (), n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return balancing_mat; +} + +template <> +aepbalance<FloatMatrix>::aepbalance (const FloatMatrix& a, bool noperm, + bool noscal) + : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal)) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + (*current_liboctave_error_handler) ("aepbalance: requires square matrix"); + + scale = FloatColumnVector (n); + + octave_idx_type info; + + F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n, + balanced_mat.fortran_vec (), n, ilo, ihi, + scale.fortran_vec (), info + F77_CHAR_ARG_LEN (1))); +} + +template <> +FloatMatrix +aepbalance<FloatMatrix>::balancing_matrix (void) const +{ + octave_idx_type n = balanced_mat.rows (); + FloatMatrix balancing_mat (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + balancing_mat.elem (i ,i) = 1.0; + + octave_idx_type info; + + char side = 'R'; + + F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, scale.data (), n, + balancing_mat.fortran_vec (), n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return balancing_mat; +} + +template <> +aepbalance<ComplexMatrix>::aepbalance (const ComplexMatrix& a, bool noperm, + bool noscal) + : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal)) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + (*current_liboctave_error_handler) ("aepbalance: requires square matrix"); + + scale = ColumnVector (n); + + octave_idx_type info; + + F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n, + balanced_mat.fortran_vec (), n, ilo, ihi, + scale.fortran_vec (), info + F77_CHAR_ARG_LEN (1))); +} + +template <> +ComplexMatrix +aepbalance<ComplexMatrix>::balancing_matrix (void) const +{ + octave_idx_type n = balanced_mat.rows (); + ComplexMatrix balancing_mat (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + balancing_mat.elem (i, i) = 1.0; + + octave_idx_type info; + + char side = 'R'; + + F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, scale.data (), n, + balancing_mat.fortran_vec (), n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return balancing_mat; +} + +template <> +aepbalance<FloatComplexMatrix>::aepbalance (const FloatComplexMatrix& a, + bool noperm, bool noscal) + : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal)) +{ + octave_idx_type n = a.cols (); + + if (a.rows () != n) + (*current_liboctave_error_handler) ("aepbalance: requires square matrix"); + + scale = FloatColumnVector (n); + + octave_idx_type info; + + F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n, + balanced_mat.fortran_vec (), n, ilo, ihi, + scale.fortran_vec (), info + F77_CHAR_ARG_LEN (1))); +} + +template <> +FloatComplexMatrix +aepbalance<FloatComplexMatrix>::balancing_matrix (void) const +{ + octave_idx_type n = balanced_mat.rows (); + FloatComplexMatrix balancing_mat (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + balancing_mat.elem (i, i) = 1.0; + + octave_idx_type info; + + char side = 'R'; + + F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, scale.data (), n, + balancing_mat.fortran_vec (), n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return balancing_mat; +} + +// Instantiations we need. + +template class aepbalance<Matrix>; + +template class aepbalance<FloatMatrix>; + +template class aepbalance<ComplexMatrix>; + +template class aepbalance<FloatComplexMatrix>;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/aepbalance.h Tue Feb 16 00:32:29 2016 -0500 @@ -0,0 +1,120 @@ +/* + +Copyright (C) 1994-2016 John W. Eaton +Copyright (C) 2008-2015 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_aepbalance_h) +#define octave_aepbalance_h 1 + +#include "octave-config.h" + +template <typename MT> +class aepbalance +{ +public: + + typedef typename MT::real_column_vector_type VT; + + aepbalance (void) : balanced_mat (), scale (), ilo (), ihi (), job () { } + + aepbalance (const MT& a, bool noperm = false, bool noscal = false); + + aepbalance (const aepbalance& a) + : balanced_mat (a.balanced_mat), scale (a.scale), + ilo(a.ilo), ihi(a.ihi), job(a.job) + { + } + + aepbalance& operator = (const aepbalance& a) + { + if (this != &a) + { + balanced_mat = a.balanced_mat; + scale = a.scale; + ilo = a.ilo; + ihi = a.ihi; + job = a.job; + } + + return *this; + } + + virtual ~aepbalance (void) { } + + MT balancing_matrix (void) const; + + MT balanced_matrix (void) const + { + return balanced_mat; + } + + VT permuting_vector (void) const + { + octave_idx_type n = balanced_mat.rows (); + + VT pv (n); + + for (octave_idx_type i = 0; i < n; i++) + pv(i) = i+1; + + for (octave_idx_type i = n-1; i >= ihi; i--) + { + octave_idx_type j = scale(i) - 1; + std::swap (pv(i), pv(j)); + } + + for (octave_idx_type i = 0; i < ilo-1; i++) + { + octave_idx_type j = scale(i) - 1; + std::swap (pv(i), pv(j)); + } + + return pv; + } + + VT scaling_vector (void) const + { + octave_idx_type n = balanced_mat.rows (); + + VT scv (n); + + for (octave_idx_type i = 0; i < ilo-1; i++) + scv(i) = 1; + + for (octave_idx_type i = ilo-1; i < ihi; i++) + scv(i) = scale(i); + + for (octave_idx_type i = ihi; i < n; i++) + scv(i) = 1; + + return scv; + } + +protected: + + MT balanced_mat; + VT scale; + octave_idx_type ilo; + octave_idx_type ihi; + char job; +}; + +#endif
--- a/liboctave/numeric/base-aepbal.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,96 +0,0 @@ -/* - -Copyright (C) 2008-2015 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_base_aepbal_h) -#define octave_base_aepbal_h 1 - -#include "octave-config.h" - -template <typename MatrixT, typename VectorT> -class base_aepbal -{ -protected: - MatrixT balanced_mat; - VectorT scale; - octave_idx_type ilo, ihi; - char job; - - base_aepbal (void) : balanced_mat (), scale (), ilo (), ihi (), job () { } - -public: - - base_aepbal (const base_aepbal& a) - : balanced_mat (a.balanced_mat), scale (a.scale), - ilo(a.ilo), ihi(a.ihi), job(a.job) - { - } - - base_aepbal& operator = (const base_aepbal& a) - { - balanced_mat = a.balanced_mat; - scale = a.scale; - ilo = a.ilo; - ihi = a.ihi; - job = a.job; - return *this; - } - - virtual ~base_aepbal (void) { } - - MatrixT balanced_matrix (void) const { return balanced_mat; } - - VectorT permuting_vector (void) const - { - octave_idx_type n = balanced_mat.rows (); - VectorT pv (n); - for (octave_idx_type i = 0; i < n; i++) - pv(i) = i+1; - for (octave_idx_type i = n-1; i >= ihi; i--) - { - octave_idx_type j = scale(i) - 1; - std::swap (pv(i), pv(j)); - } - for (octave_idx_type i = 0; i < ilo-1; i++) - { - octave_idx_type j = scale(i) - 1; - std::swap (pv(i), pv(j)); - } - - return pv; - } - - VectorT scaling_vector (void) const - { - octave_idx_type n = balanced_mat.rows (); - VectorT scv (n); - for (octave_idx_type i = 0; i < ilo-1; i++) - scv(i) = 1; - for (octave_idx_type i = ilo-1; i < ihi; i++) - scv(i) = scale(i); - for (octave_idx_type i = ihi; i < n; i++) - scv(i) = 1; - - return scv; - } -}; - -#endif
--- a/liboctave/numeric/dbleAEPBAL.cc Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,99 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -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/>. - -*/ - -#ifdef HAVE_CONFIG_H -# include <config.h> -#endif - -#include <string> - -#include "dbleAEPBAL.h" -#include "f77-fcn.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type&, - octave_idx_type&, double*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - const octave_idx_type&, const double*, - const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); -} - -AEPBALANCE::AEPBALANCE (const Matrix& a, bool noperm, bool noscal) - : base_aepbal<Matrix, ColumnVector> () -{ - octave_idx_type n = a.cols (); - - if (a.rows () != n) - (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); - - octave_idx_type info; - - scale = ColumnVector (n); - double *pscale = scale.fortran_vec (); - - balanced_mat = a; - double *p_balanced_mat = balanced_mat.fortran_vec (); - - job = noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B'); - - F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), - n, p_balanced_mat, n, ilo, ihi, pscale, info - F77_CHAR_ARG_LEN (1))); -} - -Matrix -AEPBALANCE::balancing_matrix (void) const -{ - octave_idx_type n = balanced_mat.rows (); - Matrix balancing_mat (n, n, 0.0); - for (octave_idx_type i = 0; i < n; i++) - balancing_mat.elem (i ,i) = 1.0; - - double *p_balancing_mat = balancing_mat.fortran_vec (); - const double *pscale = scale.fortran_vec (); - - octave_idx_type info; - - char side = 'R'; - - F77_XFCN (dgebak, DGEBAK, (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 balancing_mat; -}
--- a/liboctave/numeric/dbleAEPBAL.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2008-2009 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_dbleAEPBAL_h) -#define octave_dbleAEPBAL_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "base-aepbal.h" -#include "dMatrix.h" -#include "dColVector.h" - -class -OCTAVE_API -AEPBALANCE : public base_aepbal<Matrix, ColumnVector> -{ -public: - - AEPBALANCE (void) : base_aepbal<Matrix, ColumnVector> () { } - - AEPBALANCE (const Matrix& a, bool noperm = false, - bool noscal = false); - - AEPBALANCE (const AEPBALANCE& a) - : base_aepbal<Matrix, ColumnVector> (a) { } - - Matrix balancing_matrix (void) const; -}; - -#endif
--- a/liboctave/numeric/dbleGEPBAL.cc Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,122 +0,0 @@ -/* - -Copyright (C) 1994-2015 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" -#include "oct-locbuf.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"); - - if (a.dims () != b.dims ()) - err_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols()); - - 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; -}
--- a/liboctave/numeric/dbleGEPBAL.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ -/* - -Copyright (C) 1994-2015 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_dbleGEPBAL_h) -#define octave_dbleGEPBAL_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "dMatrix.h" - -class -OCTAVE_API -GEPBALANCE -{ -public: - - GEPBALANCE (void) - : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () - { } - - GEPBALANCE (const Matrix& a, const Matrix& b, const std::string& balance_job) - : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () - { - 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
--- a/liboctave/numeric/fCmplxAEPBAL.cc Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,102 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -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/>. - -*/ - -#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&, const float*, - const octave_idx_type&, FloatComplex*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); -} - -FloatComplexAEPBALANCE::FloatComplexAEPBALANCE (const FloatComplexMatrix& a, - bool noperm, bool noscal) - : base_aepbal<FloatComplexMatrix, FloatColumnVector> () -{ - octave_idx_type n = a.cols (); - - if (a.rows () != n) - (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); - - octave_idx_type info; - - scale = FloatColumnVector (n); - float *pscale = scale.fortran_vec (); - - balanced_mat = a; - FloatComplex *p_balanced_mat = balanced_mat.fortran_vec (); - - job = noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B'); - - F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), - n, p_balanced_mat, n, ilo, ihi, - pscale, info - F77_CHAR_ARG_LEN (1))); -} - -FloatComplexMatrix -FloatComplexAEPBALANCE::balancing_matrix (void) const -{ - octave_idx_type n = balanced_mat.rows (); - FloatComplexMatrix balancing_mat (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 (); - const float *pscale = scale.fortran_vec (); - - octave_idx_type info; - - 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 balancing_mat; -}
--- a/liboctave/numeric/fCmplxAEPBAL.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2008-2009 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_fCmplxAEPBAL_h) -#define octave_fCmplxAEPBAL_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "base-aepbal.h" -#include "fCMatrix.h" -#include "fColVector.h" - -class -OCTAVE_API -FloatComplexAEPBALANCE - : public base_aepbal<FloatComplexMatrix, FloatColumnVector> -{ -public: - - FloatComplexAEPBALANCE (void) - : base_aepbal<FloatComplexMatrix, FloatColumnVector> () { } - - FloatComplexAEPBALANCE (const FloatComplexMatrix& a, bool noperm = false, - bool noscal = false); - - FloatComplexAEPBALANCE (const FloatComplexAEPBALANCE& a) - : base_aepbal<FloatComplexMatrix, FloatColumnVector> (a) { } - - FloatComplexMatrix balancing_matrix (void) const; -}; - -#endif
--- a/liboctave/numeric/fCmplxGEPBAL.cc Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,126 +0,0 @@ -/* - -Copyright (C) 1994-2015 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" -#include "oct-locbuf.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 ()) - err_nonconformant ("FloatComplexGEPBALANCE", n, n, b.rows(), b.cols()); - - 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; -}
--- a/liboctave/numeric/fCmplxGEPBAL.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,94 +0,0 @@ -/* - -Copyright (C) 1994-2015 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_fCmplxGEPBAL_h) -#define octave_fCmplxGEPBAL_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "fMatrix.h" -#include "fCMatrix.h" - -class -OCTAVE_API -FloatComplexGEPBALANCE -{ -public: - - FloatComplexGEPBALANCE (void) - : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () - { } - - FloatComplexGEPBALANCE (const FloatComplexMatrix& a, - const FloatComplexMatrix& b, - const std::string& balance_job) - : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () - { - 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
--- a/liboctave/numeric/floatAEPBAL.cc Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,100 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -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/>. - -*/ - -#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&, const float*, - const octave_idx_type&, float*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); -} - -FloatAEPBALANCE::FloatAEPBALANCE (const FloatMatrix& a, - bool noperm, bool noscal) - : base_aepbal<FloatMatrix, FloatColumnVector> () -{ - octave_idx_type n = a.cols (); - - if (a.rows () != n) - (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); - - octave_idx_type info; - - scale = FloatColumnVector (n); - float *pscale = scale.fortran_vec (); - - balanced_mat = a; - float *p_balanced_mat = balanced_mat.fortran_vec (); - - job = noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B'); - - F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), - n, p_balanced_mat, n, ilo, ihi, pscale, info - F77_CHAR_ARG_LEN (1))); -} - -FloatMatrix -FloatAEPBALANCE::balancing_matrix (void) const -{ - octave_idx_type n = balanced_mat.rows (); - FloatMatrix balancing_mat (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 (); - const float *pscale = scale.fortran_vec (); - - octave_idx_type info; - - 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 balancing_mat; -}
--- a/liboctave/numeric/floatAEPBAL.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -/* - -Copyright (C) 1994-2015 John W. Eaton -Copyright (C) 2008-2009 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_floatAEPBAL_h) -#define octave_floatAEPBAL_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "base-aepbal.h" -#include "fMatrix.h" -#include "fColVector.h" - -class -OCTAVE_API -FloatAEPBALANCE : public base_aepbal<FloatMatrix, FloatColumnVector> -{ -public: - - FloatAEPBALANCE (void) : base_aepbal<FloatMatrix, FloatColumnVector> () { } - - FloatAEPBALANCE (const FloatMatrix& a, bool noperm = false, - bool noscal = false); - - FloatAEPBALANCE (const FloatAEPBALANCE& a) - : base_aepbal<FloatMatrix, FloatColumnVector> (a) { } - - FloatMatrix balancing_matrix (void) const; -}; - -#endif
--- a/liboctave/numeric/floatGEPBAL.cc Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,123 +0,0 @@ -/* - -Copyright (C) 1994-2015 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" -#include "oct-locbuf.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"); - - if (a.dims () != b.dims ()) - err_nonconformant ("FloatGEPBALANCE", n, n, b.rows(), b.cols()); - - 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; -}
--- a/liboctave/numeric/floatGEPBAL.h Mon Feb 15 22:33:45 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,88 +0,0 @@ -/* - -Copyright (C) 1994-2015 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_floatGEPBAL_h) -#define octave_floatGEPBAL_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "fMatrix.h" - -class -OCTAVE_API -FloatGEPBALANCE -{ -public: - - FloatGEPBALANCE (void) - : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () - { } - FloatGEPBALANCE (const FloatMatrix& a, const FloatMatrix& b, - const std::string& balance_job) - : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () - { - 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/gepbalance.cc Tue Feb 16 00:32:29 2016 -0500 @@ -0,0 +1,378 @@ +/* + +Copyright (C) 1994-2015 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 "Array-util.h" +#include "CMatrix.h" +#include "dMatrix.h" +#include "f77-fcn.h" +#include "fCMatrix.h" +#include "fMatrix.h" +#include "gepbalance.h" +#include "oct-locbuf.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); + + 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); + + 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 (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); +} + +template <> +octave_idx_type +gepbalance<Matrix>::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"); + + if (a.dims () != b.dims ()) + err_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols()); + + 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; +} + +template <> +octave_idx_type +gepbalance<FloatMatrix>::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"); + + if (a.dims () != b.dims ()) + err_nonconformant ("FloatGEPBALANCE", n, n, b.rows(), b.cols()); + + 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; +} + +template <> +octave_idx_type +gepbalance<ComplexMatrix>::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"); + + if (a.dims () != b.dims ()) + err_nonconformant ("ComplexGEPBALANCE", n, n, b.rows(), b.cols()); + + 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; +} + +template <> +octave_idx_type +gepbalance<FloatComplexMatrix>::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 ()) + err_nonconformant ("FloatComplexGEPBALANCE", n, n, b.rows(), b.cols()); + + 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; +} + +// Instantiations we need. + +template class gepbalance<Matrix>; + +template class gepbalance<FloatMatrix>; + +template class gepbalance<ComplexMatrix>; + +template class gepbalance<FloatComplexMatrix>;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/gepbalance.h Tue Feb 16 00:32:29 2016 -0500 @@ -0,0 +1,86 @@ +/* + +Copyright (C) 1994-2015 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 "octave-config.h" + +#include <string> + +template <typename T> +class +gepbalance +{ +public: + + typedef typename T::real_matrix_type RT; + + gepbalance (void) + : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () + { } + + gepbalance (const T& a, const T& b, const std::string& job) + : balanced_mat (), balanced_mat2 (), balancing_mat (), balancing_mat2 () + { + init (a, b, 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) { } + + T balanced_matrix (void) const { return balanced_mat; } + + T balanced_matrix2 (void) const { return balanced_mat2; } + + RT balancing_matrix (void) const { return balancing_mat; } + + RT balancing_matrix2 (void) const { return balancing_mat2; } + +private: + + T balanced_mat; + T balanced_mat2; + RT balancing_mat; + RT balancing_mat2; + + octave_idx_type init (const T& a, const T& b, const std::string& job); +}; + +#endif
--- a/liboctave/numeric/module.mk Mon Feb 15 22:33:45 2016 -0500 +++ b/liboctave/numeric/module.mk Tue Feb 16 00:32:29 2016 -0500 @@ -8,9 +8,7 @@ LIBOCTAVE_OPT_IN = $(LIBOCTAVE_OPT_INC:.h=.in) NUMERIC_INC = \ - liboctave/numeric/CmplxAEPBAL.h \ liboctave/numeric/CmplxCHOL.h \ - liboctave/numeric/CmplxGEPBAL.h \ liboctave/numeric/CmplxLU.h \ liboctave/numeric/CmplxQR.h \ liboctave/numeric/CmplxQRP.h \ @@ -31,7 +29,7 @@ liboctave/numeric/ODES.h \ liboctave/numeric/ODESFunc.h \ liboctave/numeric/Quad.h \ - liboctave/numeric/base-aepbal.h \ + liboctave/numeric/aepbalance.h \ liboctave/numeric/base-dae.h \ liboctave/numeric/base-de.h \ liboctave/numeric/base-lu.h \ @@ -39,29 +37,24 @@ liboctave/numeric/base-qr.h \ liboctave/numeric/bsxfun-decl.h \ liboctave/numeric/bsxfun.h \ - liboctave/numeric/dbleAEPBAL.h \ liboctave/numeric/dbleCHOL.h \ - liboctave/numeric/dbleGEPBAL.h \ liboctave/numeric/dbleLU.h \ liboctave/numeric/dbleQR.h \ liboctave/numeric/dbleQRP.h \ liboctave/numeric/dbleSVD.h \ liboctave/numeric/eigs-base.h \ - liboctave/numeric/fCmplxAEPBAL.h \ liboctave/numeric/fCmplxCHOL.h \ - liboctave/numeric/fCmplxGEPBAL.h \ liboctave/numeric/fCmplxLU.h \ liboctave/numeric/fCmplxQR.h \ liboctave/numeric/fCmplxQRP.h \ liboctave/numeric/fCmplxSVD.h \ liboctave/numeric/fEIG.h \ - liboctave/numeric/floatAEPBAL.h \ liboctave/numeric/floatCHOL.h \ - liboctave/numeric/floatGEPBAL.h \ liboctave/numeric/floatLU.h \ liboctave/numeric/floatQR.h \ liboctave/numeric/floatQRP.h \ liboctave/numeric/floatSVD.h \ + liboctave/numeric/gepbalance.h \ liboctave/numeric/hess.h \ liboctave/numeric/lo-mappers.h \ liboctave/numeric/lo-specfun.h \ @@ -85,9 +78,7 @@ liboctave/numeric/randpoisson.c NUMERIC_SRC = \ - liboctave/numeric/CmplxAEPBAL.cc \ liboctave/numeric/CmplxCHOL.cc \ - liboctave/numeric/CmplxGEPBAL.cc \ liboctave/numeric/CmplxLU.cc \ liboctave/numeric/CmplxQR.cc \ liboctave/numeric/CmplxQRP.cc \ @@ -100,29 +91,25 @@ liboctave/numeric/LSODE.cc \ liboctave/numeric/ODES.cc \ liboctave/numeric/Quad.cc \ - liboctave/numeric/dbleAEPBAL.cc \ + liboctave/numeric/aepbalance.cc \ liboctave/numeric/dbleCHOL.cc \ - liboctave/numeric/dbleGEPBAL.cc \ liboctave/numeric/dbleLU.cc \ liboctave/numeric/dbleQR.cc \ liboctave/numeric/dbleQRP.cc \ liboctave/numeric/dbleSVD.cc \ liboctave/numeric/eigs-base.cc \ - liboctave/numeric/fCmplxAEPBAL.cc \ liboctave/numeric/fCmplxCHOL.cc \ - liboctave/numeric/fCmplxGEPBAL.cc \ liboctave/numeric/fCmplxLU.cc \ liboctave/numeric/fCmplxQR.cc \ liboctave/numeric/fCmplxQRP.cc \ liboctave/numeric/fCmplxSVD.cc \ liboctave/numeric/fEIG.cc \ - liboctave/numeric/floatAEPBAL.cc \ liboctave/numeric/floatCHOL.cc \ - liboctave/numeric/floatGEPBAL.cc \ liboctave/numeric/floatLU.cc \ liboctave/numeric/floatQR.cc \ liboctave/numeric/floatQRP.cc \ liboctave/numeric/floatSVD.cc \ + liboctave/numeric/gepbalance.cc \ liboctave/numeric/hess.cc \ liboctave/numeric/lo-mappers.cc \ liboctave/numeric/lo-specfun.cc \
--- a/liboctave/operators/mx-defs.h Mon Feb 15 22:33:45 2016 -0500 +++ b/liboctave/operators/mx-defs.h Tue Feb 16 00:32:29 2016 -0500 @@ -58,15 +58,9 @@ class PermMatrix; -class AEPBALANCE; -class ComplexAEPBALANCE; -class FloatAEPBALANCE; -class FloatComplexAEPBALANCE; +template <typename T> class aepbalance; -class GEPBALANCE; -class ComplexGEPBALANCE; -class FloatGEPBALANCE; -class FloatComplexGEPBALANCE; +template <typename T> class gepbalance; class CHOL; class ComplexCHOL;
--- a/liboctave/operators/mx-ext.h Mon Feb 15 22:33:45 2016 -0500 +++ b/liboctave/operators/mx-ext.h Tue Feb 16 00:32:29 2016 -0500 @@ -27,8 +27,11 @@ // Result of a AEP Balance operation. -#include "dbleAEPBAL.h" -#include "CmplxAEPBAL.h" +#include "aepbalance.h" + +// Result of a GEP Balance operation. + +#include "gepbalance.h" // Result of a Determinant calculation.