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_complex_gen_matrix_fun_usd (f01fm)

## Purpose

nag_matop_complex_gen_matrix_fun_usd (f01fm) computes the matrix function, f(A)$f\left(A\right)$, of a complex n$n$ by n$n$ matrix A$A$, using analytical derivatives of f$f$ you have supplied.

## Syntax

[a, user, iflag, ifail] = f01fm(a, f, 'n', n, 'user', user)
[a, user, iflag, ifail] = nag_matop_complex_gen_matrix_fun_usd(a, f, 'n', n, 'user', user)

## Description

f(A)$f\left(A\right)$ is computed using the Schur–Parlett algorithm described in Higham (2008) and Davies and Higham (2003).
The scalar function f$f$, and the derivatives of f$f$, are returned by the function f which, given an integer m$m$, should evaluate f(m)(zi)${f}^{\left(m\right)}\left({z}_{\mathit{i}}\right)$ at a number of points zi${z}_{\mathit{i}}$, for i = 1,2,,nz$\mathit{i}=1,2,\dots ,{n}_{z}$, on the complex plane. nag_matop_complex_gen_matrix_fun_usd (f01fm) is therefore appropriate for functions that can be evaluated on the complex plane and whose derivatives, of arbitrary order, can also be evaluated on the complex plane.

## References

Davies P I and Higham N J (2003) A Schur–Parlett algorithm for computing matrix functions. SIAM J. Matrix Anal. Appl. 25(2) 464–485
Higham N J (2008) Functions of Matrices: Theory and Computation SIAM, Philadelphia, PA, USA

## Parameters

### Compulsory Input Parameters

1:     a(lda, : $:$) – complex 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$.
2:     f – function handle or string containing name of m-file
Given an integer m$m$, the function f evaluates f(m)(zi)${f}^{\left(m\right)}\left({z}_{i}\right)$ at a number of points zi${z}_{i}$.
[iflag, fz, user] = f(m, iflag, nz, z, user)

Input Parameters

1:     m – int64int32nag_int scalar
The order, m$m$, of the derivative required.
If m = 0${\mathbf{m}}=0$, f(zi)$f\left({z}_{i}\right)$ should be returned. For m > 0${\mathbf{m}}>0$, f(m)(zi)${f}^{\left(m\right)}\left({z}_{i}\right)$ should be returned.
2:     iflag – int64int32nag_int scalar
iflag will be zero.
3:     nz – int64int32nag_int scalar
nz${n}_{z}$, the number of function or derivative values required.
4:     z(nz) – complex array
The nz${n}_{z}$ points z1,z2,,znz${z}_{1},{z}_{2},\dots ,{z}_{{n}_{z}}$ at which the function f$f$ is to be evaluated.
5:     user – Any MATLAB object
f is called from nag_matop_complex_gen_matrix_fun_usd (f01fm) with the object supplied to nag_matop_complex_gen_matrix_fun_usd (f01fm).

Output Parameters

1:     iflag – int64int32nag_int scalar
iflag should either be unchanged from its entry value of zero, or may be set nonzero to indicate that there is a problem in evaluating the function f(z)$f\left(z\right)$; for instance f(zi)$f\left({z}_{i}\right)$ may not be defined for a particular zi${z}_{i}$. If iflag is returned as nonzero then nag_matop_complex_gen_matrix_fun_usd (f01fm) will terminate the computation, with ${\mathbf{ifail}}={\mathbf{2}}$.
2:     fz(nz) – complex array
The nz${n}_{z}$ function or derivative values. fz(i)${\mathbf{fz}}\left(\mathit{i}\right)$ should return the value f(m)(zi)${f}^{\left(m\right)}\left({z}_{\mathit{i}}\right)$, for i = 1,2,,nz$\mathit{i}=1,2,\dots ,{n}_{z}$.
3:     user – Any MATLAB object

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the array a.
n$n$, the order of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.
2:     user – Any MATLAB object
user is not used by nag_matop_complex_gen_matrix_fun_usd (f01fm), but is passed to f. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

lda iuser ruser

### Output Parameters

1:     a(lda, : $:$) – complex 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)$.
The n$n$ by n$n$ matrix, f(A)$f\left(A\right)$.
2:     user – Any MATLAB object
3:     iflag – int64int32nag_int scalar
iflag = 0${\mathbf{iflag}}=0$, unless iflag has been set nonzero inside f, in which case iflag will be the value set and ifail will be set to ${\mathbf{ifail}}={\mathbf{2}}$.
4:     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:
ifail = 1${\mathbf{ifail}}=1$
A Taylor series failed to converge.
ifail = 2${\mathbf{ifail}}=2$
iflag has been set nonzero by the user.
ifail = 3${\mathbf{ifail}}=3$
There was an error whilst reordering the Schur form of A$A$.
Note:  this failure should not occur and suggests that the function has been called incorrectly.
ifail = 4${\mathbf{ifail}}=4$
The function was unable to compute the Schur decomposition of A$A$.
Note:  this failure should not occur and suggests that the function has been called incorrectly.
ifail = 5${\mathbf{ifail}}=5$
ifail = 1${\mathbf{ifail}}=-1$
Input argument number _$_$ is invalid.
ifail = 3${\mathbf{ifail}}=-3$
On entry, parameter lda is invalid.
Constraint: ldan$\mathit{lda}\ge {\mathbf{n}}$.
ifail = 999${\mathbf{ifail}}=-999$
Allocation of memory failed.

