Re: Request Fortran code for Matrix exponential for a real matrix, with accuracy estimate

[Robin Vowels]

On Nov 13, 9:41 am, Zhu Wang <zhu.w...@gmail.com> wrote:
>
> I was looking for an open source Fortran code to compute Matrix
> exponential for a real matrix, with accuracy estimate, something likehttp://www.slicot.org/shared/doc/MB05OD.html. I tried to search Scilab
> but was not successful.

This might interest you. It's an edited version of W. Everling's.

/* W. Everling, "On the Evaluation of e**AT by power series", */
/* Proc. IEEE, Vol. 55, 1967, p. 413. */

EXPA: PROCEDURE (A, ERROR);
/* ERROR is relative precision. */
declare A(*,*) float (18), ERROR FLOAT (18),
(I, I1, I2, N) fixed binary;
declare ((A_H, A_I, A_J)(HBOUND(A,1),HBOUND(A,2)),
ANORM, ATEST, ETEST ) float (18);

/* Notation corresponds to that of Liou as follows: */
/* A AT */
/* I k */
/* I1 K */
/* A_H auxiliary storage for argument */
/* A_I (AT)**k/k! */
/* ATEST 10**d * |r(i,j)| */

ON ERROR SNAP SYSTEM;

N = DIM(A,1);

ANORM = RSNORM(A);
I1 = MAX(2*ANORM, 2); I2 = I1/2;
ANORM = 2*ANORM / ERROR;
A_H, A_I = A;
/* Form I + A */
DO I = 1 TO N;
A(I,I) = A(I,I) + 1;
END;

/* Loop to form A**2/2! + A**3/3! + A**4/4! + ... */
DO I = 2 BY 1;
A_J = A_H/I;
CALL MATMULT (A_I, A_J);
A = A + A_I;
IF I = I1 THEN
DO;
I1 = I1 + I2;
ATEST = RSNORM(A_I)*ANORM/(I+1);
ETEST = RSNORM(A)*ERROR**2;
IF ALL(ATEST < ABS(A) | ABS(A) < ETEST) THEN
GO TO EXPA_END;
END;
END;

EXPA_END:

END EXPA;

/* Forms the maximum_row_sum norm of A. */

RSNORM: PROCEDURE (A) RETURNS (FLOAT (18));
DECLARE A(*,*) FLOAT (18), ANORM FLOAT (18) INITIAL (0),
I FIXED BINARY;

DO I = LBOUND(A,1) TO HBOUND(A,1);
ANORM = MAX(ANORM, SUM(ABS(A(I,*))) );
END;
RETURN (ANORM);
END RSNORM;

[Arjan]

Van Loan has written an excellent book about "900 wrong ways to
calculate the matrix exponantial" (or something similar). I am afraid
that the above code is too simple to catch the subtle problems
described in that book. I'm still curious to know what the OP wants to
use the matrix exponential for! Can you give a hint? I personally use
it to estimate the decay of radionuclides and the acquired dose. With
halflives of some relevant nuclides many orders of magnitude apart,
the condition number of the matrix can be very poor. That's a yummy
for most matrix exponential routines (not)!

Arjan

[Robin Vowels]

On Nov 17, 5:29 am, Arjan <arjan.van.d...@rivm.nl> wrote:
> Van Loan has written an excellent book about "900 wrong ways to
> calculate the matrix exponantial" (or something similar). I am afraid
> that the above code is too simple to catch the subtle problems
> described in that book.

That's the well-known paper by Moler and Van Loan,
and the title is
"Nineteen dubious ways to compute the exponential of a matrix".
There's a follow-up paper, too.