Matrix exponentiation in FORTRAN
Maayan
I hope I've come to the right place - please redirect me if there's a
more appropriate forum to be asking these questions.
I am trying to figure out how to compute the exponential of a matrix
in FORTRAN - something similar to expm used in MATLAB, and MatrixExp
in Mathematica. I've gone through LAPACK's online documentation but
couldn't find an answer. Is there some way of getting this done? I'd
appreciate your kind help with this.
Maayan.
Ian Bush
Hi,
As if by magic, Maayan appeared !
Not an expert here, but one way via LAPACK is note that
if your matrix A is diagonalizable then
A=QDQ^(-1)
where Q are the evecs and D is a diagonal matrix with the evals down the
diag. From this you can show that
exp(A)=Q exp(D) Q^(-1)
and the exponential of a diagonal matrix is trivial.
Again I'm not an expert and better algorithms may exist,
Ian
Richard Maine
> As if by magic, Maayan appeared !
> > I hope I've come to the right place - please redirect me if there's a
> > more appropriate forum to be asking these questions.
> >
> > I am trying to figure out how to compute the exponential of a matrix
> > in FORTRAN - something similar to expm used in MATLAB, and MatrixExp
> > in Mathematica. I've gone through LAPACK's online documentation but
> > couldn't find an answer. Is there some way of getting this done? I'd
> > appreciate your kind help with this.
> >
> Again I'm not an expert and better algorithms may exist,
Quite a few ways exist. Pretty much all of them have "issues", depending
on the details of the needed usage. *THE* reference work on the question
is Moler's 1978 paper "Ninteen Dubious Ways to Compute the Exponential
of a Matrix". I suggest just googling the phrase "Ninteen dubious ways".
Worked for me, including finding an update paper from 2003.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
none
> Ian Bush <I.J....@ku.ca.ld> wrote:
>
>> As if by magic, Maayan appeared !
>
>> > I hope I've come to the right place - please redirect me if there's a
>> > more appropriate forum to be asking these questions.
>> >
>> > I am trying to figure out how to compute the exponential of a matrix
>> > in FORTRAN - something similar to expm used in MATLAB, and MatrixExp
>> > in Mathematica. I've gone through LAPACK's online documentation but
>> > couldn't find an answer. Is there some way of getting this done? I'd
>> > appreciate your kind help with this.
>> >
>
>> Again I'm not an expert and better algorithms may exist,
>
> Quite a few ways exist. Pretty much all of them have "issues", depending
> on the details of the needed usage. *THE* reference work on the question
> is Moler's 1978 paper "Ninteen Dubious Ways to Compute the Exponential
> of a Matrix". I suggest just googling the phrase "Ninteen dubious ways".
> Worked for me, including finding an update paper from 2003.
As well as many hits, a web search also produced
Gerry Ford
"Maayan" <maaya...@weizmann.ac.il> wrote in message
news:9b9e6cc1-21bc-41fe...@c4g2000hsg.googlegroups.com...
represent this on use net is
[[ 1, 2]
[ 3, 4]] . If what you want to do is raise this to a power, the intrinsic
matmul would suffice.
--
Gerry Ford
"The apple was really a peach."
-- Allison Dunn on the garden of eden
paul.rich...@gmail.com
> Quite a few ways exist. Pretty much all of them have "issues", depending
> on the details of the needed usage. *THE* reference work on the question
> is Moler's 1978 paper "Ninteen Dubious Ways to Compute the Exponential
> of a Matrix". I suggest just googling the phrase "Ninteen dubious ways".
> Worked for me, including finding an update paper from 2003.
Searching on "matlab expm" produces the Mathworks documentation as
first hit. This gives references
References
[1] Golub, G. H. and C. F. Van Loan, Matrix Computation, p. 384, Johns
Hopkins University Press, 1983.
[2] Moler, C. B. and C. F. Van Loan, "Nineteen Dubious Ways to Compute
the Exponential of a Matrix," SIAM Review 20, 1978, pp. 801-836.
[3] Higham, N. J., "The Scaling and Squaring Method for the Matrix
Exponential Revisited," SIAM J. Matrix Anal. Appl., 26(4) (2005), pp.
1179-1193.
Surprise, surprise, the reference that Richard provided appears :-)
but the algorithm used is referenced to [3]. Apparently it is a Pade
approximation with scaling and squaring.
The example expmdemo1.m might be a good place to start.
Cheers
Paul
widmar
Maayan wrote:
>
> I am trying to figure out how to compute the exponential of a matrix
> in FORTRAN - something similar to expm used in MATLAB, and MatrixExp
> in Mathematica. I've gone through LAPACK's online documentation but
> couldn't find an answer. Is there some way of getting this done? I'd
> appreciate your kind help with this.
see, "Nineteen dubious ways*... twentyfive years later" for a review of
algorithms. The easiest, and perhaps the least dubious, way is to
numerically integrate matrix de,
d/dt F(t) = A*F(t)
on [0,t] s.t. F(0) = I. The solution - exp(A*t) - is the matrix
exponential at t=1.
* argues that scaling/squaring (w/Pade approximation) method might be
the least dubious, however it can be as effective if combined with
numerical integration.
---
sdx - modeling, simulation.
http://www.sdynamix.com
Hans Mittelmann
Peter Nachtwey
"widmar" <j...@sdynamix.com> wrote in message
news:47AE439E...@sdynamix.com...
>
> Maayan wrote:
>>
>> I am trying to figure out how to compute the exponential of a matrix
>> in FORTRAN - something similar to expm used in MATLAB, and MatrixExp
>> in Mathematica. I've gone through LAPACK's online documentation but
>> couldn't find an answer. Is there some way of getting this done? I'd
>> appreciate your kind help with this.
>
> see, "Nineteen dubious ways*... twentyfive years later" for a review of
> algorithms. The easiest, and perhaps the least dubious, way is to
> numerically integrate matrix de,
I have read that document and implemented the most promising methods. I
worked out this
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20expm.pdf
it is a .pdf generated using Mathcad. I prototype code before coding for
real.
It isn't Fortran but it will help. I use the expm() function for doing
control simulations. In my applications the lower right hand value can be
much bigger than one. This is way causes the errors because each iteration
this number gets bigger and bigger. You can see that making the update time
T small helps reduce the size of the values of the matrix below one. This
backs up the intuitive knowledge that smaller sample periods are more
accurate as long as the round off error don't be significant.
>
> d/dt F(t) = A*F(t)
>
> on [0,t] s.t. F(0) = I. The solution - exp(A*t) - is the matrix
> exponential at t=1.
>
> * argues that scaling/squaring (w/Pade approximation) method might be
> the least dubious, however it can be as effective if combined with
> numerical integration.
Yes, and by trying many different methods I found that the scaling and
squaring part is the most important. I didn't see and significant
difference betwee Taylor's series or Pade approximant.
Peter Nachtwey
widmar
>
> > see, "Nineteen dubious ways*... twentyfive years later" for a review of
> > algorithms. The easiest, and perhaps the least dubious, way is to
> > numerically integrate matrix de,
> >
> > d/dt F(t) = A*F(t)
> >
> > on [0,t] s.t. F(0) = I. The solution - exp(A*t) - is the matrix
> > exponential at t=1.
> >
> > * argues that scaling/squaring (w/Pade approximation) method might be
> > the least dubious, however it can be as effective if combined with
> > numerical integration.
>
>
from expokit:
exp(t*A) in full using irreducible Pade, A general (dense)
apparently it does matrix exponential the same way as discussed by Moler
& Van Loan and recently by Higham (referenced earlier) with (dubious)
claim of optimal efficiency. The enduring fixation on Pade' seems odd as
there must be simpler and better ways (e.g. minimax fit) to improve
efficiency. In any case, matrix exponential implementation by Van Loan -
pade.f - is also available from netlib.
Arjan
I have used the FORTRAN-version of Expokit to estimate the
simultaneous decay
of 300 types of decaying atoms, where daughter-nuclides arise from
their mothers,
with good results. The code halts when you offer a zero-matrix, but
the author
of expokit has explained that this is a "feature". Modification of the
code was easy.
Arjan
Gerry Ford
"Arjan" <arjan.v...@rivm.nl> wrote in message
news:6997dd3e-520c-4ebc...@h25g2000hsf.googlegroups.com...
Do you mean nucleides?
--
Gerry Ford
"Er hat sich georgiert." Der Spiegel, 2008, sich auf Chimpy Eins komma null
beziehend.
Arjan
Yes, I do (=not a marriage proposal), radioactive material.
Arjan
John Harper
Arjan <arjan.v...@rivm.nl> wrote:
>> Do you mean nucleides?
>
>Yes, I do (=not a marriage proposal), radioactive material.
Eh? For those writing in English, OED says "nuclide" is the word for a
distinct species of atom or nucleus. It has no entry for "nucleide".
-- John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington 6140, New Zealand
e-mail john....@vuw.ac.nz phone (+64)(4)463 6780 fax (+64)(4)463 5045
Jan Gerrit Kootstra
It is spelled by a Dutchman. Dutch sometimes uses English like
spellings, but for this word it is a French spelling without accents.
You and Arjan are refferring to the same meaning.
Kind regards,
ir. Jan Gerrit Kootstra
Alumnus of the Dept of Mathematics, University of Groningen.
FX
> distinct species of atom or nucleus. It has no entry for "nucleide".
The IUPAC (International Union of Pure and Applied Chemistry) says in the
"Nomeclature for Isotope, Nuclear and Radioanalytical Techniques":
NUCLEIDE -- "Nuclide" is preferred.
--
FX
glen herrmannsfeldt
> The IUPAC (International Union of Pure and Applied Chemistry) says in the
> "Nomeclature for Isotope, Nuclear and Radioanalytical Techniques":
>
> NUCLEIDE -- "Nuclide" is preferred.
In either case, keep George Bush away from them
(along with the nucular weapons).
-- glen
Les
"glen herrmannsfeldt" <g...@ugcs.caltech.edu> wrote in message
news:1dednSlV5rlO_k3a...@comcast.com...
and the threat from tourism :-)
Les
widmar
>
> I have used the FORTRAN-version of Expokit to estimate the
> with good results. The code halts when you offer a zero-matrix, but
> the author of expokit has explained that this is a "feature".
Would you accept this "feature" if compilers would do the same with
their exponential function? Not likely, and neither would the language
standard. Unfortunately, it appears to be a product of convenience, as
author didn't seem to bother with returning a "pointer" to a resulting
unit matrix.
btw, Van Loan's - pade.f - doesn't take the shortcut if A=0, although it
adds a (safety?) scale bias which appears carbon copied in expokit,
iscale = int(2 + log(anorm)/log(2.))
and Higham ("matrix exponential revisited") now boasts of shaving off at
least *one* matrix multiply - what gives?
Arjan
> > the author of expokit has explained that this is a "feature".
>
> Would you accept this "feature" if compilers would do the same with
> their exponential function? Not likely, and neither would the language
> standard.
Of course not. But after my modification of the code I don't care.
Arjan
robin
news:1ibxju6.3wt1m7ga0zmmN%nos...@see.signature...
> Quite a few ways exist. Pretty much all of them have "issues", depending
> on the details of the needed usage. *THE* reference work on the question
> is Moler's 1978 paper "Ninteen Dubious Ways to Compute the Exponential
> of a Matrix". I suggest just googling the phrase "Ninteen dubious ways".
> Worked for me, including finding an update paper from 2003.
Close.
In point of fact, *THE* reference work is the 2003 survey paper
by Moler & Loan, "Nineteen Dubious Ways to Compute the
Exponential of a Matrix, Twenty-five Years Later", which includes
a revision of the 1978 paper.
Everling's algorithm in PL/I (ref. 50) is based on Liou's paper.
robin
> Everling's algorithm in PL/I (ref. 50) is based on Liou's paper.
I have transcribed this and edited it.
robin
news:6997dd3e-520c-4ebc...@h25g2000hsf.googlegroups.com...
>
Beware using this hotchpotch of Fortran 90 and F77.
The codes use the CMPLX function on double precision values,
and do not specify the kind.
This converts to a single-precision result.
As well, there are other dubious practices of using DBLE(A*B)
where apparently DPROD should have been used, or DBLE(A)*DBLE(B).
robin
Beware using this hotchpotch of Fortran 90 and F77.
The codes use the CMPLX function on double precision values,
and do not specify the kind.
This converts to a single precision result.