Mercurial > hg > octave-max

changeset 1932:682f31b20894

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-02-12 02:26:07 by jwe]
author jwe
date Mon, 12 Feb 1996 02:26:07 +0000
parents c91f81d5f72c
children 333923894848
files liboctave/CmplxHESS.cc liboctave/CmplxHESS.h liboctave/dbleHESS.cc liboctave/dbleHESS.h
diffstat 4 files changed, 120 insertions(+), 73 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CmplxHESS.cc
+++ b/liboctave/CmplxHESS.cc
@@ -60,59 +60,82 @@
 {
   int a_nr = a.rows ();
   int a_nc = a.cols ();
-   if (a_nr != a_nc)
-     {
-       (*current_liboctave_error_handler)
-	 ("ComplexHESS requires square matrix");
-       return -1;
-     }
+
+  if (a_nr != a_nc)
+    {
+      (*current_liboctave_error_handler)
+	("ComplexHESS requires square matrix");
+      return -1;
+    }
 
-   char *job = "N";
-   char *side = "R";
+  char job = 'N';
+  char side = 'R';
+
+  int n = a_nc;
+  int lwork = 32 * n;
+  int info;
+  int ilo;
+  int ihi;
 
-   int n = a_nc;
-   int lwork = 32 * n;
-   int info;
-   int ilo;
-   int ihi;
+  hess_mat = a;
+  Complex *h = hess_mat.fortran_vec ();
 
-   Complex *h = dup (a.data (), a.length ());
+  Array<double> scale (n);
+  double *pscale = scale.fortran_vec ();
+
+  F77_XFCN (zgebal, ZGEBAL, (&job, n, h, n, ilo, ihi, pscale, info,
+			     1L, 1L));
 
-   double *scale = new double [n];
-   Complex *tau = new Complex [n-1];
-   Complex *work = new Complex [lwork];
-   Complex *z = new Complex [n*n];
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
+  else
+    {
+      Array<Complex> tau (n-1);
+      Complex *ptau = tau.fortran_vec ();
 
-   F77_FCN (zgebal, ZGEBAL) (job, n, h, n, ilo, ihi, scale, info, 1L,
-			     1L);
+      Array<Complex> work (lwork);
+      Complex *pwork = work.fortran_vec ();
 
-   F77_FCN (zgehrd, ZGEHRD) (n, ilo, ihi, h, n, tau, work, lwork,
-			     info, 1L, 1L);
-
-   copy (z, h, n*n);
+      F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork,
+				 info, 1L, 1L));
 
-   F77_FCN (zunghr, ZUNGHR) (n, ilo, ihi, z, n, tau, work, lwork,
-			     info, 1L, 1L);
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in zgehrd");
+      else
+	{
+	  unitary_hess_mat = hess_mat;
+	  Complex *z = unitary_hess_mat.fortran_vec ();
 
-   F77_FCN (zgebak, ZGEBAK) (job, side, n, ilo, ihi, scale, n, z, n,
-			     info, 1L, 1L);
-
-   hess_mat = ComplexMatrix (h, n, n);
-   unitary_hess_mat = ComplexMatrix (z, n, n);
+	  F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, z, n, ptau, pwork,
+				     lwork, info, 1L, 1L));
 
-  // If someone thinks of a more graceful way of doing this (or faster
-  // for that matter :-)), please let me know!
+	  if (f77_exception_encountered)
+	    (*current_liboctave_error_handler)
+	      ("unrecoverable error in zunghr");
+	  else
+	    {
+	      F77_XFCN (zgebak, ZGEBAK, (&job, &side, n, ilo, ihi,
+					 pscale, n, z, n, info, 1L, 1L));
 
-   if (n > 2)
-     for (int j = 0; j < a_nc; j++)
-       for (int i = j+2; i < a_nr; i++)
-         hess_mat.elem (i, j) = 0;
+	      if (f77_exception_encountered)
+		(*current_liboctave_error_handler)
+		  ("unrecoverable error in zgebak");
+	      else
+		{
+		  // If someone thinks of a more graceful way of
+		  // doing this (or faster for that matter :-)),
+		  // please let me know!
 
-   delete [] work;
-   delete [] tau;
-   delete [] scale;
+		  if (n > 2)
+		    for (int j = 0; j < a_nc; j++)
+		      for (int i = j+2; i < a_nr; i++)
+			hess_mat.elem (i, j) = 0;
+		}
+	    }
+	}
+    }
 
