Mercurial > hg > octave-lyh
changeset 6958:a18c784ae599
[project @ 2007-10-04 19:21:23 by dbateman]
| author | dbateman |
|---|---|
| date | Thu, 04 Oct 2007 19:21:23 +0000 |
| parents | 768a19157591 |
| children | 47f4f4e88166 |
| files | liboctave/CMatrix.cc liboctave/ChangeLog |
| diffstat | 2 files changed, 20 insertions(+), 6 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -2627,9 +2627,6 @@ ComplexMatrix m = *this; - if (numel () == 1) - return ComplexMatrix (1, 1, exp (m(0))); - octave_idx_type nc = columns (); // Preconditioning step 1: trace normalization to reduce dynamic @@ -2644,7 +2641,11 @@ trshift /= nc; if (trshift.real () < 0.0) - trshift = trshift.imag (); + { + trshift = trshift.imag (); + if (trshift.real () > 709.0) + trshift = 709.0; + } for (octave_idx_type i = 0; i < nc; i++) m.elem (i, i) -= trshift; @@ -2722,15 +2723,23 @@ // npp, dpp: pade' approx polynomial matrices. ComplexMatrix npp (nc, nc, 0.0); + Complex *pnpp = npp.fortran_vec (); ComplexMatrix dpp = npp; + Complex *pdpp = dpp.fortran_vec (); // Now powers a^8 ... a^1. int minus_one_j = -1; for (octave_idx_type j = 7; j >= 0; j--) { - npp = m * npp + m * padec[j]; - dpp = m * dpp + m * (minus_one_j * padec[j]); + for (octave_idx_type i = 0; i < nc; i++) + { + octave_idx_type k = i * nc + i; + pnpp [k] = pnpp [k] + padec [j]; + pdpp [k] = pdpp [k] + minus_one_j * padec [j]; + } + npp = m * npp; + dpp = m * dpp; minus_one_j *= -1; }
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,8 @@ +2007-10-04 Marco Caliari <mcaliari@math.unipd.it> + + * CMatrix.cc (ComplexMatrix::expm): Limit shift to values less + than log(realmax) to avoid issues with NaN. + 2007-10-01 John W. Eaton <jwe@octave.org> * oct-time.cc (octave_strptime::init): Call mktime to propertly