expoRkit (version 0.9.4)

Rexpv: Expokit __EXPV and __PHIV wrapper.

Description

R wrapper of the Expokit Fortran subroutines __EXPV and __PHIV for sparse matrix exponentiation. In general, these routines compute the solution at time point \(t\) of the ODE $$w'(t) = x w(t) + u$$ with initial condition \(w(0) = v\).

Usage

Rexpv(a, ia, ja, n, v, t = 1, storage = "CCS", u = NULL,
  anorm = max(abs(a)), Markov = FALSE, m = 30L, tol = 0,
  itrace = 0L, mxstep = 10000L)

Arguments

a

numeric or complex non-zero entries in the \(x\)-matrix.

ia

integer index/pointer. Precise meaning depends on storage format.

ja

integer index/pointer. Precise meaning depends on storage format.

n

dimension of the (square) matrix.

v

numeric or complex vector.

t

time. Default 1.

storage

character, one of 'CCS' (Compressed Column Storage), 'CRS' (Compressed Row Storage) or 'COO' (COOrdinate list). Default 'CCS'.

u

numeric or complex vector. Default NULL.

anorm

A norm of the matrix. Default is the sup-norm.

Markov

logical, if TRUE the (transposed) matrix is taken to be a rate matrix and steps are taken to ensure that the computed result is a probability vector. Default FALSE.

m

integer, the maximum size for the Krylov basis.

tol

numeric. A value of 0 (default) means square root of machine eps.

itrace

integer, 0 (default) means no trace information from Expokit, 1 means print 'happy breakdown', and 2 means all trace information printed from Expokit.

mxstep

integer. Maximum allowable number of integration steps. The value 0 means an infinite number of steps. Default 10000.

Value

The solution, \(w\), of the ODE as a numeric or complex vector of length \(n\).

Details

The Rexpv function is the low level wrapper of the Fortran subroutines in the Expokit package. It is not intended to be used directly but rather via the expv methods. In the call the correct storage format in terms of the vectors a, ia and ja has to be specified via the storage argument. For CCS, ia contains 1-based row numbers of non-zero elements and ja contains 1-based pointers to the initial element for each column. For CRS, ja contains 1-based column numbers of non-zero elements and ia are 1-based pointers to the initial element for each row. For COO, ia and ja contain the 1-based column and row numbers, respectively, for the non-zero elements.

References

Sidje, R. B. (1998) Expokit. Software Package for Computing Matrix Exponentials. ACM Trans. Math. Softw. 24(1), 130-156.

See Also

expm, expm, expv

Examples

# NOT RUN {
### A CCS 4 by 4 real matrix. The last element in 'ja' is the number of
### non-zero elements + 1. 
a <- c(-1, 1, -2, -3, 1, 2, -1)
ia <- c(1, 3, 2, 4, 1, 2, 3) 
ja <- c(1, 3, 5, 6, 8)

v <- c(1, 1, 1, 1)
wCCS <- expoRkit:::Rexpv(a, ia, ja, 4, v = v)

### COO storage instead.
ja <- c(1, 1, 2, 2, 3, 4, 4)  
wCOO <- expoRkit:::Rexpv(a, ia, ja, 4, v = v, storage = 'COO')

### CRS storage instead.
a <- c(-1, 1, -2, 2, 1, -1, -3)
ja <- c(1, 3, 2, 4, 1, 4, 2)
ia <- c(1, 3, 5, 7, 8)
wCRS <- expoRkit:::Rexpv(a, ia, ja, 4, v = v, storage = 'CRS')

cbind(wCCS, wCOO, wCRS)

stopifnot(all.equal(wCCS, wCOO),
          all.equal(wCCS, wCRS),
          all.equal(wCRS, wCOO))

## Complex version
a <- c(-1, 1i, -2, -3i, 1, 2i, -1)
ia <- c(1, 3, 2, 4, 1, 2, 3) 
ja <- c(1, 3, 5, 6, 8)

v <- c(1, 1, 1, 1)
wCCS <- expoRkit:::Rexpv(a, ia, ja, 4, v = v)

### COO storage instead.
ja <- c(1, 1, 2, 2, 3, 4, 4)  
wCOO <- expoRkit:::Rexpv(a, ia, ja, 4, v = v, storage = 'COO')

### CRS storage instead.
a <- c(-1, 1, -2, 2i, 1i, -1, -3i)
ja <- c(1, 3, 2, 4, 1, 4, 2)
ia <- c(1, 3, 5, 7, 8)
wCRS <- expoRkit:::Rexpv(a, ia, ja, 4, v = v, storage = 'CRS')

cbind(wCCS, wCOO, wCRS)

stopifnot(all.equal(wCCS, wCOO),
          all.equal(wCCS, wCRS),
          all.equal(wCRS, wCOO))

# }