# HG changeset patch
# User jwe
# Date 941523923 0
# Node ID 13cdcb7e5066fa33576fd28aa9b9c4546ec8765a
# Parent  69136e3883bfcf73f572df756caba868ad05dd57
[project @ 1999-11-02 06:24:23 by jwe]

diff -r 69136e3883bf -r 13cdcb7e5066 liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc	Tue Nov 02 05:33:28 1999 +0000
+++ b/liboctave/CMatrix.cc	Tue Nov 02 06:25:23 1999 +0000
@@ -59,6 +59,15 @@
 
 extern "C"
 {
+  int F77_FCN (zgebal, ZGEBAL) (const char*, const int&, Complex*,
+                                const int&, int&, int&, double*, int&,
+                                long, long);
+
+  int F77_FCN (dgebak, DGEBAK) (const char*, const char*, const int&,
+                                const int&, const int&, double*,
+                                const int&, double*, const int&,
+                                int&, long, long);
+
   int F77_FCN (zgemm, ZGEMM) (const char*, const char*, const int&,
 			      const int&, const int&, const Complex&,
 			      const Complex*, const int&,
@@ -1592,18 +1601,69 @@
     m.elem (i, i) -= trshift;
 
   // Preconditioning step 2: eigenvalue balancing.
-
-  ComplexAEPBALANCE mbal (m, "B");
-  m = mbal.balanced_matrix ();
-  ComplexMatrix d = mbal.balancing_matrix ();
+  // code follows development in AEPBAL
+
+  Complex *mp = m.fortran_vec ();
+  int ilo, ihi, info;
+
+  // FIXME: should pass job as a parameter in expm
+  char job = 'B';
+
+  Array<double> scale(nc);
+  double *pscale = scale.fortran_vec ();
+
+  F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilo, ihi, pscale, info,
+			     1L, 1L));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler) ("unrecoverable error in zgebal");
+      return retval;
+    }
+
+  // construct balancing matrices dmat, dinv
+
+  Matrix dmat = Matrix (nc, nc, 0.0);
+  Matrix dinv = Matrix (nc, nc, 0.0);
+
+  for (int i = 0; i < nc; i++)
+    dmat(i,i) = dinv(i,i) = 1.0;
+
+  // pscale contains real data, so just use dgebak, dside=R => dmat := D*dmat
+  char dside = 'R';
+  F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
+			     dmat.fortran_vec(), nc, info, 1L, 1L));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler) ("unrecoverable error in dgebak");
+      return retval;
+    }
+
+  // dgebak, dside=L => dinv := dinv*D^{-1}
+  dside = 'L';
+  F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
+			     dinv.fortran_vec(), nc, info, 1L, 1L));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler) ("unrecoverable error in dgebak");
+      return retval;
+    }
 
   // Preconditioning step 3: scaling.
 
   ColumnVector work (nc);
   double inf_norm;
 
-  F77_FCN (xzlange, XZLANGE) ("I", nc, nc, m.fortran_vec (), nc,
-			      work.fortran_vec (), inf_norm);
+  F77_XFCN (xzlange, XZLANGE, ("I", nc, nc, m.fortran_vec (), nc,
+			       work.fortran_vec (), inf_norm));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler) ("unrecoverable error in zlange");
+      return retval;
+    }
 
   int sqpow = (inf_norm > 0.0
 	       ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0);
@@ -1659,14 +1719,10 @@
     }
 
   // Reverse preconditioning step 2: inverse balancing.
-  // XXX FIXME XXX -- should probably do this with Lapack calls
-  // instead of a complete matrix inversion.
-
-  retval = retval.transpose ();
-  d = d.transpose ();
-  retval = retval * d;
-  retval = d.solve (retval);
-  retval = retval.transpose ();
+  // FIXME: need to figure out how to get dgebak to do
+  // dmat*retval*dinv instead of dinv*retval*dmat
+  // This works for now
+  retval = dmat*retval*dinv;
 
   // Reverse preconditioning step 1: fix trace normalization.
 
diff -r 69136e3883bf -r 13cdcb7e5066 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Tue Nov 02 05:33:28 1999 +0000
+++ b/liboctave/ChangeLog	Tue Nov 02 06:25:23 1999 +0000
@@ -1,3 +1,9 @@
+1999-11-02  A. Scottedward Hodel <a.s.hodel@eng.auburn.edu>
+
+	* dMatrix.cc (Matrix::expm): Do balancing here instead of using
+	AEPBALANCE class.
+	* CMatrix.cc (ComplexMatrix::expm): Likewise.
+
 1999-10-29  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* oct-shlib.cc, oct-shlib.h: New files.
