expm bug fix

[A. Scottedward Hodel]

I wasn't on the octave-maintainers list until recently;
I've sent the fix below to John Eaton but thought I'd 
distribute it to the list as well.

>Actually the bug in expm can be recreated easily by doing this:
>
>octave:15> y = triu(randn(2),1); y(2,1) = 1e-32;
>octave:16> expm(y)
>ans = [](0x0)
>
>This problem also occurs in version 2.0.14.
>
>

I've solved the expm problem. I've fixed the real expm
case and will port and check the changes in the complex
expm case next.

The problem lies in the use of AEPBALANCE to construct
the balancing matrix.  

AEPBALANCE calls fortran routine dgebal which returns
balancing parameters (int) ilo, (int) ihi, and (double*) pscale.

dMatrix.cc: expm was modified to call dgebal and dgebak
 directly instead of calling AEPBALANCE.  This gives
 expm access to the balancing transformation through
 variables ilo, ihi, and pscale.
 
 Routine dgebal computes m := inv(d)*m*d where d is
        represented in variables ilo, ihi, and pscale.

 Routine dgebak can be used to directly compute the
 matrix products:
 d*a      (side = 'R')
 a*inv(d) (side = 'L')

        Unfortunately, inverse scaling requires computation of
        the matrix product

        retval = inv(d)*retval*d

        and so I chose to compute the inverse balancing via
        matrix product.  d and inv(d) are computed from dgebak as
        described above.

This is sloppy, since it's certainly possible to modify
pscale, ilo, and ihi in order to be able to get
the correct operations out of dgebak (which would be more
efficient than a matrix multiply).  However, my time's
pressed, so this is the solution I'm going with for now.

Listed below are
(1) The (simplistic) test m-file I used to check my expm fix, and
(2) the patch to dMatrix.cc and CMatrix.cc

Question: AEPBAL calls (*current_liboctave_error_handler) 
after the call to dgebal.  Is this necessary?  I did it
in dMatrix.cc but not in CMatrix.cc.  

(1) expmtest.m [taken from the user's email]
y = triu(randn(2),1)
y(2,1) = 1e-32
z = expm(y)

y = y+1e-12*sqrt(-1)
z = expm(y)

(2) Patch file to liboctave/dMatrix.cc and liboctave CMatrix.cc
===================================================================
RCS file: CDiagMatrix.cc,v
retrieving revision 1.1
diff -u -r1.1 CDiagMatrix.cc
===================================================================
RCS file: CMatrix.cc,v
retrieving revision 1.1
diff -u -r1.1 CMatrix.cc
--- 1.1 1999/06/03 14:12:42
+++ CMatrix.cc 1999/06/04 16:48:05
@@ -21,6 +21,8 @@
 
 */
 
+#undef DEBUG
+
 #if defined (__GNUG__)
 #pragma implementation
 #endif
@@ -33,6 +35,10 @@
 
 #include <iostream.h>
 
+#ifdef DEBUG
+#include <iomanip.h>
+#endif
+
 // XXX FIXME XXX
 #ifdef HAVE_SYS_TYPES_H
 #include <sys/types.h>
@@ -59,6 +65,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,10 +1607,59 @@
     m.elem (i, i) -= trshift;
 
   // Preconditioning step 2: eigenvalue balancing.
+  // code follows development in AEPBAL
 
-  ComplexAEPBALANCE mbal (m, "B");
-  m = mbal.balanced_matrix ();
-  ComplexMatrix d = mbal.balancing_matrix ();
+  #ifdef DEBUG
+  // save a copy of m for comparison later
+  ComplexMatrix m_copy = m;
+  #endif
+
+  Complex *mp = m.fortran_vec();
+  int ilo, ihi, info;
+  char job = 'B';   // FIXME: should pass job as a parameter in expm
+  Array<double> scale(nc);
+  double *pscale = scale.fortran_vec();
+
+  F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilo, ihi,
+                             pscale, info, 1L, 1L));
+
+  // construct balancing matrices dmat, dinv
+  Matrix dmat = Matrix(nc, nc, 0.0);
+  Matrix dinv = Matrix(nc, nc, 0.0);
+  int ii;
+  for( ii = 0 ; ii < nc ; ii++)
+    dmat(ii,ii) = dinv(ii,ii) = 1.0;
+
+  // pscale contains real data, so just use dgebak
+  // 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));
+
+  // 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));
+
+  #ifdef DEBUG
+  ComplexMatrix mchk = dmat*m*dinv;
+  cout << "check: balanced m, m_copy, d*m*dinv - m_copy =" << endl;
+  for( ii = 0; ii < nc ; ii++)
+    for( int jj = 0 ; jj < nc ; jj++)
+      cout 
+ << setw(3) << ii
+ << setw(3) << jj
+ << "("
+ << setw(12) << setprecision(4) << mchk(ii,jj).real() << ", "
+ << setw(12) << setprecision(4) << mchk(ii,jj).imag() << ") "
+ << "("
+ << setw(12) << setprecision(4) << m_copy(ii,jj).real() << ", "
+ << setw(12) << setprecision(4) << m_copy(ii,jj).imag() << ") "
+ << "("
+ << setw(12) << setprecision(4) << (mchk(ii,jj)-m_copy(ii,jj)).real() << ", "
+ << setw(12) << setprecision(4) << (mchk(ii,jj)-m_copy(ii,jj)).imag() << ") "
+ << endl;
+  #endif
 
   // Preconditioning step 3: scaling.
 
