[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: expm bug fix
From: |
Ross A. Lippert |
Subject: |
Re: expm bug fix |
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
> @@ -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
- expm bug fix, A. Scottedward Hodel, 1999/06/08
- Re: expm bug fix,
Ross A. Lippert <=