Mercurial > hg > octave-nkf

changeset 1930:d20ab06301e8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-02-11 22:30:18 by jwe]
author jwe
date Sun, 11 Feb 1996 22:30:39 +0000
parents 908f5b6676d7
children c91f81d5f72c
files liboctave/CmplxSCHUR.cc liboctave/CmplxSCHUR.h liboctave/CmplxSVD.cc liboctave/CmplxSVD.h liboctave/dbleSCHUR.cc liboctave/dbleSCHUR.h liboctave/dbleSVD.cc liboctave/dbleSVD.h
diffstat 8 files changed, 81 insertions(+), 58 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CmplxSCHUR.cc
+++ b/liboctave/CmplxSCHUR.cc
@@ -70,20 +70,21 @@
       return -1;
     }
 
-  char *jobvs = "V";
-  char *sense = "N";
-  char *sort = "N";
+  char jobvs = 'V';
+  char sense = 'N';
+  char sort = 'N';
 
   char ord_char = ord.empty () ? 'U' : ord[0];
 
   if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
-    sort = "S";
+    sort = 'S';
 
-  select_function selector= 0;
   if (ord_char == 'A' || ord_char == 'a')
     selector = select_ana;
   else if (ord_char == 'D' || ord_char == 'd')
     selector = select_dig;
+  else
+    selector = 0;
 
   int n = a_nc;
   int lwork = 8 * n;
@@ -116,7 +117,7 @@
 
   int *pbwork = bwork.fortran_vec ();
 
-  F77_XFCN (zgeesx, ZGEESX, (jobvs, sort, selector, sense, n, s, n,
+  F77_XFCN (zgeesx, ZGEESX, (&jobvs, &sort, selector, &sense, n, s, n,
 			     sdim, pw, q, n, rconde, rcondv, pwork,
 			     lwork, prwork, pbwork, info, 1L, 1L));
 
--- a/liboctave/CmplxSCHUR.h
+++ b/liboctave/CmplxSCHUR.h
@@ -82,6 +82,8 @@
   ComplexMatrix schur_mat;
   ComplexMatrix unitary_mat;
 
+  select_function selector;
+
   int init (const ComplexMatrix& a, const string& ord);
 };
 
--- a/liboctave/CmplxSVD.cc
+++ b/liboctave/CmplxSVD.cc
@@ -83,8 +83,8 @@
   int min_mn = m < n ? m : n;
   int max_mn = m > n ? m : n;
 
-  char *jobu = "A";
-  char *jobv = "A";
+  char jobu = 'A';
+  char jobv = 'A';
 
   int ncol_u = m;
   int nrow_vt = n;
@@ -94,12 +94,12 @@
   switch (svd_type)
     {
     case SVD::economy:
-      jobu = jobv ="S";
+      jobu = jobv = 'S';
       ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
       break;
 
     case SVD::sigma_only:
-      jobu = jobv ="N";
+      jobu = jobv = 'N';
       ncol_u = nrow_vt = 1;
       break;
 
@@ -109,34 +109,41 @@
 
   type_computed = svd_type;
 
-  Complex *u = (*jobu == 'N' ? 0 : new Complex[m * ncol_u]);
-  double *s_vec  = new double[min_mn];
-  Complex *vt = (*jobv == 'N' ? 0 : new Complex[nrow_vt * n]);
+  if (jobu != 'N')
+    left_sm.resize (m, ncol_u);
+
+  Complex *u = left_sm.fortran_vec ();
+
+  sigma.resize (nrow_s, ncol_s);
+  double *s_vec = sigma.fortran_vec ();
+
+  if (jobv != 'N')
+    right_sm.resize (nrow_vt, n);
+
+  Complex *vt = right_sm.fortran_vec ();
 
   int lwork = 2*min_mn + max_mn;
-  Complex *work = new Complex[lwork];
+
+  Array<Complex> work (lwork);
+  Complex *pwork = work.fortran_vec ();
 
   int lrwork = 5*max_mn;
-  double *rwork = new double[lrwork];
+
+  Array<double> rwork (lrwork);
+  double *prwork = rwork.fortran_vec ();
 
-  F77_FCN (zgesvd, ZGESVD) (jobu, jobv, m, n, tmp_data, m, s_vec, u,
-			    m, vt, nrow_vt, work, lwork, rwork, info,
-			    1L, 1L);
-
-  if (*jobu != 'N')
-    left_sm = ComplexMatrix (u, m, ncol_u);
+  F77_XFCN (zgesvd, ZGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec, u,
+			     m, vt, nrow_vt, pwork, lwork, prwork, info,
+			     1L, 1L));
 
-  sigma = DiagMatrix (s_vec, nrow_s, ncol_s);
-
-  if (*jobv != 'N')
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in zgesvd");
+  else
     {
-      ComplexMatrix vt_m (vt, nrow_vt, n);
-      right_sm = vt_m.hermitian ();
+      if (jobv != 'N')
+	right_sm = right_sm.hermitian ();
     }
 
-  delete [] tmp_data;
-  delete [] work;
-
   return info;
 }
 
