Re: expm bug fix

[A. Scottedward Hodel]

I agree that the solution is inelegant (about 1/2 step above a
call to solve()).  I'll need to take some time to get the format
of the pscale,ilo,ihi data to do this properly.  The scaling
data (between ilo and ihi) is easy; in octave notation this would
be
  pscale(ilo:ihi) = 1 ./ pscale(ilo:ihi).
The permutation data from 1:(ilo-1) and (ihi+1):n needs to be 
"inverted."  This probably isn't too hard, but I just need to get 
the time to see how it's done and then replace the data with
an inverse permutation.  [There's got to be an easy way to do
this, but I haven't done combinatorics like this in ages.]

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

----------
>From: "Ross A. Lippert" <address@hidden>
>To: "A. Scottedward Hodel" <address@hidden>,address@hidden
>Subject: Re: expm bug fix
>Date: Tue, Jun 8, 1999, 10:35 AM
>

>De-balancing the result by a matrix multiply and inversion
>is inelegant.  If you finish your patch to CMatrix.cc then I'll
>take what you did to expm and use dgebak to do the de-balancing
>and get new patches out in a couple of days.
>
>However, there is still the problem in dMatrix::solve() in which
>if will return a [] when non-[] data is given to it.  That needs 
>to be dealt with too, and it is a pretty fundamental problem.
>
>-r
>
>"A. Scottedward Hodel" wrote:
>> 
>> 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
>