diff -r 69136e3883bf -r 13cdcb7e5066 liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc	Tue Nov 02 05:33:28 1999 +0000
+++ b/liboctave/dMatrix.cc	Tue Nov 02 06:25:23 1999 +0000
@@ -54,6 +54,15 @@
 
 extern "C"
 {
+  int F77_FCN (dgebal, DGEBAL) (const char*, const int&, double*,
+                                const int&, int&, int&, double*,
+                                int&, long, long);
+
+  int F77_FCN (dgebak, DGEBAK) (const char*, const char*, const int&,
+                                const int&, const int&, double*,
+                                const int&, double*, const int&,
+                                int&, long, long);
+
   int F77_FCN (dgemm, DGEMM) (const char*, const char*, const int&,
 			      const int&, const int&, const double&,
 			      const double*, const int&,
@@ -1313,7 +1322,7 @@
   // range of poles, but avoid making stable eigenvalues unstable.
 
   // trace shift value
-  double trshift = 0.0;
+  volatile double trshift = 0.0;
 
   for (int i = 0; i < nc; i++)
     trshift += m.elem (i, i);
@@ -1326,45 +1335,97 @@
 	m.elem (i, i) -= trshift;
     }
 
-  // Preconditioning step 2: balancing.
-
-  AEPBALANCE mbal (m, "B");
-  m = mbal.balanced_matrix ();
-  Matrix d = mbal.balancing_matrix ();
+  // Preconditioning step 2: balancing; code follows development
+  // in AEPBAL
+
+  double *p_m = m.fortran_vec ();
+
+  Array<double> scale(nc);
+  double *pscale = scale.fortran_vec ();
+
+  int info, ilo, ihi;
+
+  // both scale and permute
+  char job = 'B';
+
+  F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilo, ihi, pscale, info,
+			     1L, 1L));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler) ("unrecoverable error in dgebal");
+      return retval;
+    }
+
+  // construct balancing matrices D, Dinv
+
+  Matrix dmat = Matrix (nc, nc, 0.0);
+  Matrix dinv = Matrix (nc, nc, 0.0);
+
+  for (int i = 0; i < nc; i++)
+    dmat(i,i) = dinv(i,i) = 1.0;
+
+  // dgebak, dside=R => dmat := D*dmat
+  char dside = 'R';
+  F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
+			     dmat.fortran_vec(), nc, info, 1L, 1L));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler) ("unrecoverable error in dgebak");
+      return retval;
+    }
+
+  // dgebak, dside=L => dinv := dinv*D^{-1}
+  dside = 'L';
+  F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
+			     dinv.fortran_vec(), nc, info, 1L, 1L));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler) ("unrecoverable error in dgebak");
+      return retval;
+    }
 
   // Preconditioning step 3: scaling.
-
+  
   ColumnVector work(nc);
   double inf_norm;
-
-  F77_FCN (xdlange, XDLANGE) ("I", nc, nc, m.fortran_vec (), nc,
-			      work.fortran_vec (), inf_norm);
+  
+  F77_XFCN (xdlange, XDLANGE, ("I", nc, nc, m.fortran_vec (), nc,
+			       work.fortran_vec (), inf_norm));
+  
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler) ("unrecoverable error in dlange");
+      return retval;
+    }
 
   int sqpow = (int) (inf_norm > 0.0
 		     ? (1.0 + log (inf_norm) / log (2.0))
 		     : 0.0);
-
+  
   // Check whether we need to square at all.
-
+  
   if (sqpow < 0)
     sqpow = 0;
-
+  
   if (sqpow > 0)
     {
       double scale_factor = 1.0;
       for (int i = 0; i < sqpow; i++)
 	scale_factor *= 2.0;
-
+  
       m = m / scale_factor;
     }
-
+  
   // npp, dpp: pade' approx polynomial matrices.
-
+  
   Matrix npp (nc, nc, 0.0);
   Matrix dpp = npp;
-
+  
   // Now powers a^8 ... a^1.
-
+  
   int minus_one_j = -1;
   for (int j = 7; j >= 0; j--)
     {
@@ -1372,38 +1433,33 @@
       dpp = m * dpp + m * (minus_one_j * padec[j]);
       minus_one_j *= -1;
     }
-
+  
   // Zero power.
-
+  
   dpp = -dpp;
   for (int j = 0; j < nc; j++)
     {
       npp.elem (j, j) += 1.0;
       dpp.elem (j, j) += 1.0;
     }
-
+  
   // Compute pade approximation = inverse (dpp) * npp.
 
-  retval = dpp.solve (npp);
-
+  retval = dpp.solve (npp, info);
+  
   // Reverse preconditioning step 3: repeated squaring.
-
+  
   while (sqpow)
     {
       retval = retval * retval;
       sqpow--;
     }
-
+  
   // Reverse preconditioning step 2: inverse balancing.
-
-  retval = retval.transpose();
-  d = d.transpose ();
-  retval = retval * d;
-  retval = d.solve (retval);
-  retval = retval.transpose ();
-
+  retval = dmat*retval*dinv;
+  
   // Reverse preconditioning step 1: fix trace normalization.
-
+  
   if (trshift > 0.0)
     retval = exp (trshift) * retval;