--- a/liboctave/CmplxSVD.h
+++ b/liboctave/CmplxSVD.h
@@ -66,6 +66,8 @@
       return *this;
     }
 
+  ~ComplexSVD (void) { }
+
   DiagMatrix singular_values (void) const { return sigma; }
 
   ComplexMatrix left_singular_matrix (void) const;
--- a/liboctave/dbleSCHUR.cc
+++ b/liboctave/dbleSCHUR.cc
@@ -71,20 +71,21 @@
       return -1;
     }
 
-  char *jobvs = "V";
-  char *sense = "N";
-  char *sort = "N";
+  char jobvs = 'V';
+  char sense = 'N';
+  char sort = 'N';
 
   char ord_char = ord.empty () ? 'U' : ord[0];
 
   if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
-    sort = "S";
+    sort = 'S';
 
-  select_function selector = 0;
   if (ord_char == 'A' || ord_char == 'a')
     selector = select_ana;
   else if (ord_char == 'D' || ord_char == 'd')
     selector = select_dig;
+  else
+    selector = 0;
 
   int n = a_nc;
   int lwork = 8 * n;
@@ -124,7 +125,7 @@
   int *piwork = iwork.fortran_vec ();
 
 
-  F77_XFCN (dgeesx, DGEESX, (jobvs, sort, selector, sense, n, s,
+  F77_XFCN (dgeesx, DGEESX, (&jobvs, &sort, selector, &sense, n, s,
 			     n, sdim, pwr, pwi, q, n, rconde, rcondv,
 			     pwork, lwork, piwork, liwork, pbwork,
 			     info, 1L, 1L));
--- a/liboctave/dbleSCHUR.h
+++ b/liboctave/dbleSCHUR.h
@@ -82,6 +82,8 @@
   Matrix schur_mat;
   Matrix unitary_mat;
 
+  select_function selector;
+
   int init (const Matrix& a, const string& ord);
 };
 
--- a/liboctave/dbleSVD.cc
+++ b/liboctave/dbleSVD.cc
@@ -78,13 +78,14 @@
   int m = a.rows ();
   int n = a.cols ();
 
-  double *tmp_data = dup (a.data (), a.length ());
+  Matrix atmp = a;
+  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";
+  char jobu = 'A';
+  char jobv = 'A';
 
   int ncol_u = m;
   int nrow_vt = n;
@@ -94,12 +95,12 @@
   switch (svd_type)
     {
     case SVD::economy:
-      jobu = jobv ="S";
+      jobu = jobv = 'S';
       ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
       break;
 
     case SVD::sigma_only:
-      jobu = jobv ="N";
+      jobu = jobv = 'N';
       ncol_u = nrow_vt = 1;
       break;
 
@@ -109,33 +110,38 @@
 
   type_computed = svd_type;
 
-  double *u = (*jobu == 'N' ? 0 : new double[m * ncol_u]);
-  double *s_vec  = new double[min_mn];
-  double *vt = (*jobv == 'N' ? 0 : new double[nrow_vt * n]);
+  if (jobu != 'N')
+    left_sm.resize (n, ncol_u);
+
+  double *u = left_sm.fortran_vec ();
+
+  sigma.resize (nrow_s, ncol_s);
+  double *s_vec  = sigma.fortran_vec ();
+
+  if (jobv != 'N')
+    right_sm.resize (nrow_vt, n);
+
+  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;
-  double *work = new double[lwork];
+
+  Array<double> work (lwork);
+  double *pwork = work.fortran_vec ();
 
-  F77_FCN (dgesvd, DGESVD) (jobu, jobv, m, n, tmp_data, m, s_vec, u,
-			    m, vt, nrow_vt, work, lwork, info, 1L,
-			    1L);
-
-  if (*jobu != 'N')
-    left_sm = Matrix (u, m, ncol_u);
+  F77_XFCN (dgesvd, DGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec,
+			     u, m, vt, nrow_vt, pwork, lwork, info,
+			     1L, 1L));
 
-  sigma = DiagMatrix (s_vec, nrow_s, ncol_s);
-
-  if (*jobv != 'N')
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in dgesvd");
+  else
     {
-      Matrix vt_m (vt, nrow_vt, n);
-      right_sm = vt_m.transpose ();
+      if (jobv != 'N')
+	right_sm = right_sm.transpose ();
     }
 
-  delete [] tmp_data;
-  delete [] work;
-
   return info;
 }
 
--- a/liboctave/dbleSVD.h
+++ b/liboctave/dbleSVD.h
@@ -69,6 +69,8 @@
       return *this;
     }
 
+  ~SVD (void) { }
+
   DiagMatrix singular_values (void) const { return sigma; }
 
   Matrix left_singular_matrix (void) const;