[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Expm giving NaN
From: |
John W. Eaton |
Subject: |
Re: Expm giving NaN |
Date: |
Mon, 01 Oct 2007 12:31:49 -0400 |
On 28-Sep-2007, Marco Caliari wrote:
| On Fri, 28 Sep 2007, David Bateman wrote:
|
| > Marco Caliari wrote:
| > > Dear maintainers,
| > >
| > > the patch against CMatrix.cc in 2.9.14 fixing the "expm giving NaN" bug.
| > > Here the result:
| > >
| > > octave:1> version
| > > ans = 2.9.14
| > > octave:2> expm(1000*i-10)
| > > ans = 2.5532e-05 + 3.7540e-05i
| > > octave:3> expm([1000*i-10 0;0 1000*i-10])
| > > ans =
| > >
| > > 2.5532e-05 + 3.7540e-05i 0.0000e+00 + 0.0000e+00i
| > > 0.0000e+00 + 0.0000e+00i 2.5532e-05 + 3.7540e-05i
| > >
| > >
| > > I also rewrote the computation of the rational approximation avoiding a
| > > matrix-scalar product and including a for loop. I'm not sure if it is
| > > better.
| > >
| > > Best regards,
| > >
| > > Marco
| > >
| >
| > Marco,
| >
| > John should probably comment on this patch before it is committed, but
| > can you resupply to patch as a diff with context. The type of patch you
| > sent is difficult to apply if the code has changed (and in fact it has).
| > To generate a diff with context, do something like
| >
| > diff -u liboctave/CMatrix.cc.orig liboctave/CMatrix.cc
|
| Enclosed.
|
| Best regards,
|
| Marco
| --- liboctave/CMatrix.orig 2007-09-24 14:20:52.000000000 +0200
| +++ liboctave/CMatrix.cc 2007-09-24 14:28:49.000000000 +0200
| @@ -2622,9 +2622,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
| @@ -2639,7 +2636,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;
| @@ -2724,8 +2725,13 @@
| 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++)
| + {
| + npp.elem (i, i) += padec[j];
| + dpp.elem (i, i) += minus_one_j * pade0c[j];
| + }
| + npp = m * npp;
| + dpp = m * dpp;
| minus_one_j *= -1;
| }
|
OK, I think the limit on trshift is OK. I see that the other part of
the change could speed things up, but will M always be a diagonal
matrix here? You could probably further improve performance by
using fortran_vec and avoiding elem.
jwe
- Re: Expm giving NaN,
John W. Eaton <=