Re: Expm giving NaN

[Marco Caliari]

On Mon, 1 Oct 2007, John W. Eaton wrote:

> 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?  

M is not supposed to be diagonal. I just used

M*P+M*a=M*(P+I*a)

> You could probably further improve performance by
> using fortran_vec and avoiding elem.

I'm sorry, but I don't know how to use it. I used elem because it was 
already used twice in expm.

Marco