Re: Expm giving NaN

[John W. Eaton]

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