Mercurial > octave-nkf
changeset 3331:13cdcb7e5066
[project @ 1999-11-02 06:24:23 by jwe]
| author | jwe |
|---|---|
| date | Tue, 02 Nov 1999 06:25:23 +0000 |
| parents | 69136e3883bf |
| children | 7c03933635c6 |
| files | liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc |
| diffstat | 3 files changed, 165 insertions(+), 47 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc Tue Nov 02 05:33:28 1999 +0000 +++ b/liboctave/CMatrix.cc Tue Nov 02 06:25:23 1999 +0000 @@ -59,6 +59,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,18 +1601,69 @@ m.elem (i, i) -= trshift; // Preconditioning step 2: eigenvalue balancing. - - ComplexAEPBALANCE mbal (m, "B"); - m = mbal.balanced_matrix (); - ComplexMatrix d = mbal.balancing_matrix (); + // code follows development in AEPBAL + + Complex *mp = m.fortran_vec (); + int ilo, ihi, info; + + // FIXME: should pass job as a parameter in expm + char job = 'B'; + + Array<double> scale(nc); + double *pscale = scale.fortran_vec (); + + F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilo, ihi, pscale, info, + 1L, 1L)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); + return retval; + } + + // construct balancing matrices dmat, dinv + + Matrix dmat = Matrix (nc, nc, 0.0); + Matrix dinv = Matrix (nc, nc, 0.0); + + for (int i = 0; i < nc; i++) + dmat(i,i) = dinv(i,i) = 1.0; + + // pscale contains real data, so just use 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)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) ("unrecoverable error in dgebak"); + return retval; + } + + // 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)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) ("unrecoverable error in dgebak"); + return retval; + } // Preconditioning step 3: scaling. ColumnVector work (nc); double inf_norm; - F77_FCN (xzlange, XZLANGE) ("I", nc, nc, m.fortran_vec (), nc, - work.fortran_vec (), inf_norm); + F77_XFCN (xzlange, XZLANGE, ("I", nc, nc, m.fortran_vec (), nc, + work.fortran_vec (), inf_norm)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) ("unrecoverable error in zlange"); + return retval; + } int sqpow = (inf_norm > 0.0 ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0); @@ -1659,14 +1719,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.
--- a/liboctave/ChangeLog Tue Nov 02 05:33:28 1999 +0000 +++ b/liboctave/ChangeLog Tue Nov 02 06:25:23 1999 +0000 @@ -1,3 +1,9 @@ +1999-11-02 A. Scottedward Hodel <a.s.hodel@eng.auburn.edu> + + * dMatrix.cc (Matrix::expm): Do balancing here instead of using + AEPBALANCE class. + * CMatrix.cc (ComplexMatrix::expm): Likewise. + 1999-10-29 John W. Eaton <jwe@bevo.che.wisc.edu> * oct-shlib.cc, oct-shlib.h: New files.
--- a/liboctave/dMatrix.cc Tue Nov 02 05:33:28 1999 +0000 +++ b/liboctave/dMatrix.cc Tue Nov 02 06:25:23 1999 +0000 @@ -54,6 +54,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&, @@ -1313,7 +1322,7 @@ // range of poles, but avoid making stable eigenvalues unstable. // trace shift value - double trshift = 0.0; + volatile double trshift = 0.0; for (int i = 0; i < nc; i++) trshift += m.elem (i, i); @@ -1326,45 +1335,97 @@ m.elem (i, i) -= trshift; } - // Preconditioning step 2: balancing. - - AEPBALANCE mbal (m, "B"); - m = mbal.balanced_matrix (); - Matrix d = mbal.balancing_matrix (); + // Preconditioning step 2: balancing; code follows development + // in AEPBAL + + double *p_m = m.fortran_vec (); + + Array<double> scale(nc); + double *pscale = scale.fortran_vec (); + + int info, ilo, ihi; + + // both scale and permute + char job = 'B'; + + 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"); + return retval; + } + + // construct balancing matrices D, Dinv + + Matrix dmat = Matrix (nc, nc, 0.0); + Matrix dinv = Matrix (nc, nc, 0.0); + + for (int i = 0; i < nc; i++) + dmat(i,i) = dinv(i,i) = 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)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) ("unrecoverable error in dgebak"); + return retval; + } + + // 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)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) ("unrecoverable error in dgebak"); + return 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); + + F77_XFCN (xdlange, XDLANGE, ("I", nc, nc, m.fortran_vec (), nc, + work.fortran_vec (), inf_norm)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) ("unrecoverable error in dlange"); + return retval; + } 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--) { @@ -1372,38 +1433,33 @@ 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); - + retval = dpp.solve (npp, info); + // Reverse preconditioning step 3: repeated squaring. - + while (sqpow) { retval = retval * retval; sqpow--; } - + // Reverse preconditioning step 2: inverse balancing. - - retval = retval.transpose(); - d = d.transpose (); - retval = retval * d; - retval = d.solve (retval); - retval = retval.transpose (); - + retval = dmat*retval*dinv; + // Reverse preconditioning step 1: fix trace normalization. - + if (trshift > 0.0) retval = exp (trshift) * retval;