Mercurial > hg > octave-lyh

changeset 537:4ecbfd3c3710

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[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)