Other questions? Check the
|Home | Download | 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
Using the compact syntax of f90, a pseudo-code for ZEXPM will read:
call ZEXPM( n, m, t, A,LDA, B,LDB, C,LDC )
C(1:n, 1:m) = exp( t * A(1:n, 1:n) ) * B(1:n, 1:m)
t is a real scalar (positive or negative)
A(LDA,*) is complex n-by-n
B(LDB,*) is complex n-by-m
C(LDC,*) is complex n-by-m
LDA, LDB, LDC are to be provided on input and represent
the so-called leading-dimension of A, B, C, respectively.
subroutine ZEXPM( n, m, t, A,LDA, B,LDB, C,LDC )
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 )
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.
if ( itrace .eq. -1 ) then
*--- example: print tk and the first component of wk
print*, t_now, w(1)
if ( itrace.gt.0 ) then
if ( itrace.ne.0 ) then
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:
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:
i j a_ij
k l a_kl
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.
i j (real_part_of_a_ij , imaginary_part_of_a_ij)
k l (real_part_of_a_kl , imaginary_part_of_a_lk)
Borrowing a bit from C/C++, a better signature would rather be:
with the corresponding
subroutine DGEXPV( n, m, t, v, w, tol, anorm, void* callbackData,
. wsp,lwsp, iwsp,liwsp, matvec, itrace,iflag )
Thus allowing to pass an arbitrary structure to get at the data needed to characterize the operation.
matvec( x, y, void* callbackData ).
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:
with a given storage format (e.g., COO) throughout the runs.
subroutine DGEXPV( n, m, t, v, w, tol, anorm, ia, ja, a,
. wsp,lwsp, iwsp,liwsp, matvec, itrace,iflag )
matvec( x, y,
double precision a(*),
integer ja(*) )
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.
Then, compile by typing
myexample: $(OBJS) myexample.o
$(FC) -o myexample myexample.o $(OBJS) $(LIBS)
this space is a tab! (to keep the makefile happy...)