Mercurial > hg > octave-max

changeset 3336:08ad797989f8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1999-11-03 21:41:34 by jwe]
author jwe
date Wed, 03 Nov 1999 21:41:35 +0000
parents f39b97e13cf2
children 9c782a34b545
files liboctave/ChangeLog liboctave/CmplxSVD.cc liboctave/dbleSVD.cc
diffstat 3 files changed, 50 insertions(+), 21 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,5 +1,8 @@
 1999-11-03  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* dbleSVD.cc (SVD::init): Let DGESVD determine work space requirement.
+	* CmplxSVD.cc (ComplexSVD::init): Likewise, for complex version.
+
 	* dbleSCHUR.cc (SCHUR::init): IWORK is always referenced by dgeesx.
 	Don't forget to pass length of third string argument to dgeesx.
 
--- a/liboctave/CmplxSVD.cc
+++ b/liboctave/CmplxSVD.cc
@@ -131,26 +131,41 @@
 
   Complex *vt = right_sm.fortran_vec ();
 
-  int lwork = 2*min_mn + max_mn;
-
-  Array<Complex> work (lwork);
-  Complex *pwork = work.fortran_vec ();
-
   int lrwork = 5*max_mn;
 
   Array<double> rwork (lrwork);
-  double *prwork = rwork.fortran_vec ();
+
+  // Ask ZGESVD what the dimension of WORK should be.
+
+  int lwork = -1;
 
-  F77_XFCN (zgesvd, ZGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec, u,
-			     m, vt, nrow_vt, pwork, lwork, prwork, info,
-			     1L, 1L));
+  Array<Complex> work (1);
+
+  F77_XFCN (zgesvd, ZGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec,
+			     u, m, vt, nrow_vt, work.fortran_vec (),
+			     lwork, rwork.fortran_vec (), info, 1L,
+			     1L));
 
   if (f77_exception_encountered)
     (*current_liboctave_error_handler) ("unrecoverable error in zgesvd");
   else
     {
-      if (! (jobv == 'N' || jobv == 'O'))
-	right_sm = right_sm.hermitian ();
+      lwork = static_cast<int> (work(0).real ());
+      work.resize (lwork);
+
+      F77_XFCN (zgesvd, ZGESVD, (&jobu, &jobv, m, n, tmp_data, m,
+				 s_vec, u, m, vt, nrow_vt,
+				 work.fortran_vec (), lwork,
+				 rwork.fortran_vec (),
+				 info, 1L, 1L));
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in zgesvd");
+      else
+	{
+	  if (! (jobv == 'N' || jobv == 'O'))
+	    right_sm = right_sm.hermitian ();
+	}
     }
 
   return info;
--- a/liboctave/dbleSVD.cc
+++ b/liboctave/dbleSVD.cc
@@ -80,7 +80,6 @@
   double *tmp_data = atmp.fortran_vec ();
 
   int min_mn = m < n ? m : n;
-  int max_mn = m > n ? m : n;
 
   char jobu = 'A';
   char jobv = 'A';
@@ -131,23 +130,35 @@
 
   double *vt = right_sm.fortran_vec ();
 
-  int tmp1 = 3*min_mn + max_mn;
-  int tmp2 = 5*min_mn - 4;
-  int lwork = tmp1 > tmp2 ? tmp1 : tmp2;
+  // Ask DGESVD what the dimension of WORK should be.
 
-  Array<double> work (lwork);
-  double *pwork = work.fortran_vec ();
+  int lwork = -1;
+
+  Array<double> work (1);
 
   F77_XFCN (dgesvd, DGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec,
-			     u, m, vt, nrow_vt, pwork, lwork, info,
-			     1L, 1L));
+			     u, m, vt, nrow_vt, work.fortran_vec (),
+			     lwork, info, 1L, 1L));
 
   if (f77_exception_encountered)
     (*current_liboctave_error_handler) ("unrecoverable error in dgesvd");
   else
     {
-      if (! (jobv == 'N' || jobv == 'O'))
-	right_sm = right_sm.transpose ();
+      lwork = static_cast<int> (work(0));
+      work.resize (lwork);
+
+      F77_XFCN (dgesvd, DGESVD, (&jobu, &jobv, m, n, tmp_data, m,
+				 s_vec, u, m, vt, nrow_vt,
+				 work.fortran_vec (), lwork, info, 1L,
+				 1L));
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in dgesvd");
+      else
+	{
+	  if (! (jobv == 'N' || jobv == 'O'))
+	    right_sm = right_sm.transpose ();
+	}
     }
 
   return info;