Mercurial > hg > octave-nkf

changeset 1933:333923894848

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-02-12 03:07:39 by jwe]
author jwe
date Mon, 12 Feb 1996 03:07:39 +0000
parents 682f31b20894
children 0e591d443ff0
files liboctave/CmplxAEPBAL.cc liboctave/CmplxAEPBAL.h liboctave/dbleAEPBAL.cc liboctave/dbleAEPBAL.h liboctave/dbleGEPBAL.cc liboctave/dbleGEPBAL.h
diffstat 6 files changed, 156 insertions(+), 112 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CmplxAEPBAL.cc
+++ b/liboctave/CmplxAEPBAL.cc
@@ -50,44 +50,47 @@
 int
 ComplexAEPBALANCE::init (const ComplexMatrix& a, const string& balance_job)
 {
-  int a_nc = a.cols ();
+  int n = a.cols ();
 
-  if (a.rows () != a_nc)
+  if (a.rows () != n)
     {
       (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix");
       return -1;
     }
 
-  int n = a_nc;
-
-  // Parameters for balance call.
-
   int info;
   int ilo;
   int ihi;
-  double *scale = new double [n];
 
-  // Copy matrix into local structure.
+  Array<double> scale (n);
+  double *pscale = scale.fortran_vec ();
 
   balanced_mat = a;
-
-  char bal_job = balance_job[0];
+  Complex *p_balanced_mat = balanced_mat.fortran_vec ();
 
-  F77_FCN (zgebal, ZGEBAL) (&bal_job, n,
-			    balanced_mat.fortran_vec (), n, ilo, ihi,
-			    scale, info, 1L, 1L);
+  char job = balance_job[0];
 
-  // Initialize balancing matrix to identity.
+  F77_XFCN (zgebal, ZGEBAL, (&job, n, p_balanced_mat, n, ilo, ihi,
+			     pscale, info, 1L, 1L));
 
-  balancing_mat = Matrix (n, n, 0.0);
-  for (int i = 0; i < n; i++)
-    balancing_mat (i, i) = 1.0;
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
+  else
+    {
+      balancing_mat = Matrix (n, n, 0.0);
+      for (int i = 0; i < n; i++)
+	balancing_mat.elem (i, i) = 1.0;
 
-  F77_FCN (zgebak, ZGEBAK) (&bal_job, "R", n, ilo, ihi, scale, n, 
-			    balancing_mat.fortran_vec (), n, info, 1L,
-			    1L);
+      Complex *p_balancing_mat = balancing_mat.fortran_vec ();
+
+      char side = 'R';
 
-  delete [] scale;
+      F77_XFCN (zgebak, ZGEBAK, (&job, &side, n, ilo, ihi, pscale, n,
+				 p_balancing_mat, n, info, 1L, 1L));
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in zgebak");
+    }
 
   return info;
 }
--- a/liboctave/CmplxAEPBAL.h
+++ b/liboctave/CmplxAEPBAL.h
@@ -59,6 +59,8 @@
       return *this;
     }
 
+  ~ComplexAEPBALANCE (void) { }
+
   ComplexMatrix balanced_matrix (void) const { return balanced_mat; }
 
   ComplexMatrix balancing_matrix (void) const { return balancing_mat; }
--- a/liboctave/dbleAEPBAL.cc
+++ b/liboctave/dbleAEPBAL.cc
@@ -49,44 +49,47 @@
 int
 AEPBALANCE::init (const Matrix& a, const string& balance_job)
 {
-  int a_nc = a.cols ();
+  int n = a.cols ();
 
-  if (a.rows () != a_nc)
+  if (a.rows () != n)
     {
       (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix");
       return -1;
     }
 
-  int n = a_nc;
-
-  // Parameters for balance call.
-
   int info;
   int ilo;
   int ihi;
-  double *scale = new double [n];
 
-  // Copy matrix into local structure.
+  Array<double> scale (n);
+  double *pscale = scale.fortran_vec ();
 
   balanced_mat = a;
-
-  char bal_job = balance_job[0];
+  double *p_balanced_mat = balanced_mat.fortran_vec ();
 
-  F77_FCN (dgebal, DGEBAL) (&bal_job, n,
-			    balanced_mat.fortran_vec (), n, ilo, ihi,
-			    scale, info, 1L, 1L);
+  char job = balance_job[0];
 
-  // Initialize balancing matrix to identity.
+  F77_XFCN (dgebal, DGEBAL, (&job, n, p_balanced_mat, n, ilo, ihi,
+			     pscale, info, 1L, 1L));
 
-  balancing_mat = Matrix (n, n, 0.0);
-  for (int i = 0; i < n; i++)
-    balancing_mat.elem (i ,i) = 1.0;
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in dgebal");
+  else
+    {
+      balancing_mat = Matrix (n, n, 0.0);
+      for (int i = 0; i < n; i++)
+	balancing_mat.elem (i ,i) = 1.0;
 
-  F77_FCN (dgebak, DGEBAK) (&bal_job, "R", n, ilo, ihi, scale, n,
-			    balancing_mat.fortran_vec (), n, info, 1L,
-			    1L);
+      double *p_balancing_mat = balancing_mat.fortran_vec ();
+
+      char side = 'R';
 
-  delete [] scale;
+      F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, pscale, n,
+				 p_balancing_mat, n, info, 1L, 1L));
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in dgebak");
+    }
 
   return info;
 }
