hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_matop_complex_gen_matrix_actexp_rcomm (f01hb)

Purpose

nag_matop_complex_gen_matrix_actexp_rcomm (f01hb) computes the action of the matrix exponential etAetA, on the matrix BB, where AA is a complex nn by nn matrix, BB is a complex nn by mm matrix and tt is a complex scalar. It uses reverse communication for evaluating matrix products, so that the matrix AA is not accessed explicitly.

Syntax

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

Description

etABetAB is computed using the algorithm described in Al–Mohy and Higham (2011) which uses a truncated Taylor series to compute the etABetAB without explicitly forming etAetA.
The algorithm does not explicity need to access the elements of AA; it only requires the result of matrix multiplications of the form AXAX or AHYAHY. 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 parameter irevcm. Between intermediate exits and re-entries, all parameters 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 00.
2:     m – int64int32nag_int scalar
The number of columns of the matrix BB.
Constraint: m0m0.
3:     b(ldb, : :) – complex array
The first dimension of the array b must be at least max (1,n)max(1,n)
The second dimension of the array must be at least mm
On initial entry: the nn by mm matrix BB.
On intermediate re-entry: must not be changed.
4:     t – complex scalar
The scalar tt.
5:     b2(ldb2, : :) – complex array
The first dimension of the array b2 must be at least max (1,n)max(1,n)
The second dimension of the array must be at least mm
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 1irevcm=1, must contain ABAB.
6:     x(ldx, : :) – complex array
The first dimension of the array x must be at least max (1,n)max(1,n)
The second dimension of the array must be at least 22
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 3irevcm=3, must contain AHYAHY.
7:     y(ldy, : :) – complex array
The first dimension of the array y must be at least max (1,n)max(1,n)
The second dimension of the array must be at least 22
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 2irevcm=2, must contain AXAX.
8:     p(n) – complex array
n, the dimension of the array, must satisfy the constraint n0n0.
On initial entry: need not be set.
n, the dimension of the array, must satisfy the constraint n0n0.
On intermediate re-entry: if irevcm = 4irevcm=4, must contain AzAz.
9:     r(n) – complex array
n, the dimension of the array, must satisfy the constraint n0n0.
On initial entry: need not be set.
n, the dimension of the array, must satisfy the constraint n0n0.
On intermediate re-entry: if irevcm = 5irevcm=5, must contain AHzAHz.
10:   z(n) – complex array
n, the dimension of the array, must satisfy the constraint n0n0.
On initial entry: need not be set.
n, the dimension of the array, must satisfy the constraint n0n0.
On intermediate re-entry: must not be changed.
11:   ccomm(n × (m + 2)n×(m+2)) – complex array
12:   comm(3 × n + 143×n+14) – double array
13:   icomm(2 × n + 402×n+40) – int64int32nag_int array

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the arrays b2, x, y, ccomm and the dimension of the arrays p, r, z. (An error is raised if these dimensions are not equal.)
nn, the order of the matrix AA.
Constraint: n0n0.
2:     tr – complex scalar
The trace of AA. If this is not available then any number can be supplied (00 is a reasonable default); however, in the trivial case, n = 1n=1, the result etrtBetrtB is immediately returned in the first row of BB. See Section [Further Comments].
Default: 00

Input Parameters Omitted from the MATLAB Interface

ldb ldb2 ldx ldy

Output Parameters

