From octave-maintainers-request at bevo dot che dot wisc dot edu Tue Jun 8 11:02:18 1999 Subject: Re: expm bug fix From: "A. Scottedward Hodel" To: octave-maintainers at bevo dot che dot wisc dot edu Date: Tue, 08 Jun 1999 11:01:32 -0500 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" >To: "A. Scottedward Hodel" ,octave-maintainers@bevo dot che dot wisc dot edu >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 >> 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 >