--- a/liboctave/dbleAEPBAL.h
+++ b/liboctave/dbleAEPBAL.h
@@ -59,6 +59,8 @@
       return *this;
     }
 
+  ~AEPBALANCE (void) { }
+
   Matrix balanced_matrix (void) const { return balanced_mat; }
 
   Matrix balancing_matrix (void) const { return balancing_mat; }
--- a/liboctave/dbleGEPBAL.cc
+++ b/liboctave/dbleGEPBAL.cc
@@ -73,14 +73,15 @@
 
   int n = a_nc;
 
-  // Parameters for balance call.
-
   int info;
   int ilo;
   int ihi;
-  double *cscale = new double [n];
-  double *cperm = new double [n];
+
+  Array<double> cscale (n);
+  double *pcscale = cscale.fortran_vec ();
+
   Matrix wk (n, 6, 0.0);
+  double *pwk = wk.fortran_vec ();
 
   // Back out the permutations:
   //
@@ -102,23 +103,18 @@
   balanced_a_mat = a;
   balanced_b_mat = b;
 
-  // Initialize balancing matrices to identity.
-
-  left_balancing_mat = Matrix (n, n, 0.0);
-  for (int i = 0; i < n; i++)
-    left_balancing_mat (i, i) = 1.0;
-
-  right_balancing_mat = left_balancing_mat;
+  double *p_balanced_a_mat = balanced_a_mat.fortran_vec ();
+  double *p_balanced_b_mat = balanced_b_mat.fortran_vec ();
 
   // Check for permutation option.
 
-  char bal_job = balance_job[0];
+  char job = balance_job[0];
 
-  if (bal_job == 'P' || bal_job == 'B')
+  if (job == 'P' || job == 'B')
     {
-      F77_FCN (reduce, REDUCE) (n, n, balanced_a_mat.fortran_vec (),
-				n, balanced_b_mat.fortran_vec (), ilo,
-				ihi, cscale, wk.fortran_vec ());
+      F77_XFCN (reduce, REDUCE, (n, n, p_balanced_a_mat, n,
+				 p_balanced_b_mat, ilo, ihi,
+				 pcscale, pwk));
     }
   else
     {
@@ -128,66 +124,102 @@
       ihi = n;
     }
 