1:     irevcm – int64int32nag_int scalar
On intermediate exit: irevcm = 1irevcm=1, 22, 33, 44 or 55. The calling program must:
(a) if irevcm = 1irevcm=1: evaluate B2 = ABB2=AB, where B2B2 is an nn by mm matrix, and store the result in b2;
if irevcm = 2irevcm=2: evaluate Y = AXY=AX, where XX and YY are nn by 22 matrices, and store the result in y;
if irevcm = 3irevcm=3: evaluate X = AHYX=AHY and store the result in x;
if irevcm = 4irevcm=4: evaluate p = Azp=Az and store the result in p;
if irevcm = 5irevcm=5: evaluate r = AHzr=AHz and store the result in r.
(b) call nag_matop_complex_gen_matrix_actexp_rcomm (f01hb) again with all other parameters unchanged.
On final exit: irevcm = 0irevcm=0.
2:     b(ldb, : :) – complex array
The first dimension of the array b will be max (1,n)max(1,n)
The second dimension of the array will be mm
ldbmax (1,n)ldbmax(1,n).
On intermediate exit: if irevcm = 1irevcm=1, contains the nn by mm matrix BB.
ldbmax (1,n)ldbmax(1,n).
On final exit: the nn by mm matrix etABetAB.
3:     b2(ldb2, : :) – complex array
The first dimension of the array b2 will be max (1,n)max(1,n)
The second dimension of the array will be mm
ldb2max (1,n)ldb2max(1,n).
On final exit: the array is undefined.
4:     x(ldx, : :) – complex array
The first dimension of the array x will be max (1,n)max(1,n)
The second dimension of the array will be 22
ldxmax (1,n)ldxmax(1,n).
On intermediate exit: if irevcm = 2irevcm=2, contains the current nn by 22 matrix XX.
ldxmax (1,n)ldxmax(1,n).
On final exit: the array is undefined.
5:     y(ldy, : :) – complex array
The first dimension of the array y will be max (1,n)max(1,n)
The second dimension of the array will be 22
ldymax (1,n)ldymax(1,n).
On intermediate exit: if irevcm = 3irevcm=3, contains the current nn by 22 matrix YY.
ldymax (1,n)ldymax(1,n).
On final exit: the array is undefined.
6:     p(n) – complex array
On final exit: the array is undefined.
7:     r(n) – complex array
On final exit: the array is undefined.
8:     z(n) – complex array
On intermediate exit: if irevcm = 4irevcm=4 or 55, contains the vector zz.
On final exit: the array is undefined.
9:     ccomm(n × (m + 2)n×(m+2)) – complex array
10:   comm(3 × n + 143×n+14) – double array
11:   icomm(2 × n + 402×n+40) – int64int32nag_int array
12:   ifail – int64int32nag_int scalar
ifail = 0ifail=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 = 1ifail=1
Note:  this failure should not occur, and suggests that the function has been called incorrectly. An unexpected internal error occurred when estimating a matrix norm.
W ifail = 2ifail=2
etABetAB has been computed using an IEEE double precision Taylor series, although the arithmetic precision is higher than IEEE double precision.
  ifail = 1ifail=-1
On initial entry.
Constraint: irevcm = 0irevcm=0.
On intermediate re-entry.
Constraint: irevcm = 1irevcm=1, 22, 33, 44 or 55.
  ifail = 2ifail=-2
On initial entry.
Constraint: n0n0.
  ifail = 3ifail=-3
On initial entry.
Constraint: m0m0.
  ifail = 5ifail=-5
On initial entry, ldb = _ldb=_ and .
Constraint: ldbmax (1,n)ldbmax(1,n).
  ifail = 9ifail=-9
On initial entry, ldb2 = _ldb2=_ and .
Constraint: ldb2max (1,n)ldb2max(1,n).
  ifail = 11ifail=-11
On initial entry, ldx = _ldx=_ and .
Constraint: ldxmax (1,n)ldxmax(1,n).
  ifail = 13ifail=-13
On initial entry, ldy = _ldy=_ and .
Constraint: ldymax (1,n)ldymax(1,n).

Accuracy

For an Hermitian matrix AA (for which AH = AAH=A) the computed matrix etABetAB 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.

Further Comments

Use of TrA

The elements of AA are not explicitly required by nag_matop_complex_gen_matrix_actexp_rcomm (f01hb). However, the trace of AA is used in the preprocessing phase of the algorithm. If Tr(A)Tr(A) is not available to the calling function then any number can be supplied (00 is recommended). This will not affect the stability of the algorithm, but it may reduce its efficiency.

When to use nag_matop_complex_gen_matrix_actexp_rcomm (f01hb)

nag_matop_complex_gen_matrix_actexp_rcomm (f01hb) is designed to be used when AA 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 etABetAB will not, in general, be sparse even if AA is sparse.
If AA is small and dense then nag_matop_complex_gen_matrix_actexp (f01ha) can be used to compute etABetAB without the use of a reverse communication interface.
The real analog of nag_matop_complex_gen_matrix_actexp_rcomm (f01hb) is nag_matop_real_gen_matrix_actexp_rcomm (f01gb).

Use in Conjunction with NAG Toolbox Functions

To compute etABetAB, 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^Hy and store in x
    case {4}
      % compute az and store in p
    case {5}
      % compute a^Hz and store in r
  end

Example

