f01 Chapter Contents
f01 Chapter Introduction
NAG Library Manual

# NAG Library Function Documentnag_matop_real_gen_matrix_actexp_rcomm (f01gbc)

## 1  Purpose

nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) computes the action of the matrix exponential ${e}^{tA}$, on the matrix $B$, where $A$ is a real $n$ by $n$ matrix, $B$ is a real $n$ by $m$ matrix and $t$ is a real scalar. It uses reverse communication for evaluating matrix products, so that the matrix $A$ is not accessed explicitly.

## 2  Specification

 #include #include
 void nag_matop_real_gen_matrix_actexp_rcomm (Integer *irevcm, Integer n, Integer m, double b[], Integer pdb, double t, double tr, double b2[], Integer pdb2, double x[], Integer pdx, double y[], Integer pdy, double p[], double r[], double z[], double comm[], Integer icomm[], NagError *fail)

## 3  Description

${e}^{tA}B$ is computed using the algorithm described in Al–Mohy and Higham (2011) which uses a truncated Taylor series to compute the ${e}^{tA}B$ without explicitly forming ${e}^{tA}$.
The algorithm does not explicity need to access the elements of $A$; it only requires the result of matrix multiplications of the form $AX$ or ${A}^{\mathrm{T}}Y$. A reverse communication interface is used, in which control is returned to the calling program whenever a matrix product is required.

## 4  References

Al–Mohy A H and Higham N J (2011) Computing the action of the matrix exponential, with an application to exponential integrators SIAM J. Sci. Statist. Comput. 33(2) 488-511
Higham N J (2008) Functions of Matrices: Theory and Computation SIAM, Philadelphia, PA, USA

## 5  Arguments

Note: this function uses reverse communication. Its use involves an initial entry, intermediate exits and re-entries, and a final exit, as indicated by the argument irevcm. Between intermediate exits and re-entries, all arguments other than b2, x, y, p and r must remain unchanged.
1:     irevcmInteger *Input/Output
On initial entry: must be set to $0$.
On intermediate exit: ${\mathbf{irevcm}}=1$, $2$, $3$, $4$ or $5$. The calling program must:
 (a) if ${\mathbf{irevcm}}=1$: evaluate ${B}_{2}=AB$, where ${B}_{2}$ is an $n$ by $m$ matrix, and store the result in b2; if ${\mathbf{irevcm}}=2$: evaluate $Y=AX$, where $X$ and $Y$ are $n$ by $2$ matrices, and store the result in y; if ${\mathbf{irevcm}}=3$: evaluate $X={A}^{\mathrm{T}}Y$ and store the result in x; if ${\mathbf{irevcm}}=4$: evaluate $p=Az$ and store the result in p; if ${\mathbf{irevcm}}=5$: evaluate $r={A}^{\mathrm{T}}z$ and store the result in r. (b) call nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) again with all other parameters unchanged.
