Mercurial > hg > octave-terminal
changeset 537:4ecbfd3c3710
[project @ 1994-07-21 22:30:34 by jwe]
author | jwe |
---|---|
date | Thu, 21 Jul 1994 22:30:47 +0000 |
parents | 83b0e891100c |
children | 8e134d3b21c9 |
files | liboctave/CmplxSVD.cc liboctave/CmplxSVD.h liboctave/dbleSVD.cc liboctave/dbleSVD.h |
diffstat | 4 files changed, 67 insertions(+), 36 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxSVD.cc +++ b/liboctave/CmplxSVD.cc @@ -42,24 +42,35 @@ } int -ComplexSVD::init (const ComplexMatrix& a) +ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type) { int info; int m = a.rows (); int n = a.cols (); - char jobu = 'A'; - char jobv = 'A'; - Complex *tmp_data = dup (a.data (), a.length ()); int min_mn = m < n ? m : n; int max_mn = m > n ? m : n; - Complex *u = new Complex[m*m]; + char jobu = 'A'; + char jobv = 'A'; + + int ncol_u = m; + int nrow_vt = n; + int nrow_s = m; + int ncol_s = n; + + if (svd_type == SVD::economy) + { + jobu = jobv = 'S'; + ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; + } + + Complex *u = new Complex[m * ncol_u]; double *s_vec = new double[min_mn]; - Complex *vt = new Complex[n*n]; + Complex *vt = new Complex[nrow_vt * n]; int lwork = 2*min_mn + max_mn; Complex *work = new Complex[lwork]; @@ -68,12 +79,12 @@ double *rwork = new double[lrwork]; F77_FCN (zgesvd) (&jobu, &jobv, &m, &n, tmp_data, &m, s_vec, u, &m, - vt, &n, work, &lwork, rwork, &info, 1L, 1L); + vt, &nrow_vt, work, &lwork, rwork, &info, 1L, 1L); - left_sm = ComplexMatrix (u, m, m); - sigma = DiagMatrix (s_vec, m, n); - ComplexMatrix vt_m (vt, n, n); - right_sm = ComplexMatrix (vt_m.hermitian ()); + left_sm = ComplexMatrix (u, m, ncol_u); + sigma = DiagMatrix (s_vec, nrow_s, ncol_s); + ComplexMatrix vt_m (vt, nrow_vt, n); + right_sm = vt_m.hermitian (); delete [] tmp_data; delete [] work;
--- a/liboctave/CmplxSVD.h +++ b/liboctave/CmplxSVD.h @@ -32,6 +32,7 @@ #include "dDiagMatrix.h" #include "CMatrix.h" +#include "dbleSVD.h" extern "C++" { @@ -43,8 +44,9 @@ ComplexSVD (void) {} - ComplexSVD (const ComplexMatrix& a); - ComplexSVD (const ComplexMatrix& a, int& info); + ComplexSVD (const ComplexMatrix& a, SVD::type svd_type = SVD::std); + ComplexSVD (const ComplexMatrix& a, int& info, + SVD::type svd_type = SVD::std); ComplexSVD (const ComplexSVD& a); @@ -58,21 +60,22 @@ private: - int init (const ComplexMatrix& a); + int init (const ComplexMatrix& a, SVD::type svd_type = SVD::std); DiagMatrix sigma; ComplexMatrix left_sm; ComplexMatrix right_sm; }; -inline ComplexSVD::ComplexSVD (const ComplexMatrix& a) +inline ComplexSVD::ComplexSVD (const ComplexMatrix& a, SVD::type svd_type) { - init (a); + init (a, svd_type); } -inline ComplexSVD::ComplexSVD (const ComplexMatrix& a, int& info) +inline ComplexSVD::ComplexSVD (const ComplexMatrix& a, int& info, + SVD::type svd_type) { - info = init (a); + info = init (a, svd_type); } inline ComplexSVD::ComplexSVD (const ComplexSVD& a)
--- a/liboctave/dbleSVD.cc +++ b/liboctave/dbleSVD.cc @@ -42,24 +42,35 @@ } int -SVD::init (const Matrix& a) +SVD::init (const Matrix& a, SVD::type svd_type) { int info; int m = a.rows (); int n = a.cols (); - char jobu = 'A'; - char jobv = 'A'; - double *tmp_data = dup (a.data (), a.length ()); int min_mn = m < n ? m : n; int max_mn = m > n ? m : n; - double *u = new double[m*m]; + char jobu = 'A'; + char jobv = 'A'; + + int ncol_u = m; + int nrow_vt = n; + int nrow_s = m; + int ncol_s = n; + + if (svd_type == SVD::economy) + { + jobu = jobv = 'S'; + ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; + } + + double *u = new double[m * ncol_u]; double *s_vec = new double[min_mn]; - double *vt = new double[n*n]; + double *vt = new double[nrow_vt * n]; int tmp1 = 3*min_mn + max_mn; int tmp2 = 5*min_mn - 4; @@ -67,12 +78,12 @@ double *work = new double[lwork]; F77_FCN (dgesvd) (&jobu, &jobv, &m, &n, tmp_data, &m, s_vec, u, &m, - vt, &n, work, &lwork, &info, 1L, 1L); + vt, &nrow_vt, work, &lwork, &info, 1L, 1L); - left_sm = Matrix (u, m, m); - sigma = DiagMatrix (s_vec, m, n); - Matrix vt_m (vt, n, n); - right_sm = Matrix (vt_m.transpose ()); + left_sm = Matrix (u, m, ncol_u); + sigma = DiagMatrix (s_vec, nrow_s, ncol_s); + Matrix vt_m (vt, nrow_vt, n); + right_sm = vt_m.transpose (); delete [] tmp_data; delete [] work;
--- a/liboctave/dbleSVD.h +++ b/liboctave/dbleSVD.h @@ -41,10 +41,16 @@ public: + enum type + { + std, + economy, + }; + SVD (void) {} - SVD (const Matrix& a); - SVD (const Matrix& a, int& info); + SVD (const Matrix& a, SVD::type svd_type = SVD::std); + SVD (const Matrix& a, int& info, SVD::type svd_type = SVD::std); SVD (const SVD& a); @@ -58,21 +64,21 @@ private: - int init (const Matrix& a); + int init (const Matrix& a, SVD::type svd_type = SVD::std); DiagMatrix sigma; Matrix left_sm; Matrix right_sm; }; -inline SVD::SVD (const Matrix& a) +inline SVD::SVD (const Matrix& a, SVD::type svd_type) { - init (a); + init (a, svd_type); } -inline SVD::SVD (const Matrix& a, int& info) +inline SVD::SVD (const Matrix& a, int& info, SVD::type svd_type) { - info = init (a); + info = init (a, svd_type); } inline SVD::SVD (const SVD& a)