Mercurial > octave
view liboctave/numeric/schur.cc @ 21274:bc536eff5eab
eliminate ostream operator << from schur classes
* liboctave/numeric/schur.h, liboctave/numeric/schur.cc:
Delete operator <<.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 16 Feb 2016 14:47:35 -0500 |
parents | e69eaee28737 |
children | 40de9f8f23a6 |
line wrap: on
line source
/* 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 "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); } // Instantiations we need. template class schur<ComplexMatrix>; template class schur<FloatComplexMatrix>; template class schur<FloatMatrix>; template class schur<Matrix>;