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

Note: the interface to this routine has changed since earlier releases of the toolbox:

At Mark 25:

m was made optional

Description

e^{t A} B

is computed using the algorithm described in Al–Mohy and Higham (2011) which uses a truncated Taylor series to compute the

e^{t A} B

without explicitly forming

e^{t A}

The algorithm does not explicity need to access the elements of

A

; it only requires the result of matrix multiplications of the form

A X

A^{T} 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 argument irevcm. Between intermediate exits and re-entries, all arguments 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$ .
2: $m$ – int64int32nag_int scalar: The number of columns of the matrix $B$ .

Constraint: $m \geq 0$ .
3: $b (ldb :)$ – double array: The first dimension of the array b must be at least $n$ .
The second dimension of the array b must be at least $m$ .

On initial entry: the $n$ by $m$ matrix $B$ .

On intermediate re-entry: must not be changed.
4: $t$ – double scalar: The scalar $t$ .
5: $b2 (ldb2 :)$ – double array: The first dimension of the array b2 must be at least $n$ .
The second dimension of the array b2 must be at least $m$ .

On initial entry: need not be set.

On intermediate re-entry: if $irevcm = 1$ , must contain $A B$ .
6: $x (ldx :)$ – double array: The first dimension of the array x must be at least $n$ .
The second dimension of the array x must be at least $2$ .

On initial entry: need not be set.

On intermediate re-entry: if $irevcm = 3$ , must contain $A^{T} Y$ .
7: $y (ldy :)$ – double array: The first dimension of the array y must be at least $n$ .
The second dimension of the array y must be at least $2$ .

On initial entry: need not be set.

On intermediate re-entry: if $irevcm = 2$ , must contain $A X$ .
8: $p (n)$ – double array: On initial entry: need not be set.

On intermediate re-entry: if $irevcm = 4$ , must contain $A z$ .
9: $r (n)$ – double array: On initial entry: need not be set.

On intermediate re-entry: if $irevcm = 5$ , must contain $A^{T} z$ .
10: $z (n)$ – double array: On initial entry: need not be set.

On intermediate re-entry: must not be changed.
11: $comm (n \times m + 3 \times n + 12)$ – double array
12: $icomm (2 \times n + 40)$ – int64int32nag_int array

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the first dimension of the arrays b, b2, x, y, comm and the dimension of the arrays icomm, p, r, z. (An error is raised if these dimensions are not equal.)
$n$ , the order of the matrix $A$ .

Constraint: $n \geq 0$ .
2: $tr$ – double scalar: Default: $0$
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^{tr t} B$ is immediately returned in the first row of $B$ . See Further Comments.

Output Parameters

1: $irevcm$ – int64int32nag_int scalar

On intermediate exit:

irevcm = 1

2

3

4

5

. The calling program must:

(a)

irevcm = 1

: evaluate

B_{2} = A B

, where

B_{2}

is an

n

m

matrix, and store the result in b2;
if

irevcm = 2

: evaluate

Y = A X

, where

X

and

Y

are

n

2

matrices, and store the result in y;
if

irevcm = 3

: evaluate

X = A^{T} Y

and store the result in x;
if

irevcm = 4

: evaluate

p = A z

and store the result in p;
if

irevcm = 5

: evaluate

r = A^{T} z

and store the result in r.

(b)

call nag_matop_real_gen_matrix_actexp_rcomm (f01gb) again with all other parameters unchanged.

On final exit:

irevcm = 0

2: $b (ldb :)$ – double array

The first dimension of the array b will be

n

The second dimension of the array b will be

m

On intermediate exit: if

irevcm = 1

, contains the

n

m

matrix

B

On final exit: the

n

m

matrix

e^{t A} B

3: $b2 (ldb2 :)$ – double array

The first dimension of the array b2 will be

n

The second dimension of the array b2 will be

m

On final exit: the array is undefined.

4: $x (ldx :)$ – double array

The first dimension of the array x will be

n

The second dimension of the array x will be

2

On intermediate exit: if

irevcm = 2

, contains the current

n

2

matrix

X

On final exit: the array is undefined.

5: $y (ldy :)$ – double array

The first dimension of the array y will be

n

The second dimension of the array y will be

2

On intermediate exit: if

irevcm = 3

, contains the current