function nag_matop_complex_gen_matrix_actexp_rcomm_example
a = [0.7+0.8i, -0.2+0.0i, 1.0+0.0i, 0.6+0.5i;
     0.3+0.7i,  0.7+0.0i, 0.9+3.0i, 1.0+0.8i;
     0.3+3.0i, -7.0+0.0i, 0.2+0.6i, 0.7+0.5i;
     0.0+0.9i,  4.0+0.0i, 0.0+0.0i, 0.2+0.0i];
b = [0.1+0.0i,  1.2+0.1i;
     1.3+0.9i, -0.2+2.0i;
     4.0+0.6i, -1.0+0.8i;
     0.4+0.0i, -0.9+0.0i];
n = 4;
m = int64(2);
t = complex(1.1);
b2 = complex(zeros(n, m));
x = complex(zeros(n, 2));
y = complex(zeros(n, 2));
p = complex(zeros(n, 1));
r = complex(zeros(n, 1));
z = complex(zeros(n, 1));
comm =  zeros(3*n+14, 1);
icomm = zeros(2*n+40, 1, 'int64');
ccomm = complex(zeros(n*(m+2), 1));
% Compute the trace of a
tr = trace(a);

% Compute exp(ta)b
irevcm = int64(0);
[irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = ...
    nag_matop_complex_gen_matrix_actexp_rcomm(irevcm, m, b, t, b2, x, y, ...
                                     p, r, z, ccomm, comm, icomm, 'tr', tr);

while irevcm ~= 0
  switch irevcm
    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^Hy and store in x
      x = a'*y;
    case {4}
      % compute az and store in p
      p = a*z;
    case {5}
      % compute a^Hz and store in r
      r = a'*z;
  end

  [irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = ...
      nag_matop_complex_gen_matrix_actexp_rcomm(irevcm, m, b, t, b2, x, y, ...
                                       p, r, z, ccomm, comm, icomm, 'tr', tr);
end

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

 exp(tA)B
 -15.3125 + 5.9123i  -4.5605 - 2.4288i
  12.3396 -50.6993i   9.2005 -10.3632i
 -65.4353 +34.3271i -17.6075 - 1.0019i
  45.6506 -28.3253i  11.3339 + 0.1127i


function f01hb_example
a = [0.7+0.8i, -0.2+0.0i, 1.0+0.0i, 0.6+0.5i;
     0.3+0.7i,  0.7+0.0i, 0.9+3.0i, 1.0+0.8i;
     0.3+3.0i, -7.0+0.0i, 0.2+0.6i, 0.7+0.5i;
     0.0+0.9i,  4.0+0.0i, 0.0+0.0i, 0.2+0.0i];
b = [0.1+0.0i,  1.2+0.1i;
     1.3+0.9i, -0.2+2.0i;
     4.0+0.6i, -1.0+0.8i;
     0.4+0.0i, -0.9+0.0i];
n = 4;
m = int64(2);
t = complex(1.1);
b2 = complex(zeros(n, m));
x = complex(zeros(n, 2));
y = complex(zeros(n, 2));
p = complex(zeros(n, 1));
r = complex(zeros(n, 1));
z = complex(zeros(n, 1));
comm =  zeros(3*n+14, 1);
icomm = zeros(2*n+40, 1, 'int64');
ccomm = complex(zeros(n*(m+2), 1));
% Compute the trace of a
tr = trace(a);

% Compute exp(ta)b
irevcm = int64(0);
[irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = ...
    f01hb(irevcm, m, b, t, b2, x, y, p, r, z, ccomm, comm, icomm, 'tr', tr);

while irevcm ~= 0
  switch irevcm
    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^Hy and store in x
      x = a'*y;
    case {4}
      % compute az and store in p
      p = a*z;
    case {5}
      % compute a^Hz and store in r
      r = a'*z;
  end

  [irevcm, b, b2, x, y, p, r, z, ccomm, comm, icomm, ifail] = ...
      f01hb(irevcm, m, b, t, b2, x, y, p, r, z, ccomm, comm, icomm, 'tr', tr);
end

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

 exp(tA)B
 -15.3125 + 5.9123i  -4.5605 - 2.4288i
  12.3396 -50.6993i   9.2005 -10.3632i
 -65.4353 +34.3271i -17.6075 - 1.0019i
  45.6506 -28.3253i  11.3339 + 0.1127i



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013