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_real_gen_matrix_pow (f01eq)

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

Purpose

nag_matop_real_gen_matrix_pow (f01eq) computes the principal real power Ap, for arbitrary p, of a real n by n matrix A.

Syntax

[a, ifail] = f01eq(a, p, 'n', n)
[a, ifail] = nag_matop_real_gen_matrix_pow(a, p, 'n', n)

Description

For a matrix A with no eigenvalues on the closed negative real line, Ap (p) can be defined as
Ap= expplogA  
where logA is the principal logarithm of A (the unique logarithm whose spectrum lies in the strip z:-π<Imz<π).
Ap is computed using the real version of the Schur–Padé algorithm described in Higham and Lin (2011) and Higham and Lin (2013).
The real number p is expressed as p=q+r where q-1,1 and r. Then Ap=AqAr. The integer power Ar is found using a combination of binary powering and, if necessary, matrix inversion. The fractional power Aq is computed, entirely in real arithmetic, using a real Schur decomposition and a Padé approximant.

References

Higham N J (2008) Functions of Matrices: Theory and Computation SIAM, Philadelphia, PA, USA
Higham N J and Lin L (2011) A Schur–Padé algorithm for fractional powers of a matrix SIAM J. Matrix Anal. Appl. 32(3) 1056–1078
Higham N J and Lin L (2013) An improved Schur–Padé algorithm for fractional powers of a matrix and their Fréchet derivatives MIMS Eprint 2013.1 Manchester Institute for Mathematical Sciences, School of Mathematics, University of Manchester http://eprints.ma.man.ac.uk/

Parameters

Compulsory Input Parameters

1:     alda: – double array
The first dimension of the array a must be at least n.
The second dimension of the array a must be at least n.
The n by n matrix A.
2:     p – double scalar
The required power of A.

Optional Input Parameters

1:     n int64int32nag_int scalar
Default: the first dimension of the array a.
n, the order of the matrix A.
Constraint: n0.

Output Parameters

1:     alda: – double array
The first dimension of the array a will be n.
The second dimension of the array a will be n.
The n by n matrix pth power, Ap.
2:     ifail int64int32nag_int scalar
ifail=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
A has eigenvalues on the negative real line. The principal pth power is not defined. nag_matop_complex_gen_matrix_pow (f01fq) can be used to find a complex, non-principal pth power.
   ifail=2
A is singular so the pth power cannot be computed.
   ifail=3
Ap has been computed using an IEEE double precision Padé approximant, although the arithmetic precision is higher than IEEE double precision.
   ifail=4
An unexpected internal error occurred. This failure should not occur and suggests that the function has been called incorrectly.
   ifail=-1
Constraint: n0.
   ifail=-3
Constraint: ldan.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   ifail=-999
Dynamic memory allocation failed.

Accuracy

For positive integer p, the algorithm reduces to a sequence of matrix multiplications. For negative integer p, the algorithm consists of a combination of matrix inversion and matrix multiplications.
For a normal matrix A (for which ATA=AAT) and non-integer p, the Schur decomposition is diagonal and the algorithm reduces to evaluating powers of the eigenvalues of A and then constructing Ap using the Schur vectors. This should give a very accurate result. In general however, no error bounds are available for the algorithm.

Further Comments

The cost of the algorithm is On3. The exact cost depends on the matrix A but if p-1,1 then the cost is independent of p. O4×n2 of real allocatable memory is required by the function.
If estimates of the condition number of Ap are required then nag_matop_real_gen_matrix_cond_pow (f01je) should be used.

Example

This example finds Ap where p=0.2 and
A = 3 3 2 1 3 1 0 2 1 1 4 3 3 0 3 1 .  
function f01eq_example


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

% Principal power p of matrix A

a = [ 3 3 2 1;
      3 1 0 2;
      1 1 4 3;
      3 0 3 1];

p = 0.2;

[pa, ifail] = f01eq(a,p);

disp('A^p:');
disp(pa);


f01eq example results

A^p:
    1.2446    0.2375    0.2172   -0.1359
    0.0925    1.1239   -0.1453    0.3731
   -0.0769    0.1972    1.3131    0.1837
    0.3985   -0.2902    0.1085    1.1560


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