f01 Chapter Contents
f01 Chapter Introduction
NAG Library Manual

# NAG Library Function Documentnag_matop_complex_gen_matrix_actexp_rcomm (f01hbc)

## 1  Purpose

nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) computes the action of the matrix exponential ${e}^{tA}$, on the matrix $B$, where $A$ is a complex $n$ by $n$ matrix, $B$ is a complex $n$ by $m$ matrix and $t$ is a complex 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_complex_gen_matrix_actexp_rcomm (Integer *irevcm, Integer n, Integer m, Complex b[], Integer pdb, Complex t, Complex tr, Complex b2[], Integer pdb2, Complex x[], Integer pdx, Complex y[], Integer pdy, Complex p[], Complex r[], Complex z[], Complex ccomm[], 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{H}}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{H}}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{H}}z$ and store the result in r. (b) call nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) 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}$]ComplexInput/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:     tComplexInput
On entry: the scalar $t$.
7:     trComplexInput
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}$]ComplexInput/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}$]ComplexInput/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{H}}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}$]ComplexInput/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]ComplexInput/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]ComplexInput/Output
On initial entry: need not be set.
On intermediate re-entry: if ${\mathbf{irevcm}}=5$, must contain ${A}^{\mathrm{H}}z$.
On final exit: the array is undefined.
16:   z[n]ComplexInput/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:   ccomm[${\mathbf{n}}×\left({\mathbf{m}}+2\right)$]ComplexCommunication Array
18:   comm[$3×{\mathbf{n}}+14$]doubleCommunication Array
19:   icomm[$2×{\mathbf{n}}+40$]IntegerCommunication Array
20:   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 an Hermitian matrix $A$ (for which ${A}^{\mathrm{H}}=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-Hermitian matrices. See Section 4 of Al–Mohy and Higham (2011) for details and further discussion.

## 8  Parallelism and Performance

Not applicable.

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

The elements of $A$ are not explicitly required by nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc). 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_complex_gen_matrix_actexp_rcomm (f01hbc)

nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) 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_complex_gen_matrix_actexp (f01hac) can be used to compute ${e}^{tA}B$ without the use of a reverse communication interface.
The real analog of nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) is nag_matop_real_gen_matrix_actexp_rcomm (f01gbc).

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

To compute ${e}^{tA}B$, the following skeleton code can normally be used:
```do {
f01hbc(&irevcm,n,m,b,tdb,t,tr,b2,tdb2,x,tdx,y,tdy,p,r,z,ccomm,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^H Y ..
}
else if (irevcm == 4){
.. Code to compute P=AZ ..
}
else if (irevcm == 5){
.. Code to compute R=A^H 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_zgemm (f16zac) can be used. If $A$ is triangular then nag_ztrmm (f16zfc) should be used. If $A$ is Hermitian, then nag_zhemm (f16zcc) should be used. If $A$ is symmetric, then nag_zsymm (f16ztc) should be used. For sparse $A$ stored in coordinate storage format nag_sparse_nherm_matvec (f11xnc) and nag_sparse_herm_matvec (f11xsc) can be used. For sparse $A$ stored in compressed column storage format (CCS) the program text of Section 10 contains the function matmul to perform matrix products.

## 10  Example

This example computes ${e}^{tA}B$ where
 $A = 0.7+0.8i -0.2+0.0i 1.0+0.0i 0.6+0.5i 0.3+0.7i 0.7+0.0i 0.9+3.0i 1.0+0.8i 0.3+3.0i -0.7+0.0i 0.2+0.6i 0.7+0.5i 0.0+0.9i 4.0+0.0i 0.0+0.0i 0.2+0.0i ,$
 $B = 0.1+0.0i 1.2+0.1i 1.3+0.9i -0.2+2.0i 4.0+0.6i -1.0+0.8i 0.4+0.0i -0.9+0.0i$
and
 $t=1.1+0.0i .$
$A$ is stored in compressed column storage format (CCS) and matrix multiplications are performed using the function matmul.

### 10.1  Program Text

Program Text (f01hbce.c)

### 10.2  Program Data

Program Data (f01hbce.d)

### 10.3  Program Results

Program Results (f01hbce.r)