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_real_gen_matrix_actexp (f01ga)

Purpose

nag_matop_real_gen_matrix_actexp (f01ga) computes the action of the matrix exponential etA${e}^{tA}$, on the matrix B$B$, where A$A$ is a real n$n$ by n$n$ matrix, B$B$ is a real n$n$ by m$m$ matrix and t$t$ is a real scalar.

Syntax

[a, b, ifail] = f01ga(m, a, b, t, 'n', n)
[a, b, ifail] = nag_matop_real_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, : $:$) – double 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, : $:$) – double 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 – double 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, : $:$) – double 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, : $:$) – double 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 symmetric matrix A$A$ (for which AT = A${A}^{\mathrm{T}}=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-symmetric 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_real_gen_matrix_exp (f01ec) and multiplying B$B$ by the result. However, experiments show that it is usually both more accurate and quicker to use nag_matop_real_gen_matrix_actexp (f01ga).
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 real allocatable memory is required by nag_matop_real_gen_matrix_actexp (f01ga).
nag_matop_complex_gen_matrix_actexp (f01ha) can be used to compute etAB${e}^{tA}B$ for complex A$A$, B$B$, and t$t$. nag_matop_real_gen_matrix_actexp_rcomm (f01gb) 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_real_gen_matrix_actexp_example
a = [0.7,-0.2, 1.0, 0.3;
0.3, 0.7, 1.2, 1.0;
0.9, 0.0, 0.2, 0.7;
2.4, 0.1, 0.0, 0.2];
b = [0.1, 1.2;
1.3, 0.2;
0.0, 1.0;
0.4, -0.9];
m = int64(2);
t = 1.2;
% Compute exp(ta)b
[a, b, ifail] = nag_matop_real_gen_matrix_actexp(m, a, b, t)
```
```

a =

0.2500   -0.8000    2.0000    0.6000
0.0750    0.2500    0.6000    0.5000
0.4500         0   -0.2500    0.7000
1.2000    0.2000         0   -0.2500

b =

0.2138    7.6756
4.9980   11.6051
0.8307    7.5468
1.2406    9.7261

ifail =

0

```
```function f01ga_example
a = [0.7,-0.2, 1.0, 0.3;
0.3, 0.7, 1.2, 1.0;
0.9, 0.0, 0.2, 0.7;
2.4, 0.1, 0.0, 0.2];
b = [0.1, 1.2;
1.3, 0.2;
0.0, 1.0;
0.4, -0.9];
m = int64(2);
t = 1.2;
% Compute exp(ta)b
[a, b, ifail] = f01ga(m, a, b, t)
```
```

a =

0.2500   -0.8000    2.0000    0.6000
0.0750    0.2500    0.6000    0.5000
0.4500         0   -0.2500    0.7000
1.2000    0.2000         0   -0.2500

b =

0.2138    7.6756
4.9980   11.6051
0.8307    7.5468
1.2406    9.7261

ifail =

0

```