Guide and Reference

SPPF, DPPF, SPOF, DPOF, CPOF, and ZPOF--Positive Definite Real Symmetric or Complex Hermitian Matrix Factorization

The SPPF and DPPF subroutines factor positive definite symmetric matrix A, stored in lower-packed storage mode, using Gaussian elimination (LDL^T) or the Cholesky factorization method. To solve a system of equations with one or more right-hand sides, follow the call to these subroutines with one or more calls to SPPS or DPPS, respectively. To find the inverse of matrix A, follow the call to these subroutines, performing Cholesky factorization, with a call to SPPICD or DPPICD, respectively.

The SPOF, DPOF, CPOF, and ZPOF subroutines factor matrix A stored in upper or lower storage mode, where:

• For SPOF and DPOF, A is a positive definite symmetric matrix.
• For CPOF and ZPOF, A is a positive definite complex Hermitian matrix.

Matrix A is factored using Cholesky factorization, (LL^T or U^TU for SPOF and DPOF and LL^H or U^HU for CPOF and ZPOF). To solve the system of equations with one or more right-hand sides, follow the call to these subroutines with a call to SPOSM, DPOSM, CPOSM, or ZPOSM. To find the inverse of matrix A, follow the call to SPOF or DPOF with a call to SPOICD or DPOICD.

Table 93. Data Types

A	Subroutine
Short-precision real	SPPF and SPOF
Long-precision real	DPPF and DPOF
Short-precision complex	CPOF
Long-precision complex	ZPOF

Note:

The output from SPPF and DPPF should be used only as input to the following subroutines for performing a solve or inverse: SPPS/SPPICD and DPPS/DPPICD, respectively. The output from SPOF, DPOF, CPOF, and ZPOF should be used only as input to the following subroutines for performing a solve or inverse: SPOSM/SPOICD, DPOSM/DPOICD, CPOSM, and ZPOSM, respectively.

Syntax:  

Fortran	CALL SPPF | DPPF (ap, n, iopt) CALL SPOF | DPOF | CPOF | ZPOF (uplo, a, lda, n)
C and C++	sppf | dppf (ap, n, iopt); spof | dpof | cpof | zpof (uplo, a, lda, n);
PL/I	CALL SPPF | DPPF (ap, n, iopt); CALL SPOF | DPOF | CPOF | ZPOF (uplo, a, lda, n);

On Entry:  

uplo

indicates whether matrix A is stored in upper or lower storage mode, where:

If uplo = 'U', A is stored in upper storage mode.

If uplo = 'L', A is stored in lower storage mode.

Specified as: a single character. It must be 'U' or 'L'.

ap

is array, referred to as AP, in which matrix A, to be factored, is stored in lower-packed storage mode.

Specified as: a one-dimensional array, containing numbers of the data type indicated in Table 93. See "Notes".

If iopt = 0, the array must have at least n(n+1)/2+n elements.

If iopt = 1, the array must have at least n(n+1)/2 elements.

a

is the positive definite matrix A, to be factored.

Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 93.

lda

is the leading dimension of the array specified for a.

Specified as: a fullword integer; lda > 0 and lda >= n.

n

is the order n of matrix A.

Specified as: a fullword integer; n >= 0.

iopt

determines the type of computation to be performed, where:

If iopt = 0, the matrix is factored using the LDL^T method.

If iopt = 1, the matrix is factored using Cholesky factorization.

Specified as: a fullword integer; iopt = 0 or 1.

On Return:  

ap

is the transformed matrix A of order n, containing the results of the factorization. See "Notes" and "Function".

Returned as: a one-dimensional array, containing numbers of the data type indicated in Table 93.

If iopt = 0, the array contains n(n+1)/2+n elements.

If iopt = 1, the array contains n(n+1)/2 elements.

a

is the transformed matrix A of order n, containing the results of the factorization. See "Function".

Returned as: a two-dimensional array, containing numbers of the data type indicated in Table 93.

Notes:  

1. All subroutines accept lowercase letters for the uplo argument.

2. In the input and output arrays specified for ap, the first n(n+1)/2 elements are matrix elements. The additional n locations, required in the array when iopt = 0, are used for working storage by this subroutine and should not be altered between calls to the factorization and solve subroutines.

3. On input, the imaginary parts of the diagonal elements of the complex Hermitian matrix A are assumed to be zero, so you do not have to set these values. On output, they are set to zero.

4. SPPF and DPPF in many cases utilize new algorithms based on recursive packed storage format. As a result, on output, the array specified for AP may be stored in this new format rather than the conventional lower packed format. (See references [52], [66], and [67]).

   The array specified for AP should not be altered between calls to the factorization and solve subroutines; otherwise unpredictable results may occur.

5. For a description of the storage modes used for the matrices, see:

   • For positive definite symmetric matrices, see "Positive Definite or Negative Definite Symmetric Matrix".

   • For positive definite complex Hermitian matrices, see "Positive Definite or Negative Definite Complex Hermitian Matrix".

Function:   The functions for these subroutines are described in the sections below.

For SPPF and DPPF:   If iopt = 0, the positive definite symmetric matrix A, stored in lower-packed storage mode, is factored using Gaussian elimination, where A is expressed as:

