|
Home | Download | Support
|
| Support |
|
Expokit is an industry-grade software designed to be very easy
to use and/or to interface with other appplications. In fact, criteria
such as user-friendliness, performance, portability, stability
and reliability were the driving focus when designing the software. Installation
instructions are provided in the package.
There have been interested users who needed to embed the
computation of
the matrix exponential in their applications written in other languages
such as C/C++,
Java, Mathematica. Some of these users have ported Expokit directly to
their native programming language of interest, while other users have
preserved the original package in Fortran, and cross-linked across
language boundaries. However, due to limited resources, we do not
distribute or support any of these derivatives. Expokit is now embedded
in several applications and packages (some commercial).
Usage / FAQs / Templates
There are certain typical usage patterns that
arise frequently in practice. We give below some hints
to assist you if you face one of these frequent cases. Join the discussion forum if
you have other
queries.
|
Scenario 1:
Compute the dense matrix
exponential applied to a set of vectors
In Matlab:
Simply use the built-in function expm.
Expokit provides another function if your matrix is sparse (see
scenario 3 below).
In Fortran:
The answer here actually involves just two executive statements
with the help of Expokit. The declaration of variables and workspace
adds
to the clutter. Specifically, to compute C = exp(t*A)*B, you can write
a
wrapper, say ZEXPM (or DEXPM if your data is real), that can be called
as:
call ZEXPM( n, m, t, A,LDA, B,LDB, C,LDC )
to compute:
C(1:n, 1:m) = exp( t * A(1:n, 1:n) ) * B(1:n, 1:m)
where
on input:
t is a real scalar (positive or negative)
A(LDA,*) is complex n-by-n
B(LDB,*) is complex n-by-m
on output:
C(LDC,*) is complex n-by-m
notes:
LDA, LDB, LDC are to be provided on input and represent
the so-called leading-dimension of A, B, C, respectively.
Using the compact syntax of f90, a pseudo-code for ZEXPM will read:
subroutine ZEXPM( n, m, t, A,LDA, B,LDB, C,LDC )
implicit none
double precision t
integer n, m, LDA, LDB, LDC
double complex A(LDA,n), B(LDA,m), C(LDC,m)
integer ideg, iflag, iexp, ns, lwsp
double complex ZERO, ONE
parameter( ideg=6, lwsp=4*n*n+ideg+1 )
parameter( ZERO=(0.0d0,0.0d0), ONE=(1.0d0,0.0d0) )
double complex, allocatable :: wsp(:)
integer, allocatable :: ipiv(:)
ALLOCATE( wsp(lwsp), ipiv(n) )
*--- compute E = exp(tA) via Expokit's Pade algorithm
call ZGPADM(ideg, n, t, A,LDA, wsp,lwsp, ipiv, iexp, ns, iflag)
*--- multiply C = E*B via BLAS
call ZGEMM('n','n', n,m,n, ONE, wsp(iexp),n, B,LDB, ZERO, C,LDC )
DEALLOCATE( ipiv, wsp )
END
Scenario 2:
Dump some entries of exp(tkA)v at a number of points tk,
e.g., for plotting purposes
This should normally come after successfully trying out your matrix (as
explained in scenario 3), but it is shorter to describe and is placed
here. When performing certain experiments, you may
wish to print some selected components of wk = exp(tkA)v
over an interval [0, t] for
0 = t0 < t1 < ... <
ts < ts+1 = t. There are different ways to
achieve this.
If your matrix isn't that large, just use DGPADM (or ZGPADM) as
described in scenario 1. Otherwise, you have a
large sparse matrix, and the easiest way to dump certain
quantities based on exp(tkA)v
is simply to locate and change the relevant
code section in the body of dgexpv.f (or zgexpv.f) during your
experiments. The red section below shows what you might add.
if ( itrace .eq. -1 ) then
*--- example: print tk and the first component of wk
print*, t_now, w(1)
endif
if ( itrace.gt.0 ) then
if ( itrace.ne.0 ) then
print*,'integration',nstep,'---------------------------------'
print*,'scale-square =',ns
print*,'step_size =',t_step
print*,'err_loc =',err_loc
print*,'next_step =',t_new
endif
With this change, you only have to call DGEXPV with its
parameter itrace set to your value (-1 in this case) for
the code to dump what you want during the step-by-step integration
process.
Scenario 3:
Compute w = exp(tA)v for my test matrix A
In Matlab:
If you are in Matlab, doing this computation is straightforward using
the provided Matlab routines of Expokit such as matlab/expv.
Just call w = expv(t, A, v).
'expv' is like our usual matrix-vector product, but does the
exponential-vector instead -- so to speak. It computes
w = expv(t, A, v) = exp(tA)*v directly.
The file matlab/expv.m
is entirely self-contained and
can be
used right after saving it (as plain text),
without having to hunt/install anything else. Once you get it, type help
expv for other details.
In Fortran:
The remainder of
this section is only concerned with Fortran. If your matrix isn't that
large, just use DGPADM (or ZGPADM) a described in scenario 1.
Otherwise, you have a large sparse matrix.
If you are able to dump your test matrix A in a file,
and are familiar with Fortran, it should not take you an overtly long
time to
try out your data. An hour or so should be enough. The following steps
summarize
what you should do:
- Setup your data: To test directly with a minimum of
effort, save your matrix in one of the formats described in the Expokit
documentation (e.g., COO - Coordinate format, or CRS - Compress Row
Storage format which is the one used
for the matrix data/c1024.crs). See
also the
included matlab/mat2coo.m, matlab/mat2crs.m which allow you
to quickly save in those formats from within Matlab. There are more
details about each format in the header of those files. The COO format
is the easiest one to use (albeit less efficient). For ease of
presentation, the sample code below assumes that you have saved your
matrix in the file named mycoo.dat. It should be
space (or tab) delimited as follows:
In the first line of the file, specify the order (n) and
the number of non-zeros (nz) entries in the matrix -- so that a
program can read just that, and allocate the sufficient amount of
memory needed to read the rest. In the subsequent lines of the file,
specify the row, column, non-zero value. That is, with the COO format,
your data file should look like:
n nz
i j a_ij
k l a_kl
etc...
Fortran (and Matlab) are not 0-based, so for seamless interfacing, dump
i as 1 for the first row and j as 1 for the first column.
If you have a real matrix, each of the a_ij, a_kl, etc, is real. If
your matrix is complex, you will have to format your data in the way
expected by Fortran. It requires that complex values be given as a pair
of real numbers enclosed in parentheses and separated by a comma.
Specifically, dump your
complex matrix as:
n nz
i j (real_part_of_a_ij , imaginary_part_of_a_ij)
k l (real_part_of_a_kl , imaginary_part_of_a_lk)
etc...
Note:
A number of routines in Expokit take an external function matvec
as input, which is to be used for matrix-vector multiplication. However
the default signature has no 'data' argument for passing state to this
function. The sample programs use static variables (passed via a
Fortran common statement) for such state, which is not always ideal.
Borrowing a bit from C/C++, a better signature would rather
be:
subroutine DGEXPV( n, m, t, v, w, tol, anorm, void* callbackData,
. wsp,lwsp, iwsp,liwsp, matvec, itrace,iflag )
with the corresponding
matvec( x, y, void* callbackData ).
Thus allowing to pass an arbitrary structure to get at the data
needed to characterize the operation.
However, once you have a positive first impression, it is very
simple to modify the code, and the preference went to simplicity at
first sight. This author used the following signatures in some of the
internal testings:
subroutine DGEXPV( n, m, t, v, w, tol, anorm, ia, ja, a,
. wsp,lwsp, iwsp,liwsp, matvec, itrace,iflag )
matvec( x, y,
integer n,
integer nz,
double precision a(*),
integer ia(*),
integer ja(*) )
with a given storage format (e.g., COO) throughout the runs.
Once you are familiar with the package and ready to make a
custom build, you might want to redefine matvec (and DGEXPV) with as
many parameters as you deem fit, depending on the setup and
storage format of your particular matrix. But for your
introductory test drive given here, simply setup your
data file as explained above, and let's continue.
- Hook to Expokit:
Create a custom driver to read your data and to call a matrix
exponential routine from Expokit. Use these ready made templates that
you can save as myexample.f in the expokit/fortran directory:
These templates come from simplifying one the default drivers included
in the Expokit bundle.
- Build:
Edit and add the following make rule in the expokit/fortran/Makefile
along side the other examples that are there:
myexample: $(OBJS) myexample.o
$(FC) -o myexample myexample.o $(OBJS) $(LIBS)
^^^^
this space is a tab! (to keep the makefile happy...)
Then, compile by typing make myexample.
- Execute: Simply type myexample > output.txt
to capture
the results in the output.txt file, and that's all. You should now
be well equipped to use Expokit in other contexts.