Home > Archive > Fortran > March 2008 > Matrix exponentiation in FORTRAN
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
Matrix exponentiation in FORTRAN
|
|
| Maayan 2008-02-07, 8:12 am |
| Hi,
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 2008-02-07, 8:12 am |
|
Hi,
As if by magic, Maayan appeared !
> Hi,
>
> 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.
>
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 2008-02-07, 10:13 pm |
| Ian Bush <I.J.Bush@ku.ca.ld> wrote:
> As if by magic, Maayan appeared !
[color=darkred]
> 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
| |
|
| On Thu, 07 Feb 2008 07:43:36 -0800, Richard Maine wrote:
> Ian Bush <I.J.Bush@ku.ca.ld> wrote:
>
>
>
>
> 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
http://www.maths.uq.edu.au/expokit/
| |
| Gerry Ford 2008-02-08, 7:18 pm |
|
"Maayan" <maayan.gal@weizmann.ac.il> wrote in message
news:9b9e6cc1-21bc-41fe-a34d-4019119edbf8@c4g2000hsg.googlegroups.com...
> Hi,
>
> 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.
Can you give an example of what you want to do with a 2 x 2? One way to
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.richard.thomas@gmail.com 2008-02-09, 4:30 am |
|
> 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 2008-02-09, 10:18 pm |
|
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 2008-02-09, 10:18 pm |
| On Feb 9, 5:23 pm, widmar <j...@sdynamix.com> wrote:
> Maayan wrote:
>
>
> 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
http://www.maths.uq.edu.au/expokit/guide.html
| |
| Peter Nachtwey 2008-02-10, 7:18 pm |
|
"widmar" <jrw@sdynamix.com> wrote in message
news:47AE439E.9285580D@sdynamix.com...
>
> Maayan wrote:
>
> 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/N...%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 2008-02-27, 10:17 pm |
| Hans Mittelmann wrote:
>
>
>
> http://www.maths.uq.edu.au/expokit/guide.html
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.
---
sdx - modeling, simulation.
http://www.sdynamix.com
| |
|
| > from expokit:
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 2008-03-05, 10:12 pm |
|
"Arjan" <arjan.van.dijk@rivm.nl> wrote in message
news:6997dd3e-520c-4ebc-bbed-9e81c7e0a325@h25g2000hsf.googlegroups.com...
>
>
> 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
>
Do you mean nucleides?
--
Gerry Ford
"Er hat sich georgiert." Der Spiegel, 2008, sich auf Chimpy Eins komma null
beziehend.
| |
|
| > Do you mean nucleides?
Yes, I do (=not a marriage proposal), radioactive material.
Arjan
| |
| John Harper 2008-03-06, 7:14 pm |
| In article <3e7d7766-d27a-4bf3-b9a0-994f306eed3c@m36g2000hse.googlegroups.com>,
Arjan <arjan.van.dijk@rivm.nl> wrote:
>
>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.harper@vuw.ac.nz phone (+64)(4)463 6780 fax (+64)(4)463 5045
| |
|
| > 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".
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 2008-03-06, 7:14 pm |
| FX wrote:
> 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
| |
|
|
"glen herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
news:1dednSlV5rlO_k3anZ2dnUVZ_oPinZ2d@co
mcast.com...
> FX wrote:
>
>
> In either case, keep George Bush away from them
> (along with the nucular weapons).
and the threat from tourism :-)
Les
| |
| widmar 2008-03-10, 7:52 pm |
| Arjan wrote:
>
> I have used the FORTRAN-version of Expokit to estimate the
> simultaneous decay...
> 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?
---
sdx - modeling, simulation.
http://www.sdynamix.com
| |
|
| > > with good results. The code halts when you offer a zero-matrix, but
>
> 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
| |
|
| "Richard Maine" <nospam@see.signature> wrote in message
news:1ibxju6.3wt1m7ga0zmmN%nospam@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.
| |
|
| "none" <none@none.net> wrote in message news:pan.2008.02.07.19.33.04.969000@none.net...
> On Thu, 07 Feb 2008 07:43:36 -0800, Richard Maine wrote:
>
>
> As well as many hits, a web search also produced
>
> http://www.maths.uq.edu.au/expokit/
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).
| |
|
| "Arjan" <arjan.van.dijk@rivm.nl> wrote in message
news:6997dd3e-520c-4ebc-bbed-9e81c7e0a325@h25g2000hsf.googlegroups.com...
>
> 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
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).
|
|
|
|
|