NAG Library Function Document

nag_matop_real_gen_matrix_actexp_rcomm (f01gbc)

void	nag_matop_real_gen_matrix_actexp_rcomm (Integer irevcm, Integer n, Integer m, double b[], Integer pdb, double t, double tr, double b2[], Integer pdb2, double x[], Integer pdx, double y[], Integer pdy, double p[], double r[], double z[], double comm[], Integer icomm[], NagError fail)

3 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.

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: irevcm – Integer *Input/Output

On initial entry: must be set to

0

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 (f01gbc) again with all other parameters unchanged.

On final exit:

irevcm = 0

2: n – IntegerInput

On entry:

n

, the order of the matrix

A

Constraint:

n \geq 0

3: m – IntegerInput

On entry: the number of columns of the matrix

B

Constraint:

m \geq 0

4: b[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array b must be at least

pdb \times m

The

(i, j)

th element of the matrix

B

is stored in

b [(j - 1) \times pdb + i - 1]

On initial entry: the

n

m

matrix

B

On intermediate exit: if

irevcm = 1

, contains the

n

m

matrix

B

On intermediate re-entry: must not be changed.

On final exit: the

n

m

matrix

e^{t A} B

5: pdb – IntegerInput

On entry: the stride separating matrix row elements in the array b.

Constraint:

pdb \geq n

6: t – doubleInput

On entry: the scalar

t

7: tr – doubleInput

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^{tr t} B

is immediately returned in the first row of

B

. See Section 9.

8: b2[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array b2 must be at least

pdb2 \times m

The

(i, j)

th element of the matrix is stored in

b2 [(j - 1) \times pdb2 + i - 1]

On initial entry: need not be set.

On intermediate re-entry: if

irevcm = 1

, must contain

A B

On final exit: the array is undefined.

9: pdb2 – IntegerInput

On entry: the stride separating matrix row elements in the array b2.

Constraint:

pdb2 \geq n

10: x[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array x must be at least

pdx \times 2

The

(i, j)

th element of the matrix

X

is stored in

x [(j - 1) \times pdx + i - 1]

On initial entry: need not be set.

On intermediate exit: if

irevcm = 2

, contains the current

n

2

matrix

X

On intermediate re-entry: if

irevcm = 3

, must contain

A^{T} Y

On final exit: the array is undefined.

11: pdx – IntegerInput

On entry: the stride separating matrix row elements in the array x.

Constraint:

pdx \geq n

12: y[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array y must be at least

pdy \times 2

The

(i, j)

th element of the matrix

Y

is stored in

y [(j - 1) \times pdy + i - 1]

On initial entry: need not be set.

On intermediate exit: if

irevcm = 3

, contains the current

n

2

matrix

Y

On intermediate re-entry: if

irevcm = 2

, must contain

A X

On final exit: the array is undefined.

13: pdy – IntegerInput

On entry: the stride separating matrix row elements in the array y.

Constraint:

pdy \geq n

14: p[n] – doubleInput/Output

On initial entry: need not be set.

On intermediate re-entry: if

irevcm = 4

, must contain

A z

On final exit: the array is undefined.

15: r[n] – doubleInput/Output

On initial entry: need not be set.

On intermediate re-entry: if

irevcm = 5

, must contain

A^{T} z

On final exit: the array is undefined.

16: z[n] – doubleInput/Output

On initial entry: need not be set.

On intermediate exit: if

irevcm = 4

5

, contains the vector

z

On intermediate re-entry: must not be changed.

On final exit: the array is undefined.

17: comm[ $n \times m + 3 \times n + 12$ ] – doubleCommunication Array

18: icomm[ $2 \times n + 40$ ] – IntegerCommunication Array

19: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_BAD_PARAM: On entry, argument $⟨value⟩$ had an illegal value.
NE_INT: On entry, $m = ⟨value⟩$ .
Constraint: $m \geq 0$ .
On entry, $n = ⟨value⟩$ .
Constraint: $n \geq 0$ .
On initial entry, $irevcm = ⟨value⟩$ .
Constraint: $irevcm = 0$ .
On intermediate re-entry, $irevcm = ⟨value⟩$ .
Constraint: $irevcm = 1$ , $2$ , $3$ , $4$ or $5$ .
NE_INT_2: On entry, $pdb = ⟨value⟩$ and $n = ⟨value⟩$ .
Constraint: $pdb \geq n$ .
On entry, $pdb2 = ⟨value⟩$ and $n = ⟨value⟩$ .
Constraint: $pdb2 \geq n$ .
On entry, $pdx = ⟨value⟩$ and $n = ⟨value⟩$ .
Constraint: $pdx \geq n$ .
On entry, $pdy = ⟨value⟩$ and $n = ⟨value⟩$ .
Constraint: $pdy \geq 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^{t A} B$ has been computed using an IEEE double precision Taylor series, although the arithmetic precision is higher than IEEE double precision.

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

8 Parallelism and Performance

nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

Please consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

9.1 Use of $T r (A)$

The elements of

A

are not explicitly required by nag_matop_real_gen_matrix_actexp_rcomm (f01gbc). 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.

9.2 When to use nag_matop_real_gen_matrix_actexp_rcomm (f01gbc)

nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) 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 (f01gac) 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 (f01gbc) is nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc).

9.3 Use in Conjunction with NAG C Library Functions

To compute

e^{t A} B

, the following skeleton code can normally be used:

do {
f01gbc(&irevcm,n,m,b,tdb,t,tr,b2,tdb2,x,tdx,y,tdy,p,r,z,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^T Y ..
  }
  else if (irevcm == 4){
    .. Code to compute P=AZ ..
  }
  else if (irevcm == 5){
    .. Code to compute R=A^T 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_dgemm (f16yac)) can be used. If

A

is triangular then nag_dtrmm (f16yfc) should be used. If

A

is symmetric, then nag_dsymm (f16ycc) should be used. For sparse

A

stored in coordinate storage format nag_sparse_nsym_matvec (f11xac) and nag_sparse_sym_matvec (f11xec) can be used. Alternatively if

A

is stored in compressed column format nag_superlu_matrix_product (f11mkc) can be used.

10 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 .

NAG Library Function Documentnag_matop_real_gen_matrix_actexp_rcomm (f01gbc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Use of TrA

9.2 When to use nag_matop_real_gen_matrix_actexp_rcomm (f01gbc)

9.3 Use in Conjunction with NAG C Library Functions

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Function Document

nag_matop_real_gen_matrix_actexp_rcomm (f01gbc)

9.1 Use of $T r (A)$