On final exit: ${\mathbf{irevcm}}=0$.
2:     nIntegerInput
On entry: $n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
3:     mIntegerInput
On entry: the number of columns of the matrix $B$.
Constraint: ${\mathbf{m}}\ge 0$.
4:     b[$\mathit{dim}$]doubleInput/Output
Note: the dimension, dim, of the array b must be at least ${\mathbf{pdb}}×{\mathbf{m}}$.
The $\left(i,j\right)$th element of the matrix $B$ is stored in ${\mathbf{b}}\left[\left(j-1\right)×{\mathbf{pdb}}+i-1\right]$.
On initial entry: the $n$ by $m$ matrix $B$.
On intermediate exit: if ${\mathbf{irevcm}}=1$, contains the $n$ by $m$ matrix $B$.
On intermediate re-entry: must not be changed.
On final exit: the $n$ by $m$ matrix ${e}^{tA}B$.
5:     pdbIntegerInput
On entry: the stride separating matrix row elements in the array b.
Constraint: ${\mathbf{pdb}}\ge {\mathbf{n}}$.
6:     tdoubleInput
On entry: the scalar $t$.
7:     trdoubleInput
On entry: the trace of $A$. If this is not available then any number can be supplied ($0$ is a reasonable default); however, in the trivial case, $n=1$, the result ${e}^{{\mathbf{tr}}t}B$ is immediately returned in the first row of $B$. See Section 9.
8:     b2[$\mathit{dim}$]doubleInput/Output
Note: the dimension, dim, of the array b2 must be at least ${\mathbf{pdb2}}×{\mathbf{m}}$.
The $\left(i,j\right)$th element of the matrix is stored in ${\mathbf{b2}}\left[\left(j-1\right)×{\mathbf{pdb2}}+i-1\right]$.
On initial entry: need not be set.
On intermediate re-entry: if ${\mathbf{irevcm}}=1$, must contain $AB$.
On final exit: the array is undefined.
9:     pdb2IntegerInput
On entry: the stride separating matrix row elements in the array b2.
Constraint: ${\mathbf{pdb2}}\ge {\mathbf{n}}$.
10:   x[$\mathit{dim}$]doubleInput/Output
Note: the dimension, dim, of the array x must be at least ${\mathbf{pdx}}×2$.
The $\left(i,j\right)$th element of the matrix $X$ is stored in ${\mathbf{x}}\left[\left(j-1\right)×{\mathbf{pdx}}+i-1\right]$.
On initial entry: need not be set.
On intermediate exit: if ${\mathbf{irevcm}}=2$, contains the current $n$ by $2$ matrix $X$.
On intermediate re-entry: if ${\mathbf{irevcm}}=3$, must contain ${A}^{\mathrm{T}}Y$.
On final exit: the array is undefined.
11:   pdxIntegerInput
On entry: the stride separating matrix row elements in the array x.
Constraint: ${\mathbf{pdx}}\ge {\mathbf{n}}$.
12:   y[$\mathit{dim}$]doubleInput/Output
Note: the dimension, dim, of the array y must be at least ${\mathbf{pdy}}×2$.
The $\left(i,j\right)$th element of the matrix $Y$ is stored in ${\mathbf{y}}\left[\left(j-1\right)×{\mathbf{pdy}}+i-1\right]$.
On initial entry: need not be set.
On intermediate exit: if ${\mathbf{irevcm}}=3$, contains the current $n$ by $2$ matrix $Y$.
On intermediate re-entry: if ${\mathbf{irevcm}}=2$, must contain $AX$.
On final exit: the array is undefined.
13:   pdyIntegerInput
On entry: the stride separating matrix row elements in the array y.
Constraint: ${\mathbf{pdy}}\ge {\mathbf{n}}$.
14:   p[n]doubleInput/Output
On initial entry: need not be set.
On intermediate re-entry: if ${\mathbf{irevcm}}=4$, must contain $Az$.
On final exit: the array is undefined.
15:   r[n]doubleInput/Output
On initial entry: need not be set.
On intermediate re-entry: if ${\mathbf{irevcm}}=5$, must contain ${A}^{\mathrm{T}}z$.
On final exit: the array is undefined.
16:   z[n]doubleInput/Output
On initial entry: need not be set.
On intermediate exit: if ${\mathbf{irevcm}}=4$ or $5$, contains the vector $z$.
On intermediate re-entry: must not be changed.
On final exit: the array is undefined.
17:   comm[${\mathbf{n}}×{\mathbf{m}}+3×{\mathbf{n}}+12$]doubleCommunication Array
18:   icomm[$2×{\mathbf{n}}+40$]IntegerCommunication Array
19:   failNagError *Input/Output
The NAG error argument (see Section 3.6 in the Essential Introduction).

## 6  Error Indicators and Warnings

