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_real_gen_matrix_fun_usd (f01em)

## Purpose

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

## Syntax

[a, user, iflag, imnorm, ifail] = f01em(a, f, 'n', n, 'user', user)
[a, user, iflag, imnorm, ifail] = nag_matop_real_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 (generally complex) points zi${z}_{\mathit{i}}$, for i = 1,2,,nz$\mathit{i}=1,2,\dots ,{n}_{z}$. For any z$z$ on the real line, f(z)$f\left(z\right)$ must also be real. nag_matop_real_gen_matrix_fun_usd (f01em) 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, : $:$) – double 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_real_gen_matrix_fun_usd (f01em) with the object supplied to nag_matop_real_gen_matrix_fun_usd (f01em).

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_real_gen_matrix_fun_usd (f01em) 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}$. If zi${z}_{i}$ lies on the real line, then so must f(m)(zi)${f}^{\left(m\right)}\left({z}_{i}\right)$.
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_real_gen_matrix_fun_usd (f01em), 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, : $:$) – double 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:     imnorm – double scalar
If A$A$ has complex eigenvalues, nag_matop_real_gen_matrix_fun_usd (f01em) will use complex arithmetic to compute f(A)$f\left(A\right)$. The imaginary part is discarded at the end of the computation, because it will theoretically vanish. imnorm contains the 1$1$-norm of the imaginary part, which should be used to check that the function has given a reliable answer.
If A$A$ has real eigenvalues, nag_matop_real_gen_matrix_fun_usd (f01em) uses real arithmetic and imnorm = 0${\mathbf{imnorm}}=0$.
5:     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 routine 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 AT A = AAT${A}^{\mathrm{T}}A=A{A}^{\mathrm{T}}$), 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.

If A$A$ has real eigenvalues then up to 6n2$6{n}^{2}$ of double allocatable memory may be required. If A$A$ has complex eigenvalues then up to 6n2$6{n}^{2}$ of complex allocatable memory may be 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_real_gen_matrix_fun_num (f01el) can be used to evaluate f(A)$f\left(A\right)$ using numerical differentiation. If A$A$ is real symmetric then it is recommended that nag_matop_real_symm_matrix_fun (f01ef) be used as it is more efficient and, in general, more accurate than nag_matop_real_gen_matrix_fun_usd (f01em).
For any z$z$ on the real line, f(z)$f\left(z\right)$ must be real. f$f$ must also be complex analytic on the spectrum of A$A$. These conditions ensure that f(A)$f\left(A\right)$ is real for real 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_real_gen_matrix_cond_usd (f01jc) should be used.
nag_matop_complex_gen_matrix_fun_usd (f01fm) can be used to find the matrix function f(A)$f\left(A\right)$ for a complex matrix A$A$.

## Example

```function nag_matop_real_gen_matrix_fun_usd_example
a =  [1,  0, -2,  1;
-1,  2,  0,  1;
2,  0,  1,  0;
1,  0, -1,  2];
% Compute f(a)
[a, user, iflag, imnorm, ifail] = nag_matop_real_gen_matrix_fun_usd(a, @f)

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

a =

-12.1880         0   -3.4747    8.3697
-13.7274   54.5982  -23.9801   82.8593
-9.7900         0  -25.4527   26.5294
-18.1597         0  -34.8991   49.2404

user =

0

iflag =

0

imnorm =

1.9974e-14

ifail =

0

```
```function f01em_example
a =  [1,  0, -2,  1;
-1,  2,  0,  1;
2,  0,  1,  0;
1,  0, -1,  2];
% Compute f(a)
[a, user, iflag, imnorm, ifail] = f01em(a, @f)

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

a =

-12.1880         0   -3.4747    8.3697
-13.7274   54.5982  -23.9801   82.8593
-9.7900         0  -25.4527   26.5294
-18.1597         0  -34.8991   49.2404

user =

0

iflag =

0

imnorm =

1.9974e-14

ifail =

0

```