Mercurial > octave-nkf
changeset 8386:a5e080076778
make balance more Matlab-compatible
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Tue, 09 Dec 2008 10:08:19 +0100 |
| parents | 6e9660cd3bf2 |
| children | 1567db1e166c |
| files | liboctave/ChangeLog liboctave/CmplxAEPBAL.cc liboctave/CmplxAEPBAL.h liboctave/Makefile.in liboctave/base-aepbal.h liboctave/dbleAEPBAL.cc liboctave/dbleAEPBAL.h liboctave/fCmplxAEPBAL.cc liboctave/fCmplxAEPBAL.h liboctave/floatAEPBAL.cc liboctave/floatAEPBAL.h src/ChangeLog src/DLD-FUNCTIONS/balance.cc |
| diffstat | 13 files changed, 277 insertions(+), 213 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/ChangeLog Tue Dec 09 10:08:19 2008 +0100 @@ -1,3 +1,11 @@ +2008-12-09 Jaroslav Hajek <highegg@gmail.com> + + * base-aepbal.h: New source. + * dbleAEPBAL.h, dbleAEPBAL.cc: Rebase AEPBAL on base_aepbal. + * floatAEPBAL.h, floatAEPBAL.cc: Rebase FloatAEPBAL on base_aepbal. + * CmplxAEPBAL.h, CmplxAEPBAL.cc: Rebase ComplexAEPBAL on base_aepbal. + * fCmplxAEPBAL.h, fCmplxAEPBAL.cc: Rebase FloatComplexAEPBAL on base_aepbal. + 2008-12-08 Jaroslav Hajek <highegg@gmail.com> * idx-vector.cc (idx_vector::idx_vector_rep::idx_vector_rep (const
--- a/liboctave/CmplxAEPBAL.cc Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/CmplxAEPBAL.cc Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -42,46 +43,52 @@ 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&, double*, + 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); } -octave_idx_type -ComplexAEPBALANCE::init (const ComplexMatrix& a, - const std::string& balance_job) +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"); - return -1; + return; } octave_idx_type info; - octave_idx_type ilo; - octave_idx_type ihi; - Array<double> scale (n); + scale = ColumnVector (n); double *pscale = scale.fortran_vec (); balanced_mat = a; Complex *p_balanced_mat = balanced_mat.fortran_vec (); - char job = balance_job[0]; + 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))); +} - balancing_mat = ComplexMatrix (n, n, 0.0); +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'; @@ -92,7 +99,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - return info; + return balancing_mat; } /*
--- a/liboctave/CmplxAEPBAL.h Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/CmplxAEPBAL.h Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -27,48 +28,25 @@ #include <iostream> #include <string> +#include "base-aepbal.h" #include "CMatrix.h" +#include "dColVector.h" class OCTAVE_API -ComplexAEPBALANCE +ComplexAEPBALANCE : public base_aepbal<ComplexMatrix, ColumnVector> { public: - ComplexAEPBALANCE (void) : balanced_mat (), balancing_mat () { } + ComplexAEPBALANCE (void) : base_aepbal<ComplexMatrix, ColumnVector> () { } - ComplexAEPBALANCE (const ComplexMatrix& a, const std::string& balance_job) - { - init (a, balance_job); - } - - ComplexAEPBALANCE (const ComplexAEPBALANCE& a) - : balanced_mat (a.balanced_mat), balancing_mat (a.balancing_mat) { } + ComplexAEPBALANCE (const ComplexMatrix& a, bool noperm = false, + bool noscal = false); - ComplexAEPBALANCE& operator = (const ComplexAEPBALANCE& a) - { - if (this != &a) - { - balanced_mat = a.balanced_mat; - balancing_mat = a.balancing_mat; - } - return *this; - } - - ~ComplexAEPBALANCE (void) { } + ComplexAEPBALANCE (const ComplexAEPBALANCE& a) + : base_aepbal<ComplexMatrix, ColumnVector> (a) { } - ComplexMatrix balanced_matrix (void) const { return balanced_mat; } - - ComplexMatrix balancing_matrix (void) const { return balancing_mat; } - - friend std::ostream& operator << (std::ostream& os, const ComplexAEPBALANCE& a); - -private: - - ComplexMatrix balanced_mat; - ComplexMatrix balancing_mat; - - octave_idx_type init (const ComplexMatrix& a, const std::string& balance_job); + ComplexMatrix balancing_matrix (void) const; }; #endif
--- a/liboctave/Makefile.in Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/Makefile.in Tue Dec 09 10:08:19 2008 +0100 @@ -43,7 +43,7 @@ MATRIX_INC := Array.h Array2.h Array3.h ArrayN.h DiagArray2.h \ Array-util.h ArrayN-idx.h MArray-defs.h \ MArray.h MArray2.h MDiagArray2.h Matrix.h MArrayN.h \ - base-lu.h dim-vector.h mx-base.h mx-op-defs.h \ + base-lu.h base-aepbal.h dim-vector.h mx-base.h mx-op-defs.h \ mx-defs.h mx-ext.h CColVector.h CDiagMatrix.h CMatrix.h \ CNDArray.h CRowVector.h CmplxAEPBAL.h CmplxCHOL.h \ CmplxGEPBAL.h CmplxHESS.h CmplxLU.h CmplxQR.h CmplxQRP.h \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/base-aepbal.h Tue Dec 09 10:08:19 2008 +0100 @@ -0,0 +1,96 @@ +/* + +Copyright (C) 2008 Jaroslav Hajek + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#if !defined (octave_base_aepbal_h) +#define octave_base_aepbal_h 1 + +#include "oct-types.h" + +template <class MatrixT, class 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; + } + + 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; + octave_idx_type k = pv(j); + pv(j) = pv(i); + pv(i) = k; + } + for (octave_idx_type i = ilo-2; i >= 0; i--) + { + octave_idx_type j = scale(i) - 1; + octave_idx_type k = pv(j); + pv(j) = pv(i); + pv(i) = k; + } + + 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/dbleAEPBAL.cc Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/dbleAEPBAL.cc Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -41,44 +42,51 @@ 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&, double*, - const octave_idx_type&, double*, const octave_idx_type&, octave_idx_type& + 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); } -octave_idx_type -AEPBALANCE::init (const Matrix& a, const std::string& balance_job) +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"); - return -1; + return; } octave_idx_type info; - octave_idx_type ilo; - octave_idx_type ihi; - Array<double> scale (n); + scale = ColumnVector (n); double *pscale = scale.fortran_vec (); balanced_mat = a; double *p_balanced_mat = balanced_mat.fortran_vec (); - char job = balance_job[0]; + 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))); +} - balancing_mat = Matrix (n, n, 0.0); +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'; @@ -89,7 +97,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - return info; + return balancing_mat; } /*
--- a/liboctave/dbleAEPBAL.h Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/dbleAEPBAL.h Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -27,48 +28,25 @@ #include <iostream> #include <string> +#include "base-aepbal.h" #include "dMatrix.h" +#include "dColVector.h" class OCTAVE_API -AEPBALANCE +AEPBALANCE : public base_aepbal<Matrix, ColumnVector> { public: - AEPBALANCE (void) : balanced_mat (), balancing_mat () { } + AEPBALANCE (void) : base_aepbal<Matrix, ColumnVector> () { } - AEPBALANCE (const Matrix& a,const std::string& balance_job) - { - init (a, balance_job); - } - - AEPBALANCE (const AEPBALANCE& a) - : balanced_mat (a.balanced_mat), balancing_mat (a.balancing_mat) { } + AEPBALANCE (const Matrix& a, bool noperm = false, + bool noscal = false); - AEPBALANCE& operator = (const AEPBALANCE& a) - { - if (this != &a) - { - balanced_mat = a.balanced_mat; - balancing_mat = a.balancing_mat; - } - return *this; - } - - ~AEPBALANCE (void) { } + AEPBALANCE (const AEPBALANCE& a) + : base_aepbal<Matrix, ColumnVector> (a) { } - Matrix balanced_matrix (void) const { return balanced_mat; } - - Matrix balancing_matrix (void) const { return balancing_mat; } - - friend std::ostream& operator << (std::ostream& os, const AEPBALANCE& a); - -private: - - Matrix balanced_mat; - Matrix balancing_mat; - - octave_idx_type init (const Matrix& a, const std::string& balance_job); + Matrix balancing_matrix (void) const; }; #endif
--- a/liboctave/fCmplxAEPBAL.cc Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/fCmplxAEPBAL.cc Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -42,46 +43,52 @@ 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&, const 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) +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"); - return -1; + return; } octave_idx_type info; - octave_idx_type ilo; - octave_idx_type ihi; - Array<float> scale (n); + scale = FloatColumnVector (n); float *pscale = scale.fortran_vec (); balanced_mat = a; FloatComplex *p_balanced_mat = balanced_mat.fortran_vec (); - char job = balance_job[0]; + 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))); +} - balancing_mat = FloatComplexMatrix (n, n, 0.0); +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'; @@ -92,7 +99,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - return info; + return balancing_mat; } /*
--- a/liboctave/fCmplxAEPBAL.h Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/fCmplxAEPBAL.h Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -27,48 +28,25 @@ #include <iostream> #include <string> +#include "base-aepbal.h" #include "fCMatrix.h" +#include "fColVector.h" class OCTAVE_API -FloatComplexAEPBALANCE +FloatComplexAEPBALANCE : public base_aepbal<FloatComplexMatrix, FloatColumnVector> { public: - FloatComplexAEPBALANCE (void) : balanced_mat (), balancing_mat () { } + FloatComplexAEPBALANCE (void) : base_aepbal<FloatComplexMatrix, FloatColumnVector> () { } - 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 (const FloatComplexMatrix& a, bool noperm = false, + bool noscal = false); - FloatComplexAEPBALANCE& operator = (const FloatComplexAEPBALANCE& a) - { - if (this != &a) - { - balanced_mat = a.balanced_mat; - balancing_mat = a.balancing_mat; - } - return *this; - } - - ~FloatComplexAEPBALANCE (void) { } + FloatComplexAEPBALANCE (const FloatComplexAEPBALANCE& a) + : base_aepbal<FloatComplexMatrix, FloatColumnVector> (a) { } - 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); + FloatComplexMatrix balancing_matrix (void) const; }; #endif
--- a/liboctave/floatAEPBAL.cc Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/floatAEPBAL.cc Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -41,44 +42,52 @@ 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& + 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); } -octave_idx_type -FloatAEPBALANCE::init (const FloatMatrix& a, const std::string& balance_job) +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"); - return -1; + return; } octave_idx_type info; - octave_idx_type ilo; - octave_idx_type ihi; - Array<float> scale (n); + scale = FloatColumnVector (n); float *pscale = scale.fortran_vec (); balanced_mat = a; float *p_balanced_mat = balanced_mat.fortran_vec (); - char job = balance_job[0]; + 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))); +} - balancing_mat = FloatMatrix (n, n, 0.0); +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'; @@ -89,7 +98,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - return info; + return balancing_mat; } /*
--- a/liboctave/floatAEPBAL.h Tue Dec 09 03:25:31 2008 +0100 +++ b/liboctave/floatAEPBAL.h Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -27,48 +28,25 @@ #include <iostream> #include <string> +#include "base-aepbal.h" #include "fMatrix.h" +#include "fColVector.h" class OCTAVE_API -FloatAEPBALANCE +FloatAEPBALANCE : public base_aepbal<FloatMatrix, FloatColumnVector> { public: - FloatAEPBALANCE (void) : balanced_mat (), balancing_mat () { } + FloatAEPBALANCE (void) : base_aepbal<FloatMatrix, FloatColumnVector> () { } - 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 (const FloatMatrix& a, bool noperm = false, + bool noscal = false); - FloatAEPBALANCE& operator = (const FloatAEPBALANCE& a) - { - if (this != &a) - { - balanced_mat = a.balanced_mat; - balancing_mat = a.balancing_mat; - } - return *this; - } - - ~FloatAEPBALANCE (void) { } + FloatAEPBALANCE (const FloatAEPBALANCE& a) + : base_aepbal<FloatMatrix, FloatColumnVector> (a) { } - 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); + FloatMatrix balancing_matrix (void) const; }; #endif
--- a/src/ChangeLog Tue Dec 09 03:25:31 2008 +0100 +++ b/src/ChangeLog Tue Dec 09 10:08:19 2008 +0100 @@ -1,3 +1,7 @@ +2008-12-09 Jaroslav Hajek <highegg@gmail.com> + + * DLD-FUNCTIONS/balance.cc (Fbalance): Exploit the new AEPBAL functionality. + 2008-12-08 Jaroslav Hajek <highegg@gmail.com> * xpow.cc ( xpow (const DiagMatrix& a, double b),
--- a/src/DLD-FUNCTIONS/balance.cc Tue Dec 09 03:25:31 2008 +0100 +++ b/src/DLD-FUNCTIONS/balance.cc Tue Dec 09 10:08:19 2008 +0100 @@ -2,6 +2,7 @@ Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006, 2007 John W. Eaton +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -50,6 +51,7 @@ "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{aa} =} balance (@var{a}, @var{opt})\n\ @deftypefnx {Loadable Function} {[@var{dd}, @var{aa}] =} balance (@var{a}, @var{opt})\n\ +@deftypefnx {Loadable Function} {[@var{p}, @var{d}, @var{aa}] =} balance (@var{a}, @var{opt})\n\ @deftypefnx {Loadable Function} {[@var{cc}, @var{dd}, @var{aa}, @var{bb}] =} balance (@var{a}, @var{b}, @var{opt})\n\ \n\ Compute @code{aa = dd \\ a * dd} in which @code{aa} is a matrix whose\n\ @@ -59,6 +61,11 @@ the equilibration to be computed without roundoff. Results of\n\ eigenvalue calculation are typically improved by balancing first.\n\ \n\ +If two output values are requested, @code{balance} returns \n\ +the permutation @code{p} and the diagonal @code{d} separately as vectors. \n\ +In this case, @code{dd = eye(n)(p,:) * diag (d)}, where @code{n} is the matrix \n\ +size. \n\ +\n\ If four output values are requested, compute @code{aa = cc*a*dd} and\n\ @code{bb = cc*b*dd)}, in which @code{aa} and @code{bb} have non-zero\n\ elements of approximately the same magnitude and @code{cc} and @code{dd}\n\ @@ -68,19 +75,11 @@ The eigenvalue balancing option @code{opt} may be one of:\n\ \n\ @table @asis\n\ -@item @code{\"N\"}, @code{\"n\"}\n\ -No balancing; arguments copied, transformation(s) set to identity.\n\ -\n\ -@item @code{\"P\"}, @code{\"p\"}\n\ -Permute argument(s) to isolate eigenvalues where possible.\n\ +@item @code{\"noperm\"}, @code{\"S\"}\n\ +Scale only; do not permute.\n\ \n\ -@item @code{\"S\"}, @code{\"s\"}\n\ -Scale to improve accuracy of computed eigenvalues.\n\ -\n\ -@item @code{\"B\"}, @code{\"b\"}\n\ -Permute and scale, in that order. Rows/columns of a (and b)\n\ -that are isolated by permutation are not scaled. This is the default\n\ -behavior.\n\ +@item @code{\"noscal\"}, @code{\"P\"}\n\ +Permute only; do not scale.\n\ @end table\n\ \n\ Algebraic eigenvalue balancing uses standard @sc{Lapack} routines.\n\ @@ -100,20 +99,11 @@ } // determine if it's AEP or GEP - int AEPcase = nargin == 1 ? 1 : args(1).is_string (); - std::string bal_job; + bool AEPcase = nargin == 1 || args(1).is_string (); // problem dimension octave_idx_type nn = args(0).rows (); - octave_idx_type arg_is_empty = empty_arg ("balance", nn, args(0).columns()); - - if (arg_is_empty < 0) - return retval; - - if (arg_is_empty > 0) - return octave_value_list (2, Matrix ()); - if (nn != args(0).columns()) { gripe_square_matrix_required ("balance"); @@ -154,75 +144,98 @@ if (AEPcase) { // Algebraic eigenvalue problem. - - if (nargin == 1) - bal_job = "B"; - else if (args(1).is_string ()) - bal_job = args(1).string_value (); - else - { - error ("balance: AEP argument 2 must be a string"); - return retval; - } + bool noperm = false, noscal = false; + if (nargin > 1) + { + std::string a1s = args(1).string_value (); + noperm = a1s == "noperm" || a1s == "S"; + noscal = a1s == "noscal" || a1s == "P"; + } // balance the AEP if (isfloat) { if (complex_case) { - FloatComplexAEPBALANCE result (fcaa, bal_job); + FloatComplexAEPBALANCE result (fcaa, noperm, noscal); if (nargout == 0 || nargout == 1) retval(0) = result.balanced_matrix (); - else + else if (nargout == 2) { retval(1) = result.balanced_matrix (); retval(0) = result.balancing_matrix (); } + else + { + retval(2) = result.balanced_matrix (); + retval(1) = result.scaling_vector (); + retval(0) = result.permuting_vector (); + } + } else { - FloatAEPBALANCE result (faa, bal_job); + FloatAEPBALANCE result (faa, noperm, noscal); if (nargout == 0 || nargout == 1) retval(0) = result.balanced_matrix (); - else + else if (nargout == 2) { retval(1) = result.balanced_matrix (); retval(0) = result.balancing_matrix (); } + else + { + retval(2) = result.balanced_matrix (); + retval(1) = result.scaling_vector (); + retval(0) = result.permuting_vector (); + } } } else { if (complex_case) { - ComplexAEPBALANCE result (caa, bal_job); + ComplexAEPBALANCE result (caa, noperm, noscal); if (nargout == 0 || nargout == 1) retval(0) = result.balanced_matrix (); - else + else if (nargout == 2) { retval(1) = result.balanced_matrix (); retval(0) = result.balancing_matrix (); } + else + { + retval(2) = result.balanced_matrix (); + retval(1) = result.scaling_vector (); + retval(0) = result.permuting_vector (); + } } else { - AEPBALANCE result (aa, bal_job); + AEPBALANCE result (aa, noperm, noscal); if (nargout == 0 || nargout == 1) retval(0) = result.balanced_matrix (); - else + else if (nargout == 2) { retval(1) = result.balanced_matrix (); retval(0) = result.balancing_matrix (); } + else + { + retval(2) = result.balanced_matrix (); + retval(1) = result.scaling_vector (); + retval(0) = result.permuting_vector (); + } } } } else { + std::string bal_job; if (nargout == 1) warning ("balance: used GEP, should have two output arguments");