Mercurial > octave
changeset 21266:e69eaee28737
make better use of templates for Schur decomposition
* liboctave/numeric/schur.h, liboctave/numeric/schur.cc:
New files generated from SCHUR.h, SCHUR.cc, CmplxSCHUR.h,
CmplxSCHUR.cc, dbleSCHUR.h, dbleSCHUR.cc, fCmplxSCHUR.h,
fCmplxSCHUR.cc, floatSCHUR.h, and floatSCHUR.cc and making them
templates.
* liboctave/numeric/module.mk: Update.
* libinterp/corefcn/schur.cc, sqrtm.cc, CMatrix.cc, dMatrix.cc,
fCMatrix.cc, fMatrix.cc, mx-defs.h, mx-ext.h: Use new template classes
and header file.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 15 Feb 2016 20:06:12 -0500 |
| parents | f780d057a3ec |
| children | f5b8c3aca5f8 |
| files | libinterp/corefcn/schur.cc libinterp/corefcn/sqrtm.cc liboctave/array/CMatrix.cc liboctave/array/dMatrix.cc liboctave/array/fCMatrix.cc liboctave/array/fMatrix.cc liboctave/numeric/CmplxSCHUR.cc liboctave/numeric/CmplxSCHUR.h liboctave/numeric/dbleSCHUR.cc liboctave/numeric/dbleSCHUR.h liboctave/numeric/fCmplxSCHUR.cc liboctave/numeric/fCmplxSCHUR.h liboctave/numeric/floatSCHUR.cc liboctave/numeric/floatSCHUR.h liboctave/numeric/module.mk liboctave/numeric/schur.cc liboctave/numeric/schur.h liboctave/operators/mx-defs.h liboctave/operators/mx-ext.h |
| diffstat | 19 files changed, 713 insertions(+), 1098 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/schur.cc Thu Jan 14 19:59:52 2016 +1100 +++ b/libinterp/corefcn/schur.cc Mon Feb 15 20:06:12 2016 -0500 @@ -26,10 +26,7 @@ #include <string> -#include "CmplxSCHUR.h" -#include "dbleSCHUR.h" -#include "fCmplxSCHUR.h" -#include "floatSCHUR.h" +#include "schur.h" #include "defun.h" #include "error.h" @@ -182,12 +179,12 @@ if (nargout <= 1) { - FloatSCHUR result (tmp, ord, false); + schur<FloatMatrix> result (tmp, ord, false); retval = ovl (result.schur_matrix ()); } else { - FloatSCHUR result (tmp, ord, true); + schur<FloatMatrix> result (tmp, ord, true); retval = ovl (result.unitary_matrix (), result.schur_matrix ()); } @@ -198,12 +195,12 @@ if (nargout <= 1) { - FloatComplexSCHUR result (ctmp, ord, false); + schur<FloatComplexMatrix> result (ctmp, ord, false); retval = ovl (mark_upper_triangular (result.schur_matrix ())); } else { - FloatComplexSCHUR result (ctmp, ord, true); + schur<FloatComplexMatrix> result (ctmp, ord, true); retval = ovl (result.unitary_matrix (), mark_upper_triangular (result.schur_matrix ())); } @@ -217,12 +214,12 @@ if (nargout <= 1) { - SCHUR result (tmp, ord, false); + schur<Matrix> result (tmp, ord, false); retval = ovl (result.schur_matrix ()); } else { - SCHUR result (tmp, ord, true); + schur<Matrix> result (tmp, ord, true); retval = ovl (result.unitary_matrix (), result.schur_matrix ()); } @@ -233,12 +230,12 @@ if (nargout <= 1) { - ComplexSCHUR result (ctmp, ord, false); + schur<ComplexMatrix> result (ctmp, ord, false); retval = ovl (mark_upper_triangular (result.schur_matrix ())); } else { - ComplexSCHUR result (ctmp, ord, true); + schur<ComplexMatrix> result (ctmp, ord, true); retval = ovl (result.unitary_matrix (), mark_upper_triangular (result.schur_matrix ())); } @@ -304,7 +301,8 @@ FloatMatrix u = args(0).float_matrix_value (); FloatMatrix t = args(1).float_matrix_value (); - FloatComplexSCHUR cs (FloatSCHUR (t, u)); + schur<FloatComplexMatrix> cs + = rsf2csf<FloatComplexMatrix, FloatMatrix> (t, u); return ovl (cs.unitary_matrix (), cs.schur_matrix ()); } @@ -313,7 +311,7 @@ Matrix u = args(0).matrix_value (); Matrix t = args(1).matrix_value (); - ComplexSCHUR cs (SCHUR (t, u)); + schur<ComplexMatrix> cs = rsf2csf<ComplexMatrix, Matrix> (t, u); return ovl (cs.unitary_matrix (), cs.schur_matrix ()); }
--- a/libinterp/corefcn/sqrtm.cc Thu Jan 14 19:59:52 2016 +1100 +++ b/libinterp/corefcn/sqrtm.cc Mon Feb 15 20:06:12 2016 -0500 @@ -27,8 +27,7 @@ #include <float.h> -#include "CmplxSCHUR.h" -#include "fCmplxSCHUR.h" +#include "schur.h" #include "lo-ieee.h" #include "lo-mappers.h" #include "oct-norm.h" @@ -177,9 +176,9 @@ do { - ComplexSCHUR schur (x, "", true); - x = schur.schur_matrix (); - u = schur.unitary_matrix (); + ComplexSCHUR schur_fact (x, "", true); + x = schur_fact.schur_matrix (); + u = schur_fact.unitary_matrix (); } while (0); // schur no longer needed. @@ -236,10 +235,10 @@ // sqrtm of a diagonal matrix is just sqrt. retval(0) = arg.sqrt (); else if (arg.is_single_type ()) - retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix, FloatComplexSCHUR> - (arg); + retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix, + schur<FloatComplexMatrix> > (arg); else if (arg.is_numeric_type ()) - retval(0) = do_sqrtm<Matrix, ComplexMatrix, ComplexSCHUR> (arg); + retval(0) = do_sqrtm<Matrix, ComplexMatrix, schur<ComplexMatrix> > (arg); if (nargout > 1) {
--- a/liboctave/array/CMatrix.cc Thu Jan 14 19:59:52 2016 +1100 +++ b/liboctave/array/CMatrix.cc Mon Feb 15 20:06:12 2016 -0500 @@ -46,7 +46,7 @@ #include "CDiagMatrix.h" #include "dDiagMatrix.h" #include "CmplxCHOL.h" -#include "CmplxSCHUR.h" +#include "schur.h" #include "CmplxSVD.h" #include "DET.h" #include "f77-fcn.h" @@ -3487,8 +3487,8 @@ // Compute Schur decompositions - ComplexSCHUR as (a, "U"); - ComplexSCHUR bs (b, "U"); + schur<ComplexMatrix> as (a, "U"); + schur<ComplexMatrix> bs (b, "U"); // Transform c to new coordinates.
--- a/liboctave/array/dMatrix.cc Thu Jan 14 19:59:52 2016 +1100 +++ b/liboctave/array/dMatrix.cc Mon Feb 15 20:06:12 2016 -0500 @@ -44,7 +44,7 @@ #include "CColVector.h" #include "PermMatrix.h" #include "DET.h" -#include "dbleSCHUR.h" +#include "schur.h" #include "dbleSVD.h" #include "dbleCHOL.h" #include "f77-fcn.h" @@ -2953,8 +2953,8 @@ // Compute Schur decompositions. - SCHUR as (a, "U"); - SCHUR bs (b, "U"); + schur<Matrix> as (a, "U"); + schur<Matrix> bs (b, "U"); // Transform c to new coordinates.
--- a/liboctave/array/fCMatrix.cc Thu Jan 14 19:59:52 2016 +1100 +++ b/liboctave/array/fCMatrix.cc Mon Feb 15 20:06:12 2016 -0500 @@ -46,7 +46,7 @@ #include "fCColVector.h" #include "fCRowVector.h" #include "fCmplxCHOL.h" -#include "fCmplxSCHUR.h" +#include "schur.h" #include "fCmplxSVD.h" #include "functor.h" #include "lo-error.h" @@ -3501,8 +3501,8 @@ // Compute Schur decompositions - FloatComplexSCHUR as (a, "U"); - FloatComplexSCHUR bs (b, "U"); + schur<FloatComplexMatrix> as (a, "U"); + schur<FloatComplexMatrix> bs (b, "U"); // Transform c to new coordinates.
--- a/liboctave/array/fMatrix.cc Thu Jan 14 19:59:52 2016 +1100 +++ b/liboctave/array/fMatrix.cc Mon Feb 15 20:06:12 2016 -0500 @@ -48,7 +48,7 @@ #include "f77-fcn.h" #include "fMatrix.h" #include "floatCHOL.h" -#include "floatSCHUR.h" +#include "schur.h" #include "floatSVD.h" #include "functor.h" #include "lo-error.h" @@ -2970,8 +2970,8 @@ // Compute Schur decompositions. - FloatSCHUR as (a, "U"); - FloatSCHUR bs (b, "U"); + schur<FloatMatrix> as (a, "U"); + schur<FloatMatrix> bs (b, "U"); // Transform c to new coordinates.
--- a/liboctave/numeric/CmplxSCHUR.cc Thu Jan 14 19:59:52 2016 +1100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,170 +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 "CmplxSCHUR.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (zgeesx, ZGEESX) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - ComplexSCHUR::select_function, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, Complex*, - const octave_idx_type&, octave_idx_type&, - Complex*, Complex*, const octave_idx_type&, - double&, double&, Complex*, - const octave_idx_type&, double*, - octave_idx_type*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (zrsf2csf, ZRSF2CSF) (const octave_idx_type&, Complex *, - Complex *, double *, double *); -} - -static octave_idx_type -select_ana (const Complex& a) -{ - return a.real () < 0.0; -} - -static octave_idx_type -select_dig (const Complex& a) -{ - return (abs (a) < 1.0); -} - -octave_idx_type -ComplexSCHUR::init (const ComplexMatrix& a, const std::string& ord, - bool calc_unitary) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (a_nr != a_nc) - (*current_liboctave_error_handler) ("ComplexSCHUR requires square matrix"); - - if (a_nr == 0) - { - schur_mat.clear (); - unitary_mat.clear (); - return 0; - } - - // Workspace requirements may need to be fixed if any of the - // following change. - - char jobvs; - char sense = 'N'; - char sort = 'N'; - - if (calc_unitary) - jobvs = 'V'; - else - jobvs = 'N'; - - char ord_char = ord.empty () ? 'U' : ord[0]; - - if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') - sort = 'S'; - - if (ord_char == 'A' || ord_char == 'a') - selector = select_ana; - else if (ord_char == 'D' || ord_char == 'd') - selector = select_dig; - else - selector = 0; - - octave_idx_type n = a_nc; - octave_idx_type lwork = 8 * n; - octave_idx_type info; - octave_idx_type sdim; - double rconde; - double rcondv; - - schur_mat = a; - if (calc_unitary) - unitary_mat.clear (n, n); - - Complex *s = schur_mat.fortran_vec (); - Complex *q = unitary_mat.fortran_vec (); - - Array<double> rwork (dim_vector (n, 1)); - double *prwork = rwork.fortran_vec (); - - Array<Complex> w (dim_vector (n, 1)); - Complex *pw = w.fortran_vec (); - - Array<Complex> work (dim_vector (lwork, 1)); - Complex *pwork = work.fortran_vec (); - - // BWORK is not referenced for non-ordered Schur. - octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; - Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); - octave_idx_type *pbwork = bwork.fortran_vec (); - - F77_XFCN (zgeesx, ZGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1), - F77_CONST_CHAR_ARG2 (&sort, 1), - selector, - F77_CONST_CHAR_ARG2 (&sense, 1), - n, s, n, sdim, pw, q, n, rconde, rcondv, - pwork, lwork, prwork, pbwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - return info; -} - -ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& s, const ComplexMatrix& u) - : schur_mat (s), unitary_mat (u), selector (0) -{ - octave_idx_type n = s.rows (); - if (s.columns () != n || u.rows () != n || u.columns () != n) - (*current_liboctave_error_handler) - ("schur: inconsistent matrix dimensions"); -} - -ComplexSCHUR::ComplexSCHUR (const SCHUR& s) - : schur_mat (s.schur_matrix ()), unitary_mat (s.unitary_matrix ()), - selector (0) -{ - octave_idx_type n = schur_mat.rows (); - if (n > 0) - { - OCTAVE_LOCAL_BUFFER (double, c, n-1); - OCTAVE_LOCAL_BUFFER (double, sx, n-1); - - F77_XFCN (zrsf2csf, ZRSF2CSF, (n, schur_mat.fortran_vec (), - unitary_mat.fortran_vec (), c, sx)); - } -}
--- a/liboctave/numeric/CmplxSCHUR.h Thu Jan 14 19:59:52 2016 +1100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,96 +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_CmplxSCHUR_h) -#define octave_CmplxSCHUR_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "CMatrix.h" -#include "dbleSCHUR.h" - -class -OCTAVE_API -ComplexSCHUR -{ -public: - - ComplexSCHUR (void) : schur_mat (), unitary_mat (), selector (0) { } - - ComplexSCHUR (const ComplexMatrix& a, const std::string& ord, - bool calc_unitary = true) - : schur_mat (), unitary_mat (), selector (0) - { - init (a, ord, calc_unitary); - } - - ComplexSCHUR (const ComplexMatrix& a, const std::string& ord, - octave_idx_type& info, - bool calc_unitary = true) - : schur_mat (), unitary_mat (), selector (0) - { - info = init (a, ord, calc_unitary); - } - - ComplexSCHUR (const ComplexSCHUR& a) - : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat), selector (0) - { } - - ComplexSCHUR (const ComplexMatrix& s, const ComplexMatrix& u); - - ComplexSCHUR (const SCHUR& s); - - ComplexSCHUR& operator = (const ComplexSCHUR& a) - { - if (this != &a) - { - schur_mat = a.schur_mat; - unitary_mat = a.unitary_mat; - } - return *this; - } - - ~ComplexSCHUR (void) { } - - ComplexMatrix schur_matrix (void) const { return schur_mat; } - - ComplexMatrix unitary_matrix (void) const { return unitary_mat; } - - friend std::ostream& operator << (std::ostream& os, const ComplexSCHUR& a); - - typedef octave_idx_type (*select_function) (const Complex&); - -private: - - ComplexMatrix schur_mat; - ComplexMatrix unitary_mat; - - select_function selector; - - octave_idx_type init (const ComplexMatrix& a, const std::string& ord, - bool calc_unitary); -}; - -#endif
--- a/liboctave/numeric/dbleSCHUR.cc Thu Jan 14 19:59:52 2016 +1100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,166 +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 <iostream> - -#include "dbleSCHUR.h" -#include "f77-fcn.h" -#include "lo-error.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (dgeesx, DGEESX) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - SCHUR::select_function, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type&, - double*, double*, double*, const octave_idx_type&, - double&, double&, double*, const octave_idx_type&, - octave_idx_type*, const octave_idx_type&, - octave_idx_type*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); -} - -static octave_idx_type -select_ana (const double& a, const double&) -{ - return (a < 0.0); -} - -static octave_idx_type -select_dig (const double& a, const double& b) -{ - return (hypot (a, b) < 1.0); -} - -octave_idx_type -SCHUR::init (const Matrix& a, const std::string& ord, bool calc_unitary) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (a_nr != a_nc) - (*current_liboctave_error_handler) ("SCHUR requires square matrix"); - - if (a_nr == 0) - { - schur_mat.clear (); - unitary_mat.clear (); - return 0; - } - - // Workspace requirements may need to be fixed if any of the - // following change. - - char jobvs; - char sense = 'N'; - char sort = 'N'; - - if (calc_unitary) - jobvs = 'V'; - else - jobvs = 'N'; - - char ord_char = ord.empty () ? 'U' : ord[0]; - - if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') - sort = 'S'; - - if (ord_char == 'A' || ord_char == 'a') - selector = select_ana; - else if (ord_char == 'D' || ord_char == 'd') - selector = select_dig; - else - selector = 0; - - octave_idx_type n = a_nc; - octave_idx_type lwork = 8 * n; - octave_idx_type liwork = 1; - octave_idx_type info; - octave_idx_type sdim; - double rconde; - double rcondv; - - schur_mat = a; - - if (calc_unitary) - unitary_mat.clear (n, n); - - double *s = schur_mat.fortran_vec (); - double *q = unitary_mat.fortran_vec (); - - Array<double> wr (dim_vector (n, 1)); - double *pwr = wr.fortran_vec (); - - Array<double> wi (dim_vector (n, 1)); - double *pwi = wi.fortran_vec (); - - Array<double> work (dim_vector (lwork, 1)); - double *pwork = work.fortran_vec (); - - // BWORK is not referenced for the non-ordered Schur routine. - octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; - Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); - octave_idx_type *pbwork = bwork.fortran_vec (); - - Array<octave_idx_type> iwork (dim_vector (liwork, 1)); - octave_idx_type *piwork = iwork.fortran_vec (); - - F77_XFCN (dgeesx, DGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1), - F77_CONST_CHAR_ARG2 (&sort, 1), - selector, - F77_CONST_CHAR_ARG2 (&sense, 1), - n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv, - pwork, lwork, piwork, liwork, pbwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - return info; -} - -std::ostream& -operator << (std::ostream& os, const SCHUR& a) -{ - os << a.schur_matrix () << "\n"; - os << a.unitary_matrix () << "\n"; - - return os; -} - -SCHUR::SCHUR (const Matrix& s, const Matrix& u) - : schur_mat (s), unitary_mat (u), selector (0) -{ - octave_idx_type n = s.rows (); - if (s.columns () != n || u.rows () != n || u.columns () != n) - (*current_liboctave_error_handler) - ("schur: inconsistent matrix dimensions"); -} -
--- a/liboctave/numeric/dbleSCHUR.h Thu Jan 14 19:59:52 2016 +1100 +++ /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_dbleSCHUR_h) -#define octave_dbleSCHUR_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "dMatrix.h" - -class -OCTAVE_API -SCHUR -{ -public: - - SCHUR (void) : schur_mat (), unitary_mat (), selector (0) { } - - SCHUR (const Matrix& a, const std::string& ord, bool calc_unitary = true) - : schur_mat (), unitary_mat (), selector (0) - { - init (a, ord, calc_unitary); - } - - SCHUR (const Matrix& a, const std::string& ord, int& info, - bool calc_unitary = true) - : schur_mat (), unitary_mat (), selector (0) - { - info = init (a, ord, calc_unitary); - } - - SCHUR (const SCHUR& a) - : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat), selector (0) - { } - - SCHUR (const Matrix& s, const Matrix& u); - - SCHUR& operator = (const SCHUR& a) - { - if (this != &a) - { - schur_mat = a.schur_mat; - unitary_mat = a.unitary_mat; - } - return *this; - } - - ~SCHUR (void) { } - - Matrix schur_matrix (void) const { return schur_mat; } - - Matrix unitary_matrix (void) const { return unitary_mat; } - - friend std::ostream& operator << (std::ostream& os, const SCHUR& a); - - typedef octave_idx_type (*select_function) (const double&, const double&); - -private: - - Matrix schur_mat; - Matrix unitary_mat; - - select_function selector; - - octave_idx_type init (const Matrix& a, const std::string& ord, - bool calc_unitary); -}; - -#endif
--- a/liboctave/numeric/fCmplxSCHUR.cc Thu Jan 14 19:59:52 2016 +1100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,171 +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 "fCmplxSCHUR.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (cgeesx, CGEESX) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - FloatComplexSCHUR::select_function, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, FloatComplex*, - const octave_idx_type&, octave_idx_type&, - FloatComplex*, FloatComplex*, - const octave_idx_type&, float&, float&, - FloatComplex*, const octave_idx_type&, - float*, octave_idx_type*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); - F77_RET_T - F77_FUNC (crsf2csf, CRSF2CSF) (const octave_idx_type&, FloatComplex *, - FloatComplex *, float *, float *); -} - -static octave_idx_type -select_ana (const FloatComplex& a) -{ - return a.real () < 0.0; -} - -static octave_idx_type -select_dig (const FloatComplex& a) -{ - return (abs (a) < 1.0); -} - -octave_idx_type -FloatComplexSCHUR::init (const FloatComplexMatrix& a, const std::string& ord, - bool calc_unitary) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (a_nr != a_nc) - (*current_liboctave_error_handler) - ("FloatComplexSCHUR requires square matrix"); - - if (a_nr == 0) - { - schur_mat.clear (); - unitary_mat.clear (); - return 0; - } - - // Workspace requirements may need to be fixed if any of the - // following change. - - char jobvs; - char sense = 'N'; - char sort = 'N'; - - if (calc_unitary) - jobvs = 'V'; - else - jobvs = 'N'; - - char ord_char = ord.empty () ? 'U' : ord[0]; - - if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') - sort = 'S'; - - if (ord_char == 'A' || ord_char == 'a') - selector = select_ana; - else if (ord_char == 'D' || ord_char == 'd') - selector = select_dig; - else - selector = 0; - - octave_idx_type n = a_nc; - octave_idx_type lwork = 8 * n; - octave_idx_type info; - octave_idx_type sdim; - float rconde; - float rcondv; - - schur_mat = a; - if (calc_unitary) - unitary_mat.clear (n, n); - - FloatComplex *s = schur_mat.fortran_vec (); - FloatComplex *q = unitary_mat.fortran_vec (); - - Array<float> rwork (dim_vector (n, 1)); - float *prwork = rwork.fortran_vec (); - - Array<FloatComplex> w (dim_vector (n, 1)); - FloatComplex *pw = w.fortran_vec (); - - Array<FloatComplex> work (dim_vector (lwork, 1)); - FloatComplex *pwork = work.fortran_vec (); - - // BWORK is not referenced for non-ordered Schur. - octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; - Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); - octave_idx_type *pbwork = bwork.fortran_vec (); - - F77_XFCN (cgeesx, CGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1), - F77_CONST_CHAR_ARG2 (&sort, 1), - selector, - F77_CONST_CHAR_ARG2 (&sense, 1), - n, s, n, sdim, pw, q, n, rconde, rcondv, - pwork, lwork, prwork, pbwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - return info; -} - -FloatComplexSCHUR::FloatComplexSCHUR (const FloatComplexMatrix& s, - const FloatComplexMatrix& u) - : schur_mat (s), unitary_mat (u), selector (0) -{ - octave_idx_type n = s.rows (); - if (s.columns () != n || u.rows () != n || u.columns () != n) - (*current_liboctave_error_handler) - ("schur: inconsistent matrix dimensions"); -} - -FloatComplexSCHUR::FloatComplexSCHUR (const FloatSCHUR& s) - : schur_mat (s.schur_matrix ()), unitary_mat (s.unitary_matrix ()), - selector (0) -{ - octave_idx_type n = schur_mat.rows (); - if (n > 0) - { - OCTAVE_LOCAL_BUFFER (float, c, n-1); - OCTAVE_LOCAL_BUFFER (float, sx, n-1); - - F77_XFCN (crsf2csf, CRSF2CSF, (n, schur_mat.fortran_vec (), - unitary_mat.fortran_vec (), c, sx)); - } -}
--- a/liboctave/numeric/fCmplxSCHUR.h Thu Jan 14 19:59:52 2016 +1100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,96 +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_fCmplxSCHUR_h) -#define octave_fCmplxSCHUR_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "fCMatrix.h" -#include "floatSCHUR.h" - -class -OCTAVE_API -FloatComplexSCHUR -{ -public: - - FloatComplexSCHUR (void) : schur_mat (), unitary_mat (), selector (0) { } - - FloatComplexSCHUR (const FloatComplexMatrix& a, const std::string& ord, - bool calc_unitary = true) - : schur_mat (), unitary_mat (), selector (0) - { - init (a, ord, calc_unitary); - } - - FloatComplexSCHUR (const FloatComplexMatrix& a, const std::string& ord, - octave_idx_type& info, bool calc_unitary = true) - : schur_mat (), unitary_mat (), selector (0) - { - info = init (a, ord, calc_unitary); - } - - FloatComplexSCHUR (const FloatComplexSCHUR& a) - : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat), selector (0) - { } - - FloatComplexSCHUR (const FloatComplexMatrix& s, const FloatComplexMatrix& u); - - FloatComplexSCHUR (const FloatSCHUR& s); - - FloatComplexSCHUR& operator = (const FloatComplexSCHUR& a) - { - if (this != &a) - { - schur_mat = a.schur_mat; - unitary_mat = a.unitary_mat; - } - return *this; - } - - ~FloatComplexSCHUR (void) { } - - FloatComplexMatrix schur_matrix (void) const { return schur_mat; } - - FloatComplexMatrix unitary_matrix (void) const { return unitary_mat; } - - friend std::ostream& operator << (std::ostream& os, - const FloatComplexSCHUR& a); - - typedef octave_idx_type (*select_function) (const FloatComplex&); - -private: - - FloatComplexMatrix schur_mat; - FloatComplexMatrix unitary_mat; - - select_function selector; - - octave_idx_type init (const FloatComplexMatrix& a, const std::string& ord, - bool calc_unitary); -}; - -#endif
--- a/liboctave/numeric/floatSCHUR.cc Thu Jan 14 19:59:52 2016 +1100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,166 +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 <iostream> - -#include "floatSCHUR.h" -#include "f77-fcn.h" -#include "lo-error.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (sgeesx, SGEESX) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - FloatSCHUR::select_function, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, float*, - const octave_idx_type&, octave_idx_type&, - float*, float*, float*, const octave_idx_type&, - float&, float&, float*, const octave_idx_type&, - octave_idx_type*, const octave_idx_type&, - octave_idx_type*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); -} - -static octave_idx_type -select_ana (const float& a, const float&) -{ - return (a < 0.0); -} - -static octave_idx_type -select_dig (const float& a, const float& b) -{ - return (hypot (a, b) < 1.0); -} - -octave_idx_type -FloatSCHUR::init (const FloatMatrix& a, const std::string& ord, - bool calc_unitary) -{ - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (a_nr != a_nc) - (*current_liboctave_error_handler) ("FloatSCHUR requires square matrix"); - - if (a_nr == 0) - { - schur_mat.clear (); - unitary_mat.clear (); - return 0; - } - - // Workspace requirements may need to be fixed if any of the - // following change. - - char jobvs; - char sense = 'N'; - char sort = 'N'; - - if (calc_unitary) - jobvs = 'V'; - else - jobvs = 'N'; - - char ord_char = ord.empty () ? 'U' : ord[0]; - - if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') - sort = 'S'; - - if (ord_char == 'A' || ord_char == 'a') - selector = select_ana; - else if (ord_char == 'D' || ord_char == 'd') - selector = select_dig; - else - selector = 0; - - octave_idx_type n = a_nc; - octave_idx_type lwork = 8 * n; - octave_idx_type liwork = 1; - octave_idx_type info; - octave_idx_type sdim; - float rconde; - float rcondv; - - schur_mat = a; - - if (calc_unitary) - unitary_mat.clear (n, n); - - float *s = schur_mat.fortran_vec (); - float *q = unitary_mat.fortran_vec (); - - Array<float> wr (dim_vector (n, 1)); - float *pwr = wr.fortran_vec (); - - Array<float> wi (dim_vector (n, 1)); - float *pwi = wi.fortran_vec (); - - Array<float> work (dim_vector (lwork, 1)); - float *pwork = work.fortran_vec (); - - // BWORK is not referenced for the non-ordered Schur routine. - octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; - Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); - octave_idx_type *pbwork = bwork.fortran_vec (); - - Array<octave_idx_type> iwork (dim_vector (liwork, 1)); - octave_idx_type *piwork = iwork.fortran_vec (); - - F77_XFCN (sgeesx, SGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1), - F77_CONST_CHAR_ARG2 (&sort, 1), - selector, - F77_CONST_CHAR_ARG2 (&sense, 1), - n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv, - pwork, lwork, piwork, liwork, pbwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - return info; -} - -FloatSCHUR::FloatSCHUR (const FloatMatrix& s, const FloatMatrix& u) - : schur_mat (s), unitary_mat (u), selector (0) -{ - octave_idx_type n = s.rows (); - if (s.columns () != n || u.rows () != n || u.columns () != n) - (*current_liboctave_error_handler) - ("schur: inconsistent matrix dimensions"); -} - -std::ostream& -operator << (std::ostream& os, const FloatSCHUR& a) -{ - os << a.schur_matrix () << "\n"; - os << a.unitary_matrix () << "\n"; - - return os; -}
--- a/liboctave/numeric/floatSCHUR.h Thu Jan 14 19:59:52 2016 +1100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,92 +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_floatSCHUR_h) -#define octave_floatSCHUR_h 1 - -#include "octave-config.h" - -#include <iosfwd> -#include <string> - -#include "fMatrix.h" - -class -OCTAVE_API -FloatSCHUR -{ -public: - - FloatSCHUR (void) : schur_mat (), unitary_mat (), selector (0) { } - - FloatSCHUR (const FloatMatrix& a, const std::string& ord, - bool calc_unitary = true) - : schur_mat (), unitary_mat (), selector (0) - { - init (a, ord, calc_unitary); - } - - FloatSCHUR (const FloatMatrix& a, const std::string& ord, int& info, - bool calc_unitary = true) - : schur_mat (), unitary_mat (), selector (0) - { - info = init (a, ord, calc_unitary); - } - - FloatSCHUR (const FloatSCHUR& a) - : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat), selector (0) - { } - - FloatSCHUR (const FloatMatrix& s, const FloatMatrix& u); - - FloatSCHUR& operator = (const FloatSCHUR& a) - { - if (this != &a) - { - schur_mat = a.schur_mat; - unitary_mat = a.unitary_mat; - } - return *this; - } - - ~FloatSCHUR (void) { } - - FloatMatrix schur_matrix (void) const { return schur_mat; } - - FloatMatrix unitary_matrix (void) const { return unitary_mat; } - - friend std::ostream& operator << (std::ostream& os, const FloatSCHUR& a); - - typedef octave_idx_type (*select_function) (const float&, const float&); - -private: - - FloatMatrix schur_mat; - FloatMatrix unitary_mat; - - select_function selector; - - octave_idx_type init (const FloatMatrix& a, const std::string& ord, - bool calc_unitary); -}; - -#endif
--- a/liboctave/numeric/module.mk Thu Jan 14 19:59:52 2016 +1100 +++ b/liboctave/numeric/module.mk Mon Feb 15 20:06:12 2016 -0500 @@ -15,7 +15,6 @@ liboctave/numeric/CmplxLU.h \ liboctave/numeric/CmplxQR.h \ liboctave/numeric/CmplxQRP.h \ - liboctave/numeric/CmplxSCHUR.h \ liboctave/numeric/CmplxSVD.h \ liboctave/numeric/CollocWt.h \ liboctave/numeric/DAE.h \ @@ -48,7 +47,6 @@ liboctave/numeric/dbleLU.h \ liboctave/numeric/dbleQR.h \ liboctave/numeric/dbleQRP.h \ - liboctave/numeric/dbleSCHUR.h \ liboctave/numeric/dbleSVD.h \ liboctave/numeric/eigs-base.h \ liboctave/numeric/fCmplxAEPBAL.h \ @@ -58,7 +56,6 @@ liboctave/numeric/fCmplxLU.h \ liboctave/numeric/fCmplxQR.h \ liboctave/numeric/fCmplxQRP.h \ - liboctave/numeric/fCmplxSCHUR.h \ liboctave/numeric/fCmplxSVD.h \ liboctave/numeric/fEIG.h \ liboctave/numeric/floatAEPBAL.h \ @@ -68,7 +65,6 @@ liboctave/numeric/floatLU.h \ liboctave/numeric/floatQR.h \ liboctave/numeric/floatQRP.h \ - liboctave/numeric/floatSCHUR.h \ liboctave/numeric/floatSVD.h \ liboctave/numeric/lo-mappers.h \ liboctave/numeric/lo-specfun.h \ @@ -80,6 +76,7 @@ liboctave/numeric/randgamma.h \ liboctave/numeric/randmtzig.h \ liboctave/numeric/randpoisson.h \ + liboctave/numeric/schur.h \ liboctave/numeric/sparse-chol.h \ liboctave/numeric/sparse-dmsolve.h \ liboctave/numeric/sparse-lu.h \ @@ -98,7 +95,6 @@ liboctave/numeric/CmplxLU.cc \ liboctave/numeric/CmplxQR.cc \ liboctave/numeric/CmplxQRP.cc \ - liboctave/numeric/CmplxSCHUR.cc \ liboctave/numeric/CmplxSVD.cc \ liboctave/numeric/CollocWt.cc \ liboctave/numeric/DASPK.cc \ @@ -115,7 +111,6 @@ liboctave/numeric/dbleLU.cc \ liboctave/numeric/dbleQR.cc \ liboctave/numeric/dbleQRP.cc \ - liboctave/numeric/dbleSCHUR.cc \ liboctave/numeric/dbleSVD.cc \ liboctave/numeric/eigs-base.cc \ liboctave/numeric/fCmplxAEPBAL.cc \ @@ -125,7 +120,6 @@ liboctave/numeric/fCmplxLU.cc \ liboctave/numeric/fCmplxQR.cc \ liboctave/numeric/fCmplxQRP.cc \ - liboctave/numeric/fCmplxSCHUR.cc \ liboctave/numeric/fCmplxSVD.cc \ liboctave/numeric/fEIG.cc \ liboctave/numeric/floatAEPBAL.cc \ @@ -135,7 +129,6 @@ liboctave/numeric/floatLU.cc \ liboctave/numeric/floatQR.cc \ liboctave/numeric/floatQRP.cc \ - liboctave/numeric/floatSCHUR.cc \ liboctave/numeric/floatSVD.cc \ liboctave/numeric/lo-mappers.cc \ liboctave/numeric/lo-specfun.cc \ @@ -144,6 +137,7 @@ liboctave/numeric/oct-norm.cc \ liboctave/numeric/oct-rand.cc \ liboctave/numeric/oct-spparms.cc \ + liboctave/numeric/schur.cc \ liboctave/numeric/sparse-chol.cc \ liboctave/numeric/sparse-dmsolve.cc \ liboctave/numeric/sparse-lu.cc \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/schur.cc Mon Feb 15 20:06:12 2016 -0500 @@ -0,0 +1,573 @@ +/* + +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 <iostream> + +#include "CMatrix.h" +#include "dMatrix.h" +#include "f77-fcn.h" +#include "fCMatrix.h" +#include "fMatrix.h" +#include "lo-error.h" +#include "schur.h" + +typedef octave_idx_type (*double_selector) (const double&, const double&); +typedef octave_idx_type (*float_selector) (const float&, const float&); +typedef octave_idx_type (*complex_selector) (const Complex&); +typedef octave_idx_type (*float_complex_selector) (const FloatComplex&); + +extern "C" +{ + F77_RET_T + F77_FUNC (dgeesx, DGEESX) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + double_selector, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type&, + double*, double*, double*, const octave_idx_type&, + double&, double&, double*, const octave_idx_type&, + octave_idx_type*, const octave_idx_type&, + octave_idx_type*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (sgeesx, SGEESX) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + float_selector, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, float*, + const octave_idx_type&, octave_idx_type&, + float*, float*, float*, const octave_idx_type&, + float&, float&, float*, const octave_idx_type&, + octave_idx_type*, const octave_idx_type&, + octave_idx_type*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zgeesx, ZGEESX) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + complex_selector, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, Complex*, + const octave_idx_type&, octave_idx_type&, + Complex*, Complex*, const octave_idx_type&, + double&, double&, Complex*, + const octave_idx_type&, double*, + octave_idx_type*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zrsf2csf, ZRSF2CSF) (const octave_idx_type&, Complex *, + Complex *, double *, double *); + F77_RET_T + F77_FUNC (cgeesx, CGEESX) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + float_complex_selector, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, FloatComplex*, + const octave_idx_type&, octave_idx_type&, + FloatComplex*, FloatComplex*, + const octave_idx_type&, float&, float&, + FloatComplex*, const octave_idx_type&, + float*, octave_idx_type*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (crsf2csf, CRSF2CSF) (const octave_idx_type&, FloatComplex *, + FloatComplex *, float *, float *); +} + +// For real types. + +template <typename T> +static octave_idx_type +select_ana (const T& a, const T&) +{ + return (a < 0.0); +} + +template <typename T> +static octave_idx_type +select_dig (const T& a, const T& b) +{ + return (hypot (a, b) < 1.0); +} + +// For complex types. + +template <typename T> +static octave_idx_type +select_ana (const T& a) +{ + return a.real () < 0.0; +} + +template <typename T> +static octave_idx_type +select_dig (const T& a) +{ + return (abs (a) < 1.0); +} + +template <> +octave_idx_type +schur<Matrix>::init (const Matrix& a, const std::string& ord, bool calc_unitary) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + if (a_nr != a_nc) + (*current_liboctave_error_handler) ("schur: requires square matrix"); + + if (a_nr == 0) + { + schur_mat.clear (); + unitary_mat.clear (); + return 0; + } + + // Workspace requirements may need to be fixed if any of the + // following change. + + char jobvs; + char sense = 'N'; + char sort = 'N'; + + if (calc_unitary) + jobvs = 'V'; + else + jobvs = 'N'; + + char ord_char = ord.empty () ? 'U' : ord[0]; + + if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') + sort = 'S'; + + double_selector selector = 0; + if (ord_char == 'A' || ord_char == 'a') + selector = select_ana<double>; + else if (ord_char == 'D' || ord_char == 'd') + selector = select_dig<double>; + else + selector = 0; + + octave_idx_type n = a_nc; + octave_idx_type lwork = 8 * n; + octave_idx_type liwork = 1; + octave_idx_type info; + octave_idx_type sdim; + double rconde; + double rcondv; + + schur_mat = a; + + if (calc_unitary) + unitary_mat.clear (n, n); + + double *s = schur_mat.fortran_vec (); + double *q = unitary_mat.fortran_vec (); + + Array<double> wr (dim_vector (n, 1)); + double *pwr = wr.fortran_vec (); + + Array<double> wi (dim_vector (n, 1)); + double *pwi = wi.fortran_vec (); + + Array<double> work (dim_vector (lwork, 1)); + double *pwork = work.fortran_vec (); + + // BWORK is not referenced for the non-ordered Schur routine. + octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; + Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); + octave_idx_type *pbwork = bwork.fortran_vec (); + + Array<octave_idx_type> iwork (dim_vector (liwork, 1)); + octave_idx_type *piwork = iwork.fortran_vec (); + + F77_XFCN (dgeesx, DGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1), + F77_CONST_CHAR_ARG2 (&sort, 1), + selector, + F77_CONST_CHAR_ARG2 (&sense, 1), + n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv, + pwork, lwork, piwork, liwork, pbwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +template <> +octave_idx_type +schur<FloatMatrix>::init (const FloatMatrix& a, const std::string& ord, + bool calc_unitary) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + if (a_nr != a_nc) + (*current_liboctave_error_handler) ("SCHUR requires square matrix"); + + if (a_nr == 0) + { + schur_mat.clear (); + unitary_mat.clear (); + return 0; + } + + // Workspace requirements may need to be fixed if any of the + // following change. + + char jobvs; + char sense = 'N'; + char sort = 'N'; + + if (calc_unitary) + jobvs = 'V'; + else + jobvs = 'N'; + + char ord_char = ord.empty () ? 'U' : ord[0]; + + if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') + sort = 'S'; + + float_selector selector = 0; + if (ord_char == 'A' || ord_char == 'a') + selector = select_ana<float>; + else if (ord_char == 'D' || ord_char == 'd') + selector = select_dig<float>; + else + selector = 0; + + octave_idx_type n = a_nc; + octave_idx_type lwork = 8 * n; + octave_idx_type liwork = 1; + octave_idx_type info; + octave_idx_type sdim; + float rconde; + float rcondv; + + schur_mat = a; + + if (calc_unitary) + unitary_mat.clear (n, n); + + float *s = schur_mat.fortran_vec (); + float *q = unitary_mat.fortran_vec (); + + Array<float> wr (dim_vector (n, 1)); + float *pwr = wr.fortran_vec (); + + Array<float> wi (dim_vector (n, 1)); + float *pwi = wi.fortran_vec (); + + Array<float> work (dim_vector (lwork, 1)); + float *pwork = work.fortran_vec (); + + // BWORK is not referenced for the non-ordered Schur routine. + octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; + Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); + octave_idx_type *pbwork = bwork.fortran_vec (); + + Array<octave_idx_type> iwork (dim_vector (liwork, 1)); + octave_idx_type *piwork = iwork.fortran_vec (); + + F77_XFCN (sgeesx, SGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1), + F77_CONST_CHAR_ARG2 (&sort, 1), + selector, + F77_CONST_CHAR_ARG2 (&sense, 1), + n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv, + pwork, lwork, piwork, liwork, pbwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +template <> +octave_idx_type +schur<ComplexMatrix>::init (const ComplexMatrix& a, const std::string& ord, + bool calc_unitary) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + if (a_nr != a_nc) + (*current_liboctave_error_handler) ("SCHUR requires square matrix"); + + if (a_nr == 0) + { + schur_mat.clear (); + unitary_mat.clear (); + return 0; + } + + // Workspace requirements may need to be fixed if any of the + // following change. + + char jobvs; + char sense = 'N'; + char sort = 'N'; + + if (calc_unitary) + jobvs = 'V'; + else + jobvs = 'N'; + + char ord_char = ord.empty () ? 'U' : ord[0]; + + if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') + sort = 'S'; + + complex_selector selector = 0; + if (ord_char == 'A' || ord_char == 'a') + selector = select_ana<Complex>; + else if (ord_char == 'D' || ord_char == 'd') + selector = select_dig<Complex>; + else + selector = 0; + + octave_idx_type n = a_nc; + octave_idx_type lwork = 8 * n; + octave_idx_type info; + octave_idx_type sdim; + double rconde; + double rcondv; + + schur_mat = a; + if (calc_unitary) + unitary_mat.clear (n, n); + + Complex *s = schur_mat.fortran_vec (); + Complex *q = unitary_mat.fortran_vec (); + + Array<double> rwork (dim_vector (n, 1)); + double *prwork = rwork.fortran_vec (); + + Array<Complex> w (dim_vector (n, 1)); + Complex *pw = w.fortran_vec (); + + Array<Complex> work (dim_vector (lwork, 1)); + Complex *pwork = work.fortran_vec (); + + // BWORK is not referenced for non-ordered Schur. + octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; + Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); + octave_idx_type *pbwork = bwork.fortran_vec (); + + F77_XFCN (zgeesx, ZGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1), + F77_CONST_CHAR_ARG2 (&sort, 1), + selector, + F77_CONST_CHAR_ARG2 (&sense, 1), + n, s, n, sdim, pw, q, n, rconde, rcondv, + pwork, lwork, prwork, pbwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +template <> +schur<ComplexMatrix> +rsf2csf<ComplexMatrix, Matrix> (const Matrix& s_arg, const Matrix& u_arg) +{ + ComplexMatrix s (s_arg); + ComplexMatrix u (u_arg); + + octave_idx_type n = s.rows (); + + if (s.columns () != n || u.rows () != n || u.columns () != n) + (*current_liboctave_error_handler) + ("rsf2csf: inconsistent matrix dimensions"); + + if (n > 0) + { + OCTAVE_LOCAL_BUFFER (double, c, n-1); + OCTAVE_LOCAL_BUFFER (double, sx, n-1); + + F77_XFCN (zrsf2csf, ZRSF2CSF, (n, s.fortran_vec (), + u.fortran_vec (), c, sx)); + } + + return schur<ComplexMatrix> (s, u); +} + +template <> +octave_idx_type +schur<FloatComplexMatrix>::init (const FloatComplexMatrix& a, + const std::string& ord, bool calc_unitary) +{ + octave_idx_type a_nr = a.rows (); + octave_idx_type a_nc = a.cols (); + + if (a_nr != a_nc) + (*current_liboctave_error_handler) ("SCHUR requires square matrix"); + + if (a_nr == 0) + { + schur_mat.clear (); + unitary_mat.clear (); + return 0; + } + + // Workspace requirements may need to be fixed if any of the + // following change. + + char jobvs; + char sense = 'N'; + char sort = 'N'; + + if (calc_unitary) + jobvs = 'V'; + else + jobvs = 'N'; + + char ord_char = ord.empty () ? 'U' : ord[0]; + + if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') + sort = 'S'; + + float_complex_selector selector = 0; + if (ord_char == 'A' || ord_char == 'a') + selector = select_ana<FloatComplex>; + else if (ord_char == 'D' || ord_char == 'd') + selector = select_dig<FloatComplex>; + else + selector = 0; + + octave_idx_type n = a_nc; + octave_idx_type lwork = 8 * n; + octave_idx_type info; + octave_idx_type sdim; + float rconde; + float rcondv; + + schur_mat = a; + if (calc_unitary) + unitary_mat.clear (n, n); + + FloatComplex *s = schur_mat.fortran_vec (); + FloatComplex *q = unitary_mat.fortran_vec (); + + Array<float> rwork (dim_vector (n, 1)); + float *prwork = rwork.fortran_vec (); + + Array<FloatComplex> w (dim_vector (n, 1)); + FloatComplex *pw = w.fortran_vec (); + + Array<FloatComplex> work (dim_vector (lwork, 1)); + FloatComplex *pwork = work.fortran_vec (); + + // BWORK is not referenced for non-ordered Schur. + octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; + Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); + octave_idx_type *pbwork = bwork.fortran_vec (); + + F77_XFCN (cgeesx, CGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1), + F77_CONST_CHAR_ARG2 (&sort, 1), + selector, + F77_CONST_CHAR_ARG2 (&sense, 1), + n, s, n, sdim, pw, q, n, rconde, rcondv, + pwork, lwork, prwork, pbwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + return info; +} + +template <> +schur<FloatComplexMatrix> +rsf2csf<FloatComplexMatrix, FloatMatrix> (const FloatMatrix& s_arg, const FloatMatrix& u_arg) +{ + FloatComplexMatrix s (s_arg); + FloatComplexMatrix u (u_arg); + + octave_idx_type n = s.rows (); + + if (s.columns () != n || u.rows () != n || u.columns () != n) + (*current_liboctave_error_handler) + ("rsf2csf: inconsistent matrix dimensions"); + + if (n > 0) + { + OCTAVE_LOCAL_BUFFER (float, c, n-1); + OCTAVE_LOCAL_BUFFER (float, sx, n-1); + + F77_XFCN (crsf2csf, CRSF2CSF, (n, s.fortran_vec (), + u.fortran_vec (), c, sx)); + } + + return schur<FloatComplexMatrix> (s, u); +} + +template <typename T> +std::ostream& +operator << (std::ostream& os, const schur<T>& a) +{ + os << a.schur_matrix () << "\n"; + os << a.unitary_matrix () << "\n"; + + return os; +} + +// Instantiations we need. + +template class schur<ComplexMatrix>; + +template class schur<FloatComplexMatrix>; + +template class schur<FloatMatrix>; + +template class schur<Matrix>; + +template <> +std::ostream& +operator << (std::ostream& os, const schur<Matrix>& a); + +template <> +std::ostream& +operator << (std::ostream& os, const schur<ComplexMatrix>& a); + +template <> +std::ostream& +operator << (std::ostream& os, const schur<FloatMatrix>& a); + +template <> +std::ostream& +operator << (std::ostream& os, const schur<FloatComplexMatrix>& a);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/schur.h Mon Feb 15 20:06:12 2016 -0500 @@ -0,0 +1,105 @@ +/* + +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_schur_h) +#define octave_schur_h 1 + +#include "octave-config.h" + +#include <iosfwd> +#include <string> + +#include "dMatrix.h" +#include "CMatrix.h" +#include "fMatrix.h" +#include "fCMatrix.h" + +template <typename T> class schur; + +template <typename T> +class +schur +{ +public: + + schur (void) : schur_mat (), unitary_mat () { } + + schur (const T& a, const std::string& ord, bool calc_unitary = true) + : schur_mat (), unitary_mat () + { + init (a, ord, calc_unitary); + } + + schur (const T& a, const std::string& ord, octave_idx_type& info, + bool calc_unitary = true) + : schur_mat (), unitary_mat () + { + info = init (a, ord, calc_unitary); + } + + // This one should really be protected or private but we need it in + // rsf2csf and I don't see how to make that function a friend of + // this class. + schur (const T& s, const T& u) : schur_mat (s), unitary_mat (u) { } + + schur (const schur& a) + + : schur_mat (a.schur_mat), unitary_mat (a.unitary_mat) + { } + + schur& operator = (const schur& a) + { + if (this != &a) + { + schur_mat = a.schur_mat; + unitary_mat = a.unitary_mat; + } + + return *this; + } + + ~schur (void) { } + + T schur_matrix (void) const { return schur_mat; } + + T unitary_matrix (void) const { return unitary_mat; } + +protected: + +private: + + T schur_mat; + T unitary_mat; + + octave_idx_type + init (const T& a, const std::string& ord, bool calc_unitary); +}; + +template <typename RT, typename AT> +extern schur<RT> +rsf2csf (const AT& s, const AT& u); + +template <typename T> +extern std::ostream& +operator << (std::ostream& os, const schur<T>& a); + +#endif
--- a/liboctave/operators/mx-defs.h Thu Jan 14 19:59:52 2016 +1100 +++ b/liboctave/operators/mx-defs.h Mon Feb 15 20:06:12 2016 -0500 @@ -80,10 +80,7 @@ class FloatHESS; class FloatComplexHESS; -class SCHUR; -class ComplexSCHUR; -class FloatSCHUR; -class FloatComplexSCHUR; +template <typename T> class schur; class SVD; class ComplexSVD;
--- a/liboctave/operators/mx-ext.h Thu Jan 14 19:59:52 2016 +1100 +++ b/liboctave/operators/mx-ext.h Mon Feb 15 20:06:12 2016 -0500 @@ -48,10 +48,7 @@ // Result of a Schur Decomposition -#include "dbleSCHUR.h" -#include "CmplxSCHUR.h" -#include "floatSCHUR.h" -#include "fCmplxSCHUR.h" +#include "schur.h" // Result of a Singular Value Decomposition.