Mercurial > hg > octave-terminal
changeset 11842:0b9c56b6bf0e release-3-0-x
partially sync Matrix::expm and ComplexMatrix::expm with development repo
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Fri, 19 Sep 2008 11:29:51 +0200 |
parents | 43fccbab412a |
children | d1b8260dbc76 |
files | liboctave/CMatrix.cc liboctave/dMatrix.cc |
diffstat | 2 files changed, 22 insertions(+), 10 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -2896,6 +2896,9 @@ if (sqpow > 0) { + if (sqpow > 1023) + sqpow = 1023; + double scale_factor = 1.0; for (octave_idx_type i = 0; i < sqpow; i++) scale_factor *= 2.0; @@ -2942,8 +2945,11 @@ // Compute pade approximation = inverse (dpp) * npp. - retval = dpp.solve (npp); - + retval = dpp.solve (npp, info); + + if (info < 0) + return retval; + // Reverse preconditioning step 3: repeated squaring. while (sqpow) @@ -2977,12 +2983,13 @@ } // construct inverse balancing permutation vector + Array<octave_idx_type> invpvec (nc); for (octave_idx_type i = 0; i < nc; i++) invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method OCTAVE_QUIT; - tmpMat = retval; + ComplexMatrix tmpMat = retval; for (octave_idx_type i = 0; i < nc; i++) for (octave_idx_type j = 0; j < nc; j++) retval(i,j) = tmpMat(invpvec(i),invpvec(j)); @@ -3002,13 +3009,12 @@ } // construct inverse balancing permutation vector - Array<octave_idx_type> invpvec (nc); for (octave_idx_type i = 0; i < nc; i++) invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method OCTAVE_QUIT; - ComplexMatrix tmpMat = retval; + tmpMat = retval; for (octave_idx_type i = 0; i < nc; i++) for (octave_idx_type j = 0; j < nc; j++) retval(i,j) = tmpMat(invpvec(i),invpvec(j));
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -2523,10 +2523,13 @@ if (sqpow > 0) { + if (sqpow > 1023) + sqpow = 1023; + double scale_factor = 1.0; for (octave_idx_type i = 0; i < sqpow; i++) scale_factor *= 2.0; - + m = m / scale_factor; } @@ -2570,7 +2573,10 @@ // Compute pade approximation = inverse (dpp) * npp. retval = dpp.solve (npp, info); - + + if (info < 0) + return retval; + // Reverse preconditioning step 3: repeated squaring. while (sqpow) @@ -2602,12 +2608,13 @@ } // construct inverse balancing permutation vector + Array<octave_idx_type> invpvec (nc); for (octave_idx_type i = 0; i < nc; i++) invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method OCTAVE_QUIT; - tmpMat = retval; + Matrix tmpMat = retval; for (octave_idx_type i = 0; i < nc; i++) for (octave_idx_type j = 0; j < nc; j++) retval(i,j) = tmpMat(invpvec(i),invpvec(j)); @@ -2627,13 +2634,12 @@ } // construct inverse balancing permutation vector - Array<octave_idx_type> invpvec (nc); for (octave_idx_type i = 0; i < nc; i++) invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method OCTAVE_QUIT; - Matrix tmpMat = retval; + tmpMat = retval; for (octave_idx_type i = 0; i < nc; i++) for (octave_idx_type j = 0; j < nc; j++) retval(i,j) = tmpMat(invpvec(i),invpvec(j));