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_linsys_real_posdef_solve_ref (f04ab)

Purpose

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 AX = BAX=B, where AA is symmetric positive definite, nag_linsys_real_posdef_solve_ref (f04ab) first computes a Cholesky factorization of AA as A = LLTA=LLT, where LL is lower triangular. An approximation to XX is found by forward and backward substitution. The residual matrix R = BAXR=B-AX is then calculated using additional precision, and a correction DD to XX is found by solving LLTD = RLLTD=R. XX is replaced by X + DX+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)max(1,n)
The second dimension of the array must be at least max (1,n)max(1,n)
The upper triangle of the nn by nn positive definite symmetric matrix AA. 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)max(1,n)
The second dimension of the array must be at least max (1,m)max(1,m)
The nn by mm right-hand side matrix BB.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the arrays a, b The second dimension of the array a.
nn, the order of the matrix AA.
Constraint: n0n0.
2:     m – int64int32nag_int scalar
Default: The second dimension of the array b.
mm, the number of right-hand sides.
Constraint: m0m0.

Input Parameters Omitted from the MATLAB Interface

lda ldb ldc wkspce ldbb

Output Parameters

1:     a(lda, : :) – double array
The first dimension of the array a will be max (1,n)max(1,n)
The second dimension of the array will be max (1,n)max(1,n)
ldamax (1,n)ldamax(1,n).
The elements of the array below the diagonal are overwritten; the upper triangle of AA is unchanged.
2:     c(ldc,mm) – double array
The first dimension of the array c will be max (1,n)max(1,n)
The second dimension of the array will be max (1,m)max(1,m)
ldcmax (1,n)ldcmax(1,n).
The nn by mm solution matrix XX.
3:     bb(ldbb,mm) – double array
The first dimension of the array bb will be max (1,n)max(1,n)
The second dimension of the array will be max (1,m)max(1,m)
ldbbmax (1,n)ldbbmax(1,n).
The final nn by mm residual matrix R = BAXR=B-AX.
4:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
The matrix AA is not positive definite, possibly due to rounding errors.
  ifail = 2ifail=2
Iterative refinement fails to improve the solution, i.e., the matrix AA is too ill-conditioned.
  ifail = 3ifail=3
On entry,n < 0n<0,
orm < 0m<0,
orlda < max (1,n)lda<max(1,n),
orldb < max (1,n)ldb<max(1,n),
orldc < max (1,n)ldc<max(1,n),
orldbb < max (1,n)ldbb<max(1,n).

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 n3n3.
If there is only one right-hand side, it is simpler to use nag_linsys_real_posdef_solve_1rhs (f04as).

Example

function nag_linsys_real_posdef_solve_ref_example
a = [5, 7, 6, 5;
     7, 10, 8, 7;
     6, 8, 10, 9;
     5, 7, 9, 10];
b = [23;
     32;
     33;
     31];
[aOut, c, bb, ifail] = nag_linsys_real_posdef_solve_ref(a, b)
 

aOut =

    5.0000    7.0000    6.0000    5.0000
    3.1305   10.0000    8.0000    7.0000
    2.6833   -0.8944   10.0000    9.0000
    2.2361         0    2.1213   10.0000


c =

     1
     1
     1
     1


bb =

     0
     0
     0
     0


ifail =

                    0


function f04ab_example
a = [5, 7, 6, 5;
     7, 10, 8, 7;
     6, 8, 10, 9;
     5, 7, 9, 10];
b = [23;
     32;
     33;
     31];
[aOut, c, bb, ifail] = f04ab(a, b)
 

aOut =

    5.0000    7.0000    6.0000    5.0000
    3.1305   10.0000    8.0000    7.0000
    2.6833   -0.8944   10.0000    9.0000
    2.2361         0    2.1213   10.0000


c =

     1
     1
     1
     1


bb =

     0
     0
     0
     0


ifail =

                    0



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