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

## Purpose

nag_matop_real_gen_matrix_fun_usd (f01em) computes the matrix function, $f\left(A\right)$, of a real $n$ by $n$ matrix $A$, using analytical derivatives of $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\left(A\right)$ is computed using the Schur–Parlett algorithm described in Higham (2008) and Davies and Higham (2003).
The scalar function $f$, and the derivatives of $f$, are returned by the function f which, given an integer $m$, should evaluate ${f}^{\left(m\right)}\left({z}_{\mathit{i}}\right)$ at a number of (generally complex) points ${z}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{n}_{z}$. For any $z$ on the real line, $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:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double array
The first dimension of the array a must be at least ${\mathbf{n}}$.
The second dimension of the array a must be at least ${\mathbf{n}}$.
The $n$ by $n$ matrix $A$.
2:     $\mathrm{f}$ – function handle or string containing name of m-file
Given an integer $m$, the function f evaluates ${f}^{\left(m\right)}\left({z}_{i}\right)$ at a number of points ${z}_{i}$.
[iflag, fz, user] = f(m, iflag, nz, z, user)

Input Parameters

1:     $\mathrm{m}$int64int32nag_int scalar
The order, $m$, of the derivative required.
If ${\mathbf{m}}=0$, $f\left({z}_{i}\right)$ should be returned. For ${\mathbf{m}}>0$, ${f}^{\left(m\right)}\left({z}_{i}\right)$ should be returned.
2:     $\mathrm{iflag}$int64int32nag_int scalar
iflag will be zero.
3:     $\mathrm{nz}$int64int32nag_int scalar
${n}_{z}$, the number of function or derivative values required.
4:     $\mathrm{z}\left({\mathbf{nz}}\right)$ – complex array
The ${n}_{z}$ points ${z}_{1},{z}_{2},\dots ,{z}_{{n}_{z}}$ at which the function $f$ is to be evaluated.
5:     $\mathrm{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:     $\mathrm{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\left(z\right)$; for instance $f\left({z}_{i}\right)$ may not be defined for a particular ${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:     $\mathrm{fz}\left({\mathbf{nz}}\right)$ – complex array
The ${n}_{z}$ function or derivative values. ${\mathbf{fz}}\left(\mathit{i}\right)$ should return the value ${f}^{\left(m\right)}\left({z}_{\mathit{i}}\right)$, for $\mathit{i}=1,2,\dots ,{n}_{z}$. If ${z}_{i}$ lies on the real line, then so must ${f}^{\left(m\right)}\left({z}_{i}\right)$.
3:     $\mathrm{user}$ – Any MATLAB object

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the array a.
$n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
2:     $\mathrm{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.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double array
The first dimension of the array a will be ${\mathbf{n}}$.
The second dimension of the array a will be ${\mathbf{n}}$.
The $n$ by $n$ matrix, $f\left(A\right)$.
2:     $\mathrm{user}$ – Any MATLAB object
3:     $\mathrm{iflag}$int64int32nag_int scalar
${\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:     $\mathrm{imnorm}$ – double scalar
If $A$ has complex eigenvalues, nag_matop_real_gen_matrix_fun_usd (f01em) will use complex arithmetic to compute $f\left(A\right)$. The imaginary part is discarded at the end of the computation, because it will theoretically vanish. imnorm contains the $1$-norm of the imaginary part, which should be used to check that the function has given a reliable answer.
If $A$ has real eigenvalues, nag_matop_real_gen_matrix_fun_usd (f01em) uses real arithmetic and ${\mathbf{imnorm}}=0$.
5:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\mathbf{ifail}}=1$
A Taylor series failed to converge.
${\mathbf{ifail}}=2$
iflag has been set nonzero by the user.
${\mathbf{ifail}}=3$
There was an error whilst reordering the Schur form of $A$.
Note:  this failure should not occur and suggests that the function has been called incorrectly.
${\mathbf{ifail}}=4$
The routine was unable to compute the Schur decomposition of $A$.
Note:  this failure should not occur and suggests that the function has been called incorrectly.
${\mathbf{ifail}}=5$
An unexpected internal error occurred. Please contact NAG.
${\mathbf{ifail}}=-1$
Input argument number $_$ is invalid.
${\mathbf{ifail}}=-3$
On entry, argument lda is invalid.
Constraint: $\mathit{lda}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-99$
An unexpected error has been triggered by this routine. Please contact NAG.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

For a normal matrix $A$ (for which ${A}^{\mathrm{T}}A=A{A}^{\mathrm{T}}$), the Schur decomposition is diagonal and the algorithm reduces to evaluating $f$ at the eigenvalues of $A$ and then constructing $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.

## Further Comments

If $A$ has real eigenvalues then up to $6{n}^{2}$ of double allocatable memory may be required. If $A$ has complex eigenvalues then up to $6{n}^{2}$ of complex allocatable memory may be required.
The cost of the Schur–Parlett algorithm depends on the spectrum of $A$, but is roughly between $28{n}^{3}$ and ${n}^{4}/3$ floating-point operations. There is an additional cost in evaluating $f$ and its derivatives. If the derivatives of $f$ are not known analytically, then nag_matop_real_gen_matrix_fun_num (f01el) can be used to evaluate $f\left(A\right)$ using numerical differentiation. If $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$ on the real line, $f\left(z\right)$ must be real. $f$ must also be complex analytic on the spectrum of $A$. These conditions ensure that $f\left(A\right)$ is real for real $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\left(A\right)$ for a complex matrix $A$.

## Example

This example finds the ${e}^{2A}$ where
 $A= 1 0 -2 1 -1 2 0 1 2 0 1 0 1 0 -1 2 .$
```function f01em_example

fprintf('f01em example results\n\n');

a =  [1,  0, -2,  1;
-1,  2,  0,  1;
2,  0,  1,  0;
1,  0, -1,  2];

% Compute f(a)
[exp2a, user, iflag, imnorm, ifail] = ...
f01em(a, @f);

disp('f(A) = exp(2A)');
disp(exp2a);

fprintf('Imnorm = %6.2f\n',imnorm);

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

f(A) = exp(2A)
-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

Imnorm =   0.00
```

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–2015