On entry, argument $⟨\mathit{\text{value}}⟩$ had an illegal value.
NE_INT
On entry, ${\mathbf{m}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{m}}\ge 0$.
On entry, ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{n}}\ge 0$.
On initial entry, ${\mathbf{irevcm}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{irevcm}}=0$.
On intermediate re-entry, ${\mathbf{irevcm}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{irevcm}}=1$, $2$, $3$, $4$ or $5$.
NE_INT_2
On entry, ${\mathbf{pdb}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdb}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{pdb2}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdb2}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{pdx}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdx}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{pdy}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdy}}\ge {\mathbf{n}}$.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
NW_SOME_PRECISION_LOSS
${e}^{tA}B$ has been computed using an IEEE double precision Taylor series, although the arithmetic precision is higher than IEEE double precision.

## 7  Accuracy

For a symmetric matrix $A$ (for which ${A}^{\mathrm{T}}=A$) the computed matrix ${e}^{tA}B$ is guaranteed to be close to the exact matrix, that is, the method is forward stable. No such guarantee can be given for non-symmetric matrices. See Section 4 of Al–Mohy and Higham (2011) for details and further discussion.

## 8  Parallelism and Performance

nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
Please consult the Users' Note for your implementation for any additional implementation-specific information.

## 9  Further Comments

### 9.1  Use of $Tr\left(A\right)$

The elements of $A$ are not explicitly required by nag_matop_real_gen_matrix_actexp_rcomm (f01gbc). However, the trace of $A$ is used in the preprocessing phase of the algorithm. If $Tr\left(A\right)$ is not available to the calling function then any number can be supplied ($0$ is recommended). This will not affect the stability of the algorithm, but it may reduce its efficiency.

### 9.2  When to use nag_matop_real_gen_matrix_actexp_rcomm (f01gbc)

nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) is designed to be used when $A$ is large and sparse. Whenever a matrix multiplication is required, the function will return control to the calling program so that the multiplication can be done in the most efficient way possible. Note that ${e}^{tA}B$ will not, in general, be sparse even if $A$ is sparse.
If $A$ is small and dense then nag_matop_real_gen_matrix_actexp (f01gac) can be used to compute ${e}^{tA}B$ without the use of a reverse communication interface.
The complex analog of nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) is nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc).

### 9.3  Use in Conjunction with NAG C Library Functions

To compute ${e}^{tA}B$, the following skeleton code can normally be used:
```do {
f01gbc(&irevcm,n,m,b,tdb,t,tr,b2,tdb2,x,tdx,y,tdy,p,r,z,comm,icomm,&fail);
if (irevcm == 1) {
.. Code to compute B2=AB ..
}
else if (irevcm == 2){
.. Code to compute Y=AX ..
}
else if (irevcm == 3){
.. Code to compute X=A^T Y ..
}
else if (irevcm == 4){
.. Code to compute P=AZ ..
}
else if (irevcm == 5){
.. Code to compute R=A^T Z ..
}
} (while irevcm !=0)
```
The code used to compute the matrix products will vary depending on the way $A$ is stored. If all the elements of $A$ are stored explicitly, then nag_dgemm (f16yac)) can be used. If $A$ is triangular then nag_dtrmm (f16yfc) should be used. If $A$ is symmetric, then nag_dsymm (f16ycc) should be used. For sparse $A$ stored in coordinate storage format nag_sparse_nsym_matvec (f11xac) and nag_sparse_sym_matvec (f11xec) can be used. Alternatively if $A$ is stored in compressed column format nag_superlu_matrix_product (f11mkc) can be used.

## 10  Example

This example computes ${e}^{tA}B$, where
 $A = 0.4 -0.2 1.3 0.6 0.3 0.8 1.0 1.0 3.0 4.8 0.2 0.7 0.5 0.0 -5.0 0.7 ,$
 $B = 0.1 1.1 1.7 -0.2 0.5 1.0 0.4 -0.2 ,$
and
 $t=-0.2 .$

### 10.1  Program Text

Program Text (f01gbce.c)

### 10.2  Program Data

Program Data (f01gbce.d)

### 10.3  Program Results

Program Results (f01gbce.r)