Expokit


Other questions? Check the
discussion forum.
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:
  1. 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.

  2. 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.

  3. 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.

  4. 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.