@@ -1659,14 +1723,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.
 
===================================================================
RCS file: boolMatrix.cc,v
retrieving revision 1.1
diff -u -r1.1 boolMatrix.cc
===================================================================
RCS file: chMatrix.cc,v
retrieving revision 1.1
diff -u -r1.1 chMatrix.cc
===================================================================
RCS file: dDiagMatrix.cc,v
retrieving revision 1.1
diff -u -r1.1 dDiagMatrix.cc
===================================================================
RCS file: dMatrix.cc,v
retrieving revision 1.1
diff -u -r1.1 dMatrix.cc
--- 1.1 1999/06/03 14:12:42
+++ dMatrix.cc 1999/06/04 16:34:21
@@ -21,6 +21,8 @@
 
 */
 
+#undef DEBUG
+
 #if defined (__GNUG__)
 #pragma implementation
 #endif
@@ -33,6 +35,10 @@
 
 #include <iostream.h>
 
+#ifdef DEBUG
+#include <iomanip.h>
+#endif
+
 #include "byte-swap.h"
 #include "dMatrix.h"
 #include "dbleAEPBAL.h"
@@ -54,6 +60,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&,
@@ -1326,86 +1341,168 @@
  m.elem (i, i) -= trshift;
     }
 
-  // Preconditioning step 2: balancing.
-
-  AEPBALANCE mbal (m, "B");
-  m = mbal.balanced_matrix ();
-  Matrix d = mbal.balancing_matrix ();
-
-  // 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);
-
-  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.
+  // Preconditioning step 2: balancing; code follows development
+  // in AEPBAL
 
-  Matrix npp (nc, nc, 0.0);
-  Matrix dpp = npp;
+  #ifdef DEBUG
+  Matrix m_copy = m;
+  cout << "dMatrix.cc: debugging code m_copy made " << endl;
+  #endif
+
+  double *p_m = m.fortran_vec();
+
+  Array<double> scale(nc);
+  double *pscale = scale.fortran_vec();
+
+  int info, ilo, ihi;
+  char job = 'B';  // both scale and permute
+
+  int ii;
+  #ifdef DEBUG
+  int jj;
+  cout << "Original m" << endl;
+  for (ii=0 ; ii < nc ; ii++)
+  {
+    for (jj=0 ; jj < nc ; jj++)
+      cout << setw(12) << setprecision(4) << m(ii,jj) << " ";
+    cout << endl;
+  }
+  #endif
 