A = LDL^T

where:

L is a unit lower triangular matrix.

L^T is the transpose of matrix L.

D is a diagonal matrix.

If iopt = 1, the positive definite symmetric matrix A is factored using Cholesky factorization, where A is expressed as:

A = LL^T

where L is a lower triangular matrix.

If n is 0, no computation is performed. See references [36] and [38].

For SPOF, DPOF, CPOF, and ZPOF:   The positive definite matrix A, stored in upper or lower storage mode, is factored using Cholesky factorization, where A is expressed as:

A = LL^T or A = U^TU    for SPOF and DPOF

A = LL^H or A = U^HU    for CPOF and ZPOF

where:

L is a lower triangular matrix.

L^T is the transpose of matrix L.

L^H is the conjugate transpose of matrix L.

U is an upper triangular matrix.

U^T is the transpose of matrix U.

U^H is the conjugate transpose of matrix U.

If n is 0, no computation is performed. See references [8], [70], and [36].

Error Conditions:  

Resource Errors:   Unable to allocate internal work area.

Computational Errors:  

1. Matrix A is not positive definite (for SPPF and DPPF when iopt = 0).
   • Processing continues to the end of the matrix.
   • One or more elements of D contain values less than or equal to 0; all elements of D are checked. The index i of the last nonpositive element encountered is identified in the computational error message.
   • The return code is set to 1.
   • i can be determined at run time by use of the ESSL error-handling facilities. To obtain this information, you must use ERRSET to change the number of allowable errors for error code 2104 in the ESSL error option table; otherwise, the default value causes your program to terminate when this error occurs. For details, see "What Can You Do about ESSL Computational Errors?".
2. Matrix A is not positive definite (for SPPF and DPPF when iopt = 1 and for SPOF, DPOF, CPOF, and ZPOF).
   • Processing stops at the first occurrence of a nonpositive definite diagonal element.
   • The order i of the first minor encountered having a nonpositive determinant is identified in the computational error message.
   • The return code is set to 1.
   • i can be determined at run time by use of the ESSL error-handling facilities. To obtain this information, you must use ERRSET to change the number of allowable errors for error code 2115 in the ESSL error option table; otherwise, the default value causes your program to terminate when this error occurs. For details, see "What Can You Do about ESSL Computational Errors?".

Input-Argument Errors:  

1. n < 0
2. iopt <> 0 or 1
3. uplo <> 'U' or 'L'
4. lda <= 0
5. n > lda

Example 1:   This example shows a factorization of positive definite symmetric matrix A of order 9, stored in lower-packed storage mode, where on input matrix A is:

             *                                             *
             | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
             | 1.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0 |
             | 1.0  2.0  3.0  3.0  3.0  3.0  3.0  3.0  3.0 |
             | 1.0  2.0  3.0  4.0  4.0  4.0  4.0  4.0  4.0 |
             | 1.0  2.0  3.0  4.0  5.0  5.0  5.0  5.0  5.0 |
             | 1.0  2.0  3.0  4.0  5.0  6.0  6.0  6.0  6.0 |
             | 1.0  2.0  3.0  4.0  5.0  6.0  7.0  7.0  7.0 |
             | 1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  8.0 |
             | 1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0 |
             *                                             *

On output, all elements of this matrix A are 1.0.

Note:

The AP arrays are formatted in a triangular arrangement for readability; however, they are stored in lower-packed storage mode.

Call Statement and Input:  

           AP  N  IOPT
           |   |   |
CALL SPPF( AP, 9,  0 )

AP = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
      3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0,
      4.0, 4.0, 4.0, 4.0, 4.0, 4.0,
      5.0, 5.0, 5.0, 5.0, 5.0,
      6.0, 6.0, 6.0, 6.0,
      7.0, 7.0, 7.0,
      8.0, 8.0,
      9.0,
       . ,  . ,  . ,  . ,  . ,  . ,  . ,  . ,  .  )

Output:  

AP = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0,
      1.0, 1.0,
      1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)

Example 2:   This example shows a factorization of the same positive definite symmetric matrix A of order 9 used in Example 1, stored in lower-packed storage mode.

Note:

The AP arrays are formatted in a triangular arrangement for readability; however, they are stored in lower-packed storage mode.

Call Statement and Input:  

           AP  N  IOPT
           |   |   |
CALL SPPF( AP, 9,  1 )

AP = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
      3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0,
      4.0, 4.0, 4.0, 4.0, 4.0, 4.0,
      5.0, 5.0, 5.0, 5.0, 5.0,
      6.0, 6.0, 6.0, 6.0,
      7.0, 7.0, 7.0,
      8.0, 8.0,
      9.0)

Output:  

AP = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0,
      1.0, 1.0,
      1.0)

Example 3:   This example shows a factorization of the same positive definite symmetric matrix A of order 9 used in Example 1, but stored in lower storage mode.

Call Statement and Input:  

           UPLO  A  LDA  N
            |    |   |   |
