nag_linsys_real_posdef_solve_ref (f04ab) calculates the accurate solution of a set of real symmetric positive definite linear equations with multiple right-hand sides, using a Cholesky factorization and iterative refinement.

Syntax

[a, c, bb, ifail] = f04ab(a, b, 'n', n, 'm', m)

[a, c, bb, ifail] = nag_linsys_real_posdef_solve_ref(a, b, 'n', n, 'm', m)

Description

Given a set of real linear equations

A X = B

, where

A

is symmetric positive definite, nag_linsys_real_posdef_solve_ref (f04ab) first computes a Cholesky factorization of

A

A = L L^{T}

, where

L

is lower triangular. An approximation to

X

is found by forward and backward substitution. The residual matrix

R = B - A X

is then calculated using additional precision, and a correction

D

X

is found by solving

L L^{T} D = R

X

is replaced by

X + D

, and this iterative refinement of the solution is repeated until full machine accuracy has been obtained.

References

Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag

Parameters

Compulsory Input Parameters

1: $a (lda :)$ – double array: The first dimension of the array a must be at least $\max (1, n)$ .
The second dimension of the array a must be at least $\max (1, n)$ .

The upper triangle of the $n$ by $n$ positive definite symmetric matrix $A$ . The elements of the array below the diagonal need not be set.
2: $b (ldb :)$ – double array: The first dimension of the array b must be at least $\max (1, n)$ .
The second dimension of the array b must be at least $\max (1, m)$ .

The $n$ by $m$ right-hand side matrix $B$ .

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the first dimension of the arrays a, b and the second dimension of the array a.
$n$ , the order of the matrix $A$ .

Constraint: $n \geq 0$ .
2: $m$ – int64int32nag_int scalar: Default: the second dimension of the array b.
$m$ , the number of right-hand sides.

Constraint: $m \geq 0$ .

Output Parameters

1: $a (lda :)$ – double array: The first dimension of the array a will be $\max (1, n)$ .
The second dimension of the array a will be $\max (1, n)$ .

The elements of the array below the diagonal are overwritten; the upper triangle of $A$ is unchanged.
2: $c (ldc, m)$ – double array: The first dimension of the array c will be $\max (1, n)$ .
The second dimension of the array c will be $\max (1, m)$ .

The $n$ by $m$ solution matrix $X$ .
3: $bb (ldbb, m)$ – double array: The first dimension of the array bb will be $\max (1, n)$ .
The second dimension of the array bb will be $\max (1, m)$ .

The final $n$ by $m$ residual matrix $R = B - A X$ .
4: $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$: The matrix $A$ is not positive definite, possibly due to rounding errors.

$ifail = 2$: Iterative refinement fails to improve the solution, i.e., the matrix $A$ is too ill-conditioned.

$ifail = 3$

On entry,	$n < 0$ ,
or	$m < 0$ ,
or	$lda < \max (1, n)$ ,
or	$ldb < \max (1, n)$ ,
or	$ldc < \max (1, n)$ ,
or	$ldbb < \max (1, n)$ .

$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

The computed solutions should be correct to full machine accuracy. For a detailed error analysis see page 39 of Wilkinson and Reinsch (1971).

Further Comments

The time taken by nag_linsys_real_posdef_solve_ref (f04ab) is approximately proportional to

n^{3}

If there is only one right-hand side, it is simpler to use nag_linsys_real_posdef_solve_1rhs (f04as).

Example

This example solves the set of linear equations

A X = B

where

A = (\begin{array}{r} 5 & 7 & 6 & 5 \\ 7 & 10 & 8 & 7 \\ 6 & 8 & 10 & 9 \\ 5 & 7 & 9 & 10 \end{array}) and B = (\begin{array}{r} 23 \\ 32 \\ 33 \\ 31 \end{array}) .

Open in the MATLAB editor: f04ab_example

function f04ab_example


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

% Solve Ax = b for positive definite A
a = [ 5,  7,  6,  5;
      7, 10,  8,  7;
      6,  8, 10,  9;
      5,  7,  9, 10];
b = [23; 32; 33; 31];

[afac, x, resid, ifail] = f04ab(a, b);

disp('Solution');
disp(x);

f04ab example results

Solution
     1
     1
     1
     1

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

Chapter Contents

Chapter Introduction

NAG Toolbox