Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_matop_complex_gen_matrix_actexp_rcomm (f01hb)

## Purpose

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

## Syntax

[irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = f01hb(irevcm, m, b, t, b2, x, y, p, r, z, ccomm, comm, icomm, 'n', n, 'tr', tr)
[irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = nag_matop_complex_gen_matrix_actexp_rcomm(irevcm, m, b, t, b2, x, y, p, r, z, ccomm, comm, icomm, 'n', n, 'tr', tr)

## Description

etAB${e}^{tA}B$ is computed using the algorithm described in Al–Mohy and Higham (2011) which uses a truncated Taylor series to compute the etAB${e}^{tA}B$ without explicitly forming etA${e}^{tA}$.
The algorithm does not explicity need to access the elements of A$A$; it only requires the result of matrix multiplications of the form AX$AX$ or AHY${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.

## 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

## Parameters

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 parameter irevcm. Between intermediate exits and re-entries, all parameters other than b2, x, y, p and r must remain unchanged.

### Compulsory Input Parameters

1:     irevcm – int64int32nag_int scalar
On initial entry: must be set to 0$0$.
2:     m – int64int32nag_int scalar
The number of columns of the matrix B$B$.
Constraint: m0${\mathbf{m}}\ge 0$.
3:     b(ldb, : $:$) – complex array
The first dimension of the array b must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array must be at least m${\mathbf{m}}$
On initial entry: the n$n$ by m$m$ matrix B$B$.
On intermediate re-entry: must not be changed.
4:     t – complex scalar
The scalar t$t$.
5:     b2(ldb2, : $:$) – complex array
The first dimension of the array b2 must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array must be at least m${\mathbf{m}}$
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 1${\mathbf{irevcm}}=1$, must contain AB$AB$.
6:     x(ldx, : $:$) – complex array
The first dimension of the array x must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array must be at least 2$2$
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 3${\mathbf{irevcm}}=3$, must contain AHY${A}^{\mathrm{H}}Y$.
7:     y(ldy, : $:$) – complex array
The first dimension of the array y must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array must be at least 2$2$
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 2${\mathbf{irevcm}}=2$, must contain AX$AX$.
8:     p(n) – complex array
n, the dimension of the array, must satisfy the constraint n0${\mathbf{n}}\ge 0$.
On initial entry: need not be set.
n, the dimension of the array, must satisfy the constraint n0${\mathbf{n}}\ge 0$.
On intermediate re-entry: if irevcm = 4${\mathbf{irevcm}}=4$, must contain Az$Az$.
9:     r(n) – complex array
n, the dimension of the array, must satisfy the constraint n0${\mathbf{n}}\ge 0$.
On initial entry: need not be set.
n, the dimension of the array, must satisfy the constraint n0${\mathbf{n}}\ge 0$.
On intermediate re-entry: if irevcm = 5${\mathbf{irevcm}}=5$, must contain AHz${A}^{\mathrm{H}}z$.
10:   z(n) – complex array
n, the dimension of the array, must satisfy the constraint n0${\mathbf{n}}\ge 0$.
On initial entry: need not be set.
n, the dimension of the array, must satisfy the constraint n0${\mathbf{n}}\ge 0$.
On intermediate re-entry: must not be changed.
11:   ccomm(n × (m + 2)${\mathbf{n}}×\left({\mathbf{m}}+2\right)$) – complex array
12:   comm(3 × n + 14$3×{\mathbf{n}}+14$) – double array
13:   icomm(2 × n + 40$2×{\mathbf{n}}+40$) – int64int32nag_int array

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the arrays b2, x, y, ccomm and the dimension of the arrays p, r, z. (An error is raised if these dimensions are not equal.)
n$n$, the order of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.
2:     tr – complex scalar
The trace of A$A$. If this is not available then any number can be supplied (0$0$ is a reasonable default); however, in the trivial case, n = 1$n=1$, the result etrtB${e}^{{\mathbf{tr}}t}B$ is immediately returned in the first row of B$B$. See Section [Further Comments].
Default: 0$0$

ldb ldb2 ldx ldy

### Output Parameters

1:     irevcm – int64int32nag_int scalar
On intermediate exit: irevcm = 1${\mathbf{irevcm}}=1$, 2$2$, 3$3$, 4$4$ or 5$5$. The calling program must:
 (a) if irevcm = 1${\mathbf{irevcm}}=1$: evaluate B2 = AB${B}_{2}=AB$, where B2${B}_{2}$ is an n$n$ by m$m$ matrix, and store the result in b2; if irevcm = 2${\mathbf{irevcm}}=2$: evaluate Y = AX$Y=AX$, where X$X$ and Y$Y$ are n$n$ by 2$2$ matrices, and store the result in y; if irevcm = 3${\mathbf{irevcm}}=3$: evaluate X = AHY$X={A}^{\mathrm{H}}Y$ and store the result in x; if irevcm = 4${\mathbf{irevcm}}=4$: evaluate p = Az$p=Az$ and store the result in p; if irevcm = 5${\mathbf{irevcm}}=5$: evaluate r = AHz$r={A}^{\mathrm{H}}z$ and store the result in r. (b) call nag_matop_complex_gen_matrix_actexp_rcomm (f01hb) again with all other parameters unchanged.
On final exit: irevcm = 0${\mathbf{irevcm}}=0$.
2:     b(ldb, : $:$) – complex array
The first dimension of the array b will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array will be m${\mathbf{m}}$
ldbmax (1,n)$\mathit{ldb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
On intermediate exit: if irevcm = 1${\mathbf{irevcm}}=1$, contains the n$n$ by m$m$ matrix B$B$.
ldbmax (1,n)$\mathit{ldb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
On final exit: the n$n$ by m$m$ matrix etAB${e}^{tA}B$.
3:     b2(ldb2, : $:$) – complex array
The first dimension of the array b2 will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array will be m${\mathbf{m}}$
ldb2max (1,n)$\mathit{ldb2}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
On final exit: the array is undefined.
4:     x(ldx, : $:$) – complex array
The first dimension of the array x will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array will be 2$2$
ldxmax (1,n)$\mathit{ldx}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
On intermediate exit: if irevcm = 2${\mathbf{irevcm}}=2$, contains the current n$n$ by 2$2$ matrix X$X$.
ldxmax (1,n)$\mathit{ldx}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
On final exit: the array is undefined.
5:     y(ldy, : $:$) – complex array
The first dimension of the array y will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array will be 2$2$
ldymax (1,n)$\mathit{ldy}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
On intermediate exit: if irevcm = 3${\mathbf{irevcm}}=3$, contains the current n$n$ by 2$2$ matrix Y$Y$.
ldymax (1,n)$\mathit{ldy}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
On final exit: the array is undefined.
6:     p(n) – complex array
On final exit: the array is undefined.
7:     r(n) – complex array
On final exit: the array is undefined.
8:     z(n) – complex array
On intermediate exit: if irevcm = 4${\mathbf{irevcm}}=4$ or 5$5$, contains the vector z$z$.
On final exit: the array is undefined.
9:     ccomm(n × (m + 2)${\mathbf{n}}×\left({\mathbf{m}}+2\right)$) – complex array
10:   comm(3 × n + 14$3×{\mathbf{n}}+14$) – double array
11:   icomm(2 × n + 40$2×{\mathbf{n}}+40$) – int64int32nag_int array
12:   ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

ifail = 1${\mathbf{ifail}}=1$
Note:  this failure should not occur, and suggests that the function has been called incorrectly. An unexpected internal error occurred when estimating a matrix norm.
W ifail = 2${\mathbf{ifail}}=2$
etAB${e}^{tA}B$ has been computed using an IEEE double precision Taylor series, although the arithmetic precision is higher than IEEE double precision.
ifail = 1${\mathbf{ifail}}=-1$
On initial entry.
Constraint: irevcm = 0${\mathbf{irevcm}}=0$.
On intermediate re-entry.
Constraint: irevcm = 1${\mathbf{irevcm}}=1$, 2$2$, 3$3$, 4$4$ or 5$5$.
ifail = 2${\mathbf{ifail}}=-2$
On initial entry.
Constraint: n0${\mathbf{n}}\ge 0$.
ifail = 3${\mathbf{ifail}}=-3$
On initial entry.
Constraint: m0${\mathbf{m}}\ge 0$.
ifail = 5${\mathbf{ifail}}=-5$
On initial entry, ldb = _$\mathit{ldb}=_$ and .
Constraint: ldbmax (1,n)$\mathit{ldb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
ifail = 9${\mathbf{ifail}}=-9$
On initial entry, ldb2 = _$\mathit{ldb2}=_$ and .
Constraint: ldb2max (1,n)$\mathit{ldb2}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
ifail = 11${\mathbf{ifail}}=-11$
On initial entry, ldx = _$\mathit{ldx}=_$ and .
Constraint: ldxmax (1,n)$\mathit{ldx}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
ifail = 13${\mathbf{ifail}}=-13$
On initial entry, ldy = _$\mathit{ldy}=_$ and .
Constraint: ldymax (1,n)$\mathit{ldy}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.

## Accuracy

For an Hermitian matrix A$A$ (for which AH = A${A}^{\mathrm{H}}=A$) the computed matrix etAB${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.

### Use of TrA

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

### When to use nag_matop_complex_gen_matrix_actexp_rcomm (f01hb)

nag_matop_complex_gen_matrix_actexp_rcomm (f01hb) is designed to be used when A$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 etAB${e}^{tA}B$ will not, in general, be sparse even if A$A$ is sparse.
If A$A$ is small and dense then nag_matop_complex_gen_matrix_actexp (f01ha) can be used to compute etAB${e}^{tA}B$ without the use of a reverse communication interface.
The real analog of nag_matop_complex_gen_matrix_actexp_rcomm (f01hb) is nag_matop_real_gen_matrix_actexp_rcomm (f01gb).

### Use in Conjunction with NAG Toolbox Functions

To compute etAB${e}^{tA}B$, the following skeleton code can normally be used:
```while irevcm ~= 0
switch irevcm
case {1}
% Compute ab and store in b2
case {2}
% Compute ax and store in y
case {3}
% Compute a^Hy and store in x
case {4}
% compute az and store in p
case {5}
% compute a^Hz and store in r
end
```

## Example

```function nag_matop_complex_gen_matrix_actexp_rcomm_example
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, -7.0+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];
n = 4;
m = int64(2);
t = complex(1.1);
b2 = complex(zeros(n, m));
x = complex(zeros(n, 2));
y = complex(zeros(n, 2));
p = complex(zeros(n, 1));
r = complex(zeros(n, 1));
z = complex(zeros(n, 1));
comm =  zeros(3*n+14, 1);
icomm = zeros(2*n+40, 1, 'int64');
ccomm = complex(zeros(n*(m+2), 1));
% Compute the trace of a
tr = trace(a);

% Compute exp(ta)b
irevcm = int64(0);
[irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = ...
nag_matop_complex_gen_matrix_actexp_rcomm(irevcm, m, b, t, b2, x, y, ...
p, r, z, ccomm, comm, icomm, 'tr', tr);

while irevcm ~= 0
switch irevcm
case {1}
% Compute ab and store in b2
b2 = a*b;
case {2}
% Compute ax and store in y
y = a*x;
case {3}
% Compute a^Hy and store in x
x = a'*y;
case {4}
% compute az and store in p
p = a*z;
case {5}
% compute a^Hz and store in r
r = a'*z;
end

[irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = ...
nag_matop_complex_gen_matrix_actexp_rcomm(irevcm, m, b, t, b2, x, y, ...
p, r, z, ccomm, comm, icomm, 'tr', tr);
end

fprintf('\n exp(tA)B\n');
disp(b);
```
```

exp(tA)B
-15.3125 + 5.9123i  -4.5605 - 2.4288i
12.3396 -50.6993i   9.2005 -10.3632i
-65.4353 +34.3271i -17.6075 - 1.0019i
45.6506 -28.3253i  11.3339 + 0.1127i

```
```function f01hb_example
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, -7.0+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];
n = 4;
m = int64(2);
t = complex(1.1);
b2 = complex(zeros(n, m));
x = complex(zeros(n, 2));
y = complex(zeros(n, 2));
p = complex(zeros(n, 1));
r = complex(zeros(n, 1));
z = complex(zeros(n, 1));
comm =  zeros(3*n+14, 1);
icomm = zeros(2*n+40, 1, 'int64');
ccomm = complex(zeros(n*(m+2), 1));
% Compute the trace of a
tr = trace(a);

% Compute exp(ta)b
irevcm = int64(0);
[irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = ...
f01hb(irevcm, m, b, t, b2, x, y, p, r, z, ccomm, comm, icomm, 'tr', tr);

while irevcm ~= 0
switch irevcm
case {1}
% Compute ab and store in b2
b2 = a*b;
case {2}
% Compute ax and store in y
y = a*x;
case {3}
% Compute a^Hy and store in x
x = a'*y;
case {4}
% compute az and store in p
p = a*z;
case {5}
% compute a^Hz and store in r
r = a'*z;
end

[irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = ...
f01hb(irevcm, m, b, t, b2, x, y, p, r, z, ccomm, comm, icomm, 'tr', tr);
end

fprintf('\n exp(tA)B\n');
disp(b);
```
```

exp(tA)B
-15.3125 + 5.9123i  -4.5605 - 2.4288i
12.3396 -50.6993i   9.2005 -10.3632i
-65.4353 +34.3271i -17.6075 - 1.0019i
45.6506 -28.3253i  11.3339 + 0.1127i

```