n

2

matrix

Y

On final exit: the array is undefined.

6: $p (n)$ – double array

On final exit: the array is undefined.

7: $r (n)$ – double array

On final exit: the array is undefined.

8: $z (n)$ – double array

On intermediate exit: if

irevcm = 4

5

, contains the vector

z

On final exit: the array is undefined.

9: $comm (n \times m + 3 \times n + 12)$ – double array

10: $icomm (2 \times n + 40)$ – int64int32nag_int array

11: $ifail$ – int64int32nag_int scalar

On final exit:

ifail = 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.

W $ifail = 2$: $e^{t A} B$ has been computed using an IEEE double precision Taylor series, although the arithmetic precision is higher than IEEE double precision.

$ifail = - 1$: On initial entry, $irevcm =_$ .
Constraint: $irevcm = 0$ .

On intermediate re-entry, $irevcm =_$ .
Constraint: $irevcm = 1$ , $2$ , $3$ , $4$ or $5$ .

$ifail = - 2$: On initial entry, $n =_$ .
Constraint: $n \geq 0$ .

$ifail = - 3$: On initial entry, $m =_$ .
Constraint: $m \geq 0$ .

$ifail = - 5$: On initial entry, $ldb =_$ and $n =_$ .
Constraint: $ldb \geq n$ .

$ifail = - 9$: On initial entry, $ldb2 =_$ and $n =_$ .
Constraint: $ldb2 \geq n$ .

$ifail = - 11$: On initial entry, $ldx =_$ and $n =_$ .
Constraint: $ldx \geq n$ .

$ifail = - 13$: On initial entry, $ldy =_$ and $n =_$ .
Constraint: $ldy \geq n$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

For a symmetric matrix

A

(for which

A^{T} = A

) the computed matrix

e^{t A} 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.

Further Comments

Use of TrA

The elements of

A

are not explicitly required by nag_matop_real_gen_matrix_actexp_rcomm (f01gb). However, the trace of

A

is used in the preprocessing phase of the algorithm. If

T r (A)

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.

When to use nag_matop_real_gen_matrix_actexp_rcomm (f01gb)

nag_matop_real_gen_matrix_actexp_rcomm (f01gb) 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^{t A} B

will not, in general, be sparse even if

A

is sparse.

A

is small and dense then nag_matop_real_gen_matrix_actexp (f01ga) can be used to compute

e^{t A} B

without the use of a reverse communication interface.

The complex analog of nag_matop_real_gen_matrix_actexp_rcomm (f01gb) is nag_matop_complex_gen_matrix_actexp_rcomm (f01hb).

Use in Conjunction with NAG Toolbox Functions

To compute

e^{t A} 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'y and store in x
    case {4}
      % compute az and store in p
    otherwise
      % compute a'z and store in r
  end

Example

This example computes

e^{t A} B

, where

A = (\begin{matrix} 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 \end{matrix}),

B = (\begin{matrix} 0.1 & 1.1 \\ 1.7 & - 0.2 \\ 0.5 & 1.0 \\ 0.4 & - 0.2 \end{matrix}),

and

t = - 0.2 .

Open in the MATLAB editor: f01gb_example

function f01gb_example


fprintf('f01gb example results\n\n');

% A and B
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];
n = 4;
m = 2;

% Initial input parameters
t = -0.2;
b2 = b;  x = b;   y = b;
p = zeros(n, 1);  r = p;  z = p;
comm =  zeros(n*m+3*n+12, 1);
icomm = zeros(2*n+40, 1, 'int64');
tr = trace(a);

% Compute exp(ta)b by reverse communication
irevcm = int64(0);
done = false;

while ~done
  [irevcm, b, b2, x, y, p, r, z, comm, icomm, ifail] = ...
  f01gb( ...
         irevcm, b, t, b2, x, y, p, r, z, comm, icomm, 'tr', tr);

  switch irevcm
    case {0}
      done = true;
    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'y and store in x
      x = a'*y;
    case {4}
      % compute az and store in p
      p = a*z;
    otherwise
      % compute a'z and store in r
      r = a'*z;
  end
end

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

f01gb example results

exp(tA)B
    0.1933    0.7812
    1.4423   -0.4055
   -1.0756    0.6686
    0.0276    0.4900

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox