Mercurial > octave
view liboctave/numeric/schur.cc @ 23450:855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
* Canvas.cc, color-picker.cc, dialog.cc, marker.cc, main-window.cc, Cell.cc,
__magick_read__.cc, __pchip_deriv__.cc, __qp__.cc, bsxfun.cc, call-stack.cc,
cellfun.cc, data.cc, defaults.cc, det.cc, eig.cc, error.cc, file-io.cc,
filter.cc, find.cc, gammainc.cc, gcd.cc, gl-render.cc, graphics.cc,
graphics.in.h, help.cc, hex2num.cc, inv.cc, load-save.cc, lookup.cc,
ls-mat4.cc, ls-mat5.cc, ls-oct-binary.cc, ls-oct-text.cc, lu.cc, max.cc,
mex.cc, oct-hist.cc, oct-map.cc, oct-procbuf.cc, oct-stream.cc, pr-output.cc,
rand.cc, regexp.cc, schur.cc, str2double.cc, strfns.cc, symtab.cc, sysdep.cc,
tril.cc, variables.cc, xdiv.cc, audiodevinfo.cc, audioread.cc, colamd.cc,
gzip.cc, ov-base-diag.cc, ov-base-mat.h, ov-base-scalar.h, ov-class.cc,
ov-classdef.cc, ov-fcn-handle.cc, ov-range.cc, ov-range.h, ov-struct.cc,
ov-usr-fcn.cc, ov.cc, octave.cc, op-class.cc, pt-eval.cc, pt-funcall.cc,
pt-idx.cc, pt-jit.cc, pt-stmt.cc, version.cc, Array-util.cc, Array.cc,
CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CSparse.cc, Range.cc, Sparse.cc,
chNDArray.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dSparse.cc,
dim-vector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fDiagMatrix.cc,
fMatrix.cc, fNDArray.cc, idx-vector.cc, idx-vector.h, Faddeeva.cc, DASPK.cc,
DASRT.cc, DASSL.cc, EIG.cc, fEIG.cc, gsvd.cc, lo-specfun.cc, oct-rand.cc,
qr.cc, qrp.cc, randgamma.cc, schur.cc, svd.cc, Sparse-diag-op-defs.h,
mx-inlines.cc, oct-env.cc, cmd-edit.cc, kpse.cc, oct-inttypes.h, oct-sort.cc:
Wrap tertiary operator in parentheses "(COND ? x : y)".
author | Rik <rik@octave.org> |
---|---|
date | Thu, 27 Apr 2017 17:33:10 -0700 |
parents | 092078913d54 |
children | d691ed308237 |
line wrap: on
line source
/* Copyright (C) 1994-2017 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 (HAVE_CONFIG_H) # include "config.h" #endif #include "CMatrix.h" #include "dMatrix.h" #include "fCMatrix.h" #include "fMatrix.h" #include "lo-error.h" #include "lo-lapack-proto.h" #include "schur.h" namespace octave { namespace math { // For real types. static F77_INT select_ana (const double& a, const double&) { return (a < 0.0); } static F77_INT select_dig (const double& a, const double& b) { return (hypot (a, b) < 1.0); } static F77_INT select_ana (const float& a, const float&) { return (a < 0.0); } static F77_INT select_dig (const float& a, const float& b) { return (hypot (a, b) < 1.0); } // For complex types. static F77_INT select_ana (const F77_DBLE_CMPLX& a_arg) { const Complex a = reinterpret_cast<const Complex&> (a_arg); return a.real () < 0.0; } static F77_INT select_dig (const F77_DBLE_CMPLX& a_arg) { const Complex& a = reinterpret_cast<const Complex&> (a_arg); return (abs (a) < 1.0); } static F77_INT select_ana (const F77_CMPLX& a_arg) { const FloatComplex& a = reinterpret_cast<const FloatComplex&> (a_arg); return a.real () < 0.0; } static F77_INT select_dig (const F77_CMPLX& a_arg) { const FloatComplex& a = reinterpret_cast<const FloatComplex&> (a_arg); return (abs (a) < 1.0); } template <> F77_INT schur<Matrix>::init (const Matrix& a, const std::string& ord, bool calc_unitary) { F77_INT a_nr = octave::to_f77_int (a.rows ()); F77_INT a_nc = octave::to_f77_int (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'; volatile double_selector selector = 0; if (ord_char == 'A' || ord_char == 'a') selector = select_ana; else if (ord_char == 'D' || ord_char == 'd') selector = select_dig; F77_INT n = a_nc; F77_INT lwork = 8 * n; F77_INT liwork = 1; F77_INT info; F77_INT 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. F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; Array<F77_INT> bwork (dim_vector (ntmp, 1)); F77_INT *pbwork = bwork.fortran_vec (); Array<F77_INT> iwork (dim_vector (liwork, 1)); F77_INT *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 <> F77_INT schur<FloatMatrix>::init (const FloatMatrix& a, const std::string& ord, bool calc_unitary) { F77_INT a_nr = octave::to_f77_int (a.rows ()); F77_INT a_nc = octave::to_f77_int (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'; volatile float_selector selector = 0; if (ord_char == 'A' || ord_char == 'a') selector = select_ana; else if (ord_char == 'D' || ord_char == 'd') selector = select_dig; F77_INT n = a_nc; F77_INT lwork = 8 * n; F77_INT liwork = 1; F77_INT info; F77_INT 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. F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; Array<F77_INT> bwork (dim_vector (ntmp, 1)); F77_INT *pbwork = bwork.fortran_vec (); Array<F77_INT> iwork (dim_vector (liwork, 1)); F77_INT *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 <> F77_INT schur<ComplexMatrix>::init (const ComplexMatrix& a, const std::string& ord, bool calc_unitary) { F77_INT a_nr = octave::to_f77_int (a.rows ()); F77_INT a_nc = octave::to_f77_int (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'; volatile complex_selector selector = 0; if (ord_char == 'A' || ord_char == 'a') selector = select_ana; else if (ord_char == 'D' || ord_char == 'd') selector = select_dig; F77_INT n = a_nc; F77_INT lwork = 8 * n; F77_INT info; F77_INT 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. F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; Array<F77_INT> bwork (dim_vector (ntmp, 1)); F77_INT *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, F77_DBLE_CMPLX_ARG (s), n, sdim, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (q), n, rconde, rcondv, F77_DBLE_CMPLX_ARG (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); F77_INT n = octave::to_f77_int (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, F77_DBLE_CMPLX_ARG (s.fortran_vec ()), F77_DBLE_CMPLX_ARG (u.fortran_vec ()), c, sx)); } return schur<ComplexMatrix> (s, u); } template <> F77_INT schur<FloatComplexMatrix>::init (const FloatComplexMatrix& a, const std::string& ord, bool calc_unitary) { F77_INT a_nr = octave::to_f77_int (a.rows ()); F77_INT a_nc = octave::to_f77_int (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'; volatile float_complex_selector selector = 0; if (ord_char == 'A' || ord_char == 'a') selector = select_ana; else if (ord_char == 'D' || ord_char == 'd') selector = select_dig; F77_INT n = a_nc; F77_INT lwork = 8 * n; F77_INT info; F77_INT 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. F77_INT ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; Array<F77_INT> bwork (dim_vector (ntmp, 1)); F77_INT *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, F77_CMPLX_ARG (s), n, sdim, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (q), n, rconde, rcondv, F77_CMPLX_ARG (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); F77_INT n = octave::to_f77_int (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, F77_CMPLX_ARG (s.fortran_vec ()), F77_CMPLX_ARG (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>; } }