-   return info;
+  return info;
 }
 
 /*
--- a/liboctave/CmplxHESS.h
+++ b/liboctave/CmplxHESS.h
@@ -56,6 +56,8 @@
       return *this;
     }
 
+  ~ComplexHESS (void) { }
+
   ComplexMatrix hess_matrix (void) const { return hess_mat; }
 
   ComplexMatrix unitary_hess_matrix (void) const
--- a/liboctave/dbleHESS.cc
+++ b/liboctave/dbleHESS.cc
@@ -59,14 +59,15 @@
 {
   int a_nr = a.rows ();
   int a_nc = a.cols ();
+
   if (a_nr != a_nc)
     {
       (*current_liboctave_error_handler) ("HESS requires square matrix");
       return -1;
     }
 
-  char *jobbal = "N";
-  char *side = "R";
+  char job = 'N';
+  char side = 'R';
 
   int n = a_nc;
   int lwork = 32 * n;
@@ -74,44 +75,63 @@
   int ilo;
   int ihi;
 
-  double *h = dup (a.data (), a.length ());
+  hess_mat = a;
+  double *h = hess_mat.fortran_vec ();
 
-  double *tau = new double [n+1];
-  double *scale = new double [n];
-  double *z = new double [n*n];
-  double *work = new double [lwork];
+  Array<double> scale (n);
+  double *pscale = scale.fortran_vec ();
+
+  F77_XFCN (dgebal, DGEBAL, (&job, n, h, n, ilo, ihi, pscale, info,
+			     1L, 1L));
 
-  F77_FCN (dgebal, DGEBAL) (jobbal, n, h, n, ilo, ihi, scale, info,
-			    1L, 1L);
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in dgebal");
+  else
+    {
+      Array<double> tau (n-1);
+      double *ptau = tau.fortran_vec ();
+
+      Array<double> work (lwork);
+      double *pwork = work.fortran_vec ();
 
-  F77_FCN (dgehrd, DGEHRD) (n, ilo, ihi, h, n, tau, work, lwork, info,
-			    1L, 1L);
+      F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
+				 lwork, info, 1L, 1L));
 
-  copy (z, h, n*n);
-
-  F77_FCN (dorghr, DORGHR) (n, ilo, ihi, z, n, tau, work, lwork, info,
-			    1L, 1L);
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in dgehrd");
+      else
+	{
+	  unitary_hess_mat = hess_mat;
+	  double *z = unitary_hess_mat.fortran_vec ();
 
-  F77_FCN (dgebak, DGEBAK) (jobbal, side, n, ilo, ihi, scale, n, z, n,
-			    info, 1L, 1L);
+	  F77_XFCN (dorghr, DORGHR, (n, ilo, ihi, z, n, ptau, pwork,
+				     lwork, info, 1L, 1L));
 
-  // We need to clear out all of the area below the sub-diagonal which
-  // was used to store the unitary matrix.
-
-  hess_mat = Matrix (h, n, n);
-  unitary_hess_mat = Matrix (z, n, n);
+	  if (f77_exception_encountered)
+	    (*current_liboctave_error_handler)
+	      ("unrecoverable error in dorghr");
+	  else
+	    {
+	      F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi,
+					 pscale, n, z, n, info, 1L, 1L));
 
-  // If someone thinks of a more graceful way of doing this (or faster
-  // for that matter :-)), please let me know! 
+	      if (f77_exception_encountered)
+		(*current_liboctave_error_handler)
+		  ("unrecoverable error in dgebak");
+	      else
+		{
+		  // If someone thinks of a more graceful way of doing
+		  // this (or faster for that matter :-)), please let
+		  // me know!
 
-  if (n > 2)
-    for (int j = 0; j < a_nc; j++)
-      for (int i = j+2; i < a_nr; i++)
-        hess_mat.elem (i, j) = 0;
-
-  delete [] tau;
-  delete [] work;
-  delete [] scale;
+		  if (n > 2)
+		    for (int j = 0; j < a_nc; j++)
+		      for (int i = j+2; i < a_nr; i++)
+			hess_mat.elem (i, j) = 0;
+		}
+	    }
+	}
+    }
 
   return info;
 }
--- a/liboctave/dbleHESS.h
+++ b/liboctave/dbleHESS.h
@@ -56,6 +56,8 @@
       return *this;
     }
 
+  ~HESS (void) { }
+
   Matrix hess_matrix (void) const { return hess_mat; }
 
   Matrix unitary_hess_matrix (void) const { return unitary_hess_mat; }