nag_lapack_dpftri (f07wj) computes the inverse of a real symmetric positive definite matrix using the Cholesky factorization computed by nag_lapack_dpftrf (f07wd) stored in Rectangular Full Packed (RFP) format.

Syntax

[ar, info] = f07wj(transr, uplo, n, ar)

[ar, info] = nag_lapack_dpftri(transr, uplo, n, ar)

Description

nag_lapack_dpftri (f07wj) is used to compute the inverse of a real symmetric positive definite matrix

A

, stored in RFP format. The RFP storage format is described in Rectangular Full Packed (RFP) Storage in the F07 Chapter Introduction. The function must be preceded by a call to nag_lapack_dpftrf (f07wd), which computes the Cholesky factorization of

A

uplo ='U'

A = U^{T} U

and

A^{- 1}

is computed by first inverting

U

and then forming

(U^{- 1}) U^{- T}

uplo ='L'

A = L L^{T}

and

A^{- 1}

is computed by first inverting

L

and then forming

L^{- T} (L^{- 1})

References

Du Croz J J and Higham N J (1992) Stability of methods for matrix inversion IMA J. Numer. Anal. 12 1–19

Gustavson F G, Waśniewski J, Dongarra J J and Langou J (2010) Rectangular full packed format for Cholesky's algorithm: factorization, solution, and inversion ACM Trans. Math. Software 37, 2

Parameters

Compulsory Input Parameters

1: $transr$ – string (length ≥ 1)

Specifies whether the RFP representation of

A

is normal or transposed.

$transr ='N'$: The matrix $A$ is stored in normal RFP format.
$transr ='T'$: The matrix $A$ is stored in transposed RFP format.

Constraint:

transr ='N'

'T'

2: $uplo$ – string (length ≥ 1)

Specifies how

A

has been factorized.

$uplo ='U'$: $A = U^{T} U$ , where $U$ is upper triangular.
$uplo ='L'$: $A = L L^{T}$ , where $L$ is lower triangular.

Constraint:

uplo ='U'

'L'

3: $n$ – int64int32nag_int scalar

n

, the order of the matrix

A

Constraint:

n \geq 0

4: $ar (n \times (n + 1) / 2)$ – double array

The Cholesky factorization of

A

stored in RFP format, as returned by nag_lapack_dpftrf (f07wd).

Optional Input Parameters

None.

Output Parameters

1: $ar (n \times (n + 1) / 2)$ – double array: The factorization stores the $n$ by $n$ matrix $A^{- 1}$ stored in RFP format.
2: $info$ – int64int32nag_int scalar: $info = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

$info < 0$: If $info = - i$ , argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.

W $info > 0$: The leading minor of order $_$ is not positive definite and the factorization could not be completed. Hence $A$ itself is not positive definite. This may indicate an error in forming the matrix $A$ . There is no function specifically designed to invert a symmetric matrix stored in RFP format which is not positive definite; the matrix must be treated as a full symmetric matrix, by calling nag_lapack_dsytri (f07mj).

Accuracy

The computed inverse

X

satisfies

{‖X A - I‖}_{2} \leq c (n) ε κ_{2} (A) and {‖A X - I‖}_{2} \leq c (n) ε κ_{2} (A),

where

c (n)

is a modest function of

n

ε

is the machine precision and

κ_{2} (A)

is the condition number of

A

defined by

κ_{2} (A) = {‖A‖}_{2} {‖A^{- 1}‖}_{2} .

Further Comments

The total number of floating-point operations is approximately

\frac{2}{3} n^{3}

The complex analogue of this function is nag_lapack_zpftri (f07ww).

Example

This example computes the inverse of the matrix

A

, where

A = (\begin{matrix} 4.16 & - 3.12 & 0.56 & - 0.10 \\ - 3.12 & 5.03 & - 0.83 & 1.18 \\ 0.56 & - 0.83 & 0.76 & 0.34 \\ - 0.10 & 1.18 & 0.34 & 1.18 \end{matrix}) .

Here

A

is symmetric positive definite, stored in RFP format, and must first be factorized by nag_lapack_dpftrf (f07wd).

Open in the MATLAB editor: f07wj_example

function f07wj_example


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

% Symmetric matrix in RFP format
transr = 'n';
uplo   = 'l';
ar = [ 0.76   0.34; 
       4.16   1.18;
      -3.12   5.03;
       0.56  -0.83;
      -0.10   1.18];
n  = int64(4);
n2 = (n*(n+1))/2;
ar  = reshape(ar,[n2,1]);

% Factorize a
[ar, info] = f07wd(transr, uplo, n, ar);

if info == 0
  % Compute inverse of a
  [ar, info] = f07wj(transr, uplo, n, ar);
  % Convert inverse to full array form, and print it
  [a, info] = f01vg(transr, uplo, n, ar);
  fprintf('\n');
  [ifail] = x04ca(uplo, 'n', a, 'Inverse');
else
  fprintf('\na is not positive definite.\n');
end

f07wj example results


 Inverse
             1          2          3          4
 1      0.6995
 2      0.7769     1.4239
 3      0.7508     1.8255     4.0688
 4     -0.9340    -1.8841    -2.9342     3.4978

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox