NAG Library Function Document

nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc)

nag_matop_complex_gen_matrix_actexp_rcomm (Integer *irevcm, Integer n, Integer m, Complex b[], Integer pdb, Complex t, Complex tr, Complex b2[], Integer pdb2, Complex x[], Integer pdx, Complex y[], Integer pdy, Complex p[], Complex r[], Complex z[], Complex ccomm[], 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^{H} 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^{H} 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^{H} z

and store the result in r.

(b)

call nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) 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$ ] – ComplexInput/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 – ComplexInput

On entry: the scalar

t

7: tr – ComplexInput

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$ ] – ComplexInput/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$ ] – ComplexInput/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^{H} 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$ ] – ComplexInput/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] – ComplexInput/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] – ComplexInput/Output

On initial entry: need not be set.

On intermediate re-entry: if

irevcm = 5

, must contain

A^{H} z

On final exit: the array is undefined.

16: z[n] – ComplexInput/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: ccomm[ $n \times (m + 2)$ ] – ComplexCommunication Array

18: comm[ $3 \times n + 14$ ] – doubleCommunication Array

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

20: 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 an Hermitian matrix

A

(for which

A^{H} = 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-Hermitian matrices. See Section 4 of Al–Mohy and Higham (2011) for details and further discussion.

8 Parallelism and Performance

Not applicable.

9 Further Comments

9.1 Use of $T r (A)$

The elements of

A

are not explicitly required by nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc). 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_complex_gen_matrix_actexp_rcomm (f01hbc)

nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) 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_complex_gen_matrix_actexp (f01hac) can be used to compute

e^{t A} B

without the use of a reverse communication interface.

The real analog of nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) is nag_matop_real_gen_matrix_actexp_rcomm (f01gbc).

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

A

is triangular then nag_ztrmm (f16zfc) should be used. If

A

is Hermitian, then nag_zhemm (f16zcc) should be used. If

A

is symmetric, then nag_zsymm (f16ztc) should be used. For sparse

A

stored in coordinate storage format nag_sparse_nherm_matvec (f11xnc) and nag_sparse_herm_matvec (f11xsc) can be used. For sparse

A

stored in compressed column storage format (CCS) the program text of Section 10 contains the function matmul to perform matrix products.

10 Example

This example computes

e^{t A} B

where

A = (\begin{array}{r} 0.7 + 0.8 i & - 0.2 + 0.0 i & 1.0 + 0.0 i & 0.6 + 0.5 i \\ 0.3 + 0.7 i & 0.7 + 0.0 i & 0.9 + 3.0 i & 1.0 + 0.8 i \\ 0.3 + 3.0 i & - 0.7 + 0.0 i & 0.2 + 0.6 i & 0.7 + 0.5 i \\ 0.0 + 0.9 i & 4.0 + 0.0 i & 0.0 + 0.0 i & 0.2 + 0.0 i \end{array}),

B = (\begin{matrix} 0.1 + 0.0 i & 1.2 + 0.1 i \\ 1.3 + 0.9 i & - 0.2 + 2.0 i \\ 4.0 + 0.6 i & - 1.0 + 0.8 i \\ 0.4 + 0.0 i & - 0.9 + 0.0 i \end{matrix})

and

t = 1.1 + 0.0 i .

A

is stored in compressed column storage format (CCS) and matrix multiplications are performed using the function matmul.

NAG Library Function Documentnag_matop_complex_gen_matrix_actexp_rcomm (f01hbc)

+− 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_complex_gen_matrix_actexp_rcomm (f01hbc)

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_complex_gen_matrix_actexp_rcomm (f01hbc)

9.1 Use of $T r (A)$