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_lapack_zpstrf (f07kr)

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

Purpose

nag_lapack_zpstrf (f07kr) computes the Cholesky factorization with complete pivoting of a complex Hermitian positive semidefinite matrix.

Syntax

[a, piv, rank, info] = f07kr(uplo, a, 'n', n, 'tol', tol)
[a, piv, rank, info] = nag_lapack_zpstrf(uplo, a, 'n', n, 'tol', tol)

Description

nag_lapack_zpstrf (f07kr) forms the Cholesky factorization of a complex Hermitian positive semidefinite matrix A either as PTAP=UHU if uplo='U' or PTAP=LLH if uplo='L', where P is a permutation matrix, U is an upper triangular matrix and L is lower triangular.
This algorithm does not attempt to check that A is positive semidefinite.

References

Higham N J (2002) Accuracy and Stability of Numerical Algorithms (2nd Edition) SIAM, Philadelphia
Lucas C (2004) LAPACK-style codes for Level 2 and 3 pivoted Cholesky factorizations LAPACK Working Note No. 161. Technical Report CS-04-522 Department of Computer Science, University of Tennessee, 107 Ayres Hall, Knoxville, TN 37996-1301, USA http://www.netlib.org/lapack/lawnspdf/lawn161.pdf

Parameters

Compulsory Input Parameters

1:     uplo – string (length ≥ 1)
Specifies whether the upper or lower triangular part of A is stored and how A is to be factorized.
uplo='U'
The upper triangular part of A is stored and A is factorized as UHU, where U is upper triangular.
uplo='L'
The lower triangular part of A is stored and A is factorized as LLH, where L is lower triangular.
Constraint: uplo='U' or 'L'.
2:     alda: – complex array
The first dimension of the array a must be at least max1,n.
The second dimension of the array a must be at least max1,n.
The n by n Hermitian positive semidefinite matrix A.
  • If uplo='U', the upper triangular part of a must be stored and the elements of the array below the diagonal are not referenced.
  • If uplo='L', the lower triangular part of a must be stored and the elements of the array above the diagonal are not referenced.

Optional Input Parameters

1:     n int64int32nag_int scalar
Default: the first dimension of the array a and the second dimension of the array a.
n, the order of the matrix A.
Constraint: n0.
2:     tol – double scalar
Default: -1
User defined tolerance. If tol<0, then n×maxk=1,nAkk×machine precision will be used. The algorithm terminates at the rth step if the r+1th step pivot <tol.

Output Parameters

1:     alda: – complex array
The first dimension of the array a will be max1,n.
The second dimension of the array a will be max1,n.
If uplo='U', the first rank rows of the upper triangle of A are overwritten with the nonzero elements of the Cholesky factor U, and the remaining rows of the triangle are destroyed.
If uplo='L', the first rank columns of the lower triangle of A are overwritten with the nonzero elements of the Cholesky factor L, and the remaining columns of the triangle are destroyed.
2:     pivn int64int32nag_int array
piv is such that the nonzero entries of P are P pivk,k =1, for k=1,2,,n.
3:     rank int64int32nag_int scalar
The computed rank of A given by the number of steps the algorithm completed.
4:     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=1
The matrix A is not positive definite. It is either positive semidefinite with computed rank as returned in rank and less than n, or it may be indefinite, see Further Comments.

Accuracy

If uplo='L' and rank=r, the computed Cholesky factor L and permutation matrix P satisfy the following upper bound
A - PLLHPT 2 A2 2r cr ε W 2 + 1 2 + Oε2 ,  
where
W = L 11 -1 L12 ,   L = L11 0 L12 0 ,   L11 r×r ,  
cr is a modest linear function of r, ε is machine precision, and
W2 13 n-r 4r-1 .  
So there is no guarantee of stability of the algorithm for large n and r, although W2 is generally small in practice.

Further Comments

The total number of real floating-point operations is approximately 4nr2-8/3r3, where r is the computed rank of A.
This algorithm does not attempt to check that A is positive semidefinite, and in particular the rank detection criterion in the algorithm is based on A being positive semidefinite. If there is doubt over semidefiniteness then you should use the indefinite factorization nag_lapack_zhetrf (f07mr). See Lucas (2004) for further information.
The real analogue of this function is nag_lapack_dpstrf (f07kd).

Example

This example computes the Cholesky factorization of the matrix A, where
A= 12.40+0.00i 2.39+0.00i 5.50+0.05i 4.47+0.00i 11.89+0.00i 2.39+0.00i 1.63+0.00i 1.04+0.10i 1.14+0.00i 1.81+0.00i 5.50+0.05i 1.04+0.10i 2.45+0.00i 1.98-0.03i 5.28-0.02i 4.47+0.00i 1.14+0.00i 1.98-0.03i 1.71+0.00i 4.14+0.00i 11.89+0.00i 1.81+0.00i 5.28-0.02i 4.14+0.00i 11.63+0.00i .  
function f07kr_example


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

% Semidefinite matrix A
a = [12.40,       2.39,       5.50+0.05i, 4.47,       11.89;
      2.39,       1.63,       1.04+0.10i, 1.14,        1.81;
      5.50+0.05i, 1.04+0.10i, 2.45,       1.98-0.03i,  5.28-0.02i;
      4.47,       1.14,       1.98-0.03i, 1.71,        4.14;
     11.89,       1.81,       5.28-0.02i, 4.14,       11.63];

% Catch warnings about rank defficient matrix ifail=1
wstat = warning();
warning('OFF');

% Factorize a
uplo = 'l';
[afac, piv, rank, info] = f07kr( ...
                                 uplo, a);

fprintf('\nComputed rank: %d\n\n', rank);
% Zero out columns rank+1 onwards and display
afac(:, rank+1:5) = 0;

ncol   = int64(80);
indent = int64(0);
[ifail] = x04db( ...
                 uplo, 'non-unit', afac, 'bracketed', 'f5.2', 'Factor', ...
                 'Integer', 'Integer', ncol, indent);

fprintf('\n piv:\n   ');
fprintf('%14d', piv);
fprintf('\n');

warning(wstat);


f07kr example results


Computed rank: 3

 Factor
                1             2             3             4             5
 1  ( 3.52, 0.00)
 2  ( 0.68, 0.00) ( 1.08, 0.00)
 3  ( 1.27, 0.00) ( 0.26, 0.00) ( 0.18, 0.00)
 4  ( 1.56, 0.01) (-0.02, 0.08) ( 0.01,-0.05) ( 0.00, 0.00)
 5  ( 3.38, 0.00) (-0.45, 0.00) (-0.17, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)

 piv:
                1             2             4             3             5

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