Re: Expm giving NaN

[John W. Eaton]

On  2-Oct-2007, David Bateman wrote:

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

OK, will you please check this in?

Thanks,

jwe