# HG changeset patch
# User jwe
# Date 745014485 0
# Node ID 2cd2476fb32d813239051b9363924e44f4591a58
# Parent  1b59f5c6c3709949dd6e27ab3a90a7b728555163
[project @ 1993-08-10 20:28:05 by jwe]
Add init functions for balancing classes.

diff -r 1b59f5c6c370 -r 2cd2476fb32d liboctave/Matrix-ext.cc
--- a/liboctave/Matrix-ext.cc	Tue Aug 10 20:18:30 1993 +0000
+++ b/liboctave/Matrix-ext.cc	Tue Aug 10 20:28:05 1993 +0000
@@ -29,6 +29,206 @@
 #include "mx-inlines.cc"
 
 /*
+ * AEPBALANCE operations
+ */
+
+int
+AEPBALANCE::init (const Matrix& a, const char *balance_job)
+{
+  if (a.nr != a.nc)
+    FAIL;
+
+  int n = a.nc;
+
+// Parameters for balance call.
+
+  int info;
+  int ilo;
+  int ihi;
+  double *scale = new double [n];
+
+// Copy matrix into local structure.
+
+  balanced_mat = a;
+
+  F77_FCN (dgebal) (balance_job, &n, balanced_mat.fortran_vec (), 
+		    &n, &ilo, &ihi, scale, &info, 1L, 1L);
+
+// Initialize balancing matrix to identity.
+
+  balancing_mat = Matrix (n, n, 0.0);
+  for (int i = 0; i < n; i++)
+    balancing_mat.elem (i ,i) = 1.0;
+
+  F77_FCN (dgebak) (balance_job, "R", &n, &ilo, &ihi, scale, &n, 
+		    balancing_mat.fortran_vec (), &n, &info, 1L, 1L);
+
+  delete [] scale;
+
+  return info;
+}
+
+int
+ComplexAEPBALANCE::init (const ComplexMatrix& a, const char *balance_job)
+{
+
+  int n = a.nc;
+
+// Parameters for balance call.
+
+  int info;
+  int ilo;
+  int ihi;
+  double *scale = new double [n];
+
+// Copy matrix into local structure.
+
+  balanced_mat = a;
+
+  F77_FCN (zgebal) (balance_job, &n, balanced_mat.fortran_vec (),
+		    &n, &ilo, &ihi, scale, &info, 1L, 1L);
+
+// Initialize balancing matrix to identity.
+
+  balancing_mat = Matrix (n, n, 0.0);
+  for (int i = 0; i < n; i++)
+    balancing_mat (i, i) = 1.0;
+
+  F77_FCN (zgebak) (balance_job, "R", &n, &ilo, &ihi, scale, &n, 
+		    balancing_mat.fortran_vec(), &n, &info, 1L, 1L);
+
+  delete [] scale;
+
+  return info;
+}
+
+/*
+ * GEPBALANCE operations
+ */
+
+int
+GEPBALANCE::init (const Matrix& a, const Matrix& b, const char *balance_job)
+{
+  if (a.nr != a.nc || a.nr != b.nr || b.nr != b.nc)
+    FAIL;
+
+  int n = a.nc;
+
+// Parameters for balance call.
+
+  int info;
+  int ilo;
+  int ihi;
+  double *cscale = new double [n];
+  double *cperm = new double [n];
+  Matrix wk (n, 6, 0.0);
+
+// Back out the permutations:
+//
+// cscale contains the exponents of the column scaling factors in its 
+// ilo through ihi locations and the reducing column permutations in 
+// its first ilo-1 and its ihi+1 through n locations.
+//
+// cperm contains the column permutations applied in grading the a and b 
+// submatrices in its ilo through ihi locations.
+//
+// wk contains the exponents of the row scaling factors in its ilo 
+// through ihi locations, the reducing row permutations in its first 
+// ilo-1 and its ihi+1 through n locations, and the row permutations
+// applied in grading the a and b submatrices in its n+ilo through 
+// n+ihi locations.
+  
+// Copy matrices into local structure.
+
+  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;
+
+// Check for permutation option.
+
+  if (*balance_job == 'P' || *balance_job == 'B')
+    {
+      F77_FCN(reduce)(&n, &n, balanced_a_mat.fortran_vec (),
+		      &n, balanced_b_mat.fortran_vec (), &ilo, &ihi,
+		      cscale, wk.fortran_vec ());
+    }
+  else
+    {
+
+// Set up for scaling later.
+
+      ilo = 1;
+      ihi = n;
+    }
+
+// Check for scaling option.
+
+  if ((*balance_job == 'S' || *balance_job == 'B') && ilo != ihi)
+    {
+      F77_FCN(scaleg)(&n, &n, balanced_a_mat.fortran_vec (), 
+		      &n, balanced_b_mat.fortran_vec (), &ilo, &ihi,
+		      cscale, cperm, wk.fortran_vec ());
+    }
+  else
+    {
+
+// Set scaling data to 0's.
+
+      for (int tmp = ilo-1; tmp < ihi; tmp++)
+	{
+	  cscale[tmp] = 0.0;
+	  wk.elem(tmp,0) = 0.0;
+	}
+    }
+
+// 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) (balance_job, "R", &n, &ilo, &ihi, cscale, &n, 
+		    right_balancing_mat.fortran_vec (), &n, &info, 1L,
+		    1L);
+    
+// Row permutations/scaling.
+
+  F77_FCN (dgebak) (balance_job, "L", &n, &ilo, &ihi, &wk.elem (0, 0), &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 ((*balance_job == 'P' || *balance_job == 'B') && ilo != ihi)
+    {
+      F77_FCN (gradeq) (&n, &n, balanced_a_mat.fortran_vec (),
+			&n, balanced_b_mat.fortran_vec (), &ilo, &ihi,
+			cperm, &wk.elem (0, 1));
+    }
+#endif
+
+// Transpose for aa = cc*a*dd convention...
+  left_balancing_mat = left_balancing_mat.transpose ();
+
+  delete [] cscale;
+  delete [] cperm;
+
+  return info;
+}
+
+/*
  * HESS stuff
  */
 
@@ -38,7 +238,7 @@
   if (a.nr != a.nc)
     FAIL;
 
-  char jobbal = 'S';
+  char jobbal = 'N';
   char side = 'R';
 
   int n = a.nc;
@@ -96,7 +296,7 @@
    if (a.nr != a.nc)
      FAIL;
 
-   char job = 'S';
+   char job = 'N';
    char side = 'R';
 
    int n = a.nc;