
Home  Download  Support 
Support  
Expokit is an industrygrade software designed to be very easy
to use and/or to interface with other appplications. In fact, criteria
such as userfriendliness, 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 crosslinked 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 pseudocode for ZEXPM will read:
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 nbyn
B(LDB,*) is complex nbym
on output:
C(LDC,*) is complex nbym
notes:
LDA, LDB, LDC are to be provided on input and represent
the socalled leadingdimension of A, B, C, respectively.
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
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 stepbystep integration process.
if ( itrace .eq. 1 ) then
* example: print t_{k} and the first component of w_{k}
print*, t_now, w(1)
endif
if ( itrace.gt.0 ) then
if ( itrace.ne.0 ) then
print*,'integration',nstep,''
print*,'scalesquare =',ns
print*,'step_size =',t_step
print*,'err_loc =',err_loc
print*,'next_step =',t_new
endif
In the first line of the file, specify the order (n) and the number of nonzeros (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, nonzero value. That is, with the COO format, your data file should look like:
Fortran (and Matlab) are not 0based, 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 a_ij
k l a_kl
etc...
Note: A number of routines in Expokit take an external function matvec as input, which is to be used for matrixvector 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.
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...
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,
integer n,
integer nz,
double precision a(*),
integer ia(*),
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...)