octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Expm giving NaN


From: David Bateman
Subject: Re: Expm giving NaN
Date: Tue, 02 Oct 2007 10:33:34 +0200
User-agent: Thunderbird 1.5.0.7 (X11/20060921)

Marco Caliari wrote:
>> 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.
>   
Assuming the algorithm itself is ok, what John means is something like
the attached patch. Note there appears to have been a typo in your patch
as "pade0c" should read "padec".. I tested this against a few examples
and compared with matlab2007a and it appears to be correct.

The test I ran was

n = 4; m = 20; for i = 1:m, a{i} = 1000*randn(n,n) + 10*randn(n,n); end;
save("-v7","expm_test.mat","a")

then I ran

load expm_test.mat
for i=1:length(a), expm(a{i}), end

In both Octave 2.9.14+ and matlab2007a. Where numerical results were
returned they were nearly identical.. In some cases matlab2007a returned
a matrix of NaN values where Octave returned a matrix of Inf values..
For the numerical values there is a relative difference between between
matlab and Octave on the order of 1e3 * EPS, which I consider not bad
for numerical values on the order of 10e270.

Regards
David

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

*** ./liboctave/CMatrix.cc.orig15       2007-10-02 09:59:18.626311084 +0200
--- ./liboctave/CMatrix.cc      2007-10-02 10:10:21.665670245 +0200
***************
*** 2627,2635 ****
  
    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
--- 2627,2632 ----
***************
*** 2644,2650 ****
    trshift /= nc;
  
    if (trshift.real () < 0.0)
!     trshift = trshift.imag ();
  
    for (octave_idx_type i = 0; i < nc; i++)
      m.elem (i, i) -= trshift;
--- 2641,2651 ----
    trshift /= nc;
  
    if (trshift.real () < 0.0)
!     {
!       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,2736 ****
    // npp, dpp: pade' approx polynomial matrices.
  
    ComplexMatrix npp (nc, nc, 0.0);
    ComplexMatrix dpp = npp;
  
    // 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]);
        minus_one_j *= -1;
      }
  
--- 2723,2745 ----
    // 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--)
      {
!       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;
      }
  

reply via email to

[Prev in Thread] Current Thread [Next in Thread]