EXPOKIT
Software for Computing the Matrix Exponential.

EXPOKIT is a FORTRAN77 library which handles various forms of the matrix exponential problem, by Roger Sidje.

Formally, for a square matrix A and scalar t, the matrix exponential exp(A*t) can be defined as the sum:

exp(A*t) = sum ( 0 <= i < oo ) A^i t^i / i!

The simplest form of the matrix exponential problem asks for the value when t = 1. Even for this simple case, and for a matrix of small order, it can be quite difficult to compute the matrix exponential accurately.

EXPOKIT provides a set of routines aimed at the accurate computation of the matrix exponential itself, or of computing the product of the matrix exponential with a given vector, or the solution of a system of linear ODE's of the form w' = A*w(t) + B*u(t).

The software can handle real or complex matrices, and has special interfaces for symmetric or Hermitian matrices.

The version of EXPOKIT exhibited here has been modified from the original. Persons interested in using EXPOKIT are urged to consult the original source, available at "http://www.maths.uq.edu.au/expokit"

Licensing:

Permission to use, copy, modify, and distribute EXPOKIT and its supporting documentation for non-commercial purposes, is hereby granted without fee, provided that this permission message and copyright notice appear in all copies. Approval must be sought for commercial purposes.

Neither the Institution (University of Queensland) nor the Authors make any representations about the suitability of this software for any purpose. This software is provided ``as is'' without express or implied warranty.

Languages:

EXPOKIT is available in a FORTRAN77 version.

Related Data and Programs:

ARPACK, a FORTRAN90 library which computes eigenvalues for large matrices;

DLAP, a data directory which contains examples of DLAP files, a sparse matrix file format used by DLAP;

DSP, a data directory which contains examples of DSP files, a sparse matrix file format, storing just (I,J,A(I,J)), and using one-based indexing.

HB, a data directory which contains examples of HB files, the Harwell Boeing sparse matrix file format;

HBSMC, a dataset directory which contains the Harwell Boeing Sparse Matrix Collection;

MM, a data directory which contains examples of MM files, the Matrix Market Sparse Matrix File Format;

SPARSE_CC, a data directory which contains examples of SPARSE CC files, a three-file format for the sparse compressed column matrix format;

SPARSE_CR, a data directory which contains examples of SPARSE CR files, a three-file format for the sparse compressed row matrix format;

ST, a data directory which contains examples of ST files, the Sparse Triplet format, a sparse matrix file format, storing just (I,J,A(I,J)), and using zero-based indexing.

TEST_MAT, a FORTRAN77 library which defines test matrices.

Author:

Roger Sidje,
Department of Mathematics,
University of Queensland,
Brisbane, Australia.

Reference:

1. Cleve Moler, Charles VanLoan,
   Nineteen Dubious Ways to Compute the Exponential of a Matrix, SIAM Review,
   Volume 20, Number 4, October 1978, pages 801-836.
2. Cleve Moler, Charles VanLoan,
   Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later,
   SIAM Review,
   Volume 45, Number 1, March 2003, pages 3-49.
3. Roger Sidje,
   EXPOKIT: Software Package for Computing Matrix Exponentials,
   ACM Transactions on Mathematical Software,
   Volume 24, Number 1, March 1998, pages 130-156.

Source Code:

• expokit.f, the source code.
• expokit.sh, BASH commands to compile the source code.

Examples and Tests:

• expokit_prb.f, a sample calling program.
• expokit_prb.sh, BASH commands to compile and run the sample program.
• expokit_prb_output.txt, the output file.

List of Routines:

• CLOCK returns a reading from the real time clock.
• DASUM sums the absolute value of the entries of a vector.
• DAXPY adds a multiple of one vector to another.
• DCABS1 computes the L1 norm of a complex number.
• DCMPAC compacts the array IX and sorts IXX and XX.
• DCOPY copies a vector.
• DDOT computes the dot product of two vectors.
• DGCCSV computes y = A * x for A in compressed column storage (CCS).
• DGCHBV computes y=exp(t*H)*y, general matrix H, using rational Chebyshev.
• DGCNVR converts between sparse matrix storage formats.
• DGCOOV computes y = A * x, with A stored in "COO" format.
• DGCRSV computes y = A * x, for A stored in Compressed Row Storage (CRS) format.
• DGEFA factors a double precision matrix by gaussian elimination.
• DGEMM performs matrix-matrix operations of the form C=alpha*A*B+beta*C.
• DGEMV performs a matrix-vector operation of the form y=alpha*A*x+beta*y.
• DGESL solves a factored linear system of the form A*x=b.
• DGESV solves multiple general linear systems of the form A*X=B.
• DGEXPV computes w=exp(t*A)*v for a general matrix A.
• DGPADM computes exp(t*H) using an irreducible rational Pade approximation.
• DGPHIV computes the solution of the ODE w'=A*w+u, w(0)=v.
• DLARAN returns a random real number uniformly selected from (0,1).
• DMEXPV computes w = exp(t*A)*v, customized for Markov chains.
• DNCHBV computes y = exp(t*H)*y using a partial fraction expansion.
• DNRM2 computes the Euclidean norm of a vector.
• DSCAL scales a vector by a constant.
• DSCHBV computes y = exp(t*H)*y using a partial fraction expansion.
• DSEXPV computes w = exp(t*A)*v, for a symmetric matrix A.
• DSIFA factors a symmetric matrix.
• DSISL solves a factored symmetric linear system.
• DSPADM computes exp(t*H), the matrix exponential of a symmetric matrix.
• DSPHIV solves the ODE w' = Aw + u, w(0) = v, for A symmetric.
• DSWAP interchanges two vectors.
• DSYSV solves multiple symmetric linear systems of the form A*X=B.
• DZNRM2 computes the unitary norm of a complex vector.
• GETPAT loads a Harwell-Boeing pattern.
• IDAMAX indexes the vector entry of maximum absolute value.
• IDSRT1 is an indirect sort; sorts IX and carries XX along.
• IDSRT2 is an indirect sort: sorts IX and carries IXX and XX along.
• IZAMAX indexes the complex vector entry of maximum magnitude.
• IZSRT1 is an indirect sort; sorts IX and carries XX along.
• IZSRT2 is an indirect sort: sorts IX and carries IXX and XX along.
• LOADHB loads a Harwell-Boeing sparse matrix from a file.
• LSAME compares two characters for identity.
• TIMESTAMP prints out the current YMDHMS date as a timestamp.
• TNSPOS transposes a CRS matrix.
• XERBLA is an error handler for the level 2 BLAS routines.
• ZAXPY adds a multiple of one vector to another.
• ZCMPAC compacts the array IX and sorts IXX and XX.
• ZCOPY copies a vector.
• ZDOTC computes the conjugated dot product of complex vectors.
• ZDOTU computes the unconjugated dot product of complex vectors.
• ZDSCAL scales a complex vector by a real constant.
• ZGCCSV computes y = A * x, for A in compressed column storage (CCS).
• ZGCHBV computes y = exp(t*H)*y using a partial fraction expansion.
• ZGCNVR transforms a sparse storage format into another sparse storage format.
• ZGCOOV computes y = A * x, with A stored in COO storage format.
• ZGCRSV computes y = A * x, with A stored in compressed row storage (CRS).
• ZGEFA factors a complex general linear system.c
• ZGEMM performs a matrix-matrix operation of the form C=alpha*A*B+beta*C.
• ZGEMV performs a matrix-vector operation of the form y=alpha*A*x+beta*y.
• ZGESL solves a factored complex general linear system.
• ZGESV solves multiple general linear systems of the form A*X=B.
• ZGEXPV computes w = exp(t*A)*v.
• ZGPADM computes exp(t*H), the matrix exponential of a general complex matrix.
• ZGPHIV solves the ODE w' = Aw + u, w(0) = v.
• ZHESV solves multiple hermitian linear systems of the form A*X=B.
• ZHEXPV computes w = exp(t*A)*v for a Hermitian matrix A.
• ZHIFA factors a complex Hermitian linear system.
• ZHISL solves a factored complex Hermitian linear system.
• ZHPADM computes exp(t*H), the matrix exponential of a Hermitian matrix.
• ZHPHIV solves the linear ODE w' = Aw + u, w(0) = v for A hermitian.
• ZNCHBV computes y = exp(t*H)*y using a partial fraction expansion.
• ZSCAL scales a vector by a constant.
• ZSIFA factors a complex symmetric matrix.
• ZSISL solves a factored complex symmetric linear system.
• ZSWAP interchanges two vectors.
• ZSYSV solves multiple symmetric linear systems of the form A*X=B.

You can go up one level to the FORTRAN77 source codes.