Mercurial > hg > octave-lyh

changeset 1543:d6e96e0bc681

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