-  // Now powers a^8 ... a^1.
-
-  int minus_one_j = -1;
-  for (int j = 7; j >= 0; j--)
-    {
-      npp = m * npp + m * padec[j];
-      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);
-
-  // Reverse preconditioning step 3: repeated squaring.
-
-  while (sqpow)
-    {
-      retval = retval * retval;
-      sqpow--;
+  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");
+  else
+  {
+    #ifdef DEBUG
+    cout << "Balanced m" << endl;
+    for (ii=0 ; ii < nc ; ii++)
+    {
+      for (jj=0 ; jj < nc ; jj++)
+        cout << setw(12) << setprecision(4) << m(ii,jj) << " ";
+      cout << endl;
+    }
+    cout << "Check: trying to construct m_copy from scale, and m" << endl;
+    Matrix bal_m_copy = m;
+    #endif
+
+    // construct balancing matrices D, Dinv; initialize dmat, dinv to identity
+    Matrix dmat = Matrix(nc,nc,0.0);
+    Matrix dinv = Matrix(nc,nc,0.0);
+    for(ii=0 ; ii < nc ; ii++)
+      dmat(ii,ii) = dinv(ii,ii) = 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));
+
+    // 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));
+
+    #ifdef DEBUG
+    cout << "d,dinv,d*dinv" << endl;
+    Matrix d_dinv = dmat*dinv;
+    for (ii=0 ; ii < nc ; ii++)
+    {
+      for (jj = 0; jj < nc ; jj++)
+        cout  
+ << setw(3) << ii 
+ << setw(3) << jj 
+ << setw(12) << setprecision(4) << dmat(ii,jj) << " "
+ << setw(12) << setprecision(4) << dinv(ii,jj) << " "
+ << setw(12) << setprecision(4) << d_dinv(ii,jj) << " "
+ << endl;
+    }
+    
+    bal_m_copy = dmat*bal_m_copy*dinv;
+    
+    // print out difference
+    cout << "dMatrix:expm: dinv*m*d - orig_m = " << endl;
+    for (ii=0 ; ii < nc ; ii++)
+    {
+      for (jj = 0; jj < nc ; jj++)
+ cout << setw(3) << ii 
+ << setw(3) << jj 
+ << setw(12) << setprecision(4) << bal_m_copy(ii,jj) << " "
+ << setw(12) << setprecision(4) << m_copy(ii,jj) << " "
+ << setw(12) << setprecision(4) << bal_m_copy(ii,jj) - m_copy(ii,jj)
+ << endl;
     }
-
-  // Reverse preconditioning step 2: inverse balancing.
+    #endif
 
-  retval = retval.transpose();
-  d = d.transpose ();
-  retval = retval * d;
-  retval = d.solve (retval);
-  retval = retval.transpose ();
-
-  // Reverse preconditioning step 1: fix trace normalization.
-
-  if (trshift > 0.0)
-    retval = exp (trshift) * 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);
+  
+    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--)
+      {
+        npp = m * npp + m * padec[j];
+        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);
+  
+    // Reverse preconditioning step 3: repeated squaring.
+  
+    while (sqpow)
+      {
+        retval = retval * retval;
+        sqpow--;
+      }
+  
+    // Reverse preconditioning step 2: inverse balancing.
+    retval = dmat*retval*dinv;
+  
+    // Reverse preconditioning step 1: fix trace normalization.
+  
+    if (trshift > 0.0)
+      retval = exp (trshift) * retval;
+  } // else of dgebal exception test
 
   return retval;
 }


A S Hodel Assoc. Prof. Dept Elect Eng, Auburn Univ,AL  36849-5201
On leave at NASA Marshall Space Flight Center (256) 544-1426
Address until 15 Mar 2000:Mail Code TD-55, MSFC, Alabama, 35812
http://www.eng.auburn.edu/~scotte