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 (f01ha)

## Purpose

nag_matop_complex_gen_matrix_actexp (f01ha) 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.

## Syntax

[a, b, ifail] = f01ha(m, a, b, t, 'n', n)
[a, b, ifail] = nag_matop_complex_gen_matrix_actexp(m, a, b, t, 'n', n)

## 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 product etAB${e}^{tA}B$ without explicitly forming etA${e}^{tA}$.

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

### Compulsory Input Parameters

1:     m – int64int32nag_int scalar
m$m$, the number of columns of the matrix B$B$.
Constraint: m0${\mathbf{m}}\ge 0$.
2:     a(lda, : $:$) – complex array
The first dimension of the array a 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 n${\mathbf{n}}$
The n$n$ by n$n$ matrix A$A$.
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}}$
The n$n$ by m$m$ matrix B$B$.
4:     t – complex scalar
The scalar t$t$.

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the arrays a, b.
n$n$, the order of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.

lda ldb

### Output Parameters

1:     a(lda, : $:$) – complex array
The first dimension of the array a 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 n${\mathbf{n}}$
ldamax (1,n)$\mathit{lda}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
A$A$ is overwritten during the computation.
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)$.
The n$n$ by m$m$ matrix etAB${e}^{tA}B$.
3:     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 trying to balance the matrix A$A$.
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$
Constraint: n0${\mathbf{n}}\ge 0$.
ifail = 2${\mathbf{ifail}}=-2$
Constraint: m0${\mathbf{m}}\ge 0$.
ifail = 4${\mathbf{ifail}}=-4$
Constraint: ldamax (1,n)$\mathit{lda}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
ifail = 6${\mathbf{ifail}}=-6$
Constraint: ldbmax (1,n)$\mathit{ldb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
ifail = 999${\mathbf{ifail}}=-999$
Allocation of memory failed.

## Accuracy

For a 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.

The matrix etAB${e}^{tA}B$ could be computed by explicitly forming etA${e}^{tA}$ using nag_matop_complex_gen_matrix_exp (f01fc) and multiplying B$B$ by the result. However, experiments show that it is usually both more accurate and quicker to use nag_matop_complex_gen_matrix_actexp (f01ha).
The cost of the algorithm is O(n2m)$\mathit{O}\left({n}^{2}m\right)$. The precise cost depends on A$A$ since a combination of balancing, shifting and scaling is used prior to the Taylor series evaluation.
Approximately n2 + (2m + 8) n ${n}^{2}+\left(2m+8\right)n$ of complex allocatable memory is required by nag_matop_complex_gen_matrix_actexp (f01ha).
nag_matop_real_gen_matrix_actexp (f01ga) can be used to compute etAB${e}^{tA}B$ for real A$A$, B$B$, and t$t$. nag_matop_complex_gen_matrix_actexp_rcomm (f01hb) provides an implementation of the algorithm with a reverse communication interface, which returns control to the user when matrix multiplications are required. This should be used if A$A$ is large and sparse.

## Example

```function nag_matop_complex_gen_matrix_actexp_example
a = [0.5+0.0i, -0.2+0.0i, 1.0+0.1i, 0.0+0.4i;
0.3+0.0i,  0.5+1.2i, 3.1+0.0i, 1.0+0.2i;
0.0+2.0i,  0.1+0.0i, 1.2+0.2i, 0.5+0.0i;
1.0+0.3i,  0.0+0.2i, 0.0+0.9i, 0.5+0.0i];
b = [0.4+0.0i,  1.2+0.0i;
1.3+0.0i, -0.2+0.1i;
0.0+0.3i,  2.1+0.0i;
0.4+0.0i, -0.9+0.0i];
m = int64(2);
t = complex(-0.5);
% Compute exp(ta)b
[a, b, ifail] = nag_matop_complex_gen_matrix_actexp(m, a, b, t)
```
```

a =

-0.1750 - 0.3500i  -0.8000 + 0.0000i   1.0000 + 0.1000i   0.0000 + 0.8000i
0.0750 + 0.0000i  -0.1750 + 0.8500i   0.7750 + 0.0000i   0.5000 + 0.1000i
0.0000 + 2.0000i   0.4000 + 0.0000i   0.5250 - 0.1500i   1.0000 + 0.0000i
0.5000 + 0.1500i   0.0000 + 0.4000i   0.0000 + 0.4500i  -0.1750 - 0.3500i

b =

0.4251 - 0.1061i  -0.0220 + 0.3289i
0.7229 - 0.5940i  -1.7931 + 1.4952i
-0.1394 - 0.1151i   1.4781 - 0.4514i
0.1054 - 0.0786i  -1.0059 - 0.7079i

ifail =

0

```
```function f01ha_example
a = [0.5+0.0i, -0.2+0.0i, 1.0+0.1i, 0.0+0.4i;
0.3+0.0i,  0.5+1.2i, 3.1+0.0i, 1.0+0.2i;
0.0+2.0i,  0.1+0.0i, 1.2+0.2i, 0.5+0.0i;
1.0+0.3i,  0.0+0.2i, 0.0+0.9i, 0.5+0.0i];
b = [0.4+0.0i,  1.2+0.0i;
1.3+0.0i, -0.2+0.1i;
0.0+0.3i,  2.1+0.0i;
0.4+0.0i, -0.9+0.0i];
m = int64(2);
t = complex(-0.5);
% Compute exp(ta)b
[a, b, ifail] = f01ha(m, a, b, t)
```
```

a =

-0.1750 - 0.3500i  -0.8000 + 0.0000i   1.0000 + 0.1000i   0.0000 + 0.8000i
0.0750 + 0.0000i  -0.1750 + 0.8500i   0.7750 + 0.0000i   0.5000 + 0.1000i
0.0000 + 2.0000i   0.4000 + 0.0000i   0.5250 - 0.1500i   1.0000 + 0.0000i
0.5000 + 0.1500i   0.0000 + 0.4000i   0.0000 + 0.4500i  -0.1750 - 0.3500i

b =

0.4251 - 0.1061i  -0.0220 + 0.3289i
0.7229 - 0.5940i  -1.7931 + 1.4952i
-0.1394 - 0.1151i   1.4781 - 0.4514i
0.1054 - 0.0786i  -1.0059 - 0.7079i

ifail =

0

```