Mercurial > octave-nkf
changeset 10601:3ce0c530a9c9
implement svd_driver
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Mon, 03 May 2010 13:21:35 +0200 |
| parents | 036bdc2d0af0 |
| children | 38eae0c3a003 |
| files | liboctave/ChangeLog liboctave/CmplxSVD.cc liboctave/CmplxSVD.h liboctave/dbleSVD.cc liboctave/dbleSVD.h liboctave/fCmplxSVD.cc liboctave/fCmplxSVD.h liboctave/floatSVD.cc liboctave/floatSVD.h src/ChangeLog src/DLD-FUNCTIONS/svd.cc |
| diffstat | 11 files changed, 307 insertions(+), 86 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/ChangeLog Mon May 03 13:21:35 2010 +0200 @@ -1,3 +1,15 @@ +2010-05-03 Jaroslav Hajek <highegg@gmail.com> + + * dbleSVD.h (SVD::driver): New enum. + (SVD::SVD, SVD::init): Add driver option. + * floatSVD.h (FloatSVD::FloatSVD, FloatSVD::init): Add driver option. + * CmplxSVD.h (ComplexSVD::ComplexSVD, ComplexSVD::init): Add driver option. + * fCmplxSVD.h (FloatComplexSVD::FloatComplexSVD, FloatComplexSVD::init): Add driver option. + * dbleSVD.cc (SVD::init): Optionally use xGESDD driver. + * floatSVD.cc (FloatSVD::init): Ditto. + * CmplxSVD.cc (ComplexSVD::init): Ditto. + * fCmplxSVD.cc (FloatComplexSVD::init): Ditto. + 2010-04-28 John W. Eaton <jwe@octave.org> * dim-vector.h (dim_vector (const octave_idx_type *, size_t)): Delete.
--- a/liboctave/CmplxSVD.cc Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/CmplxSVD.cc Mon May 03 13:21:35 2010 +0200 @@ -28,6 +28,7 @@ #include "CmplxSVD.h" #include "f77-fcn.h" #include "lo-error.h" +#include "oct-locbuf.h" extern "C" { @@ -40,6 +41,14 @@ 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 @@ -69,7 +78,7 @@ } octave_idx_type -ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type) +ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type, SVD::driver svd_driver) { octave_idx_type info; @@ -144,24 +153,50 @@ octave_idx_type one = 1; octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one); - 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))); + if (svd_driver == SVD::GESVD) + { + 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 (lwork, 1); - lwork = static_cast<octave_idx_type> (work(0).real ()); - work.resize (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_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); - 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))); + 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 (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 ();
--- a/liboctave/CmplxSVD.h Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/CmplxSVD.h Mon May 03 13:21:35 2010 +0200 @@ -38,15 +38,16 @@ ComplexSVD (void) { } - ComplexSVD (const ComplexMatrix& a, SVD::type svd_type = SVD::std) + ComplexSVD (const ComplexMatrix& a, + SVD::type svd_type = SVD::std, SVD::driver svd_driver = SVD::GESVD) { - init (a, svd_type); + { init (a, svd_type, svd_driver); } } ComplexSVD (const ComplexMatrix& a, octave_idx_type& info, - SVD::type svd_type = SVD::std) + SVD::type svd_type = SVD::std, SVD::driver svd_driver = SVD::GESVD) { - info = init (a, svd_type); + info = init (a, svd_type, svd_driver); } ComplexSVD (const ComplexSVD& a) @@ -83,7 +84,9 @@ ComplexMatrix left_sm; ComplexMatrix right_sm; - octave_idx_type init (const ComplexMatrix& a, SVD::type svd_type = SVD::std); + octave_idx_type init (const ComplexMatrix& a, + SVD::type svd_type = SVD::std, + SVD::driver svd_driver = SVD::GESVD); }; #endif
--- a/liboctave/dbleSVD.cc Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/dbleSVD.cc Mon May 03 13:21:35 2010 +0200 @@ -29,6 +29,7 @@ #include "dbleSVD.h" #include "f77-fcn.h" +#include "oct-locbuf.h" extern "C" { @@ -41,6 +42,14 @@ 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 @@ -70,7 +79,7 @@ } octave_idx_type -SVD::init (const Matrix& a, SVD::type svd_type) +SVD::init (const Matrix& a, SVD::type svd_type, SVD::driver svd_driver) { octave_idx_type info; @@ -140,22 +149,48 @@ octave_idx_type one = 1; octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one); - 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))); + 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 (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))); - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (lwork, 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))); - 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 (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 ();
--- a/liboctave/dbleSVD.h Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/dbleSVD.h Mon May 03 13:21:35 2010 +0200 @@ -42,13 +42,22 @@ sigma_only }; + enum driver + { + GESVD, + GESDD + }; + SVD (void) : sigma (), left_sm (), right_sm () { } - SVD (const Matrix& a, type svd_type = SVD::std) { init (a, svd_type); } + SVD (const Matrix& a, + type svd_type = SVD::std, driver svd_driver = SVD::GESVD) + { init (a, svd_type, svd_driver); } - SVD (const Matrix& a, octave_idx_type& info, type svd_type = SVD::std) + SVD (const Matrix& a, octave_idx_type& info, + type svd_type = SVD::std, driver svd_driver = SVD::GESVD) { - info = init (a, svd_type); + info = init (a, svd_type, svd_driver); } SVD (const SVD& a) @@ -86,7 +95,8 @@ Matrix left_sm; Matrix right_sm; - octave_idx_type init (const Matrix& a, type svd_type = std); + octave_idx_type init (const Matrix& a, + type svd_type = std, driver svd_driver = GESVD); }; #endif
--- a/liboctave/fCmplxSVD.cc Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/fCmplxSVD.cc Mon May 03 13:21:35 2010 +0200 @@ -28,6 +28,7 @@ #include "fCmplxSVD.h" #include "f77-fcn.h" #include "lo-error.h" +#include "oct-locbuf.h" extern "C" { @@ -40,6 +41,14 @@ 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 @@ -69,7 +78,7 @@ } octave_idx_type -FloatComplexSVD::init (const FloatComplexMatrix& a, SVD::type svd_type) +FloatComplexSVD::init (const FloatComplexMatrix& a, SVD::type svd_type, SVD::driver svd_driver) { octave_idx_type info; @@ -144,24 +153,50 @@ octave_idx_type one = 1; octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one); - 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))); + if (svd_driver == SVD::GESVD) + { + 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 (lwork, 1); - lwork = static_cast<octave_idx_type> (work(0).real ()); - work.resize (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_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); - 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))); + 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 (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 ();
--- a/liboctave/fCmplxSVD.h Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/fCmplxSVD.h Mon May 03 13:21:35 2010 +0200 @@ -38,15 +38,16 @@ FloatComplexSVD (void) { } - FloatComplexSVD (const FloatComplexMatrix& a, SVD::type svd_type = SVD::std) + FloatComplexSVD (const FloatComplexMatrix& a, + SVD::type svd_type = SVD::std, SVD::driver svd_driver = SVD::GESVD) { - init (a, svd_type); + init (a, svd_type, svd_driver); } FloatComplexSVD (const FloatComplexMatrix& a, octave_idx_type& info, - SVD::type svd_type = SVD::std) + SVD::type svd_type = SVD::std, SVD::driver svd_driver = SVD::GESVD) { - info = init (a, svd_type); + info = init (a, svd_type, svd_driver); } FloatComplexSVD (const FloatComplexSVD& a) @@ -83,7 +84,9 @@ FloatComplexMatrix left_sm; FloatComplexMatrix right_sm; - octave_idx_type init (const FloatComplexMatrix& a, SVD::type svd_type = SVD::std); + octave_idx_type init (const FloatComplexMatrix& a, + SVD::type svd_type = SVD::std, + SVD::driver svd_driver = SVD::GESVD); }; #endif
--- a/liboctave/floatSVD.cc Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/floatSVD.cc Mon May 03 13:21:35 2010 +0200 @@ -29,6 +29,7 @@ #include "floatSVD.h" #include "f77-fcn.h" +#include "oct-locbuf.h" extern "C" { @@ -41,6 +42,14 @@ 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 @@ -70,7 +79,7 @@ } octave_idx_type -FloatSVD::init (const FloatMatrix& a, SVD::type svd_type) +FloatSVD::init (const FloatMatrix& a, SVD::type svd_type, SVD::driver svd_driver) { octave_idx_type info; @@ -140,22 +149,48 @@ octave_idx_type one = 1; octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one); - 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))); + 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 (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))); - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (lwork, 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))); - 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 (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 ();
--- a/liboctave/floatSVD.h Sun May 02 22:05:41 2010 -0700 +++ b/liboctave/floatSVD.h Mon May 03 13:21:35 2010 +0200 @@ -38,11 +38,14 @@ FloatSVD (void) : sigma (), left_sm (), right_sm () { } - FloatSVD (const FloatMatrix& a, SVD::type svd_type = SVD::std) { init (a, svd_type); } + FloatSVD (const FloatMatrix& a, + SVD::type svd_type = SVD::std, SVD::driver svd_driver = SVD::GESVD) + { init (a, svd_type, svd_driver); } - FloatSVD (const FloatMatrix& a, octave_idx_type& info, SVD::type svd_type = SVD::std) + FloatSVD (const FloatMatrix& a, octave_idx_type& info, + SVD::type svd_type = SVD::std, SVD::driver svd_driver = SVD::GESVD) { - info = init (a, svd_type); + info = init (a, svd_type, svd_driver); } FloatSVD (const FloatSVD& a) @@ -80,7 +83,9 @@ FloatMatrix left_sm; FloatMatrix right_sm; - octave_idx_type init (const FloatMatrix& a, SVD::type svd_type = SVD::std); + octave_idx_type init (const FloatMatrix& a, + SVD::type svd_type = SVD::std, + SVD::driver svd_driver = SVD::GESVD); }; #endif
--- a/src/ChangeLog Sun May 02 22:05:41 2010 -0700 +++ b/src/ChangeLog Mon May 03 13:21:35 2010 +0200 @@ -1,3 +1,9 @@ +2010-05-03 Jaroslav Hajek <highegg@gmail.com> + + * DLD-FUNCTIONS/svd.cc (driver): New static var. + (Fsvd): Use it to select a driver. + (Fsvd_driver): New DEFUN. + 2010-05-02 Rik <octave@nomad.inbox5.com> * ov-cell.cc: Wrap documentation line to prevent overfull hbox
--- a/src/DLD-FUNCTIONS/svd.cc Sun May 02 22:05:41 2010 -0700 +++ b/src/DLD-FUNCTIONS/svd.cc Mon May 03 13:21:35 2010 +0200 @@ -37,6 +37,8 @@ #include "pr-output.h" #include "utils.h" +static SVD::driver driver = SVD::GESVD; + DEFUN_DLD (svd, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{s} =} svd (@var{a})\n\ @@ -203,7 +205,7 @@ return retval; } - FloatSVD result (tmp, type); + FloatSVD result (tmp, type, driver); FloatDiagMatrix sigma = result.singular_values (); @@ -231,7 +233,7 @@ return retval; } - FloatComplexSVD result (ctmp, type); + FloatComplexSVD result (ctmp, type, driver); FloatDiagMatrix sigma = result.singular_values (); @@ -262,7 +264,7 @@ return retval; } - SVD result (tmp, type); + SVD result (tmp, type, driver); DiagMatrix sigma = result.singular_values (); @@ -290,7 +292,7 @@ return retval; } - ComplexSVD result (ctmp, type); + ComplexSVD result (ctmp, type, driver); DiagMatrix sigma = result.singular_values (); @@ -396,3 +398,43 @@ %!error <Invalid call to svd.*> [u, v] = svd ([1, 2; 3, 4]); */ + +DEFUN_DLD (svd_driver, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{old} =} svd_driver (@var{new})\n\ +Sets or queries the underlying LAPACK driver used by svd.\n\ +Currently recognized are \"gesvd\" and \"gesdd\". Default is \"gesvd\".\n\ +@seealso{svd}\n\ +@end deftypefn") +{ + octave_value retval; + int nargin = args.length (); + + if (nargin == 0 || (nargin == 1 && args(0).is_string ())) + { + if (driver == SVD::GESVD) + retval = "gesvd"; + else if (driver == SVD::GESDD) + retval = "gesdd"; + else + panic_impossible (); + + + if (nargin == 1) + { + std::string opt = args(0).xtolower ().string_value (); + if (opt == "gesvd") + driver = SVD::GESVD; + else if (opt == "gesdd") + driver = SVD::GESDD; + else + error ("svd_driver: invalid driver spec: %s", opt.c_str ()); + } + } + else + print_usage (); + + return retval; +} + +