[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;
}