Mercurial > hg > octave-terminal
changeset 1543:d6e96e0bc681
[project @ 1995-10-06 05:49:21 by jwe]
author | jwe |
---|---|
date | Fri, 06 Oct 1995 05:54:07 +0000 |
parents | 5bd8c07faacf |
children | f1931fc63ce9 |
files | liboctave/CmplxSVD.cc liboctave/CmplxSVD.h liboctave/dbleSVD.cc liboctave/dbleSVD.h src/svd.cc |
diffstat | 5 files changed, 104 insertions(+), 18 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxSVD.cc +++ b/liboctave/CmplxSVD.cc @@ -31,6 +31,7 @@ #include "CmplxSVD.h" #include "f77-uscore.h" +#include "lo-error.h" #include "mx-inlines.cc" extern "C" @@ -43,6 +44,32 @@ long); } +ComplexMatrix +ComplexSVD::left_singular_matrix (void) const +{ + if (type == SVD::sigma_only) + { + (*current_liboctave_error_handler) + ("ComplexSVD: U not computed because type == SVD::sigma_only"); + return ComplexMatrix (); + } + else + return left_sm; +} + +ComplexMatrix +ComplexSVD::right_singular_matrix (void) const +{ + if (type == SVD::sigma_only) + { + (*current_liboctave_error_handler) + ("ComplexSVD: V not computed because type == SVD::sigma_only"); + return ComplexMatrix (); + } + else + return right_sm; +} + int ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type) { @@ -64,15 +91,25 @@ int nrow_s = m; int ncol_s = n; - if (svd_type == SVD::economy) + switch (svd_type) { - jobu = jobv = "S"; + case SVD::economy: + jobu = jobv ="S"; ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; + break; + + case SVD::sigma_only: + jobu = jobv ="N"; + ncol_u = nrow_vt = 0; + break; + + default: + break; } - Complex *u = new Complex[m * ncol_u]; + Complex *u = (ncol_u > 0 ? new Complex[m * ncol_u] : 0); double *s_vec = new double[min_mn]; - Complex *vt = new Complex[nrow_vt * n]; + Complex *vt = (nrow_vt > 0 ? new Complex[nrow_vt * n] : 0); int lwork = 2*min_mn + max_mn; Complex *work = new Complex[lwork]; @@ -84,10 +121,16 @@ m, vt, nrow_vt, work, lwork, rwork, info, 1L, 1L); - left_sm = ComplexMatrix (u, m, ncol_u); + if (ncol_u > 0) + 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 (); + + if (nrow_vt > 0) + { + 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 @@ -72,9 +72,9 @@ DiagMatrix singular_values (void) const { return sigma; } - ComplexMatrix left_singular_matrix (void) const { return left_sm; } + ComplexMatrix left_singular_matrix (void); - ComplexMatrix right_singular_matrix (void) const { return right_sm; } + ComplexMatrix right_singular_matrix (void); friend ostream& operator << (ostream& os, const ComplexSVD& a);
--- a/liboctave/dbleSVD.cc +++ b/liboctave/dbleSVD.cc @@ -42,6 +42,32 @@ long, long); } +Matrix +left_singular_matrix (void) const +{ + if (type == SVD::sigma_only) + { + (*current_liboctave_error_handler) + ("ComplexSVD: U not computed because type == SVD::sigma_only"); + return Matrix (); + } + else + return left_sm; +} + +Matrix +right_singular_matrix (void) const +{ + if (type == SVD::sigma_only) + { + (*current_liboctave_error_handler) + ("ComplexSVD: V not computed because type == SVD::sigma_only"); + return Matrix (); + } + else + return right_sm; +} + int SVD::init (const Matrix& a, SVD::type svd_type) { @@ -63,15 +89,25 @@ int nrow_s = m; int ncol_s = n; - if (svd_type == SVD::economy) + switch (svd_type) { + case SVD::economy: jobu = jobv ="S"; ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; + break; + + case SVD::sigma_only: + jobu = jobv ="N"; + ncol_u = nrow_vt = 0; + break; + + default: + break; } - double *u = new double[m * ncol_u]; + double *u = (ncol_u > 0 ? new double[m * ncol_u] : 0); double *s_vec = new double[min_mn]; - double *vt = new double[nrow_vt * n]; + double *vt = (ncol_vt > 0 ? new double[nrow_vt * n] : 0); int tmp1 = 3*min_mn + max_mn; int tmp2 = 5*min_mn - 4; @@ -82,10 +118,16 @@ m, vt, nrow_vt, work, lwork, info, 1L, 1L); - left_sm = Matrix (u, m, ncol_u); + if (ncol_u > 0) + 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 (); + + if (nrow_vt > 0) + { + 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 @@ -74,9 +74,9 @@ DiagMatrix singular_values (void) const { return sigma; } - Matrix left_singular_matrix (void) const { return left_sm; } + Matrix left_singular_matrix (void); - Matrix right_singular_matrix (void) const { return right_sm; } + Matrix right_singular_matrix (void); friend ostream& operator << (ostream& os, const SVD& a);
--- a/src/svd.cc +++ b/src/svd.cc @@ -66,7 +66,8 @@ else if (arg_is_empty > 0) return Octave_object (3, Matrix ()); - SVD::type type = (nargin == 2) ? SVD::economy : SVD::std; + SVD::type type = ((nargout == 1) ? SVD::sigma_only + : (nargin == 2) ? SVD::economy : SVD::std); if (arg.is_real_type ()) {