Mercurial > octave
diff liboctave/numeric/aepbalance.cc @ 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 | liboctave/numeric/dbleAEPBAL.cc@f7121e111991 |
| children | 40de9f8f23a6 |
line wrap: on
line diff
--- /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>;