CALL SPOF( 'L' , A , 9 , 9 )

        *                                             *
        | 1.0   .    .    .    .    .    .    .    .  |
        | 1.0  2.0   .    .    .    .    .    .    .  |
        | 1.0  2.0  3.0   .    .    .    .    .    .  |
        | 1.0  2.0  3.0  4.0   .    .    .    .    .  |
A    =  | 1.0  2.0  3.0  4.0  5.0   .    .    .    .  |
        | 1.0  2.0  3.0  4.0  5.0  6.0   .    .    .  |
        | 1.0  2.0  3.0  4.0  5.0  6.0  7.0   .    .  |
        | 1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0   .  |
        | 1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0 |
        *                                             *

Output:  

        *                                             *
        | 1.0   .    .    .    .    .    .    .    .  |
        | 1.0  1.0   .    .    .    .    .    .    .  |
        | 1.0  1.0  1.0   .    .    .    .    .    .  |
        | 1.0  1.0  1.0  1.0   .    .    .    .    .  |
A    =  | 1.0  1.0  1.0  1.0  1.0   .    .    .    .  |
        | 1.0  1.0  1.0  1.0  1.0  1.0   .    .    .  |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0   .    .  |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0   .  |
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        *                                             *

Example 4:   This example shows a factorization of the same positive definite symmetric matrix A of order 9 used in Example 1, but stored in upper storage mode.

Call Statement and Input:  

           UPLO  A  LDA  N
            |    |   |   |
CALL SPOF( 'U' , A , 9 , 9 )

        *                                             *
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        |  .   2.0  2.0  2.0  2.0  2.0  2.0  2.0  2.0 |
        |  .    .   3.0  3.0  3.0  3.0  3.0  3.0  3.0 |
        |  .    .    .   4.0  4.0  4.0  4.0  4.0  4.0 |
A    =  |  .    .    .    .   5.0  5.0  5.0  5.0  5.0 |
        |  .    .    .    .    .   6.0  6.0  6.0  6.0 |
        |  .    .    .    .    .    .   7.0  7.0  7.0 |
        |  .    .    .    .    .    .    .   8.0  8.0 |
        |  .    .    .    .    .    .    .    .   9.0 |
        *                                             *

Output:  

        *                                             *
        | 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        |  .   1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        |  .    .   1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
        |  .    .    .   1.0  1.0  1.0  1.0  1.0  1.0 |
A    =  |  .    .    .    .   1.0  1.0  1.0  1.0  1.0 |
        |  .    .    .    .    .   1.0  1.0  1.0  1.0 |
        |  .    .    .    .    .    .   1.0  1.0  1.0 |
        |  .    .    .    .    .    .    .   1.0  1.0 |
        |  .    .    .    .    .    .    .    .   1.0 |
        *                                             *

Example 5:   This example shows a factorization of positive definite complex Hermitian matrix A of order 3, stored in lower storage mode, where on input matrix A is:

              *                                         *
              |  (25.0, 0.0)  (-5.0, -5.0)  (10.0, 5.0) |
              |  (-5.0, 5.0)   (51.0, 0.0)  (4.0, -6.0) |
              | (10.0, -5.0)    (4.0, 6.0)  (71.0, 0.0) |
              *                                         *

Note:

On input, the imaginary parts of the diagonal elements of the complex Hermitian matrix A are assumed to be zero, so you do not have to set these values. On output, they are set to zero.

Call Statement and Input:  

           UPLO  A  LDA  N
            |    |   |   |
CALL CPOF( 'L' , A , 3 , 3 )

        *                                     *
        |   (25.0, . )      .           .     |
A    =  |  (-5.0, 5.0) (51.0,  . )      .     |
        | (10.0, -5.0)  (4.0, 6.0) (71.0, . ) |
        *                                     *

Output:  

        *                                     *
        |  (5.0, 0.0)      .           .      |
A    =  | (-1.0, 1.0)  (7.0, 0.0)      .      |
        | (2.0, -1.0)  (1.0, 1.0)  (8.0, 0.0) |
        *                                     *

Example 6:   This example shows a factorization of positive definite complex Hermitian matrix A of order 3, stored in upper storage mode, where on input matrix A is:

               *                                     *
               |  (9.0, 0.0)  (3.0, 3.0) (3.0, -3.0) |
               | (3.0, -3.0) (18.0, 0.0) (8.0, -6.0) |
               |  (3.0, 3.0)  (8.0, 6.0) (43.0, 0.0) |
               *                                     *

Note:

On input, the imaginary parts of the diagonal elements of the complex Hermitian matrix A are assumed to be zero, so you do not have to set these values. On output, they are set to zero.

Call Statement and Input:  

           UPLO  A  LDA  N
            |    |   |   |
CALL CPOF( 'U' , A , 3 , 3 )

        *                                     *
        | (9.0,  . )   (3.0,3.0)   (3.0,-3.0) |
A    =  |     .       (18.0, . )   (8.0,-6.0) |
        |     .            .      (43.0,  . ) |
        *                                     *

Output:  

        *                                       *
        | (3.0, 0.0)    (1.0, 1.0)  (1.0, -1.0) |
A    =  |     .         (4.0, 0.0)  (2.0, -1.0) |
        |     .            .         (6.0, 0.0) |
        *                                       *