-  // Check for scaling option.
-
-  if ((bal_job == 'S' || bal_job == 'B') && ilo != ihi)
-    {
-      F77_FCN (scaleg, SCALEG) (n, n, balanced_a_mat.fortran_vec (), 
-				n, balanced_b_mat.fortran_vec (), ilo,
-				ihi, cscale, cperm, wk.fortran_vec ());
-    }
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in reduce");
   else
     {
-      // Set scaling data to 0's.
+      Array<double> cperm (n);
+      double *pcperm = cperm.fortran_vec ();
+
+      // Check for scaling option.
+
+      if ((job == 'S' || job == 'B') && ilo != ihi)
+	{
+	  F77_XFCN (scaleg, SCALEG, (n, n, p_balanced_a_mat, n,
+				     p_balanced_b_mat, ilo, ihi,
+				     pcscale, pcperm, pwk));
+	}
+      else
+	{
+	  // Set scaling data to 0's.
+
+	  for (int i = ilo-1; i < ihi; i++)
+	    {
+	      cscale.elem (i) = 0.0;
+	      wk.elem (i, 0) = 0.0;
+	    }
+	}
+
+      if (f77_exception_encountered)
+	(*current_liboctave_error_handler) ("unrecoverable error in scaleg");
+      else
+	{
+	  // Scaleg returns exponents, not values, so...
+
+	  for (int i = ilo-1; i < ihi; i++)
+	    {
+	      cscale.elem (i) = pow (2.0, cscale.elem (i));
+	      wk.elem (i, 0) = pow (2.0, -wk.elem (i, 0));
+	    }
+
+	  // Initialize balancing matrices to identity.
+
+	  left_balancing_mat = Matrix (n, n, 0.0);
+	  for (int i = 0; i < n; i++)
+	    left_balancing_mat (i, i) = 1.0;
+
+	  right_balancing_mat = left_balancing_mat;
+
+	  double *p_right_balancing_mat = right_balancing_mat.fortran_vec ();
+	  double *p_left_balancing_mat = left_balancing_mat.fortran_vec ();
 
-      for (int tmp = ilo-1; tmp < ihi; tmp++)
-	{
-	  cscale[tmp] = 0.0;
-	  wk.elem (tmp, 0) = 0.0;
+	  // Column permutations/scaling.
+
+	  char side = 'R';
+
+	  F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, pcscale,
+				     n, p_right_balancing_mat, n, info,
+				     1L, 1L));
+    
+	  if (f77_exception_encountered)
+	    (*current_liboctave_error_handler)
+	      ("unrecoverable error in dgebak");
+	  else
+	    {
+	      // Row permutations/scaling.
+
+	      side = 'L';
+
+	      F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, pwk,
+					 n, p_left_balancing_mat, n,
+					 info, 1L, 1L));
+
+#if 0
+	      // XXX FIXME XXX --- these four lines need to be added and
+	      // debugged.  GEPBALANCE::init will work without them, though, so
+	      // here they are.
+
+	      if ((job == 'P' || job == 'B') && ilo != ihi)
+		{
+		  F77_XFCN (gradeq, GRADEQ, (n, n, p_balanced_a_mat, n,
+					     p_balanced_b_mat, ilo, ihi,
+					     pcperm, pwk));
+		}
+#endif
+
+	      if (f77_exception_encountered)
+		(*current_liboctave_error_handler)
+		  ("unrecoverable error in dgebak");
+	      else
+		{
+		  // Transpose for aa = cc*a*dd convention...
+
+		  left_balancing_mat = left_balancing_mat.transpose ();
+		}
+	    }
 	}
     }
 
-  // Scaleg returns exponents, not values, so...
-
-  for (int tmp = ilo-1; tmp < ihi; tmp++)
-    {
-      cscale[tmp] = pow (2.0, cscale[tmp]);
-      wk.elem (tmp, 0) = pow (2.0, -wk.elem (tmp, 0));
-    }
-
-  // Column permutations/scaling.
-
-  F77_FCN (dgebak, DGEBAK) (&bal_job, "R", n, ilo, ihi, cscale, n, 
-			    right_balancing_mat.fortran_vec (), n,
-			    info, 1L, 1L);
-    
-  // Row permutations/scaling.
-
-  F77_FCN (dgebak, DGEBAK) (&bal_job, "L", n, ilo, ihi,
-			    wk.fortran_vec (), n,
-			    left_balancing_mat.fortran_vec (), n,
-			    info, 1L, 1L);
-
-  // XXX FIXME XXX --- these four lines need to be added and
-  // debugged.  GEPBALANCE::init will work without them, though, so
-  // here they are.
-
-#if 0
-  if ((bal_job == 'P' || bal_job == 'B') && ilo != ihi)
-    {
-      F77_FCN (gradeq, GRADEQ) (n, n, balanced_a_mat.fortran_vec (),
-				n, balanced_b_mat.fortran_vec (), ilo,
-				ihi, cperm, wk.fortran_vec ());
-    }
-#endif
-
-  // Transpose for aa = cc*a*dd convention...
-
-  left_balancing_mat = left_balancing_mat.transpose ();
-
-  delete [] cscale;
-  delete [] cperm;
-
   return info;
 }
 
--- a/liboctave/dbleGEPBAL.h
+++ b/liboctave/dbleGEPBAL.h
@@ -66,6 +66,8 @@
       return *this;
     }
 
+  ~GEPBALANCE (void) { }
+
   Matrix balanced_a_matrix (void) const { return balanced_a_mat; }
   Matrix balanced_b_matrix (void) const { return balanced_b_mat; }