Mercurial > octave
changeset 21273:cbced1c09916
better use of templates for svd classes
* liboctave/numeric/svd.h, liboctave/numeric/svd.cc: New files for svd
classes generated from CmplxSVD.cc, CmplxSVD.h, dbleSVD.cc, dbleSVD.h,
fCmplxSVD.cc, fCmplxSVD.h, floatSVD.cc, and floatSVD.h and converted
to templates.
* liboctave/numeric/module.mk: Update.
* __qp__.cc, svd.cc, CMatrix.cc, CMatrix.h, dDiagMatrix.h, dMatrix.cc,
dMatrix.h, fCMatrix.cc, fCMatrix.h, fDiagMatrix.h, fMatrix.cc,
fMatrix.h, oct-norm.cc, mx-defs.h, mx-ext.h: Use new classes.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 16 Feb 2016 14:41:06 -0500 |
| parents | 987c1a79d33f |
| children | bc536eff5eab |
| files | libinterp/corefcn/__qp__.cc libinterp/corefcn/svd.cc liboctave/array/CMatrix.cc liboctave/array/CMatrix.h liboctave/array/dDiagMatrix.h liboctave/array/dMatrix.cc liboctave/array/dMatrix.h liboctave/array/fCMatrix.cc liboctave/array/fCMatrix.h liboctave/array/fDiagMatrix.h liboctave/array/fMatrix.cc liboctave/array/fMatrix.h liboctave/numeric/CmplxSVD.cc liboctave/numeric/CmplxSVD.h liboctave/numeric/dbleSVD.cc liboctave/numeric/dbleSVD.h liboctave/numeric/fCmplxSVD.cc liboctave/numeric/fCmplxSVD.h liboctave/numeric/floatSVD.cc liboctave/numeric/floatSVD.h liboctave/numeric/module.mk liboctave/numeric/oct-norm.cc liboctave/numeric/svd.cc liboctave/numeric/svd.h liboctave/operators/mx-defs.h liboctave/operators/mx-ext.h |
| diffstat | 26 files changed, 969 insertions(+), 1423 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/__qp__.cc Tue Feb 16 14:39:52 2016 -0500 +++ b/libinterp/corefcn/__qp__.cc Tue Feb 16 14:41:06 2016 -0500 @@ -27,7 +27,7 @@ #include <cfloat> #include "chol.h" -#include "dbleSVD.h" +#include "svd.h" #include "mx-m-dm.h" #include "EIG.h" @@ -47,7 +47,7 @@ if (! A.is_empty ()) { - SVD A_svd (A); + svd<Matrix> A_svd (A); DiagMatrix S = A_svd.singular_values ();
--- a/libinterp/corefcn/svd.cc Tue Feb 16 14:39:52 2016 -0500 +++ b/libinterp/corefcn/svd.cc Tue Feb 16 14:41:06 2016 -0500 @@ -24,10 +24,7 @@ # include <config.h> #endif -#include "CmplxSVD.h" -#include "dbleSVD.h" -#include "fCmplxSVD.h" -#include "floatSVD.h" +#include "svd.h" #include "defun.h" #include "error.h" @@ -37,7 +34,23 @@ #include "utils.h" #include "variables.h" -static int Vsvd_driver = SVD::GESVD; +static std::string Vsvd_driver = "gesvd"; + +template <typename T> +static typename svd<T>::type +svd_type (int nargin, int nargout) +{ + return ((nargout == 0 || nargout == 1) + ? svd<T>::sigma_only + : (nargin == 2) ? svd<T>::economy : svd<T>::std); +} + +template <typename T> +static typename svd<T>::driver +svd_driver (void) +{ + return Vsvd_driver == "gesvd" ? svd<T>::GESVD : svd<T>::GESDD; +} DEFUN (svd, args, nargout, "-*- texinfo -*-\n\ @@ -138,142 +151,93 @@ bool isfloat = arg.is_single_type (); - SVD::type type = ((nargout == 0 || nargout == 1) - ? SVD::sigma_only - : (nargin == 2) ? SVD::economy : SVD::std); - - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); - - SVD::driver driver = static_cast<SVD::driver> (Vsvd_driver); + if (isfloat) + { + if (arg.is_real_type ()) + { + FloatMatrix tmp = arg.float_matrix_value (); - if (nr == 0 || nc == 0) - { - if (isfloat) - { - switch (type) - { - case SVD::std: - retval = ovl (FloatDiagMatrix (nr, nr, 1.0f), - FloatMatrix (nr, nc), - FloatDiagMatrix (nc, nc, 1.0f)); - break; + if (tmp.any_element_is_inf_or_nan ()) + error ("svd: cannot take SVD of matrix containing Inf or NaN values"); + + svd<FloatMatrix> result (tmp, + svd_type<FloatMatrix> (nargin, nargout), + svd_driver<FloatMatrix> ()); + + FloatDiagMatrix sigma = result.singular_values (); - case SVD::economy: - retval = ovl (FloatDiagMatrix (nr, 0, 1.0f), - FloatMatrix (0, 0), - FloatDiagMatrix (0, nc, 1.0f)); - break; - - case SVD::sigma_only: default: - retval(0) = FloatMatrix (0, 1); - break; - } + if (nargout == 0 || nargout == 1) + retval(0) = sigma.extract_diag (); + else + retval = ovl (result.left_singular_matrix (), + sigma, + result.right_singular_matrix ()); } - else + else if (arg.is_complex_type ()) { - switch (type) - { - case SVD::std: - retval = ovl (DiagMatrix (nr, nr, 1.0), - Matrix (nr, nc), - DiagMatrix (nc, nc, 1.0)); - break; + FloatComplexMatrix ctmp = arg.float_complex_matrix_value (); + + if (ctmp.any_element_is_inf_or_nan ()) + error ("svd: cannot take SVD of matrix containing Inf or NaN values"); + + svd<FloatComplexMatrix> result (ctmp, + svd_type<FloatComplexMatrix> (nargin, nargout), + svd_driver<FloatComplexMatrix> ()); - case SVD::economy: - retval = ovl (DiagMatrix (nr, 0, 1.0), - Matrix (0, 0), - DiagMatrix (0, nc, 1.0)); - break; + FloatDiagMatrix sigma = result.singular_values (); - case SVD::sigma_only: default: - retval(0) = Matrix (0, 1); - break; - } + if (nargout == 0 || nargout == 1) + retval(0) = sigma.extract_diag (); + else + retval = ovl (result.left_singular_matrix (), + sigma, + result.right_singular_matrix ()); } } else { - if (isfloat) + if (arg.is_real_type ()) { - if (arg.is_real_type ()) - { - FloatMatrix tmp = arg.float_matrix_value (); + Matrix tmp = arg.matrix_value (); - if (tmp.any_element_is_inf_or_nan ()) - error ("svd: cannot take SVD of matrix containing Inf or NaN values"); + if (tmp.any_element_is_inf_or_nan ()) + error ("svd: cannot take SVD of matrix containing Inf or NaN values"); - FloatSVD result (tmp, type, driver); + svd<Matrix> result (tmp, + svd_type<Matrix> (nargin, nargout), + svd_driver<Matrix> ()); - FloatDiagMatrix sigma = result.singular_values (); + DiagMatrix sigma = result.singular_values (); - if (nargout == 0 || nargout == 1) - retval(0) = sigma.extract_diag (); - else - retval = ovl (result.left_singular_matrix (), - sigma, - result.right_singular_matrix ()); - } - else if (arg.is_complex_type ()) - { - FloatComplexMatrix ctmp = arg.float_complex_matrix_value (); + if (nargout == 0 || nargout == 1) + retval(0) = sigma.extract_diag (); + else + retval = ovl (result.left_singular_matrix (), + sigma, + result.right_singular_matrix ()); + } + else if (arg.is_complex_type ()) + { + ComplexMatrix ctmp = arg.complex_matrix_value (); - if (ctmp.any_element_is_inf_or_nan ()) - error ("svd: cannot take SVD of matrix containing Inf or NaN values"); + if (ctmp.any_element_is_inf_or_nan ()) + error ("svd: cannot take SVD of matrix containing Inf or NaN values"); - FloatComplexSVD result (ctmp, type, driver); - - FloatDiagMatrix sigma = result.singular_values (); + svd<ComplexMatrix> result (ctmp, + svd_type<ComplexMatrix> (nargin, nargout), + svd_driver<ComplexMatrix> ()); - if (nargout == 0 || nargout == 1) - retval(0) = sigma.extract_diag (); - else - retval = ovl (result.left_singular_matrix (), - sigma, - result.right_singular_matrix ()); - } + DiagMatrix sigma = result.singular_values (); + + if (nargout == 0 || nargout == 1) + retval(0) = sigma.extract_diag (); + else + retval = ovl (result.left_singular_matrix (), + sigma, + result.right_singular_matrix ()); } else - { - if (arg.is_real_type ()) - { - Matrix tmp = arg.matrix_value (); - - if (tmp.any_element_is_inf_or_nan ()) - error ("svd: cannot take SVD of matrix containing Inf or NaN values"); - - SVD result (tmp, type, driver); - - DiagMatrix sigma = result.singular_values (); - - if (nargout == 0 || nargout == 1) - retval(0) = sigma.extract_diag (); - else - retval = ovl (result.left_singular_matrix (), - sigma, - result.right_singular_matrix ()); - } - else if (arg.is_complex_type ()) - { - ComplexMatrix ctmp = arg.complex_matrix_value (); - - if (ctmp.any_element_is_inf_or_nan ()) - error ("svd: cannot take SVD of matrix containing Inf or NaN values"); - - ComplexSVD result (ctmp, type, driver); - - DiagMatrix sigma = result.singular_values (); - - if (nargout == 0 || nargout == 1) - retval(0) = sigma.extract_diag (); - else - retval = ovl (result.left_singular_matrix (), - sigma, - result.right_singular_matrix ()); - } - else - err_wrong_type_arg ("svd", arg); - } + err_wrong_type_arg ("svd", arg); } return retval;
--- a/liboctave/array/CMatrix.cc Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/CMatrix.cc Tue Feb 16 14:41:06 2016 -0500 @@ -47,7 +47,7 @@ #include "CDiagMatrix.h" #include "dDiagMatrix.h" #include "schur.h" -#include "CmplxSVD.h" +#include "svd.h" #include "DET.h" #include "f77-fcn.h" #include "functor.h" @@ -1137,7 +1137,7 @@ { ComplexMatrix retval; - ComplexSVD result (*this, SVD::economy); + svd<ComplexMatrix> result (*this, svd<ComplexMatrix>::economy); DiagMatrix S = result.singular_values (); ComplexMatrix U = result.left_singular_matrix ();
--- a/liboctave/array/CMatrix.h Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/CMatrix.h Tue Feb 16 14:41:06 2016 -0500 @@ -50,6 +50,9 @@ typedef Matrix real_matrix_type; typedef ComplexMatrix complex_matrix_type; + typedef DiagMatrix real_diag_matrix_type; + typedef ComplexDiagMatrix complex_diag_matrix_type; + typedef double real_elt_type; typedef Complex complex_elt_type;
--- a/liboctave/array/dDiagMatrix.h Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/dDiagMatrix.h Tue Feb 16 14:41:06 2016 -0500 @@ -37,9 +37,6 @@ OCTAVE_API DiagMatrix : public MDiagArray2<double> { - friend class SVD; - friend class ComplexSVD; - public: DiagMatrix (void) : MDiagArray2<double> () { }
--- a/liboctave/array/dMatrix.cc Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/dMatrix.cc Tue Feb 16 14:41:06 2016 -0500 @@ -46,7 +46,7 @@ #include "PermMatrix.h" #include "DET.h" #include "schur.h" -#include "dbleSVD.h" +#include "svd.h" #include "f77-fcn.h" #include "functor.h" #include "lo-error.h" @@ -825,7 +825,7 @@ Matrix Matrix::pseudo_inverse (double tol) const { - SVD result (*this, SVD::economy); + svd<Matrix> result (*this, svd<Matrix>::economy); DiagMatrix S = result.singular_values (); Matrix U = result.left_singular_matrix ();
--- a/liboctave/array/dMatrix.h Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/dMatrix.h Tue Feb 16 14:41:06 2016 -0500 @@ -49,6 +49,9 @@ typedef Matrix real_matrix_type; typedef ComplexMatrix complex_matrix_type; + typedef DiagMatrix real_diag_matrix_type; + typedef ComplexDiagMatrix complex_diag_matrix_type; + typedef double real_elt_type; typedef Complex complex_elt_type;
--- a/liboctave/array/fCMatrix.cc Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/fCMatrix.cc Tue Feb 16 14:41:06 2016 -0500 @@ -47,7 +47,7 @@ #include "fCColVector.h" #include "fCRowVector.h" #include "schur.h" -#include "fCmplxSVD.h" +#include "svd.h" #include "functor.h" #include "lo-error.h" #include "lo-ieee.h" @@ -1141,7 +1141,7 @@ { FloatComplexMatrix retval; - FloatComplexSVD result (*this, SVD::economy); + svd<FloatComplexMatrix> result (*this, svd<FloatComplexMatrix>::economy); FloatDiagMatrix S = result.singular_values (); FloatComplexMatrix U = result.left_singular_matrix ();
--- a/liboctave/array/fCMatrix.h Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/fCMatrix.h Tue Feb 16 14:41:06 2016 -0500 @@ -50,6 +50,9 @@ typedef FloatMatrix real_matrix_type; typedef FloatComplexMatrix complex_matrix_type; + typedef FloatDiagMatrix real_diag_matrix_type; + typedef FloatComplexDiagMatrix complex_diag_matrix_type; + typedef float real_elt_type; typedef FloatComplex complex_elt_type;
--- a/liboctave/array/fDiagMatrix.h Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/fDiagMatrix.h Tue Feb 16 14:41:06 2016 -0500 @@ -37,9 +37,6 @@ OCTAVE_API FloatDiagMatrix : public MDiagArray2<float> { - friend class FloatSVD; - friend class FloatComplexSVD; - public: FloatDiagMatrix (void) : MDiagArray2<float> () { }
--- a/liboctave/array/fMatrix.cc Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/fMatrix.cc Tue Feb 16 14:41:06 2016 -0500 @@ -49,7 +49,7 @@ #include "f77-fcn.h" #include "fMatrix.h" #include "schur.h" -#include "floatSVD.h" +#include "svd.h" #include "functor.h" #include "lo-error.h" #include "lo-ieee.h" @@ -832,7 +832,7 @@ FloatMatrix FloatMatrix::pseudo_inverse (float tol) const { - FloatSVD result (*this, SVD::economy); + svd<FloatMatrix> result (*this, svd<FloatMatrix>::economy); FloatDiagMatrix S = result.singular_values (); FloatMatrix U = result.left_singular_matrix ();
--- a/liboctave/array/fMatrix.h Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/array/fMatrix.h Tue Feb 16 14:41:06 2016 -0500 @@ -49,6 +49,9 @@ typedef FloatMatrix real_matrix_type; typedef FloatComplexMatrix complex_matrix_type; + typedef FloatDiagMatrix real_diag_matrix_type; + typedef FloatComplexDiagMatrix complex_diag_matrix_type; + typedef float real_elt_type; typedef FloatComplex complex_elt_type;
--- a/liboctave/numeric/CmplxSVD.cc Tue Feb 16 14:39:52 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,209 +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 "CmplxSVD.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - Complex*, const octave_idx_type&, - double*, Complex*, const octave_idx_type&, - Complex*, const octave_idx_type&, Complex*, - const octave_idx_type&, double*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (zgesdd, ZGESDD) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - Complex*, const octave_idx_type&, - double*, Complex*, const octave_idx_type&, - Complex*, const octave_idx_type&, Complex*, - const octave_idx_type&, double*, - octave_idx_type *, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); -} - -ComplexMatrix -ComplexSVD::left_singular_matrix (void) const -{ - if (type_computed == SVD::sigma_only) - (*current_liboctave_error_handler) - ("ComplexSVD: U not computed because type == SVD::sigma_only"); - - return left_sm; -} - -ComplexMatrix -ComplexSVD::right_singular_matrix (void) const -{ - if (type_computed == SVD::sigma_only) - (*current_liboctave_error_handler) - ("ComplexSVD: V not computed because type == SVD::sigma_only"); - - return right_sm; -} - -octave_idx_type -ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type, - SVD::driver svd_driver) -{ - octave_idx_type info; - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - - ComplexMatrix atmp = a; - Complex *tmp_data = atmp.fortran_vec (); - - octave_idx_type min_mn = m < n ? m : n; - octave_idx_type max_mn = m > n ? m : n; - - char jobu = 'A'; - char jobv = 'A'; - - octave_idx_type ncol_u = m; - octave_idx_type nrow_vt = n; - octave_idx_type nrow_s = m; - octave_idx_type ncol_s = n; - - switch (svd_type) - { - case SVD::economy: - jobu = jobv = 'S'; - ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; - break; - - case SVD::sigma_only: - - // Note: for this case, both jobu and jobv should be 'N', but - // there seems to be a bug in dgesvd from Lapack V2.0. To - // demonstrate the bug, set both jobu and jobv to 'N' and find - // the singular values of [eye(3), eye(3)]. The result is - // [-sqrt(2), -sqrt(2), -sqrt(2)]. - // - // For Lapack 3.0, this problem seems to be fixed. - - jobu = jobv = 'N'; - ncol_u = nrow_vt = 1; - break; - - default: - break; - } - - type_computed = svd_type; - - if (! (jobu == 'N' || jobu == 'O')) - left_sm.resize (m, ncol_u); - - Complex *u = left_sm.fortran_vec (); - - sigma.resize (nrow_s, ncol_s); - double *s_vec = sigma.fortran_vec (); - - if (! (jobv == 'N' || jobv == 'O')) - right_sm.resize (nrow_vt, n); - - Complex *vt = right_sm.fortran_vec (); - - // Query ZGESVD for the correct dimension of WORK. - - octave_idx_type lwork = -1; - - Array<Complex> work (dim_vector (1, 1)); - - octave_idx_type one = 1; - octave_idx_type m1 = std::max (m, one); - octave_idx_type nrow_vt1 = std::max (nrow_vt, one); - - if (svd_driver == SVD::GESVD) - { - octave_idx_type lrwork = 5*max_mn; - Array<double> rwork (dim_vector (lrwork, 1)); - - F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, - rwork.fortran_vec (), info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - lwork = static_cast<octave_idx_type> (work(0).real ()); - work.resize (dim_vector (lwork, 1)); - - F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, - rwork.fortran_vec (), info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - } - else if (svd_driver == SVD::GESDD) - { - assert (jobu == jobv); - char jobz = jobu; - - octave_idx_type lrwork; - if (jobz == 'N') - lrwork = 7*min_mn; - else - lrwork = 5*min_mn*min_mn + 5*min_mn; - Array<double> rwork (dim_vector (lrwork, 1)); - - OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); - - F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, - rwork.fortran_vec (), iwork, info - F77_CHAR_ARG_LEN (1))); - - lwork = static_cast<octave_idx_type> (work(0).real ()); - work.resize (dim_vector (lwork, 1)); - - F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, - rwork.fortran_vec (), iwork, info - F77_CHAR_ARG_LEN (1))); - } - else - assert (0); // impossible - - if (! (jobv == 'N' || jobv == 'O')) - right_sm = right_sm.hermitian (); - - return info; -}
--- a/liboctave/numeric/CmplxSVD.h Tue Feb 16 14:39:52 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,99 +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_CmplxSVD_h) -#define octave_CmplxSVD_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "dDiagMatrix.h" -#include "CMatrix.h" -#include "dbleSVD.h" - -class -OCTAVE_API -ComplexSVD -{ -public: - - ComplexSVD (void) - : type_computed (), sigma (), left_sm (), right_sm () - { } - - ComplexSVD (const ComplexMatrix& a, SVD::type svd_type = SVD::std, - SVD::driver svd_driver = SVD::GESVD) - : type_computed (), sigma (), left_sm (), right_sm () - { - init (a, svd_type, svd_driver); - } - - ComplexSVD (const ComplexMatrix& a, octave_idx_type& info, - SVD::type svd_type = SVD::std, - SVD::driver svd_driver = SVD::GESVD) - : type_computed (), sigma (), left_sm (), right_sm () - { - info = init (a, svd_type, svd_driver); - } - - ComplexSVD (const ComplexSVD& a) - : type_computed (a.type_computed), sigma (a.sigma), - left_sm (a.left_sm), right_sm (a.right_sm) - { } - - ComplexSVD& operator = (const ComplexSVD& a) - { - if (this != &a) - { - type_computed = a.type_computed; - sigma = a.sigma; - left_sm = a.left_sm; - right_sm = a.right_sm; - } - return *this; - } - - ~ComplexSVD (void) { } - - DiagMatrix singular_values (void) const { return sigma; } - - ComplexMatrix left_singular_matrix (void) const; - - ComplexMatrix right_singular_matrix (void) const; - - friend std::ostream& operator << (std::ostream& os, const ComplexSVD& a); - -private: - - SVD::type type_computed; - - DiagMatrix sigma; - ComplexMatrix left_sm; - ComplexMatrix right_sm; - - octave_idx_type init (const ComplexMatrix& a, - SVD::type svd_type = SVD::std, - SVD::driver svd_driver = SVD::GESVD); -}; - -#endif
--- a/liboctave/numeric/dbleSVD.cc Tue Feb 16 14:39:52 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,205 +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 "dbleSVD.h" -#include "f77-fcn.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - double*, const octave_idx_type&, double*, - double*, const octave_idx_type&, 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 (dgesdd, DGESDD) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - double*, const octave_idx_type&, double*, - double*, const octave_idx_type&, double*, - const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type *, - octave_idx_type& - F77_CHAR_ARG_LEN_DECL); -} - -Matrix -SVD::left_singular_matrix (void) const -{ - if (type_computed == SVD::sigma_only) - (*current_liboctave_error_handler) - ("SVD: U not computed because type == SVD::sigma_only"); - - return left_sm; -} - -Matrix -SVD::right_singular_matrix (void) const -{ - if (type_computed == SVD::sigma_only) - (*current_liboctave_error_handler) - ("SVD: V not computed because type == SVD::sigma_only"); - - return right_sm; -} - -octave_idx_type -SVD::init (const Matrix& a, SVD::type svd_type, SVD::driver svd_driver) -{ - octave_idx_type info; - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - - Matrix atmp = a; - double *tmp_data = atmp.fortran_vec (); - - octave_idx_type min_mn = m < n ? m : n; - - char jobu = 'A'; - char jobv = 'A'; - - octave_idx_type ncol_u = m; - octave_idx_type nrow_vt = n; - octave_idx_type nrow_s = m; - octave_idx_type ncol_s = n; - - switch (svd_type) - { - case SVD::economy: - jobu = jobv = 'S'; - ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; - break; - - case SVD::sigma_only: - - // Note: for this case, both jobu and jobv should be 'N', but - // there seems to be a bug in dgesvd from Lapack V2.0. To - // demonstrate the bug, set both jobu and jobv to 'N' and find - // the singular values of [eye(3), eye(3)]. The result is - // [-sqrt(2), -sqrt(2), -sqrt(2)]. - // - // For Lapack 3.0, this problem seems to be fixed. - - jobu = jobv = 'N'; - ncol_u = nrow_vt = 1; - break; - - default: - break; - } - - type_computed = svd_type; - - if (! (jobu == 'N' || jobu == 'O')) - left_sm.resize (m, ncol_u); - - double *u = left_sm.fortran_vec (); - - sigma.resize (nrow_s, ncol_s); - double *s_vec = sigma.fortran_vec (); - - if (! (jobv == 'N' || jobv == 'O')) - right_sm.resize (nrow_vt, n); - - double *vt = right_sm.fortran_vec (); - - // Query DGESVD for the correct dimension of WORK. - - octave_idx_type lwork = -1; - - Array<double> work (dim_vector (1, 1)); - - octave_idx_type one = 1; - octave_idx_type m1 = std::max (m, one); - octave_idx_type nrow_vt1 = std::max (nrow_vt, one); - - if (svd_driver == SVD::GESVD) - { - F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (dim_vector (lwork, 1)); - - F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - } - else if (svd_driver == SVD::GESDD) - { - assert (jobu == jobv); - char jobz = jobu; - OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); - - F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, - work.fortran_vec (), lwork, iwork, info - F77_CHAR_ARG_LEN (1))); - - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (dim_vector (lwork, 1)); - - F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, - work.fortran_vec (), lwork, iwork, info - F77_CHAR_ARG_LEN (1))); - - } - else - assert (0); // impossible - - if (! (jobv == 'N' || jobv == 'O')) - right_sm = right_sm.transpose (); - - return info; -} - -std::ostream& -operator << (std::ostream& os, const SVD& a) -{ - os << a.left_singular_matrix () << "\n"; - os << a.singular_values () << "\n"; - os << a.right_singular_matrix () << "\n"; - - return os; -}
--- a/liboctave/numeric/dbleSVD.h Tue Feb 16 14:39:52 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,108 +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_dbleSVD_h) -#define octave_dbleSVD_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "dDiagMatrix.h" -#include "dMatrix.h" - -class -OCTAVE_API -SVD -{ -public: - - enum type - { - std, - economy, - sigma_only - }; - - enum driver - { - GESVD, - GESDD - }; - - SVD (void) : type_computed (), sigma (), left_sm (), right_sm () { } - - SVD (const Matrix& a, - type svd_type = SVD::std, driver svd_driver = SVD::GESVD) - : type_computed (), sigma (), left_sm (), right_sm () - { - init (a, svd_type, svd_driver); - } - - SVD (const Matrix& a, octave_idx_type& info, - type svd_type = SVD::std, driver svd_driver = SVD::GESVD) - : type_computed (), sigma (), left_sm (), right_sm () - { - info = init (a, svd_type, svd_driver); - } - - SVD (const SVD& a) - : type_computed (a.type_computed), sigma (a.sigma), - left_sm (a.left_sm), right_sm (a.right_sm) - { } - - SVD& operator = (const SVD& a) - { - if (this != &a) - { - type_computed = a.type_computed; - sigma = a.sigma; - left_sm = a.left_sm; - right_sm = a.right_sm; - } - - return *this; - } - - ~SVD (void) { } - - DiagMatrix singular_values (void) const { return sigma; } - - Matrix left_singular_matrix (void) const; - - Matrix right_singular_matrix (void) const; - - friend std::ostream& operator << (std::ostream& os, const SVD& a); - -private: - - SVD::type type_computed; - - DiagMatrix sigma; - Matrix left_sm; - Matrix right_sm; - - octave_idx_type init (const Matrix& a, - type svd_type = std, driver svd_driver = GESVD); -}; - -#endif
--- a/liboctave/numeric/fCmplxSVD.cc Tue Feb 16 14:39:52 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,216 +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 "fCmplxSVD.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, float*, - FloatComplex*, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, - float*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (cgesdd, CGESDD) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, float*, - FloatComplex*, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, - float*, octave_idx_type *, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); -} - -FloatComplexMatrix -FloatComplexSVD::left_singular_matrix (void) const -{ - if (type_computed == SVD::sigma_only) - { - (*current_liboctave_error_handler) - ("FloatComplexSVD: U not computed because type == SVD::sigma_only"); - return FloatComplexMatrix (); - } - else - return left_sm; -} - -FloatComplexMatrix -FloatComplexSVD::right_singular_matrix (void) const -{ - if (type_computed == SVD::sigma_only) - { - (*current_liboctave_error_handler) - ("FloatComplexSVD: V not computed because type == SVD::sigma_only"); - return FloatComplexMatrix (); - } - else - return right_sm; -} - -octave_idx_type -FloatComplexSVD::init (const FloatComplexMatrix& a, SVD::type svd_type, - SVD::driver svd_driver) -{ - octave_idx_type info; - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - - FloatComplexMatrix atmp = a; - FloatComplex *tmp_data = atmp.fortran_vec (); - - octave_idx_type min_mn = m < n ? m : n; - octave_idx_type max_mn = m > n ? m : n; - - char jobu = 'A'; - char jobv = 'A'; - - octave_idx_type ncol_u = m; - octave_idx_type nrow_vt = n; - octave_idx_type nrow_s = m; - octave_idx_type ncol_s = n; - - switch (svd_type) - { - case SVD::economy: - jobu = jobv = 'S'; - ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; - break; - - case SVD::sigma_only: - - // Note: for this case, both jobu and jobv should be 'N', but - // there seems to be a bug in dgesvd from Lapack V2.0. To - // demonstrate the bug, set both jobu and jobv to 'N' and find - // the singular values of [eye(3), eye(3)]. The result is - // [-sqrt(2), -sqrt(2), -sqrt(2)]. - // - // For Lapack 3.0, this problem seems to be fixed. - - jobu = jobv = 'N'; - ncol_u = nrow_vt = 1; - break; - - default: - break; - } - - type_computed = svd_type; - - if (! (jobu == 'N' || jobu == 'O')) - left_sm.resize (m, ncol_u); - - FloatComplex *u = left_sm.fortran_vec (); - - sigma.resize (nrow_s, ncol_s); - float *s_vec = sigma.fortran_vec (); - - if (! (jobv == 'N' || jobv == 'O')) - right_sm.resize (nrow_vt, n); - - FloatComplex *vt = right_sm.fortran_vec (); - - // Query CGESVD for the correct dimension of WORK. - - octave_idx_type lwork = -1; - - Array<FloatComplex> work (dim_vector (1, 1)); - - octave_idx_type one = 1; - octave_idx_type m1 = std::max (m, one); - octave_idx_type nrow_vt1 = std::max (nrow_vt, one); - - if (svd_driver == SVD::GESVD) - { - octave_idx_type lrwork = 5*max_mn; - Array<float> rwork (dim_vector (lrwork, 1)); - - F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, - rwork.fortran_vec (), info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - lwork = static_cast<octave_idx_type> (work(0).real ()); - work.resize (dim_vector (lwork, 1)); - - F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, - rwork.fortran_vec (), info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - } - else if (svd_driver == SVD::GESDD) - { - assert (jobu == jobv); - char jobz = jobu; - - octave_idx_type lrwork; - if (jobz == 'N') - lrwork = 5*min_mn; - else - lrwork = min_mn * std::max (5*min_mn+7, 2*max_mn+2*min_mn+1); - Array<float> rwork (dim_vector (lrwork, 1)); - - OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); - - F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, - rwork.fortran_vec (), iwork, info - F77_CHAR_ARG_LEN (1))); - - lwork = static_cast<octave_idx_type> (work(0).real ()); - work.resize (dim_vector (lwork, 1)); - - F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, - rwork.fortran_vec (), iwork, info - F77_CHAR_ARG_LEN (1))); - } - else - assert (0); // impossible - - if (! (jobv == 'N' || jobv == 'O')) - right_sm = right_sm.hermitian (); - - return info; -}
--- a/liboctave/numeric/fCmplxSVD.h Tue Feb 16 14:39:52 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,101 +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_fCmplxSVD_h) -#define octave_fCmplxSVD_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "fDiagMatrix.h" -#include "fCMatrix.h" -#include "dbleSVD.h" - -class -OCTAVE_API -FloatComplexSVD -{ -public: - - FloatComplexSVD (void) - : type_computed (), sigma (), left_sm (), right_sm () - { } - - FloatComplexSVD (const FloatComplexMatrix& a, - SVD::type svd_type = SVD::std, - SVD::driver svd_driver = SVD::GESVD) - : type_computed (), sigma (), left_sm (), right_sm () - { - init (a, svd_type, svd_driver); - } - - FloatComplexSVD (const FloatComplexMatrix& a, octave_idx_type& info, - SVD::type svd_type = SVD::std, - SVD::driver svd_driver = SVD::GESVD) - : type_computed (), sigma (), left_sm (), right_sm () - { - info = init (a, svd_type, svd_driver); - } - - FloatComplexSVD (const FloatComplexSVD& a) - : type_computed (a.type_computed), sigma (a.sigma), - left_sm (a.left_sm), right_sm (a.right_sm) - { } - - FloatComplexSVD& operator = (const FloatComplexSVD& a) - { - if (this != &a) - { - type_computed = a.type_computed; - sigma = a.sigma; - left_sm = a.left_sm; - right_sm = a.right_sm; - } - return *this; - } - - ~FloatComplexSVD (void) { } - - FloatDiagMatrix singular_values (void) const { return sigma; } - - FloatComplexMatrix left_singular_matrix (void) const; - - FloatComplexMatrix right_singular_matrix (void) const; - - friend std::ostream& operator << (std::ostream& os, - const FloatComplexSVD& a); - -private: - - SVD::type type_computed; - - FloatDiagMatrix sigma; - FloatComplexMatrix left_sm; - FloatComplexMatrix right_sm; - - octave_idx_type init (const FloatComplexMatrix& a, - SVD::type svd_type = SVD::std, - SVD::driver svd_driver = SVD::GESVD); -}; - -#endif
--- a/liboctave/numeric/floatSVD.cc Tue Feb 16 14:39:52 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,206 +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 "floatSVD.h" -#include "f77-fcn.h" -#include "oct-locbuf.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (sgesvd, SGESVD) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - float*, const octave_idx_type&, float*, - float*, const octave_idx_type&, 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 (sgesdd, SGESDD) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - float*, const octave_idx_type&, float*, - float*, const octave_idx_type&, float*, - const octave_idx_type&, float*, - const octave_idx_type&, octave_idx_type *, - octave_idx_type& - F77_CHAR_ARG_LEN_DECL); -} - -FloatMatrix -FloatSVD::left_singular_matrix (void) const -{ - if (type_computed == SVD::sigma_only) - (*current_liboctave_error_handler) - ("FloatSVD: U not computed because type == SVD::sigma_only"); - - return left_sm; -} - -FloatMatrix -FloatSVD::right_singular_matrix (void) const -{ - if (type_computed == SVD::sigma_only) - (*current_liboctave_error_handler) - ("FloatSVD: V not computed because type == SVD::sigma_only"); - - return right_sm; -} - -octave_idx_type -FloatSVD::init (const FloatMatrix& a, SVD::type svd_type, - SVD::driver svd_driver) -{ - octave_idx_type info; - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - - FloatMatrix atmp = a; - float *tmp_data = atmp.fortran_vec (); - - octave_idx_type min_mn = m < n ? m : n; - - char jobu = 'A'; - char jobv = 'A'; - - octave_idx_type ncol_u = m; - octave_idx_type nrow_vt = n; - octave_idx_type nrow_s = m; - octave_idx_type ncol_s = n; - - switch (svd_type) - { - case SVD::economy: - jobu = jobv = 'S'; - ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; - break; - - case SVD::sigma_only: - - // Note: for this case, both jobu and jobv should be 'N', but - // there seems to be a bug in dgesvd from Lapack V2.0. To - // demonstrate the bug, set both jobu and jobv to 'N' and find - // the singular values of [eye(3), eye(3)]. The result is - // [-sqrt(2), -sqrt(2), -sqrt(2)]. - // - // For Lapack 3.0, this problem seems to be fixed. - - jobu = jobv = 'N'; - ncol_u = nrow_vt = 1; - break; - - default: - break; - } - - type_computed = svd_type; - - if (! (jobu == 'N' || jobu == 'O')) - left_sm.resize (m, ncol_u); - - float *u = left_sm.fortran_vec (); - - sigma.resize (nrow_s, ncol_s); - float *s_vec = sigma.fortran_vec (); - - if (! (jobv == 'N' || jobv == 'O')) - right_sm.resize (nrow_vt, n); - - float *vt = right_sm.fortran_vec (); - - // Query SGESVD for the correct dimension of WORK. - - octave_idx_type lwork = -1; - - Array<float> work (dim_vector (1, 1)); - - octave_idx_type one = 1; - octave_idx_type m1 = std::max (m, one); - octave_idx_type nrow_vt1 = std::max (nrow_vt, one); - - if (svd_driver == SVD::GESVD) - { - F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (dim_vector (lwork, 1)); - - F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, - nrow_vt1, work.fortran_vec (), lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - } - else if (svd_driver == SVD::GESDD) - { - assert (jobu == jobv); - char jobz = jobu; - OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); - - F77_XFCN (sgesdd, SGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, - work.fortran_vec (), lwork, iwork, info - F77_CHAR_ARG_LEN (1))); - - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (dim_vector (lwork, 1)); - - F77_XFCN (sgesdd, SGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), - m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, - work.fortran_vec (), lwork, iwork, info - F77_CHAR_ARG_LEN (1))); - - } - else - assert (0); // impossible - - if (! (jobv == 'N' || jobv == 'O')) - right_sm = right_sm.transpose (); - - return info; -} - -std::ostream& -operator << (std::ostream& os, const FloatSVD& a) -{ - os << a.left_singular_matrix () << "\n"; - os << a.singular_values () << "\n"; - os << a.right_singular_matrix () << "\n"; - - return os; -}
--- a/liboctave/numeric/floatSVD.h Tue Feb 16 14:39:52 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,98 +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_floatSVD_h) -#define octave_floatSVD_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "fDiagMatrix.h" -#include "fMatrix.h" -#include "dbleSVD.h" - -class -OCTAVE_API -FloatSVD -{ -public: - - FloatSVD (void) : type_computed (), sigma (), left_sm (), right_sm () { } - - FloatSVD (const FloatMatrix& a, - SVD::type svd_type = SVD::std, SVD::driver svd_driver = SVD::GESVD) - : type_computed (), sigma (), left_sm (), right_sm () - { - init (a, svd_type, svd_driver); - } - - FloatSVD (const FloatMatrix& a, octave_idx_type& info, - SVD::type svd_type = SVD::std, - SVD::driver svd_driver = SVD::GESVD) - : type_computed (), sigma (), left_sm (), right_sm () - { - info = init (a, svd_type, svd_driver); - } - - FloatSVD (const FloatSVD& a) - : type_computed (a.type_computed), sigma (a.sigma), - left_sm (a.left_sm), right_sm (a.right_sm) - { } - - FloatSVD& operator = (const FloatSVD& a) - { - if (this != &a) - { - type_computed = a.type_computed; - sigma = a.sigma; - left_sm = a.left_sm; - right_sm = a.right_sm; - } - - return *this; - } - - ~FloatSVD (void) { } - - FloatDiagMatrix singular_values (void) const { return sigma; } - - FloatMatrix left_singular_matrix (void) const; - - FloatMatrix right_singular_matrix (void) const; - - friend std::ostream& operator << (std::ostream& os, const FloatSVD& a); - -private: - - SVD::type type_computed; - - FloatDiagMatrix sigma; - FloatMatrix left_sm; - FloatMatrix right_sm; - - octave_idx_type init (const FloatMatrix& a, - SVD::type svd_type = SVD::std, - SVD::driver svd_driver = SVD::GESVD); -}; - -#endif
--- a/liboctave/numeric/module.mk Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/numeric/module.mk Tue Feb 16 14:41:06 2016 -0500 @@ -10,7 +10,6 @@ NUMERIC_INC = \ liboctave/numeric/CmplxQR.h \ liboctave/numeric/CmplxQRP.h \ - liboctave/numeric/CmplxSVD.h \ liboctave/numeric/CollocWt.h \ liboctave/numeric/DAE.h \ liboctave/numeric/DAEFunc.h \ @@ -37,15 +36,12 @@ liboctave/numeric/chol.h \ liboctave/numeric/dbleQR.h \ liboctave/numeric/dbleQRP.h \ - liboctave/numeric/dbleSVD.h \ liboctave/numeric/eigs-base.h \ liboctave/numeric/fCmplxQR.h \ liboctave/numeric/fCmplxQRP.h \ - liboctave/numeric/fCmplxSVD.h \ liboctave/numeric/fEIG.h \ liboctave/numeric/floatQR.h \ liboctave/numeric/floatQRP.h \ - liboctave/numeric/floatSVD.h \ liboctave/numeric/gepbalance.h \ liboctave/numeric/hess.h \ liboctave/numeric/lo-mappers.h \ @@ -63,7 +59,8 @@ liboctave/numeric/sparse-chol.h \ liboctave/numeric/sparse-dmsolve.h \ liboctave/numeric/sparse-lu.h \ - liboctave/numeric/sparse-qr.h + liboctave/numeric/sparse-qr.h \ + liboctave/numeric/svd.h NUMERIC_C_SRC = \ liboctave/numeric/randgamma.c \ @@ -73,7 +70,6 @@ NUMERIC_SRC = \ liboctave/numeric/CmplxQR.cc \ liboctave/numeric/CmplxQRP.cc \ - liboctave/numeric/CmplxSVD.cc \ liboctave/numeric/CollocWt.cc \ liboctave/numeric/DASPK.cc \ liboctave/numeric/DASRT.cc \ @@ -86,15 +82,12 @@ liboctave/numeric/chol.cc \ liboctave/numeric/dbleQR.cc \ liboctave/numeric/dbleQRP.cc \ - liboctave/numeric/dbleSVD.cc \ liboctave/numeric/eigs-base.cc \ liboctave/numeric/fCmplxQR.cc \ liboctave/numeric/fCmplxQRP.cc \ - liboctave/numeric/fCmplxSVD.cc \ liboctave/numeric/fEIG.cc \ liboctave/numeric/floatQR.cc \ liboctave/numeric/floatQRP.cc \ - liboctave/numeric/floatSVD.cc \ liboctave/numeric/gepbalance.cc \ liboctave/numeric/hess.cc \ liboctave/numeric/lo-mappers.cc \ @@ -110,6 +103,7 @@ liboctave/numeric/sparse-dmsolve.cc \ liboctave/numeric/sparse-lu.cc \ liboctave/numeric/sparse-qr.cc \ + liboctave/numeric/svd.cc \ $(NUMERIC_C_SRC) LIBOCTAVE_TEMPLATE_SRC += \
--- a/liboctave/numeric/oct-norm.cc Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/numeric/oct-norm.cc Tue Feb 16 14:41:06 2016 -0500 @@ -33,33 +33,32 @@ #include <iostream> #include <vector> -#include "oct-cmplx.h" +#include "Array-util.h" +#include "Array.h" +#include "CColVector.h" +#include "CMatrix.h" +#include "CRowVector.h" +#include "CSparse.h" +#include "dColVector.h" +#include "dDiagMatrix.h" +#include "dMatrix.h" +#include "dRowVector.h" +#include "dSparse.h" +#include "fCColVector.h" +#include "fCMatrix.h" +#include "fCRowVector.h" +#include "fColVector.h" +#include "fDiagMatrix.h" +#include "fMatrix.h" +#include "fRowVector.h" #include "lo-error.h" #include "lo-ieee.h" #include "mx-cm-s.h" -#include "mx-s-cm.h" #include "mx-fcm-fs.h" #include "mx-fs-fcm.h" -#include "Array.h" -#include "Array-util.h" -#include "CMatrix.h" -#include "dMatrix.h" -#include "fCMatrix.h" -#include "fMatrix.h" -#include "CColVector.h" -#include "dColVector.h" -#include "CRowVector.h" -#include "dRowVector.h" -#include "fCColVector.h" -#include "fColVector.h" -#include "fCRowVector.h" -#include "fRowVector.h" -#include "CSparse.h" -#include "dSparse.h" -#include "dbleSVD.h" -#include "CmplxSVD.h" -#include "floatSVD.h" -#include "fCmplxSVD.h" +#include "mx-s-cm.h" +#include "oct-cmplx.h" +#include "svd.h" // Theory: norm accumulator is an object that has an accum method able // to handle both real and complex element, and a cast operator @@ -475,14 +474,14 @@ static int max_norm_iter = 100; // version with SVD for dense matrices -template <typename MatrixT, typename VectorT, typename SVDT, typename R> -R matrix_norm (const MatrixT& m, R p, VectorT, SVDT) +template <typename MatrixT, typename VectorT, typename R> +R svd_matrix_norm (const MatrixT& m, R p, VectorT) { R res = 0; if (p == 2) { - SVDT svd (m, SVD::sigma_only); - res = svd.singular_values () (0,0); + svd<MatrixT> fact (m, svd<MatrixT>::sigma_only); + res = fact.singular_values () (0,0); } else if (p == 1) res = xcolnorms (m, 1).max (); @@ -529,7 +528,7 @@ OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \ { return vector_norm (x, p); } \ OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \ - { return matrix_norm (x, p, PREFIX##Matrix (), PREFIX##SVD ()); } \ + { return svd_matrix_norm (x, p, PREFIX##Matrix ()); } \ OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \ { return vector_norm (x, static_cast<RTYPE> (2)); }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/svd.cc Tue Feb 16 14:41:06 2016 -0500 @@ -0,0 +1,722 @@ +/* + +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 <cassert> + +#include "CMatrix.h" +#include "dDiagMatrix.h" +#include "fDiagMatrix.h" +#include "dMatrix.h" +#include "f77-fcn.h" +#include "fCMatrix.h" +#include "fMatrix.h" +#include "lo-error.h" +#include "oct-locbuf.h" +#include "svd.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + double*, const octave_idx_type&, double*, + double*, const octave_idx_type&, 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 (dgesdd, DGESDD) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + double*, const octave_idx_type&, double*, + double*, const octave_idx_type&, double*, + const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type *, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (sgesvd, SGESVD) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + float*, const octave_idx_type&, float*, + float*, const octave_idx_type&, 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 (sgesdd, SGESDD) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + float*, const octave_idx_type&, float*, + float*, const octave_idx_type&, float*, + const octave_idx_type&, float*, + const octave_idx_type&, octave_idx_type *, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + Complex*, const octave_idx_type&, + double*, Complex*, const octave_idx_type&, + Complex*, const octave_idx_type&, Complex*, + const octave_idx_type&, double*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zgesdd, ZGESDD) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + Complex*, const octave_idx_type&, + double*, Complex*, const octave_idx_type&, + Complex*, const octave_idx_type&, Complex*, + const octave_idx_type&, double*, + octave_idx_type *, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + F77_RET_T + F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL, + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, float*, + FloatComplex*, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + float*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (cgesdd, CGESDD) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, float*, + FloatComplex*, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + FloatComplex*, const octave_idx_type&, + float*, octave_idx_type *, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); +} + +template <typename T> +T +svd<T>::left_singular_matrix (void) const +{ + if (type_computed == svd::sigma_only) + (*current_liboctave_error_handler) + ("svd: U not computed because type == svd::sigma_only"); + + return left_sm; +} + +template <typename T> +T +svd<T>::right_singular_matrix (void) const +{ + if (type_computed == svd::sigma_only) + (*current_liboctave_error_handler) + ("svd: V not computed because type == svd::sigma_only"); + + return right_sm; +} + +template <typename T> +octave_idx_type +svd<T>::empty_init (octave_idx_type nr, octave_idx_type nc, svd::type svd_type) +{ + assert (nr == 0 || nc == 0); + + static typename T::element_type zero (0); + static typename T::element_type one (1); + + switch (svd_type) + { + case svd::std: + left_sm = T (nr, nr, zero); + for (octave_idx_type i = 0; i < nr; i++) + left_sm.xelem (i, i) = one; + sigma = DM_T (nr, nc); + right_sm = T (nc, nc, zero); + for (octave_idx_type i = 0; i < nc; i++) + right_sm.xelem (i, i) = one; + break; + + case svd::economy: + left_sm = T (nr, 0, zero); + sigma = DM_T (0, 0); + right_sm = T (0, nc, zero); + break; + + case svd::sigma_only: + default: + sigma = DM_T (0, 1); + break; + } + + return 0; +} + +// Specializations. + +template <> +octave_idx_type +svd<Matrix>::init (const Matrix& a, svd::type svd_type, svd::driver svd_driver) +{ + octave_idx_type info = 0; + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + if (m == 0 || n == 0) + return empty_init (m, n, svd_type); + + Matrix atmp = a; + double *tmp_data = atmp.fortran_vec (); + + octave_idx_type min_mn = m < n ? m : n; + + char jobu = 'A'; + char jobv = 'A'; + + octave_idx_type ncol_u = m; + octave_idx_type nrow_vt = n; + octave_idx_type nrow_s = m; + octave_idx_type ncol_s = n; + + switch (svd_type) + { + case svd::economy: + jobu = jobv = 'S'; + ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; + break; + + case svd::sigma_only: + + // Note: for this case, both jobu and jobv should be 'N', but + // there seems to be a bug in dgesvd from Lapack V2.0. To + // demonstrate the bug, set both jobu and jobv to 'N' and find + // the singular values of [eye(3), eye(3)]. The result is + // [-sqrt(2), -sqrt(2), -sqrt(2)]. + // + // For Lapack 3.0, this problem seems to be fixed. + + jobu = jobv = 'N'; + ncol_u = nrow_vt = 1; + break; + + default: + break; + } + + type_computed = svd_type; + + if (! (jobu == 'N' || jobu == 'O')) + left_sm.resize (m, ncol_u); + + double *u = left_sm.fortran_vec (); + + sigma.resize (nrow_s, ncol_s); + double *s_vec = sigma.fortran_vec (); + + if (! (jobv == 'N' || jobv == 'O')) + right_sm.resize (nrow_vt, n); + + double *vt = right_sm.fortran_vec (); + + // Query DGESVD for the correct dimension of WORK. + + octave_idx_type lwork = -1; + + Array<double> work (dim_vector (1, 1)); + + octave_idx_type one = 1; + octave_idx_type m1 = std::max (m, one); + octave_idx_type nrow_vt1 = std::max (nrow_vt, one); + + if (svd_driver == svd::GESVD) + { + F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + lwork = static_cast<octave_idx_type> (work(0)); + work.resize (dim_vector (lwork, 1)); + + F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + } + else if (svd_driver == svd::GESDD) + { + assert (jobu == jobv); + char jobz = jobu; + OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); + + F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, + work.fortran_vec (), lwork, iwork, info + F77_CHAR_ARG_LEN (1))); + + lwork = static_cast<octave_idx_type> (work(0)); + work.resize (dim_vector (lwork, 1)); + + F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, + work.fortran_vec (), lwork, iwork, info + F77_CHAR_ARG_LEN (1))); + + } + else + abort (); + + if (! (jobv == 'N' || jobv == 'O')) + right_sm = right_sm.transpose (); + + return info; +} + +template <> +octave_idx_type +svd<FloatMatrix>::init (const FloatMatrix& a, svd::type svd_type, + svd::driver svd_driver) +{ + octave_idx_type info; + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + if (m == 0 || n == 0) + return empty_init (m, n, svd_type); + + FloatMatrix atmp = a; + float *tmp_data = atmp.fortran_vec (); + + octave_idx_type min_mn = m < n ? m : n; + + char jobu = 'A'; + char jobv = 'A'; + + octave_idx_type ncol_u = m; + octave_idx_type nrow_vt = n; + octave_idx_type nrow_s = m; + octave_idx_type ncol_s = n; + + switch (svd_type) + { + case svd::economy: + jobu = jobv = 'S'; + ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; + break; + + case svd::sigma_only: + + // Note: for this case, both jobu and jobv should be 'N', but + // there seems to be a bug in dgesvd from Lapack V2.0. To + // demonstrate the bug, set both jobu and jobv to 'N' and find + // the singular values of [eye(3), eye(3)]. The result is + // [-sqrt(2), -sqrt(2), -sqrt(2)]. + // + // For Lapack 3.0, this problem seems to be fixed. + + jobu = jobv = 'N'; + ncol_u = nrow_vt = 1; + break; + + default: + break; + } + + type_computed = svd_type; + + if (! (jobu == 'N' || jobu == 'O')) + left_sm.resize (m, ncol_u); + + float *u = left_sm.fortran_vec (); + + sigma.resize (nrow_s, ncol_s); + float *s_vec = sigma.fortran_vec (); + + if (! (jobv == 'N' || jobv == 'O')) + right_sm.resize (nrow_vt, n); + + float *vt = right_sm.fortran_vec (); + + // Query SGESVD for the correct dimension of WORK. + + octave_idx_type lwork = -1; + + Array<float> work (dim_vector (1, 1)); + + octave_idx_type one = 1; + octave_idx_type m1 = std::max (m, one); + octave_idx_type nrow_vt1 = std::max (nrow_vt, one); + + if (svd_driver == svd::GESVD) + { + F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + lwork = static_cast<octave_idx_type> (work(0)); + work.resize (dim_vector (lwork, 1)); + + F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + } + else if (svd_driver == svd::GESDD) + { + assert (jobu == jobv); + char jobz = jobu; + OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); + + F77_XFCN (sgesdd, SGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, + work.fortran_vec (), lwork, iwork, info + F77_CHAR_ARG_LEN (1))); + + lwork = static_cast<octave_idx_type> (work(0)); + work.resize (dim_vector (lwork, 1)); + + F77_XFCN (sgesdd, SGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, + work.fortran_vec (), lwork, iwork, info + F77_CHAR_ARG_LEN (1))); + + } + else + abort (); + + if (! (jobv == 'N' || jobv == 'O')) + right_sm = right_sm.transpose (); + + return info; +} + +template <> +octave_idx_type +svd<ComplexMatrix>::init (const ComplexMatrix& a, svd::type svd_type, + svd::driver svd_driver) +{ + octave_idx_type info; + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + if (m == 0 || n == 0) + return empty_init (m, n, svd_type); + + ComplexMatrix atmp = a; + Complex *tmp_data = atmp.fortran_vec (); + + octave_idx_type min_mn = m < n ? m : n; + octave_idx_type max_mn = m > n ? m : n; + + char jobu = 'A'; + char jobv = 'A'; + + octave_idx_type ncol_u = m; + octave_idx_type nrow_vt = n; + octave_idx_type nrow_s = m; + octave_idx_type ncol_s = n; + + switch (svd_type) + { + case svd::economy: + jobu = jobv = 'S'; + ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; + break; + + case svd::sigma_only: + + // Note: for this case, both jobu and jobv should be 'N', but + // there seems to be a bug in dgesvd from Lapack V2.0. To + // demonstrate the bug, set both jobu and jobv to 'N' and find + // the singular values of [eye(3), eye(3)]. The result is + // [-sqrt(2), -sqrt(2), -sqrt(2)]. + // + // For Lapack 3.0, this problem seems to be fixed. + + jobu = jobv = 'N'; + ncol_u = nrow_vt = 1; + break; + + default: + break; + } + + type_computed = svd_type; + + if (! (jobu == 'N' || jobu == 'O')) + left_sm.resize (m, ncol_u); + + Complex *u = left_sm.fortran_vec (); + + sigma.resize (nrow_s, ncol_s); + double *s_vec = sigma.fortran_vec (); + + if (! (jobv == 'N' || jobv == 'O')) + right_sm.resize (nrow_vt, n); + + Complex *vt = right_sm.fortran_vec (); + + // Query ZGESVD for the correct dimension of WORK. + + octave_idx_type lwork = -1; + + Array<Complex> work (dim_vector (1, 1)); + + octave_idx_type one = 1; + octave_idx_type m1 = std::max (m, one); + octave_idx_type nrow_vt1 = std::max (nrow_vt, one); + + if (svd_driver == svd::GESVD) + { + octave_idx_type lrwork = 5*max_mn; + Array<double> rwork (dim_vector (lrwork, 1)); + + F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, + rwork.fortran_vec (), info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + lwork = static_cast<octave_idx_type> (work(0).real ()); + work.resize (dim_vector (lwork, 1)); + + F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, + rwork.fortran_vec (), info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + else if (svd_driver == svd::GESDD) + { + assert (jobu == jobv); + char jobz = jobu; + + octave_idx_type lrwork; + if (jobz == 'N') + lrwork = 7*min_mn; + else + lrwork = 5*min_mn*min_mn + 5*min_mn; + Array<double> rwork (dim_vector (lrwork, 1)); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); + + F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, + rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1))); + + lwork = static_cast<octave_idx_type> (work(0).real ()); + work.resize (dim_vector (lwork, 1)); + + F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, + rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1))); + } + else + abort (); + + if (! (jobv == 'N' || jobv == 'O')) + right_sm = right_sm.hermitian (); + + return info; +} + +template <> +octave_idx_type +svd<FloatComplexMatrix>::init (const FloatComplexMatrix& a, svd::type svd_type, + svd::driver svd_driver) +{ + octave_idx_type info; + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + + if (m == 0 || n == 0) + return empty_init (m, n, svd_type); + + FloatComplexMatrix atmp = a; + FloatComplex *tmp_data = atmp.fortran_vec (); + + octave_idx_type min_mn = m < n ? m : n; + octave_idx_type max_mn = m > n ? m : n; + + char jobu = 'A'; + char jobv = 'A'; + + octave_idx_type ncol_u = m; + octave_idx_type nrow_vt = n; + octave_idx_type nrow_s = m; + octave_idx_type ncol_s = n; + + switch (svd_type) + { + case svd::economy: + jobu = jobv = 'S'; + ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; + break; + + case svd::sigma_only: + + // Note: for this case, both jobu and jobv should be 'N', but + // there seems to be a bug in dgesvd from Lapack V2.0. To + // demonstrate the bug, set both jobu and jobv to 'N' and find + // the singular values of [eye(3), eye(3)]. The result is + // [-sqrt(2), -sqrt(2), -sqrt(2)]. + // + // For Lapack 3.0, this problem seems to be fixed. + + jobu = jobv = 'N'; + ncol_u = nrow_vt = 1; + break; + + default: + break; + } + + type_computed = svd_type; + + if (! (jobu == 'N' || jobu == 'O')) + left_sm.resize (m, ncol_u); + + FloatComplex *u = left_sm.fortran_vec (); + + sigma.resize (nrow_s, ncol_s); + float *s_vec = sigma.fortran_vec (); + + if (! (jobv == 'N' || jobv == 'O')) + right_sm.resize (nrow_vt, n); + + FloatComplex *vt = right_sm.fortran_vec (); + + // Query CGESVD for the correct dimension of WORK. + + octave_idx_type lwork = -1; + + Array<FloatComplex> work (dim_vector (1, 1)); + + octave_idx_type one = 1; + octave_idx_type m1 = std::max (m, one); + octave_idx_type nrow_vt1 = std::max (nrow_vt, one); + + if (svd_driver == svd::GESVD) + { + octave_idx_type lrwork = 5*max_mn; + Array<float> rwork (dim_vector (lrwork, 1)); + + F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, + rwork.fortran_vec (), info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + lwork = static_cast<octave_idx_type> (work(0).real ()); + work.resize (dim_vector (lwork, 1)); + + F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, + rwork.fortran_vec (), info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + else if (svd_driver == svd::GESDD) + { + assert (jobu == jobv); + char jobz = jobu; + + octave_idx_type lrwork; + if (jobz == 'N') + lrwork = 5*min_mn; + else + lrwork = min_mn * std::max (5*min_mn+7, 2*max_mn+2*min_mn+1); + Array<float> rwork (dim_vector (lrwork, 1)); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); + + F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, + rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1))); + + lwork = static_cast<octave_idx_type> (work(0).real ()); + work.resize (dim_vector (lwork, 1)); + + F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, + rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1))); + } + else + abort (); + + if (! (jobv == 'N' || jobv == 'O')) + right_sm = right_sm.hermitian (); + + return info; +} + +// Instantiations we need. + +template class svd<Matrix>; + +template class svd<FloatMatrix>; + +template class svd<ComplexMatrix>; + +template class svd<FloatComplexMatrix>;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/svd.h Tue Feb 16 14:41:06 2016 -0500 @@ -0,0 +1,109 @@ +/* + +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_svd_h) +#define octave_svd_h 1 + +#include "octave-config.h" + +#include <iosfwd> + +template <typename T> +class +svd +{ +public: + + typedef typename T::real_diag_matrix_type DM_T; + + enum type + { + std, + economy, + sigma_only + }; + + enum driver + { + GESVD, + GESDD + }; + + svd (void) + : type_computed (), left_sm (), sigma (), right_sm () + { } + + svd (const T& a, type svd_type = svd::std, driver svd_driver = svd::GESVD) + : type_computed (), left_sm (), sigma (), right_sm () + { + init (a, svd_type, svd_driver); + } + + svd (const T& a, octave_idx_type& info, type svd_type = svd::std, + driver svd_driver = svd::GESVD) + : type_computed (), left_sm (), sigma (), right_sm () + { + info = init (a, svd_type, svd_driver); + } + + svd (const svd& a) + : type_computed (a.type_computed), left_sm (a.left_sm), + sigma (a.sigma), right_sm (a.right_sm) + { } + + svd& operator = (const svd& a) + { + if (this != &a) + { + type_computed = a.type_computed; + left_sm = a.left_sm; + sigma = a.sigma; + right_sm = a.right_sm; + } + + return *this; + } + + ~svd (void) { } + + T left_singular_matrix (void) const; + + DM_T singular_values (void) const { return sigma; } + + T right_singular_matrix (void) const; + +private: + + svd::type type_computed; + + T left_sm; + DM_T sigma; + T right_sm; + + octave_idx_type + init (const T& a, type svd_type, driver svd_driver); + + octave_idx_type + empty_init (octave_idx_type nr, octave_idx_type nc, type svd_type); +}; + +#endif
--- a/liboctave/operators/mx-defs.h Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/operators/mx-defs.h Tue Feb 16 14:41:06 2016 -0500 @@ -70,10 +70,7 @@ template <typename T> class schur; -class SVD; -class ComplexSVD; -class FloatSVD; -class FloatComplexSVD; +template <typename T> class svd; template <typename T> class lu;
--- a/liboctave/operators/mx-ext.h Tue Feb 16 14:39:52 2016 -0500 +++ b/liboctave/operators/mx-ext.h Tue Feb 16 14:41:06 2016 -0500 @@ -51,10 +51,7 @@ // Result of a Singular Value Decomposition. -#include "dbleSVD.h" -#include "CmplxSVD.h" -#include "floatSVD.h" -#include "fCmplxSVD.h" +#include "svd.h" // Result of an Eigenvalue computation.