From octave-maintainers-request at bevo dot che dot wisc dot edu Tue Jun 8 10:36:29 1999 Subject: Re: expm bug fix From: "Ross A. Lippert" To: "A. Scottedward Hodel" , octave-maintainers@bevo.che.wisc.edu Date: Tue, 08 Jun 1999 09:35:53 -0600 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 > at @ -21,6 +21,8 @@ > > */ > > +#undef DEBUG > + > #if defined (__GNUG__) > #pragma implementation > #endif > at @ -33,6 +35,10 @@ > > #include > > +#ifdef DEBUG > +#include > +#endif > + > // XXX FIXME XXX > #ifdef HAVE_SYS_TYPES_H > #include > at @ -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&, > at @ -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 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. > > at @ -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 > at @ -21,6 +21,8 @@ > > */ > > +#undef DEBUG > + > #if defined (__GNUG__) > #pragma implementation > #endif > at @ -33,6 +35,10 @@ > > #include > > +#ifdef DEBUG > +#include > +#endif > + > #include "byte-swap.h" > #include "dMatrix.h" > #include "dbleAEPBAL.h" > at @ -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&, > at @ -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 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