Engineering and Scientific Subroutine Library for AIX Version 3 Release 3: Guide and Reference

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

|These subroutines factor matrix A as explained |below: |

|SPPF and DPPF

|The SPPF and DPPF subroutines factor positive definite |real 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.

|SPOF, DPOF, CPOF, ZPOF, SPOTRF, DPOTRF, CPOTRF, and ZPOTRF

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

• |For SPOF, DPOF, SPOTRF, and DPOTRF, A is a positive definite |real symmetric matrix.
• |For CPOF, ZPOF, CPOTRF, and ZPOTRF, A is a positive definite |complex Hermitian matrix. |

|Matrix A is factored using Cholesky factorization as |follows: |

• |For SPOF, DPOF, SPOTRF, and DPOTRF, using LL^T or |U^TU.
• |For CPOF, ZPOF, CPOTRF, and ZPOTRF, using LL^H or |U^HU. |

|To solve the system of equations with one or more right-hand sides, follow |the call to SPOF, DPOF, CPOF, ZPOF, SPOTRF, DPOTRF, CPOTRF, or ZPOTRF with a |call to SPOSM, DPOSM, CPOSM, ZPOSM, SPOTRS, DPOTRS, CPOTRS, or ZPOTRS, |respectively.

|To find the inverse of matrix A, follow the call to SPOF, DPOF, |SPOTRF, or DPOTRF with a call to SPOICD, DPOICD, SPOTRI, or DPOTRI. |

Table 93. Data Types

Note:

|The output from each of these subroutines should be used only as |input for specific other subroutines, as shown in the table |below.
|

Output from this subroutine:	Should be used only as input to the following subroutine(s) for performing a solve or inverse:
SPPF	SPPS, SPPICD
DPPF	DPPS, DPPICD
SPOF	SPOSM, SPOICD
DPOF	DPOSM, DPOICD
CPOF	CPOSM
ZPOF	ZPOSM
SPOTRF	SPOTRS, SPOTRI
DPOTRF	DPOTRS, DPOTRI
CPOTRF	CPOTRS
ZPOTRF	ZPOTRS

|

Syntax

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 |or 10, the array must have at least n(n+1)/2+n elements.

If iopt = 1 |or 11, 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, |and the output is stored in an internal format.

If iopt = 1, the matrix is factored using Cholesky factorization, |and the output is stored in an internal format.

|If iopt = 10, the matrix is factored using the |LDL^T method, and the output is stored in lower-packed |storage mode.

|If iopt = 11, the matrix is factored using Cholesky |factorization, and the output is stored in lower-packed storage |mode.

Specified as: a fullword integer; iopt = |0, 1, 10, or 11.

| info

|See On Return.

On Return

ap

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

|If iopt is 0 or 1, the transformed matrix is stored in an |internal format and should only be used as input to the corresponding solve or |inverse subroutine.

|If iopt is 10 or 11, the transformed matrix is stored in |lower-packed storage mode.

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

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

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

|See Notes and see Function.

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.

| info

|has the following meaning:

|If info = 0, the factorization completed |successfully.

|If info > 0, info is set equal to the order |i of the first minor encountered having a nonpositive |determinant.

|Specified as: a fullword integer; |info >= 0.

Notes

1. |In your C program, argument info must be passed by |reference.
2. All subroutines accept lowercase letters for the uplo argument.
3. 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 |or 10, are used for working storage by this subroutine and should not be altered between calls to the factorization and solve subroutines.
4. 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.
5. |If iopt is 0 or 1, SPPF and DPPF in some cases utilize |algorithms based on recursive packed storage format. As a result, on output, |if iopt is 0 or 1, 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.

6. |The way _POTRF subroutines handle computational errors differs from |LAPACK. These subroutines use the info argument to provide |information about the computational error, like LAPACK, but also provide an |error message.
7. |On both input and output, matrix A conforms to LAPACK |format.
8. 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 |or 10, 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 |or 11, 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, ZPOF, SPOTRF, DPOTRF, CPOTRF, and ZPOTRF

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, DPOF, SPOTRF, and DPOTRF

|A = LL^H or |A = U^HU    for |CPOF, ZPOF, CPOTRF, and ZPOTRF

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 |or 10).
   • 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 |or 11 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?.
3. |Matrix A is not positive definite (for SPOTRF, DPOTRF, |CPOTRF, and ZPOTRF). |
   • |The order i of the first minor encountered having a |nonpositive determinant is identified in the computational error |message.
   • |The computational error message may occur multiple times with processing |continuing after each error, because the default for the number of allowable |errors for error code 2148 is set to be unlimited in the ESSL error option |table. |

Input-Argument Errors

1. n < 0
2. |iopt <> 0, 1, 10, or 11
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 )
 

|or

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

        *                                             *
        | 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 |
        *                                             *

|INFO  =   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 )

|or

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

        *                                             *
        | 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 |
        *                                             *

|INFO  =   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 )

|or

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

        *                                     *
        | (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) |
        *                                      *

|INFO =      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 )

|or

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

        *                                     *
        | (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) |
        *                                       *

|INFO =     0