## Accuracy

For a normal matrix A$A$ (for which AH A = AAH${A}^{\mathrm{H}}A=A{A}^{\mathrm{H}}$), the Schur decomposition is diagonal and the algorithm reduces to evaluating f$f$ at the eigenvalues of A$A$ and then constructing f(A)$f\left(A\right)$ using the Schur vectors. This should give a very accurate result. In general, however, no error bounds are available for the algorithm. See Section 9.4 of Higham (2008) for further discussion of the Schur–Parlett algorithm.

Up to 6n2$6{n}^{2}$ of complex allocatable memory is required.
The cost of the Schur–Parlett algorithm depends on the spectrum of A$A$, but is roughly between 28n3$28{n}^{3}$ and n4 / 3${n}^{4}/3$ floating point operations. There is an additional cost in evaluating f$f$ and its derivatives. If the derivatives of f$f$ are not known analytically, then nag_matop_complex_gen_matrix_fun_num (f01fl) can be used to evaluate f(A)$f\left(A\right)$ using numerical differentiation. If A$A$ is complex Hermitian then it is recommended that nag_matop_complex_herm_matrix_fun (f01ff) be used as it is more efficient and, in general, more accurate than nag_matop_complex_gen_matrix_fun_usd (f01fm).
Note that f$f$ must be analytic in the region of the complex plane containing the spectrum of A$A$.
For further information on matrix functions, see Higham (2008).
If estimates of the condition number of the matrix function are required then nag_matop_complex_gen_matrix_cond_usd (f01kc) should be used.
nag_matop_real_gen_matrix_fun_usd (f01em) can be used to find the matrix function f(A)$f\left(A\right)$ for a real matrix A$A$.

## Example

```function nag_matop_complex_gen_matrix_fun_usd_example
a = [ 1.0+0.0i, 0.0+0.0i,  1.0+0.0i, 0.0+2.0i;
0.0+1.0i, 1.0+0.0i, -1.0+0.0i, 1.0+0.0i;
-1.0+0.0i, 0.0+1.0i,  0.0+1.0i, 0.0+1.0i;
1.0+1.0i, 0.0+2.0i, -1.0+0.0i, 0.0+1.0i];
% Compute exp(3*a)
[a, user, iflag, ifail] = nag_matop_complex_gen_matrix_fun_usd(a, @f)

function [iflag, fz, user] = f(m, iflag, nz, z, user)
fz = double(3^m)*exp(3*z);
```
```

a =

-10.3264 +14.8082i  -1.4883 +74.3369i -12.1206 -47.0956i  41.5622 +32.2927i
63.3909 -40.5336i -21.0117 -62.7073i  16.5106 +35.2787i  -5.1725 +17.9413i
-6.3954 +56.4708i  25.4246 +13.8034i -14.4937 - 9.2397i -20.3167 + 2.8647i
31.4957 +23.2757i  28.6003 +21.4573i -23.8034 -11.6547i  23.9841 +18.7737i

user =

0

iflag =

0

ifail =

0

```
```function f01fm_example
a = [ 1.0+0.0i, 0.0+0.0i,  1.0+0.0i, 0.0+2.0i;
0.0+1.0i, 1.0+0.0i, -1.0+0.0i, 1.0+0.0i;
-1.0+0.0i, 0.0+1.0i,  0.0+1.0i, 0.0+1.0i;
1.0+1.0i, 0.0+2.0i, -1.0+0.0i, 0.0+1.0i];
% Compute exp(3*a)
[a, user, iflag, ifail] = f01fm(a, @f)

function [iflag, fz, user] = f(m, iflag, nz, z, user)
fz = double(3^m)*exp(3*z);
```
```

a =

-10.3264 +14.8082i  -1.4883 +74.3369i -12.1206 -47.0956i  41.5622 +32.2927i
63.3909 -40.5336i -21.0117 -62.7073i  16.5106 +35.2787i  -5.1725 +17.9413i
-6.3954 +56.4708i  25.4246 +13.8034i -14.4937 - 9.2397i -20.3167 + 2.8647i
31.4957 +23.2757i  28.6003 +21.4573i -23.8034 -11.6547i  23.9841 +18.7737i

user =

0

iflag =

0

ifail =

0

```