Guide and Reference

Linear Least Squares Considerations

This section provides some key points about using the linear least squares subroutines.

Use Considerations

If you want to use a singular value decomposition method to compute the minimal norm linear least squares solution of AX  is congruent to  B, calls to SGESVF or DGESVF should be followed by calls to SGESVS or DGESVS, respectively.

Performance and Accuracy Considerations

1. Least squares solutions obtained by using a singular value decomposition require more storage and run time than those obtained using a QR decomposition with column pivoting. The singular value decomposition method, however, is a more reliable way to handle rank deficiency.

2. The short-precision subroutines provide increased accuracy by accumulating intermediate results in long precision. Occasionally, for performance reasons, these intermediate results are stored.

3. The accuracy of the resulting singular values and singular vectors varies between the short- and long-precision versions of each subroutine. The degree of difference depends on the size and conditioning of the matrix computation.

4. There are ESSL-specific rules that apply to the results of computations on the workstation processors using the ANSI/IEEE standards. For details, see "What Data Type Standards Are Used by ESSL, and What Exceptions Should You Know About?".

Dense Linear Algebraic Equation Subroutines

This section contains the dense linear algebraic equation subroutine descriptions.

SGEF, DGEF, CGEF, and ZGEF--General Matrix Factorization

This subroutine factors a square general matrix A using Gaussian elimination with partial pivoting. To solve the system of equations with one or more right-hand sides, follow the call to these subroutines with one or more calls to SGES/SGESM, DGES/DGESM, CGES/CGESM, or ZGES/ZGESM, respectively. To compute the inverse of matrix A, follow the call to these subroutines with a call to SGEICD or DGEICD, respectively.

Table 87. Data Types

A	Subroutine
Short-precision real	SGEF
Long-precision real	DGEF
Short-precision complex	CGEF
Long-precision complex	ZGEF

Note:

The output from these factorization subroutines should be used only as input to the following subroutines for performing a solve or inverse: SGES/SGESM/SGEICD, DGES/DGESM/DGEICD, CGES/CGESM, and ZGES/ZGESM, respectively.

Syntax

On Entry

a

is the n by n general matrix A to be factored. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 87.

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 of matrix A. Specified as: a fullword integer; 0 <= n <= lda.

ipvt

See 'On Return'.

On Return

a

is the n by n transformed matrix A, containing the results of the factorization. See "Function". Returned as: an lda by (at least) n array, containing numbers of the data type indicated in Table 87.

ipvt

is the integer vector ipvt of length n, containing the pivot information necessary to construct matrix L from the information contained in the output array a. Returned as: a one-dimensional array of (at least) length n, containing fullword integers.

Notes

1. ipvt is not a permutation vector in the strict sense. It is used to record row interchanges in L due to partial pivoting.

2. Calling SGEFCD or DGEFCD with iopt = 0 is equivalent to calling SGEF or DGEF.

Function

The matrix A is factored using Gaussian elimination with partial pivoting to compute the LU factorization of A, where:

ipvt is the vector containing the pivoting information.

L is a unit lower triangular matrix.

U is an upper triangular matrix.

The transformed matrix A contains U in the upper triangle. In its strict lower triangle, it contains the multipliers necessary to construct, with the help of ipvt, a matrix L, such that A = LU.

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

Error Conditions

Resource Errors

Unable to allocate internal work area.

Computational Errors

Matrix A is singular.

• One or more columns of L and the corresponding diagonal of U contain all zeros (all columns of L are checked). The first column, i, of L with a corresponding U = 0 diagonal element 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 2103 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. lda <= 0
2. n < 0
3. n > lda

Example 1

This example shows a factorization of a real general matrix A of order 9.

Call Statement and Input

           A  LDA  N   IPVT
           |   |   |    |
CALL SGEF( A , 9 , 9 , IPVT )

        *                                                *
        | 1.0  1.0  1.0  1.0  0.0  0.0   0.0   0.0   0.0 |
        | 1.0  1.0  1.0  1.0  1.0  0.0   0.0   0.0   0.0 |
        | 4.0  1.0  1.0  1.0  1.0  1.0   0.0   0.0   0.0 |
        | 0.0  5.0  1.0  1.0  1.0  1.0   1.0   0.0   0.0 |
A    =  | 0.0  0.0  6.0  1.0  1.0  1.0   1.0   1.0   0.0 |
        | 0.0  0.0  0.0  7.0  1.0  1.0   1.0   1.0   1.0 |
        | 0.0  0.0  0.0  0.0  8.0  1.0   1.0   1.0   1.0 |
        | 0.0  0.0  0.0  0.0  0.0  9.0   1.0   1.0   1.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  10.0  11.0  12.0 |
        *                                                *

Output

        *                                                                             *
        | 4.0000  1.0000  1.0000  1.0000   1.0000   1.0000   0.0000   0.0000   0.0000 |
        | 0.0000  5.0000  1.0000  1.0000   1.0000   1.0000   1.0000   0.0000   0.0000 |
        | 0.0000  0.0000  6.0000  1.0000   1.0000   1.0000   1.0000   1.0000   0.0000 |
        | 0.0000  0.0000  0.0000  7.0000   1.0000   1.0000   1.0000   1.0000   1.0000 |
A    =  | 0.0000  0.0000  0.0000  0.0000   8.0000   1.0000   1.0000   1.0000   1.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   9.0000   1.0000   1.0000   1.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000  10.0000  11.0000  12.0000 |
        | 0.2500  0.1500  0.1000  0.0714   0.0536  -0.0694  -0.0306   0.1806   0.3111 |
        | 0.2500  0.1500  0.1000  0.0714  -0.0714  -0.0556  -0.0194   0.9385  -0.0031 |
        *                                                                             *

IPVT     =  (3, 4, 5, 6, 7, 8, 9, 8, 9)

Example 2

This example shows a factorization of a complex general matrix A of order 4.

Call Statement and Input

           A  LDA  N   IPVT
           |   |   |    |
CALL CGEF( A , 4 , 4 , IPVT )
 
        *                                             *
        | (1.0, 2.0) (1.0, 7.0) (2.0, 4.0) (3.0, 1.0) |
A    =  | (2.0, 0.0) (1.0, 3.0) (4.0, 4.0) (2.0, 3.0) |
        | (2.0, 1.0) (5.0, 0.0) (3.0, 6.0) (0.0, 0.0) |
        | (8.0, 5.0) (1.0, 9.0) (6.0, 6.0) (8.0, 1.0) |
        *                                             *

Output

        *                                                                          *
        |  (8.0000, 5.0000)   (1.0000, 9.0000)  (6.0000, 6.0000)  (8.0000, 1.0000) |
A    =  |  (0.2022, 0.1236)   (1.9101, 5.0562)  (1.5281, 2.0449) (1.5056, -0.1910) |
        | (0.2360, -0.0225) (-0.0654, -0.9269) (-0.3462, 6.2692) (-1.6346, 1.3269) |
        | (0.1798, -0.1124)   (0.2462, 0.1308) (0.4412, -0.3655)  (0.2900, 2.3864) |
        *                                                                          *

IPVT     =  (4, 4, 3, 4)

SGES, DGES, CGES, and ZGES--General Matrix, Its Transpose, or Its Conjugate Transpose Solve

These subroutines solve the system Ax = b for x, where A is a general matrix and x and b are vectors. Using the iopt argument, they can also solve the real system A^Tx = b or the complex system A^Hx = b for x. These subroutines use the results of the factorization of matrix A, produced by a preceding call to SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, or ZGEF, respectively.

Table 88. Data Types

A, b, x	Subroutine
Short-precision real	SGES
Long-precision real	DGES
Short-precision complex	CGES
Long-precision complex	ZGES

Note:

The input to these solve subroutines must be the output from the factorization subroutines SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, and ZGEF, respectively.

Syntax

On Entry

a

is the factorization of matrix A, produced by a preceding call to SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, or ZGEF, respectively. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 88.

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 of matrix A. Specified as: a fullword integer; 0 <= n <= lda.

ipvt

is the integer vector ipvt of length n, produced by a preceding call to SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, or ZGEF, respectively. It contains the pivot information necessary to construct matrix L from the information contained in the array specified for a.

Specified as: a one-dimensional array of (at least) length n, containing fullword integers.

bx

is the vector b of length n, containing the right-hand side of the system. Specified as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 88.

iopt

determines the type of computation to be performed, where:

If iopt = 0, A is used in the computation.

If iopt = 1, A^T is used in SGES and DGES. A^H is used in CGES and ZGES.

Note:

No data should be moved to form AT or AH; that is, the matrix A should always be stored in its untransposed form.

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

On Return

bx

is the solution vector x of length n, containing the results of the computation. Returned as: a one-dimensional array, containing numbers of the data type indicated in Table 88.

Notes

1. The scalar data specified for input arguments lda and n for these subroutines must be the same as the corresponding input arguments specified for SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, and ZGEF, respectively.

2. The array data specified for input arguments a and ipvt for these subroutines must be the same as the corresponding output arguments for SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, and ZGEF, respectively.

3. The vectors and matrices used in this computation must have no common elements; otherwise, results are unpredictable. See "Concepts".

Function

The system Ax = b is solved for x, where A is a general matrix and x and b are vectors. Using the iopt argument, this subroutine can also solve the real system A^Tx = b or the complex system A^Hx = b for x. These subroutines use the results of the factorization of matrix A, produced by a preceding call to SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, or ZGEF, respectively. The transformed matrix A consists of the upper triangular matrix U and the multipliers necessary to construct L using ipvt, as defined in "Function". For a description of how A is factored, see SGEF, DGEF, CGEF, and ZGEF--General Matrix Factorization.

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

Error Conditions

Computational Errors

None

Note:

If the factorization performed by SGEF, DGEF, CGEF, ZGEF, SGEFCD, DGEFCD, or DGEFP failed because a pivot element is zero, the results returned by this subroutine are unpredictable, and there may be a divide-by-zero program exception message.

Input-Argument Errors

1. lda <= 0
2. n < 0
3. n > lda
4. iopt <> 0 or 1

Example 1

Part 1

This part of the example shows how to solve the system Ax = b, where matrix A is the same matrix factored in the "Example 1" for SGEF and DGEF.

Call Statement and Input

           A  LDA  N   IPVT   BX  IOPT
           |   |   |    |     |    |
CALL SGES( A , 9 , 9 , IPVT , BX , 0  )

IPVT     =  (3, 4, 5, 6, 7, 8, 9, 8, 9)
BX       =  (4.0, 5.0, 9.0, 10.0, 11.0, 12.0, 12.0, 12.0, 33.0)
A        =(same as output A in
"Example 1")

Output

BX       =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)

Part 2

This part of the example shows how to solve the system A^Tx = b, where matrix A is the input matrix factored in "Example 1" for SGEF and DGEF. Most of the input is the same in Part 2 as in Part 1.

Call Statement and Input

           A  LDA  N   IPVT   BX  IOPT
           |   |   |    |     |    |
CALL SGES( A , 9 , 9 , IPVT , BX , 1  )

IPVT     =  (3, 4, 5, 6, 7, 8, 9, 8, 9)
BX       =  (6.0, 8.0, 10.0, 12.0, 13.0, 14.0, 15.0, 15.0, 15.0)
A        =(same as output A in
"Example 1")

Output

BX       =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)

Example 2

Part 1

This part of the example shows how to solve the system Ax = b, where matrix A is the same matrix factored in the "Example 2" for CGEF and ZGEF.

Call Statement and Input

           A  LDA  N   IPVT   BX  IOPT
           |   |   |    |     |    |
CALL CGES( A , 4 , 4 , IPVT , BX , 0  )

IPVT     =  (4, 4, 3, 4)
BX       =  ((-10.0, 85.0), (-6.0, 61.0), (10.0, 38.0),
             (58.0, 168.0))
A        =(same as output A in
"Example 1")

Output

BX       =  ((9.0, 0.0), (5.0, 1.0), (1.0, 6.0), (3.0, 4.0))

Part 2

This part of the example shows how to solve the system A^Hx = b, where matrix A is the input matrix factored in "Example 2" for CGEF and ZGEF. Most of the input is the same in Part 2 as in Part 1.

Call Statement and Input

           A  LDA  N   IPVT   BX  IOPT
           |   |   |    |     |    |
CALL CGES( A , 4 , 4 , IPVT , BX , 1  )

IPVT     =  (4, 4, 3, 4)
BX       =  ((71.0, 12.0), (61.0, -70.0), (123.0, -34.0),
             (68.0, 7.0))
A        =(same as output A in
"Example 1")

Output

BX       =  ((9.0, 0.0), (5.0, 1.0), (1.0, 6.0), (3.0, 4.0))

SGESM, DGESM, CGESM, and ZGESM--General Matrix, Its Transpose, or Its Conjugate Transpose Multiple Right-Hand Side Solve

These subroutines solve the following systems of equations for multiple right-hand sides, where A, X, and B are general matrices. SGESM and DGESM solve one of the following:

1. AX = B

2. A^TX = B

CGESM and ZGESM solve one of the following:

1. AX = B

2. A^TX = B

3. A^HX = B

These subroutines use the results of the factorization of matrix A, produced by a preceding call to SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, or ZGEF, respectively.

Table 89. Data Types

A, B, X	Subroutine
Short-precision real	SGESM
Long-precision real	DGESM
Short-precision complex	CGESM
Long-precision complex	ZGESM

Note:

The input to these solve subroutines must be the output from the factorization subroutines SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, and ZGEF, respectively.

Syntax

On Entry

trans

indicates the form of matrix A to use in the computation, where:

If transa = 'N', A is used in the computation, resulting in equation 1.

If transa = 'T', A^T is used in the computation, resulting in equation 2.

If transa = 'C', A^H is used in the computation, resulting in equation 3.

Specified as: a single character. It must be 'N', 'T', or 'C'.

a

is the factorization of matrix A, produced by a preceding call to SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, or ZGEF, respectively. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 89.

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 of matrix A. Specified as: a fullword integer; 0 <= n <= lda.

ipvt

is the integer vector ipvt of length n, produced by a preceding call to SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, or ZGEF, respectively. It contains the pivot information necessary to construct matrix L from the information contained in the array specified for a.

Specified as: a one-dimensional array of (at least) length n, containing fullword integers.

b

is the matrix B, containing the nrhs right-hand sides of the system. The right-hand sides, each of length n, reside in the columns of matrix B. Specified as: an ldb by (at least) nrhs array, containing numbers of the data type indicated in Table 89.

ldb

is the leading dimension of the array specified for b. Specified as: a fullword integer; ldb > 0 and ldb >= n.

nrhs

is the number of right-hand sides in the system to be solved. Specified as: a fullword integer; nrhs >= 0.

On Return

b

is the matrix B, containing the nrhs solutions to the system in the columns of B. Specified as: an ldb by (at least) nrhs array, containing numbers of the data type indicated in Table 89.

Notes

1. For SGESM and DGESM, if you specify 'C' for the trans argument, it is interpreted as though you specified 'T'.

2. The scalar data specified for input arguments lda and n for these subroutines must be the same as the corresponding input arguments specified for SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, and ZGEF, respectively.

3. The array data specified for input arguments a and ipvt for these subroutines must be the same as the corresponding output arguments for SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, and ZGEF, respectively.

4. The vectors and matrices used in this computation must have no common elements; otherwise, results are unpredictable. See "Concepts".

Function

One of the following systems of equations is solved for multiple right-hand sides:

1. AX = B

2. A^TX = B

3. A^HX = B (only for CGESM and ZGESM)

where A, B, and X are general matrices. These subroutines use the results of the factorization of matrix A, produced by a preceding call to SGEF/SGEFCD, DGEF/DGEFP/DGEFCD, CGEF, or ZGEF, respectively. The transformed matrix A consists of the upper triangular matrix U and the multipliers necessary to construct L using ipvt, as defined in "Function". For a description of how A is factored, see SGEF, DGEF, CGEF, and ZGEF--General Matrix Factorization.

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

Error Conditions

Computational Errors

None

Note:

If the factorization performed by SGEF, DGEF, CGEF, ZGEF, SGEFCD, DGEFCD, or DGEFP failed because a pivot element is zero, the results returned by this subroutine are unpredictable, and there may be a divide-by-zero program exception message.

Input-Argument Errors

1. trans <> 'N', 'T', or 'C'
2. lda, ldb <= 0
3. n < 0
4. n > lda, ldb
5. nrhs < 0

Example 1

Part 1

This part of the example shows how to solve the system AX = B for two right-hand sides, where matrix A is the same matrix factored in the "Example 1" for SGEF and DGEF.

Call Statement and Input

           TRANS  A  LDA  N   IPVT   B  LDB  NRHS
             |    |   |   |    |     |   |    |
CALL SGESM( 'N' , A , 9 , 9 , IPVT , B , 9 ,  2  )

IPVT     =  (3, 4, 5, 6, 7, 8, 9, 8, 9)
A        =(same as output A in
"Example 1")

        *             *
        |  4.0   10.0 |
        |  5.0   15.0 |
        |  9.0   24.0 |
        | 10.0   35.0 |
B    =  | 11.0   48.0 |
        | 12.0   63.0 |
        | 12.0   70.0 |
        | 12.0   78.0 |
        | 33.0  266.0 |
        *             *

Output

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

Part 2

This part of the example shows how to solve the system A^TX = B for two right-hand sides, where matrix A is the input matrix factored in "Example 1" for SGEF and DGEF.

Call Statement and Input

           TRANS  A  LDA  N   IPVT   B  LDB  NRHS
             |    |   |   |    |     |   |    |
CALL SGESM( 'T' , A , 9 , 9 , IPVT , B , 9 ,  2  )

IPVT     =  (3, 4, 5, 6, 7, 8, 9, 8, 9)
A        =(same as output A in
"Example 1")

        *             *
        |  6.0   15.0 |
        |  8.0   26.0 |
        | 10.0   40.0 |
        | 12.0   57.0 |
B    =  | 13.0   76.0 |
        | 14.0   97.0 |
        | 15.0  120.0 |
        | 15.0  125.0 |
        | 15.0  129.0 |
        *             *

Output

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

Example 2

Part 1

This part of the example shows how to solve the system AX = B for two right-hand sides, where matrix A is the same matrix factored in the "Example 2" for CGEF and ZGEF.

Call Statement and Input

           TRANS  A  LDA  N   IPVT   B  LDB  NRHS
             |    |   |   |    |     |   |    |
CALL CGESM( 'N' , A , 4 , 4 , IPVT , B , 4 ,  2  )

IPVT     =  (4, 4, 3, 4)
A        =(same as output A in
"Example 2")

        *                              *
        | (-10.0, 85.0)  (-11.0, 53.0) |
B    =  |  (-6.0, 61.0)   (-6.0, 54.0) |
        |  (10.0, 38.0)    (2.0, 40.0) |
        | (58.0, 168.0)  (15.0, 105.0) |
        *                              *

Output

        *                        *
        | (9.0, 0.0)  (1.0, 1.0) |
B    =  | (5.0, 1.0)  (2.0, 2.0) |
        | (1.0, 6.0)  (3.0, 3.0) |
        | (3.0, 4.0)  (4.0, 4.0) |
        *                        *

Part 2

This part of the example shows how to solve the system A^TX = B for two right-hand sides, where matrix A is the input matrix factored in "Example 2" for CGEF and ZGEF.

Call Statement and Input

           TRANS  A  LDA  N   IPVT   B  LDB  NRHS
             |    |   |   |    |     |   |    |
CALL CGESM( 'T' , A , 4 , 4 , IPVT , B , 4 ,  2  )

IPVT     =  (4, 4, 3, 4)
A        =(same as output A in
"Example 2")

        *                               *
        |   (71.0, 12.0)   (18.0, 68.0) |
B    =  |  (61.0, -70.0)  (-27.0, 71.0) |
        | (123.0, -34.0)  (-11.0, 97.0) |
        |    (68.0, 7.0)   (28.0, 50.0) |
        *                               *

Output

        *                        *
        | (9.0, 0.0)  (1.0, 1.0) |
B    =  | (5.0, 1.0)  (2.0, 2.0) |
        | (1.0, 6.0)  (3.0, 3.0) |
        | (3.0, 4.0)  (4.0, 4.0) |
        *                        *

Part 3

This part of the example shows how to solve the system A^HX = B for two right-hand sides, where matrix A is the input matrix factored in "Example 2" for CGEF and ZGEF.

Call Statement and Input

           TRANS  A  LDA  N   IPVT   B  LDB  NRHS
             |    |   |   |    |     |   |    |
CALL CGESM( 'C' , A , 4 , 4 , IPVT , B , 4 ,  2  )

IPVT     =  (4, 4, 3, 4)
A        =(same as output A in
"Example 2")

        *                              *
        |  (58.0, -3.0)   (45.0, 20.0) |
B    =  | (68.0, -31.0)  (83.0, -20.0) |
        | (89.0, -22.0)    (98.0, 1.0) |
        |  (53.0, 15.0)   (45.0, 25.0) |
        *                              *

Output

        *                        *
        | (1.0, 4.0)  (4.0, 5.0) |
B    =  | (2.0, 3.0)  (3.0, 4.0) |
        | (3.0, 2.0)  (2.0, 3.0) |
        | (4.0, 1.0)  (1.0, 2.0) |
        *                        *

SGETRF, DGETRF, CGETRF and ZGETRF--General Matrix Factorization

These subroutines factor general matrix A using Gaussian elimination with partial pivoting. To solve the system of equations with one or more right-hand sides, follow the call to these subroutines with one or more calls to SGETRS, DGETRS CGETRS, or ZGETRS, respectively. To compute the inverse of matrix A, follow the call to these subroutines with a call to SGEICD or DGEICD, respectively.

Table 90. Data Types

A	Subroutine
Short-precision real	SGETRF
Long-precision real	DGETRF
Short-precision complex	CGETRF
Long-precision complex	ZGETRF

Note:

The output from these factorization subroutines should be used only as input to the following subroutines for performing a solve or inverse: SGETRS, DGETRS, CGETRS, ZGETRS, SGEICD or DGEICD respectively.

Syntax

On Entry

m

the number of rows in general matrix A used in the computation. Specified as: a fullword integer; 0 <= m <= lda.

n

the number of columns in general matrix A used in the computation. Specified as: a fullword integer; n >= 0.

a

is the m by n general matrix A to be factored. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 90.

lda

is the leading dimension of matrix A. Specified as: a fullword integer; lda > 0 and lda >= m.

ipvt

See 'On Return'.

info

See 'On Return'.

On Return

a

is the m by n transformed matrix A, containing the results of the factorization. See "Function". Returned as: an lda by (at least) n array, containing numbers of the data type indicated in Table 90.

ipvt

is the integer vector ipvt of length min(m,n), containing the pivot information necessary to construct matrix L from the information contained in the output array A. Returned as: a one-dimensional array of (at least) length min(m,n), containing fullword integers,where 1 <= ipvt(i) <= m.

info

has the following meaning:

If info = 0, the factorization of general matrix A completed successfully.

If info > 0, info is set equal to the first i, where U_ii is singular and its inverse could not be computed.

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

Notes

1. In your C program, argument info must be passed by reference.

2. The matrix A and vector ipvt must have no common elements; otherwise results are unpredictable.

3. The way these subroutines handle singularity differs from LAPACK. These subroutines us the info argument to provide information about the singularity of A, like LAPACK, but also provide an error message.

4. On both input and output, matrix A conforms to LAPACK format.

Function

The matrix A is factored using Gaussian elimination with partial pivoting to compute the LU factorization of A, where:

ipvt is the vector containing the pivoting information.

L is a unit lower triangular matrix.

U is an upper triangular matrix.

On output, the transformed matrix A contains U in the upper triangle (if m >= n) or upper trapezoid (if m < n). In its strict lower triangle (if m <= n) or lower trapezoid (if m > n), it contains the multipliers necessary to construct, with the help of ipvt, a matrix L, such that A = LU.

If m or n is 0, no computation is performed and the subroutine returns after doing some parameter checking. See references [36] and [59].

Error Conditions

Resource Errors

Unable to allocate internal work area.

Computational Errors

Matrix A is singular.

• The first column, i, of L with a corresponding U_ii = 0 diagonal element 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 2146 is set to be unlimited in the ESSL error option table.

Input-Argument Errors

1. m < 0
2. n < 0
3. m > lda
4. lda <= 0

Example 1

This example shows a factorization of a real general matrix A of order 9.

Call Statement and Input

             M   N   A  LDA  IPVT  INFO
             |   |   |   |     |     |
CALL DGETRF( 9 , 9 , A,  9 , IPVT, INFO )

        *                                             *
        | 1.0  1.2  1.4  1.6  1.8  2.0  2.2  2.4  2.6 |
        | 1.2  1.0  1.2  1.4  1.6  1.8  2.0  2.2  2.4 |
        | 1.4  1.2  1.0  1.2  1.4  1.6  1.8  2.0  2.2 |
        | 1.6  1.4  1.2  1.0  1.2  1.4  1.6  1.8  2.0 |
A    =  | 1.8  1.6  1.4  1.2  1.0  1.2  1.4  1.6  1.8 |
        | 2.0  1.8  1.6  1.4  1.2  1.0  1.2  1.4  1.6 |
        | 2.2  2.0  1.8  1.6  1.4  1.2  1.0  1.2  1.4 |
        | 2.4  2.2  2.0  1.8  1.6  1.4  1.2  1.0  1.2 |
        | 2.6  2.4  2.2  2.0  1.8  1.6  1.4  1.2  1.0 |
        *                                             *

Output

        *                                              *
        | 2.6   2.4  2.2  2.0  1.8  1.6  1.4  1.2  1.0 |
        | 0.4   0.3  0.6  0.8  1.1  1.4  1.7  1.9  2.2 |
        | 0.5  -0.4  0.4  0.8  1.2  1.6  2.0  2.4  2.8 |
        | 0.5  -0.3  0.0  0.4  0.8  1.2  1.6  2.0  2.4 |
A    =  | 0.6  -0.3  0.0  0.0  0.4  0.8  1.2  1.6  2.0 |
        | 0.7  -0.2  0.0  0.0  0.0  0.4  0.8  1.2  1.6 |
        | 0.8  -0.2  0.0  0.0  0.0  0.0  0.4  0.8  1.2 |
        | 0.8  -0.1  0.0  0.0  0.0  0.0  0.0  0.4  0.8 |
        | 0.9  -0.1  0.0  0.0  0.0  0.0  0.0  0.0  0.4 |
        *                                              *

IPVT     =  (9, 9, 9, 9, 9, 9, 9, 9, 9)
INFO     =  0

Example 2

This example shows a factorization of a complex general matrix A of order 9.

Call Statement and Input

             M   N   A  LDA  IPVT  INFO
             |   |   |   |     |     |
CALL ZGETRF( 9 , 9 , A,  9 , IPVT, INFO )

 
        *                                                                                                     *
        | (2.0, 1.0) (2.4,-1.0) (2.8,-1.0) (3.2,-1.0)  (3.6,-1.0) (4.0,-1.0) (4.4,-1.0) (4.8,-1.0) (5.2,-1.0) |
        | (2.4, 1.0) (2.0, 1.0) (2.4,-1.0) (2.8,-1.0)  (3.2,-1.0) (3.6,-1.0) (4.0,-1.0) (4.4,-1.0) (4.8,-1.0) |
        | (2.8, 1.0) (2.4, 1.0) (2.0, 1.0) (2.4,-1.0)  (2.8,-1.0) (3.2,-1.0) (3.6,-1.0) (4.0,-1.0) (4.4,-1.0) |
        | (3.2, 1.0) (2.8, 1.0) (2.4, 1.0) (2.0, 1.0)  (2.4,-1.0) (2.8,-1.0) (3.2,-1.0) (3.6,-1.0) (4.0,-1.0) |
A    =  | (3.6, 1.0) (3.2, 1.0) (2.8, 1.0) (2.4, 1.0)  (2.0, 1.0) (2.4,-1.0) (2.8,-1.0) (3.2,-1.0) (3.6,-1.0) |
        | (4.0, 1.0) (3.6, 1.0) (3.2, 1.0) (2.8, 1.0)  (2.4, 1.0) (2.0, 1.0) (2.4,-1.0) (2.8,-1.0) (3.2,-1.0) |
        | (4.4, 1.0) (4.0, 1.0) (3.6, 1.0) (3.2, 1.0)  (2.8, 1.0) (2.4, 1.0) (2.0, 1.0) (2.4,-1.0) (2.8,-1.0) |
        | (4.8, 1.0) (4.4, 1.0) (4.0, 1.0) (3.6, 1.0)  (3.2, 1.0) (2.8, 1.0) (2.4, 1.0) (2.0, 1.0) (2.4,-1.0) |
        | (5.2, 1.0) (4.8, 1.0) (4.4, 1.0) (4.0, 1.0)  (3.6, 1.0) (3.2, 1.0) (2.8, 1.0) (2.4, 1.0) (2.0, 1.0) |
        *                                                                                                     *

Output

        *                                                                                                         *
        | (5.2, 1.0) (4.8, 1.0) (4.4, 1.0)   (4.0, 1.0)   (3.6, 1.0)  (3.2, 1.0) (2.8, 1.0) (2.4, 1.0) (2.0, 1.0) |
        | (0.4, 0.1) (0.6,-2.0) (1.1,-1.9)   (1.7,-1.9)   (2.3,-1.8)  (2.8,-1.8) (3.4,-1.7) (3.9,-1.7) (4.5,-1.6) |
        | (0.5, 0.1) (0.0,-0.1) (0.6,-1.9)   (1.2,-1.8)   (1.8,-1.7)  (2.5,-1.6) (3.1,-1.5) (3.7,-1.4) (4.3,-1.3) |
        | (0.6, 0.1) (0.0,-0.1) (-0.1,-0.1)  (0.7,-1.9)   (1.3,-1.7)  (2.0,-1.6) (2.7,-1.5) (3.4,-1.4) (4.0,-1.2) |
A    =  | (0.6, 0.1) (0.0,-0.1) (-0.1,-0.1) (-0.1, 0.0)   (0.7,-1.9)  (1.5,-1.7) (2.2,-1.6) (2.9,-1.5) (3.7,-1.3) |
        | (0.7, 0.1) (0.0,-0.1)  (0.0, 0.0) (-0.1, 0.0)  (-0.1, 0.0)  (0.8,-1.9) (1.6,-1.8) (2.4,-1.6) (3.2,-1.5) |
        | (0.8, 0.0) (0.0, 0.0)  (0.0, 0.0)  (0.0, 0.0)   (0.0, 0.0)  (0.0, 0.0) (0.8,-1.9) (1.7,-1.8) (2.5,-1.8) |
        | (0.9, 0.0) (0.0, 0.0)  (0.0, 0.0)  (0.0, 0.0)   (0.0, 0.0)  (0.0, 0.0) (0.0, 0.0) (0.8,-2.0) (1.7,-1.9) |
        | (0.9, 0.0) (0.0, 0.0)  (0.0, 0.0)  (0.0, 0.0)   (0.0, 0.0)  (0.0, 0.0) (0.0, 0.0) (0.0, 0.0) (0.8,-2.0) |
        *                                                                                                         *

IPVT     =  (9, 9, 9, 9, 9, 9, 9, 9, 9)
INFO     =  0

SGETRS, DGETRS, CGETRS, and ZGETRS--General Matrix Multiple Right-Hand Side Solve

SGETRS and DGETRS solve one of the following systems of equations for multiple right-hand sides:

1. AX = B

2. A^TX = B

CGETRS and ZGETRS solve one of the following systems of equations for multiple right-hand sides:

1. AX = B

2. A^TX = B

3. A^HX = B

In the formulas above:

A represents the general matrix A containing the LU factorization.

B represents the general matrix B containing the right-hand sides in its columns.

X represents the general matrix B containing the solution vectors in its columns.

These subroutines use the results of the factorization of matrix A, produced by a preceding call to SGETRF, DGETRF, CGETRF, or ZGETRF, respectively.

Table 91. Data Types

A, B	Subroutine
Short-precision real	SGETRS
Long-precision real	DGETRS
Short-precision complex	CGETRS
Long-precision complex	ZGETRS

Note:

The input to these solve subroutines must be the output from the factorization subroutines SGETRF, DGETRF, CGETRF and ZGETRF, respectively.

Syntax

On Entry

transa

indicates the form of matrix A to use in the computation, where:

If transa = 'N', A is used in the computation, resulting in solution 1.

If transa = 'T', A^T is used in the computation, resulting in solution 2.

If transa = 'C', A^H is used in the computation, resulting in solution 3.

Specified as: a single character; transa = 'N', 'T', or 'C'.

n

is the order of factored matrix A and the number of rows in matrix B. Specified as: a fullword integer; n >= 0.

nrhs

the number of right-hand sides--that is, the number of columns in matrix B used in the computation. Specified as: a fullword integer; nrhs >= 0.

a

is the factorization of matrix A, produced by a preceding call to SGETRF, DGETRF, CGETRF, or ZGETRF, respectively. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 91.

lda

is the leading dimension of the array specified for a. Specified as: a fullword integer; lda > 0 and lda >= n.

ipvt

is the integer vector ipvt of length n, produced by a preceding call to SGETRF, DGETRF, CGETRF or ZGETRF, respectively. It contains the pivot information necessary to construct matrix L from the information contained in the array specified for a.

Specified as: a one-dimensional array of (at least) length n, containing fullword integers, where 1 <= ipvt(i) <= n.

bx

is the general matrix B containing the right-hand side of the system. Specified as: an ldb by (at least) nrhs array, containing numbers of the data type indicated in Table 91.

ldb

is the leading dimension of the array specified for b. Specified as: a fullword integer; ldb > 0 and ldb >= n.

info

See 'On Return'.

On Return

bx

is the solution X containing the results of the computation. Returned as: an ldb by (at least) nrhs array, containing numbers of the data type indicated in Table 91.

info

info has the following meaning:

If info = 0, the solve of general matrix A completed successfully.

Notes

1. In your C program, argument info must be passed by reference.

2. These subroutines accept lower case letters for the transa argument.

3. For SGETRS and DGETRS, if you specify 'C' for the transa argument, it is interpreted as though you specified 'T'.

4. The scalar data specified for input argument n must be the same for both _GETRF and _GETRS. In addition, the scalar data specified for input argument m in _GETRF must be the same as input argument n in both _GETRF and _GETRS.

   If, however, you do not plan to call _GETRS after calling _GETRF, then input arguments m and n in _GETRF do not need to be equal.

5. The array data specified for input arguments a and ipvt for these subroutines must be the same as the corresponding output arguments for SGETRF, DGETRF, CGETRF, and ZGETRF, respectively.

6. The matrices and vector used in this computation must have no common elements; otherwise, results are unpredictable. See "Concepts".

7. On both input and output, matrices A and B conform to LAPACK format.

Function

One of the following systems of equations is solved for multiple right-hand sides:

1. AX = B

2. A^TX = B

3. A^HX = B (only for CGETRS and ZGETRS)

where A, B, and X are general matrices. These subroutines uses the results of the factorization of matrix A, produced by a preceding call to SGETRF, DGETRF, CGETRF or ZGETRF, respectively. On input, the transformed matrix A consists of the upper triangular matrix U and the multipliers necessary to construct L using ipvt, as defined in "Function". For details on the factorization, see SGETRF, DGETRF, CGETRF and ZGETRF--General Matrix Factorization.

If n = 0 or nrhs = 0, no computation is performed and the subroutine returns after doing some parameter checking. See references [36] and [59].

Error Conditions

Computational Errors

None

Note:

If the factorization performed by SGETRF, DGETRF, CGETRF or ZGETRF failed because a pivot element is zero, the results returned by this subroutine are unpredictable, and there may be a divide-by-zero program exception message.

Input-Argument Errors

1. transa <> 'N', 'T', or 'C'
2. n < 0
3. nrhs < 0
4. n > lda
5. n > ldb
6. lda <= 0
7. ldb <= 0

Example 1

This example shows how to solve the system AX = B, where matrix A is the same matrix factored in the "Example 1" for DGETRF.

Call Statement and Input

           TRANSA  N  NRHS  A  LDA  IPIV  BX LDB  INFO
             |     |    |   |   |     |   |   |     |
CALL DGETRS('N' ,  9 ,  5 , A , 9 , IPIV, B , 9 , INFO)

IPVT     =  (9, 9, 9, 9, 9, 9, 9, 9, 9)
A      =(same as output A in
"Example 1")

        *                                  *
        | 93.0  186.0  279.0  372.0  465.0 |
        | 84.4  168.8  253.2  337.6  422.0 |
        | 76.6  153.2  229.8  306.4  383.0 |
        | 70.0  140.0  210.0  280.0  350.0 |
B    =  | 65.0  130.0  195.0  260.0  325.0 |
        | 62.0  124.0  186.0  248.0  310.0 |
        | 61.4  122.8  184.2  245.6  307.0 |
        | 63.6  127.2  190.8  254.4  318.0 |
        | 69.0  138.0  207.0  276.0  345.0 |
        *                                  *

Output

        *                             *
        | 1.0   2.0   3.0   4.0   5.0 |
        | 2.0   4.0   6.0   8.0  10.0 |
        | 3.0   6.0   9.0  12.0  15.0 |
        | 4.0   8.0  12.0  16.0  20.0 |
B    =  | 5.0  10.0  15.0  20.0  25.0 |
        | 6.0  12.0  18.0  24.0  30.0 |
        | 7.0  14.0  21.0  28.0  35.0 |
        | 8.0  16.0  24.0  32.0  40.0 |
        | 9.0  18.0  27.0  36.0  45.0 |
        *                             *

INFO     =  0

Example 2

This example shows how to solve the system AX = b, where matrix A is the same matrix factored in the "Example 2" for ZGETRF.

Call Statement and Input

           TRANS   N  NRHS  A  LDA  IPIV  B  LDB  INFO
             |     |    |   |   |     |   |   |     |
CALL ZGETRS('N' ,  9 ,  5 , A , 9 , IPIV, B , 9 , INFO)

IPVT     =  (9, 9, 9, 9, 9, 9, 9, 9, 9)
A        =(same as output A in
"Example 2")

        *                                                                           *
        | (193.0,-10.6)  (200.0, 21.8)  (207.0, 54.2)  (214.0, 86.6)  (221.0,119.0) |
        | (173.8, -9.4)  (178.8, 20.2)  (183.8, 49.8)  (188.8, 79.4)  (193.8,109.0) |
        | (156.2, -5.4)  (159.2, 22.2)  (162.2, 49.8)  (165.2, 77.4)  (168.2,105.0) |
        | (141.0,  1.4)  (142.0, 27.8)  (143.0, 54.2)  (144.0, 80.6)  (145.0,107.0) |
B    =  | (129.0, 11.0)  (128.0, 37.0)  (127.0, 63.0)  (126.0, 89.0)  (125.0,115.0) |
        | (121.0, 23.4)  (118.0, 49.8)  (115.0, 76.2)  (112.0,102.6)  (109.0,129.0) |
        | (117.8, 38.6)  (112.8, 66.2)  (107.8, 93.8)  (102.8,121.4)   (97.8,149.0) |
        | (120.2, 56.6)  (113.2, 86.2)  (106.2,115.8)   (99.2,145.4)   (92.2,175.0) |
        | (129.0, 77.4)  (120.0,109.8)  (111.0,142.2). (102.0,174.6)   (93.0,207.0) |
        *                                                                           *

Output

        *                                                     *
        | (1.0,1.0)  (1.0,2.0)  (1.0,3.0) (1.0,4.0) (1.0,5.0) |
        | (2.0,1.0)  (2.0,2.0)  (2.0,3.0) (2.0,4.0) (2.0,5.0) |
        | (3.0,1.0)  (3.0,2.0)  (3.0,3.0) (3.0,4.0) (3.0,5.0) |
        | (4.0,1.0)  (4.0,2.0)  (4.0,3.0) (4.0,4.0) (4.0,5.0) |
B     = | (5.0,1.0)  (5.0,2.0)  (5.0,3.0) (5.0,4.0) (5.0,5.0) |
        | (6.0,1.0)  (6.0,2.0)  (6.0,3.0) (6.0,4.0) (6.0,5.0) |
        | (7.0,1.0)  (7.0,2.0)  (7.0,3.0) (7.0,4.0) (7.0,5.0) |
        | (8.0,1.0)  (8.0,2.0)  (8.0,3.0) (8.0,4.0) (8.0,5.0) |
        | (9.0,1.0)  (9.0,2.0)  (9.0,3.0) (9.0,4.0) (9.0,5.0) |
        *                                                     *

INFO     =  0

SGEFCD and DGEFCD--General Matrix Factorization, Condition Number Reciprocal, and Determinant

These subroutines factor general matrix A using Gaussian elimination. An estimate of the reciprocal of the condition number and the determinant of matrix A can also be computed. 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 SGES/SGESM or DGES/DGESM, respectively. To compute the inverse of matrix A, follow the call to these subroutines with a call to SGEICD and DGEICD, respectively.

Table 92. Data Types

A, aux, rcond, det	Subroutine
Short-precision real	SGEFCD
Long-precision real	DGEFCD

Note:

The output from these factorization subroutines should be used only as input to the following subroutines for performing a solve or inverse: SGES/SGESM/SGEICD and DGES/DGESM/DGEICD, respectively.

Syntax

On Entry

a

is a general matrix A of order n, whose factorization, reciprocal of condition number, and determinant are computed. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 92.

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 of matrix A. Specified as: a fullword integer; 0 <= n <= lda.

ipvt

See 'On Return'.

iopt

indicates the type of computation to be performed, where:

If iopt = 0, the matrix is factored.

If iopt = 1, the matrix is factored, and the reciprocal of the condition number is computed.

If iopt = 2, the matrix is factored, and the determinant is computed.

If iopt = 3, the matrix is factored, and the reciprocal of the condition number and the determinant are computed.

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

rcond

See 'On Return'.

det

See 'On Return'.

aux

has the following meaning:

If naux = 0 and error 2015 is unrecoverable, aux is ignored.

Otherwise, it is is a storage work area used by this subroutine. Its size is specified by naux.

Specified as: an area of storage, containing numbers of the data type indicated in Table 92.

naux

is the size of the work area specified by aux--that is, the number of elements in aux. Specified as: a fullword integer, where:

If naux = 0 and error 2015 is unrecoverable, SGEFCD and DGEFCD dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, naux >= n.

On Return

a

is the transformed matrix A of order n, containing the results of the factorization. See "Function". Returned as: an lda by (at least) n array, containing numbers of the data type indicated in Table 92.

ipvt

is the integer vector ipvt of length n, containing the pivot information necessary to construct matrix L from the information contained in output matrix A. Returned as: a one-dimensional array of (at least) length n, containing fullword integers.

rcond

is an estimate of the reciprocal of the condition number, rcond, of matrix A. Returned as: a number of the data type indicated in Table 92; rcond >= 0.

det

is the vector det, containing the two components, det₁ and det₂, of the determinant of matrix A. The determinant is:

where 1 <= det₁ < 10. Returned as: an array of length 2, containing numbers of the data type indicated in Table 92.

Notes

1. In your C program, argument rcond must be passed by reference.

2. When iopt = 0, these subroutines provide the same function as a call to SGEF or DGEF, respectively.

3. You have the option of having the minimum required value for naux dynamically returned to your program. For details, see "Using Auxiliary Storage in ESSL".

Function

Matrix A is factored using Gaussian elimination with partial pivoting to compute the LU factorization of A, where:

ipvt is the vector containing the pivoting information.

L is a unit lower triangular matrix.

U is an upper triangular matrix.

The transformed matrix A contains U in the upper triangle. In its strict lower triangle, it contains the multipliers necessary to construct, with the help of ipvt, a matrix L, such that A = LU.

An estimate of the reciprocal of the condition number, rcond, and the determinant, det, can also be computed by this subroutine. The estimate of the condition number uses an enhanced version of the algorithm described in references [63] and [64].

If n is 0, no computation is performed. See reference [36].

These subroutines call SGEF and DGEF, respectively, to perform the factorization. ipvt is an output vector of SGEF and DGEF. It is returned for use by SGES/SGESM and DGES/DGESM, the solve subroutines.

Error Conditions

Resource Errors

Error 2015 is unrecoverable, naux = 0, and unable to allocate work area.

Computational Errors

Matrix A is singular.

• If your program is not terminated by SGEF and DGEF, then SGEFCD and DGEFCD, respectively, return 0 for rcond and det.
• One or more columns of L and the corresponding diagonal of U contain all zeros (all columns of L are checked). The first column, i, of L with a corresponding U = 0 diagonal element is identified in the computational error message, issued by SGEF or DGEF, respectively.
• i can be determined at run time by using the ESSL error-handling facilities. To obtain this information, you must use ERRSET to change the number of allowable errors for error code 2103 in the ESSL error option table; otherwise, the default value causes your program to be terminated by SGEF or DGEF, respectively, when this error occurs. If your program is not terminated by SGEF or DGEF, respectively, the return code is set to 2. For details, see "What Can You Do about ESSL Computational Errors?".

Input-Argument Errors

1. lda <= 0
2. n < 0
3. n > lda
4. iopt <> 0, 1, 2, or 3
5. Error 2015 is recoverable or naux<>0, and naux is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.

Example

This example shows a factorization of matrix A of order 9. The input is the same as used in SGEF and DGEF. See "Example 1". The reciprocal of the condition number and the determinant of matrix A are also computed. The values used to estimate the reciprocal of the condition number in this example are obtained with the following values:

||A||₁ = max(6.0, 8.0, 10.0, 12.0, 13.0, 14.0, 15.0, 15.0, 15.0) = 15.0

Estimate of ||A^-1||₁ = 1091.87

This estimate is equal to the actual rcond of 5.436(10^-5), which is computed by SGEICD and DGEICD. (See "Example 1".) On output, the value in det, |A|, is equal to 336.

Call Statement and Input

             A  LDA  N   IPVT  IOPT  RCOND   DET   AUX  NAUX
             |   |   |    |     |     |       |     |    |
CALL DGEFCD( A , 9 , 9 , IPVT , 3  , RCOND , DET , AUX , 9  )

A        =(same as input A in
"Example 1")

Output

A        =(same as output A in
"Example 1")
IPVT     =  (3, 4, 5, 6, 7, 8, 9, 8, 9)
RCOND    =  0.00005436
DET      =  (3.36, 2.00)

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

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. 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], [64], 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) |
              *                                         *

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

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

SPPS and DPPS--Positive Definite Real Symmetric Matrix Solve

These subroutines solve the system Ax = b for x, where A is a positive definite symmetric matrix, and x and b are vectors. The subroutines use the results of the factorization of matrix A, produced by a preceding call to SPPF/SPPFCD or DPPF/DPPFP/DPPFCD, respectively.

Table 94. Data Types

A, b, x	Subroutine
Short-precision real	SPPS
Long-precision real	DPPS

Note:

The input to these solve subroutines must be the output from the factorization subroutines SPPF/SPPFCD and DPPF/DPPFP/DPPFCD, respectively.

Syntax

On Entry

ap

is the factorization of matrix A, produced by a preceding call to SPPF/SPPFCD or DPPF/DPPFP/DPPFCD, respectively. Specified as: a one-dimensional array, containing numbers of the data type indicated in Table 94, where:

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

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

n

is the order of matrix A used in the factorization, and the lengths of vectors b and x. Specified as: a fullword integer; n >= 0.

bx

is the vector b of length n, containing the right-hand side of the system. Specified as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 94.

iopt

indicates the type of factorization that was performed on matrix A, where:

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

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

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

On Return

bx

is the solution vector x of length n, containing the results of the computation. Specified as: a one-dimensional array, containing numbers of the data type indicated in Table 94.

Notes

1. The array data specified for input argument ap for these subroutines must be the same as the corresponding output argument for SPPF/SPPFCD and DPPF/DPPFP/DPPFCD, respectively.

2. The scalar data specified for input argument n for these subroutines must be the same as that specified for SPPF/SPPFCD and DPPF/DPPFP/DPPFCD, respectively.

3. When you call these subroutines after calling SPPF or DPPF, the value of input argument iopt must be the same as that specified for SPPF and DPPF.

4. When you call these subroutines after calling SPPFCD or DPPFCD, the value of input argument iopt must be 0.

5. When you call these subroutines after calling DPPFP, the value of input argument iopt must be 1.

6. In the input array 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.

7. The vectors and matrices used in this computation must have no common elements; otherwise, results are unpredictable. See "Concepts".

8. For a description of how a positive definite symmetric matrix is stored in lower-packed storage mode in an array, see "Symmetric Matrix".

Function

The system Ax = b is solved for x, where A is a positive definite symmetric matrix, stored in lower-packed storage mode in array AP, and x and b are vectors. These subroutines use the results of the factorization of matrix A, produced by a preceding call to SPPF/SPPFCD or DPPF/DPPFP/DPPFCD, respectively.

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

Error Conditions

Computational Errors

None

Note:

If a call to SPPF, DPPF, SPPFCD, DPPFCD, or DPPFP resulted in a nonpositive definite matrix, error 2104 or 2115, SPPS or DPPS results may be unpredictable or numerically unstable.

Input-Argument Errors

1. n < 0
2. iopt <> 0 or 1

Example 1

This example shows how to solve the system Ax = b, where matrix A is the same matrix factored in the "Example 1" for SPPF and DPPF.

Call Statement and Input

            AP   N   BX   IOPT
            |    |   |     |
CALL SPPS ( AP , 9 , BX ,  0 )

AP       =(same as output AP in
"Example 1" for SPPF and DPPF)
BX       =  (9.0, 17.0, 24.0, 30.0, 35.0, 39.0, 42.0, 44.0, 45.0)

Output

BX       =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)

Example 2

This example shows how to solve the same system as in Example 1, where matrix A is the same matrix factored in the "Example 2" for SPPF and DPPF.

Call Statement and Input

           AP   N   BX  IOPT
           |    |   |    |
CALL SPPS( AP , 9 , BX , 1 )

AP       =(same as output AP in
"Example 2" for SPPF and DPPF)
BX       =  (9.0, 17.0, 24.0, 30.0, 35.0, 39.0, 42.0, 44.0, 45.0)

Output

BX       =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)

SPOSM, DPOSM, CPOSM, and ZPOSM--Positive Definite Real Symmetric or Complex Hermitian Matrix Multiple Right-Hand Side Solve

These subroutines solve the system AX = B for X, using multiple right-hand sides, where X and B are general matrices and:

• For SPOSM and DPOSM, A is a positive definite symmetric matrix.
• For CPOSM and ZPOSM, A is a positive definite complex Hermitian matrix.

These subroutines use the results of the factorization of matrix A, produced by a preceding call to SPOF/SPOFCD, DPOF/DPOFCD, CPOF, or ZPOF, respectively.

Table 95. Data Types

A, B, X	Subroutine
Short-precision real	SPOSM
Long-precision real	DPOSM
Short-precision complex	CPOSM
Long-precision complex	ZPOSM

Note:

The input to these solve subroutines must be the output from the factorization subroutines SPOF/SPOFCD, DPOF/DPOFCD, CPOF, and ZPOF, respectively.

Syntax

On Entry

uplo

indicates whether the original 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'.

a

is the factorization of positive definite matrix A, produced by a preceding call to SPOF/SPOFCD, DPOF/DPOFCD, CPOF, or ZPOF, respectively. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 95.

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 of matrix A. Specified as: a fullword integer; 0 <= n <= lda.

b

is the matrix B, containing the nrhs right-hand sides of the system. The right-hand sides, each of length n, reside in the columns of matrix B. Specified as: an ldb by (at least) nrhs array, containing numbers of the data type indicated in Table 95.

ldb

is the leading dimension of the array specified for b. Specified as: a fullword integer; ldb > 0 and ldb >= n.

nrhs

is the number of right-hand sides in the system to be solved. Specified as: a fullword integer; nrhs >= 0.

On Return

b

is the matrix B, containing the nrhs solutions to the system in the columns of B. Specified as: an ldb by (at least) nrhs array, containing numbers of the data type indicated in Table 95.

Notes

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

2. The scalar data specified for input arguments uplo, lda, and n for these subroutines must be the same as the corresponding input arguments specified for SPOF/SPOFCD, DPOF/DPOFCD, CPOF, and ZPOF, respectively.

3. The array data specified for input argument a for these subroutines must be the same as the corresponding output arguments for SPOF/SPOFCD, DPOF/DPOFCD, CPOF, and ZPOF, respectively.

4. The vectors and matrices used in this computation must have no common elements; otherwise, results are unpredictable. See "Concepts".

5. For a description of how the matrices are stored:

   • 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 system AX = B is solved for X, using multiple right-hand sides, where X and B are general matrices, and A is a positive definite symmetric matrix for SPOSM and DPOSM and a positive definite complex Hermitian matrix for CPOSM and ZPOSM. These subroutines use the results of the factorization of matrix A, produced by a preceding call to SPOF/SPOFCD, DPOF/DPOFCD, CPOF, or ZPOF, respectively. For a description of how A is factored, see SPPF, DPPF, SPOF, DPOF, CPOF, and ZPOF--Positive Definite Real Symmetric or Complex Hermitian Matrix Factorization.

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

Error Conditions

Computational Errors

None

Note:

If the factorization performed by SPOF, DPOF, CPOF, ZPOF, SPOFCD, or DPOFCD failed because matrix A was not positive definite, the results returned by this subroutine are unpredictable, and there may be a divide-by-zero program exception message.

Input-Argument Errors

1. uplo <> 'U' or 'L'
2. lda, ldb <= 0
3. n < 0
4. n > lda
5. n > ldb
6. nrhs < 0

Example 1

This example shows how to solve the system AX = B for two right-hand sides, where matrix A is the same matrix factored in the "Example 3" for SPOF.

Call Statement and Input

            UPLO  A  LDA  N   B  LDB  NRHS
             |    |   |   |   |   |    |
CALL SPOSM( 'L' , A , 9 , 9 , B , 9 ,  2  )

A        =(same as output A in
"Example 3")

        *             *
        |  9.0   45.0 |
        | 17.0   89.0 |
        | 24.0  131.0 |
        | 30.0  170.0 |
B    =  | 35.0  205.0 |
        | 39.0  235.0 |
        | 42.0  259.0 |
        | 44.0  276.0 |
        | 45.0  285.0 |
        *             *

Output

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

Example 2

This example shows how to solve the system A^TX = B for two right-hand sides, where matrix A is the input matrix factored in "Example 4" for SPOF.

Call Statement and Input

            UPLO  A  LDA  N   B  LDB  NRHS
             |    |   |   |   |   |    |
CALL SPOSM( 'U' , A , 9 , 9 , B , 9 ,  2  )

A        =(same as output A in
"Example 4")

        *             *
        |  9.0   45.0 |
        | 17.0   89.0 |
        | 24.0  131.0 |
        | 30.0  170.0 |
B    =  | 35.0  205.0 |
        | 39.0  235.0 |
        | 42.0  259.0 |
        | 44.0  276.0 |
        | 45.0  285.0 |
        *             *

Output

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

Example 3

This example shows how to solve the system AX = B for two right-hand sides, where matrix A is the same matrix factored in the "Example 5" for CPOF.

Call Statement and Input

            UPLO  A  LDA  N   B  LDB  NRHS
             |    |   |   |   |   |    |
CALL CPOSM( 'L' , A , 3 , 3 , B , 3 ,  2  )

A        =(same as output A in
"Example 5")

        *                                *
        |  (60.0, -55.0)    (70.0, 10.0) |
B    =  |   (34.0, 58.0)  (-51.0, 110.0) |
        | (13.0, -152.0)    (75.0, 63.0) |
        *                                *

Output

        *                          *
        | (2.0, -1.0)   (2.0, 0.0) |
B    =  |  (1.0, 1.0)  (-1.0, 2.0) |
        | (0.0, -2.0)   (1.0, 1.0) |
        *                          *

Example 4

This example shows how to solve the system AX = B for two right-hand sides, where matrix A is the input matrix factored in "Example 6" for CPOF.

Call Statement and Input

            UPLO  A  LDA  N   B  LDB  NRHS
             |    |   |   |   |   |    |
CALL CPOSM( 'U' , A , 3 , 3 , B , 3 ,  2  )

A        =(same as output A in
"Example 6")

        *                              *
        | (33.0, -18.0)   (15.0, -3.0) |
B    =  | (45.0, -45.0)    (8.0, -2.0) |
        |  (152.0, 1.0)  (43.0, -29.0) |
        *                              *

Output

        *                          *
        | (2.0, -1.0)   (2.0, 0.0) |
B    =  | (1.0, -1.0)   (0.0, 1.0) |
        |  (3.0, 0.0)  (1.0, -1.0) |
        *                          *

SPPFCD, DPPFCD, SPOFCD, and DPOFCD--Positive Definite Real Symmetric Matrix Factorization, Condition Number Reciprocal, and Determinant

The SPPFCD and DPPFCD subroutines factor positive definite symmetric matrix A, stored in lower-packed storage mode, using Gaussian elimination (LDL^T). The reciprocal of the condition number and the determinant of matrix A can also be computed. To solve the 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.

The SPOFCD and DPOFCD subroutines factor positive definite symmetric matrix A, stored in upper or lower storage mode, using Cholesky factorization (LL^T or U^TU). The reciprocal of the condition number and the determinant of matrix A can also be computed. To solve the system of equations with one or more right-hand sides, follow the call to these subroutines with a call to SPOSM or DPOSM, respectively. To find the inverse of matrix A, follow the call to these subroutines with a call to SPOICD or DPOICD, respectively.

Table 96. Data Types

A, aux, rcond, det	Subroutine
Short-precision real	SPPFCD and SPOFCD
Long-precision real	DPPFCD and DPOFCD

Note:

The output factorization from SPPFCD and DPPFCD should be used only as input to the solve subroutines SPPS and DPPS, respectively. The output from SPOFCD and DPOFCD should be used only as input to the following subroutines for performing a solve or inverse: SPOSM/SPOICD and DPOSM/DPOICD, respectively.

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 the array, referred to as AP, in which the matrix A, to be factored, is stored in lower-packed storage mode. Specified as: a one-dimensional array of (at least) length n(n+1)/2+n, containing numbers of the data type indicated in Table 96.

a

is the positive definite symmetric matrix A, to be factored. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 96.

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, where:

For SPPFCD and DPPFCD, n >= 0.

For SPOFCD and DPOFCD, 0 <= n <= lda.

iopt

indicates the type of computation to be performed, where:

If iopt = 0, the matrix is factored.

If iopt = 1, the matrix is factored, and the reciprocal of the condition number is computed.

If iopt = 2, the matrix is factored, and the determinant is computed.

If iopt = 3, the matrix is factored and the reciprocal of the condition number and the determinant are computed.

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

rcond

See 'On Return'.

det

See 'On Return'.

aux

has the following meaning:

If naux = 0 and error 2015 is unrecoverable, aux is ignored.

Otherwise, is the storage work area used by these subroutines. Its size is specified by naux. Specified as: an area of storage, containing numbers of the data type indicated in Table 96.

naux

is the size of the work area specified by aux--that is, the number of elements in aux. Specified as: a fullword integer, where:

If naux = 0 and error 2015 is unrecoverable, SPPFCD, DPPFCD, SPOFCD, and DPOFCD dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, naux >= n.

On Return

ap

is the transformed matrix A of order n, containing the results of the factorization. See "Function". Returned as: a one-dimensional array of (at least) length n(n+1)/2+n, containing numbers of the data type indicated in Table 96.

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 96.

rcond

is the estimate of the reciprocal of the condition number, rcond, of matrix A. Returned as: a number of the data type indicated in Table 96; rcond >= 0.

det

is the vector det, containing the two components det₁ and det₂ of the determinant of matrix A. The determinant is:

where 1 <= det₁ < 10. Returned as: an array of length 2, containing numbers of the data type indicated in Table 96.

Notes

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

2. In your C program, argument rcond must be passed by reference.

3. When iopt = 0, SPPFCD and DPPFCD provide the same function as a call to SPPF or DPPF, respectively. When iopt = 0, SPOFCD and DPOFCD provide the same function as a call to SPOF or DPOF, respectively.

4. See "Notes" for information on specifying a value for iopt in the SPPS and DPPS subroutines after calling SPPFCD and DPPFCD, respectively.

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

6. For a description of how a positive definite symmetric matrix is stored in lower-packed storage mode in an array, see "Symmetric Matrix". For a description of how a positive definite symmetric matrix is stored in upper or lower storage mode, see "Positive Definite or Negative Definite Symmetric Matrix".

7. You have the option of having the minimum required value for naux dynamically returned to your program. For details, see "Using Auxiliary Storage in ESSL".

Function

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

For SPPFCD and DPPFCD

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.

An estimate of the reciprocal of the condition number, rcond, and the determinant, det, can also be computed by this subroutine. The estimate of the condition number uses an enhanced version of the algorithm described in references [63] and [64].

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

These subroutines call SPPF and DPPF, respectively, to perform the factorization using Gaussian elimination (LDL^T). If you want to use the Cholesky factorization method, you must call SPPF and DPPF directly.

For SPOFCD and DPOFCD

The positive definite symmetric 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

where:

L is a lower triangular matrix.

L^T is the transpose of matrix L.

U is an upper triangular matrix.

U^T is the transpose of matrix U.

If specified, the estimate of the reciprocal of the condition number and the determinant can also be computed. The estimate of the condition number uses an enhanced version of the algorithm described in references [63] and [64].

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

Error Conditions

Resource Errors

Error 2015 is unrecoverable, naux = 0, and unable to allocate work area.

Computational Errors

1. Matrix A is not positive definite (for SPPFCD and DPPFCD).
   • If matrix A is singular (at least one of the diagonal elements are 0), then rcond and det, if you requested them, are set to 0.
   • If matrix A is nonsingular and nonpositive definite (none of the diagonal elements are 0 and at least one diagonal element is negative), then rcond and det, if you requested them, are computed.
   • 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, issued by SPPF or DPPF, respectively.
   • i can be determined at run time by using 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 be terminated by SPPF or DPPF, respectively, when this error occurs. If your program is not terminated by SPPF or DPPF, respectively, the return code is set to 2. For details, see "What Can You Do about ESSL Computational Errors?".
2. Matrix A is not positive definite (for SPOFCD and DPOFCD).
   • If matrix A is singular (at least one of the diagonal elements are 0), then rcond and det, if you requested them, are set to 0.
   • If matrix A is nonsingular and nonpositive definite (none of the diagonal elements are 0 and at least one diagonal element is negative), then rcond and det, if you requested them, are computed.
   • 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.
   • i can be determined at run time by using 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 be terminated by SPPF or DPPF, respectively, when this error occurs. If your program is not terminated by SPPF or DPPF, respectively, the return code is set to 2. For details, see "What Can You Do about ESSL Computational Errors?".

Input-Argument Errors

1. uplo <> 'U' or 'L'
2. lda <= 0
3. lda < n
4. n < 0
5. iopt <> 0, 1, 2, or 3
6. Error 2015 is recoverable or naux<>0, and naux is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.

Example 1

This example computes the factorization, reciprocal of the condition number, and determinant of matrix A. The input is the same as used in "Example 1" for SPPF.

The values used to estimate the reciprocal of the condition number are obtained with the following values:

||A||₁ = max(9.0, 17.0, 24.0, 30.0, 35.0, 39.0, 42.0, 44.0, 45.0) = 45.0

Estimate of ||A|| = 4.0

On output, the value in det, |A|, is equal to 1.

Call Statement and Input

             AP   N  IOPT  RCOND   DET   AUX  NAUX
             |    |   |      |      |     |    |
CALL DPPFCD( AP , 9 , 3  , RCOND , DET , AUX , 9  )

AP       =(same as input AP in
"Example 1")

Output

AP       =(same as output AP in
"Example 1")
RCOND    =  0.0055555
DET      =  (1.0, 0.0)

Example 2

This example computes the factorization, reciprocal of the condition number, and determinant of matrix A. The input is the same as used in "Example 3" for SPOF.

The values used to estimate the reciprocal of the condition number are obtained with the following values:

||A||₁ = max(9.0, 17.0, 24.0, 30.0, 35.0, 39.0, 42.0, 44.0, 45.0) = 45.0

Estimate of ||A|| = 4.0

On output, the value in det, |A|, is equal to 1.

Call Statement and Input

             UPLO A  LDA  N IOPT  RCOND   DET   AUX  NAUX
              |   |   |   |   |     |      |     |    |
CALL SPOFCD( 'L', A , 9 , 9 , 3 , RCOND , DET , AUX , 9  )

A        =(same as input A in
"Example 3")

Output

A        =(same as output A in
"Example 3")
RCOND    =  0.0055555
DET      =  (1.0, 0.0)

Example 3

This example computes the factorization, reciprocal of the condition number, and determinant of matrix A. The input is the same as used in "Example 4" for SPOF.

The values used to estimate the reciprocal of the condition number are obtained with the following values:

||A||₁ = max(9.0, 17.0, 24.0, 30.0, 35.0, 39.0, 42.0, 44.0, 45.0) = 45.0

Estimate of ||A|| = 4.0

On output, the value in det, |A|, is equal to 1.

Call Statement and Input

             UPLO A  LDA  N IOPT  RCOND   DET   AUX  NAUX
              |   |   |   |   |     |      |     |    |
CALL SPOFCD( 'U', A , 9 , 9 , 3 , RCOND , DET , AUX , 9  )

A        =(same as input A in
"Example 4")

Output

A        =(same as output A in
"Example 4")
RCOND    =  0.0055555
DET      =  (1.0, 0.0)

SGEICD and DGEICD--General Matrix Inverse, Condition Number Reciprocal, and Determinant

These subroutines find the inverse, the reciprocal of the condition number, and the determinant of matrix A.

Table 97. Data Types

A, aux, rcond, det	Subroutine
Short-precision real	SGEICD
Long-precision real	DGEICD

Note:

If you call these subroutines with iopt = 0, 1, 2, or 3 the input must be the output from the factorization subroutines SGEF/SGEFCD/SGETRF or DGEF/DGEFCD/DGEFP/DGETRF, respectively.

Syntax

On Entry

a

has the following meaning, where:

If iopt = 0, 1, 2, or 3, it is matrix A of order n, whose inverse, reciprocal of condition number, and determinant are computed.

If iopt = 4, it is the transformed matrix A of order n, resulting from the factorization performed in a previous call to SGEF/SGEFCD or DGEF/DGEFCD/DGEFP, respectively, whose inverse is computed.

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

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 of matrix A. Specified as: a fullword integer; 0 <= n <= lda.

iopt

indicates the type of computation to be performed, where:

If iopt = 0, the inverse is computed for matrix A.

If iopt = 1, the inverse and the reciprocal of the condition number are computed for matrix A.

If iopt = 2, the inverse and the determinant are computed for matrix A.

If iopt = 3, the inverse, the reciprocal of the condition number, and the determinant are computed for matrix A.

If iopt = 4, the inverse is computed using the factored matrix A.

Specified as: a fullword integer; iopt = 0, 1, 2, 3, 4.

rcond

See 'On Return'.

det

See 'On Return'.

aux

has the following meaning, and its size is specified by naux:

If iopt = 0, 1, 2, or 3, then if naux = 0 and error 2015 is unrecoverable, aux is ignored. Otherwise, it is the storage work area used by this subroutine.

If iopt = 4, aux has the following meaning:

• For SGEICD, the first n locations in aux must contain the ipvt integer vector of length n, resulting from a previous call to SGEF, SGETRF, or SGEFCD.

• For DGEICD, the first ceiling(n/2) locations in aux must contain the ipvt integer vector of length n, resulting from a previous call to DGEF, DGETRF, DGEFCD, or DGEFP.

Specified as: an area of storage, containing numbers of the data type indicated in Table 97.

naux

is the size of the work area specified by aux--that is, the number of elements in aux. Specified as: a fullword integer, where:

If iopt <> 4, then if naux = 0 and error 2015 is unrecoverable, SGEICD and DGEICD dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise naux must have the following value:

For the RS/6000 POWER or PowerPC processors, naux >= 100n.

For the RS/6000 POWER2 processors, naux >= 200n.

Note:

naux values specified for releases prior to ESSL Version 2 Release 2 will still work, but you may not achieve optimal performance.

On Return

a

is the resulting inverse of matrix A of order n. Returned as: an lda by (at least) n array, containing numbers of the data type indicated in Table 97.

rcond

is the reciprocal of the condition number, rcond, of matrix A. Returned as: a real number of the data type indicated in Table 97; rcond >= 0.

det

is the vector det, containing the two components det₁ and det₂ of the determinant of matrix A. The determinant is:

where 1 <= det₁ < 10. Returned as: an array of length 2, containing numbers of the data type indicated in Table 97.

Notes

1. In your C program, argument rcond must be passed by reference.

2. If iopt = 4, the following input arguments for SGEICD and DGEICD must be set to the same values in the previous call to SGEF/SGEFCD or DGEF/DGEFCD/DGEFP, respectively:
   

   For _GEF_	For _GEICD
   Input arguments n and lda	Input arguments n and lda
   Output arguments a and ipvt	Input arguments a and aux

3. You have the option of having the value for naux dynamically returned to your program. For details, see "Using Auxiliary Storage in ESSL".

Function

The inverse, the reciprocal of the condition number, and the determinant of a general square matrix A are computed using partial pivoting to preserve accuracy, where:

• A^-1 is the inverse of matrix A, where AA^-1 = A^-1A = I, and I is the identity matrix.

• 1/(||A||₁)(||A^-1||₁) is the reciprocal of the condition number, where ||A||₁ is the one-norm of matrix A.

• |A| is the determinant of matrix A, where |A| is expressed as:

  
  
  
  
  

The iopt argument is used to determine the combination of output items produced by this subroutine: the inverse, the reciprocal of the condition number, and the determinant.

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

Error Conditions

Resource Errors

If iopt = 0, 1, 2, or 3, then error 2015 is unrecoverable, naux = 0, and unable to allocate work area.

Computational Errors

Matrix A is singular or nearly singular.

• The index i of the first pivot element having a value equal to 0, is identified in the computational error message.
• These subroutines return 0 for rcond and det, if you requested them.
• The return code is set to 2.
• 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 2105 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. lda <= 0
2. n < 0
3. n > lda
4. iopt <> 0, 1, 2, 3, or 4
5. Error 2015 is recoverable or naux<>0, and naux is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.

Example 1

This example computes the inverse, the reciprocal of the condition number, and the determinant of matrix A. The values used to compute the reciprocal of the condition number in this example are obtained with the following values:

||A||₁ = max(6.0, 8.0, 10.0, 12.0, 13.0, 14.0, 15.0, 15.0, 15.0) = 15.0

||A^-1||₁ = 1226.33

On output, the value in det, |A|, is equal to 336.

Call Statement and Input

             A  LDA  N  IOPT  RCOND   DET   AUX   NAUX
             |   |   |   |      |      |     |     |
CALL DGEICD( A , 9 , 9 , 3  , RCOND , DET , AUX , 293 )
 
        *                                                *
        | 1.0  1.0  1.0  1.0  0.0  0.0   0.0   0.0   0.0 |
        | 1.0  1.0  1.0  1.0  1.0  0.0   0.0   0.0   0.0 |
        | 4.0  1.0  1.0  1.0  1.0  1.0   0.0   0.0   0.0 |
        | 0.0  5.0  1.0  1.0  1.0  1.0   1.0   0.0   0.0 |
A    =  | 0.0  0.0  6.0  1.0  1.0  1.0   1.0   1.0   0.0 |
        | 0.0  0.0  0.0  7.0  1.0  1.0   1.0   1.0   1.0 |
        | 0.0  0.0  0.0  0.0  8.0  1.0   1.0   1.0   1.0 |
        | 0.0  0.0  0.0  0.0  0.0  9.0   1.0   1.0   1.0 |
        | 0.0  0.0  0.0  0.0  0.0  0.0  10.0  11.0  12.0 |
        *                                                *

Output

        *                                                                      *
        |    0.333   -0.667   0.333  0.000  0.000  0.000   0.042 -0.042  0.000 |
        |   56.833  -52.167  -1.167 -0.500 -0.500 -0.357   6.836 -0.479 -0.500 |
        |  -55.167   51.833   0.833  0.500  0.500  0.214  -6.735  0.521  0.500 |
        |   -1.000    1.000   0.000  0.000  0.000  0.143  -0.143  0.000  0.000 |
A    =  |   -1.000    1.000   0.000  0.000  0.000  0.000   0.000  0.000  0.000 |
        |   -1.000    1.000   0.000  0.000  0.000  0.000  -0.125  0.125  0.000 |
        | -226.000  206.000   5.000  3.000  2.000  1.429 -27.179  1.750  2.000 |
        |  560.000 -520.000 -10.000 -6.000 -4.000 -2.857  67.857 -5.000 -5.000 |
        | -325.000  305.000   5.000  3.000  2.000  1.429 -39.554  3.125  3.000 |
        *                                                                      *
 
RCOND    =  0.00005436
DET      =  (3.36, 2.00)

Example 2

This example computes the inverse of matrix A, where iopt = 4 and matrix A is the transformed matrix factored in "Example 1" by SGEF. The input contents of AUX, shown here, is the same as the output contents of IPVT in that example.

Call Statement and Input

             A  LDA  N  IOPT  RCOND   DET   AUX   NAUX
             |   |   |   |      |      |     |     |
CALL SGEICD( A , 9 , 9 , 4  , RCOND , DET , AUX , 300 )

A        =(same as output A in
"Example 1")
AUX      =  (3, 4, 5, 6, 7, 8, 9, 8, 9)

Output

        *                                                                      *
        |    0.333   -0.667   0.333  0.000  0.000  0.000   0.042 -0.042  0.000 |
        |   56.833  -52.167  -1.167 -0.500 -0.500 -0.357   6.836 -0.479 -0.500 |
        |  -55.167   51.833   0.833  0.500  0.500  0.214  -6.735  0.521  0.500 |
        |   -1.000    1.000   0.000  0.000  0.000  0.143  -0.143  0.000  0.000 |
A    =  |   -1.000    1.000   0.000  0.000  0.000  0.000   0.000  0.000  0.000 |
        |   -1.000    1.000   0.000  0.000  0.000  0.000  -0.125  0.125  0.000 |
        | -226.000  206.000   5.000  3.000  2.000  1.429 -27.179  1.750  2.000 |
        |  560.000 -520.000 -10.000 -6.000 -4.000 -2.857  67.857 -5.000 -5.000 |
        | -325.000  305.000   5.000  3.000  2.000  1.429 -39.554  3.125  3.000 |
        *                                                                      *

SPPICD, DPPICD, SPOICD, and DPOICD--Positive Definite Real Symmetric Matrix Inverse, Condition Number Reciprocal, and Determinant

These subroutines find the inverse, the reciprocal of the condition number, and the determinant of positive definite symmetric matrix A using Cholesky factorization, where:

• For SPPICD and DPPICD, A is stored in lower-packed storage mode.
• For SPOICD and DPOICD, A is stored in upper or lower storage mode.

Table 98. Data Types

A, aux, rcond, det	Subroutine
Short-precision real	SPPICD and SPOICD
Long-precision real	DPPICD and DPOICD

Note:

If you call these subroutines with iopt = 4, the input must be the output from the factorization subroutines SPPF, DPPF, SPOF/SPOFCD, or DPOF/DPOFCD, respectively, where Cholesky factorization was performed.

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 the array, referred to as AP, where:

If iopt = 0, 1, 2, or 3, then AP contains the positive definite real symmetric matrix A, whose inverse, condition number reciprocal, and determinant are computed, where matrix A is stored in lower-packed storage mode.

If iopt = 4, then AP contains the transformed matrix A of order n, resulting from the Cholesky factorization performed in a previous call to SPPF or DPPF, respectively, whose inverse is computed.

Specified as: a one-dimensional array of (at least) length n(n+1)/2, containing numbers of the data type indicated in Table 98.

a

has the following meaning, where:

If iopt = 0, 1, 2, or 3, it is the positive definite real symmetric matrix A, whose inverse, condition number reciprocal, and determinant are computed, where matrix A is stored in upper or lower storage mode.

If iopt = 4, it is the transformed matrix A of order n, containing results of the factorization from a previous call to SPOF/SPOFCD or DPOF/DPOFCD, respectively, whose inverse is computed.

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

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

indicates the type of computation to be performed, where:

If iopt = 0, the inverse is computed for matrix A.

If iopt = 1, the inverse and the reciprocal of the condition number are computed for matrix A.

If iopt = 2, the inverse and the determinant are computed for matrix A.

If iopt = 3, the inverse, the reciprocal of the condition number, and the determinant are computed for matrix A.

If iopt = 4, the inverse is computed for the (Cholesky) factored matrix A.

Specified as: a fullword integer; iopt = 0, 1, 2, 3, or 4.

rcond

See 'On Return'.

det

See 'On Return'.

aux

has the following meaning:

If naux = 0 and error 2015 is unrecoverable, aux is ignored.

Otherwise, it is the storage work area used by this subroutine. Its size is specified by naux. Specified as: an area of storage, containing numbers of the data type indicated in Table 98.

naux

is the size of the work area specified by aux--that is, the number of elements in aux. Specified as: a fullword integer, where:

If naux = 0 and error 2015 is unrecoverable, SPPICD, DPPICD, SPOICD, AND DPOICD dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, naux >= n.

On Return

ap

is the resulting array, referred to as AP, containing the inverse of the matrix in lower-packed storage mode. Returned as: a one-dimensional array of (at least) length n(n+1)/2, containing numbers of the data type indicated in Table 98.

a

is the transformed matrix A of order n, containing the inverse of the matrix in upper or lower storage mode. Returned as: a two-dimensional array, containing numbers of the data type indicated in Table 98.

rcond

is the reciprocal of the condition number, rcond, of matrix A. Returned as: a real number of the data type indicated in Table 98; rcond >= 0.

det

is the vector det, containing the two components det₁ and det₂ of the determinant of matrix A. The determinant is:

where 1 <= det₁ < 10. Returned as: an array of length 2, containing numbers of the data type indicated in Table 98.

Notes

1. For these subroutines, when you specify iopt = 4, you must do the following:

   • For SPPICD and DPPICD, use Cholesky factorization in the previous call to SPPF and DPPF, respectively.

   • For SPOICD and DPOICD, specify the same storage mode for matrix A that was specified in the previous call to SPOF/SPOFCD and DPOF/DPOFCD, respectively.

   • The scalar data specified for input arguments uplo, lda, and n for these subroutines must be the same as the corresponding input arguments specified for SPOF/SPOFCD and DPOF/DPOFCD, respectively.

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

3. In your C program, argument rcond must be passed by reference.

4. For a description of how a positive definite symmetric matrix is stored in lower-packed storage mode in an array, see "Symmetric Matrix". For a description of how a positive definite symmetric matrix is stored in upper or lower storage mode, see "Positive Definite or Negative Definite Symmetric Matrix".

5. You have the option of having the minimum required value for naux dynamically returned to your program. For details, see "Using Auxiliary Storage in ESSL".

Function

These subroutines find the inverse, the reciprocal of the condition number, and the determinant of positive definite symmetric matrix A using Cholesky factorization, where:

• A^-1 is the inverse of matrix A, where AA^-1 = A^-1A = I, and I is the identity matrix.

• 1/(||A||₁)(||A^-1||₁) is the reciprocal of the condition number, where ||A||₁ is the one-norm of matrix A.

• |A| is the determinant of matrix A, where |A| is expressed as:

  
  
  
  
  

The iopt argument is used to determine the combination of output items produced by this subroutine: the inverse, the reciprocal of the condition number, and the determinant.

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

Error Conditions

Resource Errors

• Error 2015 is unrecoverable, naux = 0, and unable to allocate work area.
• Unable to allocate internal work area.

Computational Errors

Matrix A is not positive definite.

• These subroutines do not perform the inverse, determinant, and reciprocal of the condition number computations.
• For iopt = 1, 2, or 3, the leading minor of order i has a nonpositive determinant. The order i is identified in the computational error message, issued by SPPF, DPPF, SPOF, or DPOF, respectively.

  For iopt = 4 for SPPICD and DPPICD, if the Cholesky factorization performed by SPPF or DPPF, respectively, failed due to a nonpositive definite matrix A, the results from STPI or DTPI, respectively, are unpredictable, and a computational error message may be issued.

  For iopt = 4 for SPOICD and DPOICD, if the factorization performed by SPOF/SPOFCD or DPOF/DPOFCD, respectively, failed due to a nonpositive definite matrix A, the results from STRI or DTRI, respectively, are unpredictable, and a computational error message may be issued.

• i can be determined at run time by using 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 be terminated by SPPF, DPPF, SPOF, or DPOF, respectively, when this error occurs. If your program is not terminated by SPPF, DPPF, SPOF, or DPOF, respectively, the return code is set to 2. For details, see "What Can You Do about ESSL Computational Errors?".

Input-Argument Errors

1. uplo <> 'U' or 'L'
2. n < 0
3. lda <= 0
4. lda < n
5. iopt <> 0, 1, 2, 3, or 4
6. Error 2015 is recoverable or naux<>0, and naux is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.

Example 1

This example uses SPPICD to compute the inverse, reciprocal of the condition number, and determinant of matrix A. Where 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 |
             *                                             *

The values used to compute the reciprocal of the condition number in this example are obtained with the following values:

||A||₁ = max(9.0, 17.0, 24.0, 30.0, 35.0, 39.0, 42.0, 44.0, 45.0) = 45.0

||A^-1||₁ = 4.0

On output, the value in det, |A|, is equal to 1, and RCOND = 1/180.

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 RCOND   DET   AUX  NAUX
             |    |   |     |      |     |    |
CALL SPPICD( AP , 9 , 3 , RCOND , DET , AUX , 9  )
 
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   =   (2.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0,
          2.0, -1.0, 0.0,
          2.0, -1.0,
          1.0)
 
RCOND    =  0.005556
DET      =  (1.0, 0.0)

Example 2

This example uses SPPICD to compute the inverse of matrix A, where iopt = 4, and matrix A is the transformed matrix factored in "Example 1" by SPPF.

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 RCOND   DET   AUX NAUX
            |    |   |     |      |     |    |
CALL SPPICD(AP , 9 , 4 , RCOND , DET , AUX , 9)

AP       =(same as output AP in
"Example 2" for SPPF)

Output

AP   =   (2.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0, 0.0,
          2.0, -1.0, 0.0, 0.0,
          2.0, -1.0, 0.0,
          2.0, -1.0,
          1.0)

Example 3

This example uses SPOICD to compute the inverse, reciprocal of the condition number, and determinant of the same matrix A used in Example 1; however, matrix A is stored in upper storage mode in this example.

The values used to compute the reciprocal of the condition number in this example are obtained with the following values:

||A||₁ = max(9.0, 17.0, 24.0, 30.0, 35.0, 39.0, 42.0, 44.0, 45.0) = 45.0

||A^-1||₁ = 4.0

On output, the value in det, |A|, is equal to 1, and RCOND = 1/180.

Call Statement and Input

             UPLO  A  LDA   N  IOPT RCOND   DET   AUX  NAUX
              |    |   |    |   |     |      |     |    |
CALL SPOICD( 'U' , A , 9 ,  9 , 3 , RCOND , DET , AUX , 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

        *                                                     *
        | 2.0  -1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0 |
        |  .    2.0  -1.0   0.0   0.0   0.0   0.0   0.0   0.0 |
        |  .     .    2.0  -1.0   0.0   0.0   0.0   0.0   0.0 |
        |  .     .     .    2.0  -1.0   0.0   0.0   0.0   0.0 |
A    =  |  .     .     .     .    2.0  -1.0   0.0   0.0   0.0 |
        |  .     .     .     .     .    2.0  -1.0   0.0   0.0 |
        |  .     .     .     .     .     .    2.0  -1.0   0.0 |
        |  .     .     .     .     .     .     .    2.0  -1.0 |
        |  .     .     .     .     .     .     .     .    1.0 |
        *                                                     *
 
RCOND    =  0.005555556
DET      =  (1.0, 0.0)

Example 4

This example uses SPOICD to compute the inverse of matrix A, where iopt = 4, and matrix A is the transformed matrix factored in "Example 1" by SPOF.

Call Statement and Input

             UPLO  A  LDA   N  IOPT RCOND   DET   AUX  NAUX
              |    |   |    |   |     |      |     |    |
CALL SPOICD( 'U' , A , 9 ,  9 , 4 , RCOND , DET , AUX , 9  )

A        =(same as output A in
"Example 4" for SPOF)

Output

        *                                                     *
        | 2.0  -1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0 |
        |  .    2.0  -1.0   0.0   0.0   0.0   0.0   0.0   0.0 |
        |  .     .    2.0  -1.0   0.0   0.0   0.0   0.0   0.0 |
        |  .     .     .    2.0  -1.0   0.0   0.0   0.0   0.0 |
A    =  |  .     .     .     .    2.0  -1.0   0.0   0.0   0.0 |
        |  .     .     .     .     .    2.0  -1.0   0.0   0.0 |
        |  .     .     .     .     .     .    2.0  -1.0   0.0 |
        |  .     .     .     .     .     .     .    2.0  -1.0 |
        |  .     .     .     .     .     .     .     .    1.0 |
        *                                                     *

STRSV, DTRSV, CTRSV, ZTRSV, STPSV, DTPSV, CTPSV, and ZTPSV--Solution of a Triangular System of Equations with a Single Right-Hand Side

STRSV, DTRSV, STPSV, and DTPSV perform one of the following solves for a triangular system of equations with a single right-hand side, using the vector x and triangular matrix A or its transpose:

Solution	Equation
1. x <-- A-1x	Ax = b
2. x <-- A-Tx	ATx = b

CTRSV, ZTRSV, CTPSV, and ZTPSV perform one of the following solves for a triangular system of equations with a single right-hand side, using the vector x and and triangular matrix A, its transpose, or its conjugate transpose:

Solution	Equation
1. x <-- A-1x	Ax = b
2. x <-- A-Tx	ATx = b
3. x <-- A-Tx	AHx = b

Matrix A can be either upper or lower triangular, where:

• For the _TRSV subroutines, it is stored in upper- or lower-triangular storage mode, respectively.

• For the _TPSV subroutines, it is stored in upper- or lower-triangular-packed storage mode, respectively.

Note:

The term b used in the systems of equations listed above represents the right-hand side of the system. It is important to note that in these subroutines the right-hand side of the equation is actually provided in the input-output argument x.

Table 99. Data Types

A, x	Subroutine
Short-precision real	STRSV and STPSV
Long-precision real	DTRSV and DTPSV
Short-precision complex	CTRSV and CTPSV
Long-precision complex	ZTRSV and ZTPSV

Syntax

On Entry

uplo

indicates whether matrix A is an upper or lower triangular matrix, where:

If uplo = 'U', A is an upper triangular matrix.

If uplo = 'L', A is a lower triangular matrix.

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

transa

indicates the form of matrix A used in the system of equations, where:

If transa = 'N', A is used, resulting in solution 1.

If transa = 'T', A^T is used, resulting in solution 2.

If transa = 'C', A^H is used, resulting in solution 3.

Specified as: a single character. It must be 'N', 'T', or 'C'.

diag

indicates the characteristics of the diagonal of matrix A, where:

If diag = 'U', A is a unit triangular matrix.

If diag = 'N', A is not a unit triangular matrix.

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

n

is the order of triangular matrix A. Specified as: a fullword integer; n >= 0 and n <= lda.

a

is the upper or lower triangular matrix A of order n, stored in upper- or lower-triangular storage mode, respectively. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 99.

lda

is the leading dimension of the array specified for a. Specified as: a fullword integer; lda > 0 and lda >= n.

ap

is the upper or lower triangular matrix A of order n, stored in upper- or lower-triangular-packed storage mode, respectively. Specified as: a one-dimensional array of (at least) length n(n+1)/2, containing numbers of the data type indicated in Table 99.

x

is the vector x of length n, containing the right-hand side of the triangular system to be solved. Specified as: a one-dimensional array of (at least) length 1+(n-1)|incx|, containing numbers of the data type indicated in Table 99.

incx

is the stride for vector x. Specified as: a fullword integer; incx > 0 or incx < 0.

On Return

x

is the solution vector x of length n, containing the results of the computation. Returned as: a one-dimensional array, containing numbers of the data type indicated in Table 99.

Notes

1. These subroutines accept lowercase letters for the uplo, transa, and diag arguments.

2. For STRSV, DTRSV, STPSV, and DTPSV, if you specify 'C' for the transa argument, it is interpreted as though you specified 'T'.

3. Matrix A and vector x must have no common elements; otherwise, results are unpredictable.

4. ESSL assumes certain values in your array for parts of a triangular matrix. As a result, you do not have to set these values. For unit diagonal matrices, the elements of the diagonal are assumed to be 1.0 for real matrices and (1.0, 0.0) for complex matrices. When using upper- or lower-triangular storage, the unreferenced elements in the lower and upper triangular part, respectively, are assumed to be zero.

5. For a description of triangular matrices and how they are stored in upper- and lower-triangular storage mode and in upper- and lower-triangular-packed storage mode, see "Triangular Matrix".

Function

These subroutines solve a triangular system of equations with a single right-hand side. The solution x may be any of the following, where triangular matrix A, its transpose, or its conjugate transpose is used, and where A can be either upper- or lower-triangular:

1. x <-- A^-1x

2. x <-- A^-Tx

3. x <-- A^-Tx (only for CTRSV, ZTRSV, CTPSV, and ZTPSV)

where:

x is a vector of length n.

A is an upper or lower triangular matrix of order n. For _TRSV, it is stored in upper- or lower-triangular storage mode, respectively. For _TPSV, it is stored in upper- or lower-triangular-packed storage mode, respectively.

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

Error Conditions

Computational Errors

None

Input-Argument Errors

1. uplo <> 'L' or 'U'
2. transa <> 'T', 'N', or 'C'
3. diag <> 'N' or 'U'
4. n < 0
5. lda <= 0
6. lda < n
7. incx = 0

Example 1

This example shows the solution x <-- A^-1x. Matrix A is a real 4 by 4 lower unit triangular matrix, stored in lower-triangular storage mode. Vector x is a vector of length 4.

Note:

Because matrix A is unit triangular, the diagonal elements are not referenced. ESSL assumes a value of 1.0 for the diagonal elements.

Call Statement and Input

            UPLO TRANSA DIAG  N   A  LDA  X  INCX
             |     |     |    |   |   |   |   |
CALL STRSV( 'L' , 'N' , 'U' , 4 , A , 4 , X , 1  )
 
        *                  *
        |  .    .    .   . |
        | 1.0   .    .   . |
A    =  | 2.0  3.0   .   . |
        | 3.0  4.0  3.0  . |
        *                  *
 
X        =  (1.0, 3.0, 11.0, 24.0)

Output

X        =  (1.0, 2.0, 3.0, 4.0)

Example 2

This example shows the solution x <-- A^-Tx. Matrix A is a real 4 by 4 upper nonunit triangular matrix, stored in upper-triangular storage mode. Vector x is a vector of length 4.

Call Statement and Input

            UPLO TRANSA DIAG  N   A  LDA  X  INCX
             |     |     |    |   |   |   |   |
CALL STRSV( 'U' , 'T' , 'N' , 4 , A , 4 , X , 1  )
 
        *                    *
        | 1.0  2.0  3.0  2.0 |
A    =  |  .   2.0  2.0  5.0 |
        |  .    .   3.0  3.0 |
        |  .    .    .   1.0 |
        *                    *
 
X        =  (5.0, 18.0, 32.0, 41.0)

Output

X        =  (5.0, 4.0, 3.0, 2.0)

Example 3

This example shows the solution x <-- A^-Tx. Matrix A is a complex 4 by 4 upper unit triangular matrix, stored in upper-triangular storage mode. Vector x is a vector of length 4.

Note:

Because matrix A is unit triangular, the diagonal elements are not referenced. ESSL assumes a value of (1.0, 0.0) for the diagonal elements.

Call Statement and Input

            UPLO TRANSA DIAG  N   A  LDA  X  INCX
             |     |     |    |   |   |   |   |
CALL CTRSV( 'U' , 'C' , 'U' , 4 , A , 4 , X , 1  )
 
        *                                     *
        | . (2.0, 2.0) (3.0,  3.0) (2.0, 2.0) |
A    =  | .     .      (2.0,  2.0) (5.0, 5.0) |
        | .     .          .       (3.0, 3.0) |
        | .     .          .           .      |
        *                                     *
 
X        =  ((5.0, 5.0), (24.0, 4.0), (49.0, 3.0), (80.0, 2.0))

Output

X        =  ((5.0, 5.0), (4.0, 4.0), (3.0, 3.0), (2.0, 2.0))

Example 4

This example shows the solution x <-- A^-1x. Matrix A is a real 4 by 4 lower unit triangular matrix, stored in lower-triangular-packed storage mode. Vector x is a vector of length 4. Matrix A is:

                          *                    *
                          | 1.0   .    .    .  |
                          | 1.0  1.0   .    .  |
                          | 2.0  3.0  1.0   .  |
                          | 3.0  4.0  3.0  1.0 |
                          *                    *

Call Statement and Input

            UPLO TRANSA DIAG  N   AP   X  INCX
             |     |     |    |   |    |   |
CALL STPSV( 'L' , 'N' , 'U' , 4 , AP , X , 1  )
 
AP       =  ( . , 1.0, 2.0, 3.0, . , 3.0, 4.0, . , 3.0, . )
X        =  (1.0, 3.0, 11.0, 24.0)

Output

X        =  (1.0, 2.0, 3.0, 4.0)

Example 5

This example shows the solution x <-- A^-Tx. Matrix A is a real 4 by 4 upper nonunit triangular matrix, stored in upper-triangular-packed storage mode. Vector x is a vector of length 4. Matrix A is:

                         *                    *
                         | 1.0  2.0  3.0  2.0 |
                         |  .   2.0  2.0  5.0 |
                         |  .    .   3.0  3.0 |
                         |  .    .    .   1.0 |
                         *                    *

Call Statement and Input

            UPLO TRANSA DIAG  N   AP   X  INCX
             |     |     |    |   |    |   |
CALL STPSV( 'U' , 'T' , 'N' , 4 , AP , X , 1  )
 
AP       =  (1.0, 2.0, 2.0, 3.0, 2.0, 3.0, 2.0, 5.0, 3.0, 1.0)
X        =  (5.0, 18.0, 32.0, 41.0)

Output

X        =  (5.0, 4.0, 3.0, 2.0)

Example 6

This example shows the solution x <-- A^-Tx. Matrix A is a complex 4 by 4 upper unit triangular matrix, stored in upper-triangular-packed storage mode. Vector x is a vector of length 4. Matrix A is:

                *                                                *
                | (1.0, 0.0)  (2.0, 2.0)  (3.0, 3.0)  (2.0, 2.0) |
                |     .       (1.0, 0.0)  (2.0, 2.0)  (5.0, 5.0) |
                |     .           .       (1.0, 0.0)  (3.0, 3.0) |
                |     .           .           .       (1.0, 0.0) |
                *                                                *

Call Statement and Input

            UPLO TRANSA DIAG  N   AP   X  INCX
             |     |     |    |   |    |   |
CALL CTPSV( 'U' , 'C' , 'U' , 4 , AP , X , 1  )
 
AP       =  ( . , (2.0, 2.0), . , (3.0, 3.0), (2.0, 2.0), . ,
             (2.0, 2.0), (5.0, 5.0), (3.0, 3.0), . )
X        =  ((5.0, 5.0), (24.0, 4.0), (49.0, 3.0), (80.0, 2.0))

Output

X        =  ((5.0, 5.0), (4.0, 4.0), (3.0, 3.0), (2.0, 2.0))

STRSM, DTRSM, CTRSM, and ZTRSM--Solution of Triangular Systems of Equations with Multiple Right-Hand Sides

STRSM and DTRSM perform one of the following solves for a triangular system of equations with multiple right-hand sides, using scalar alpha, rectangular matrix B, and triangular matrix A or its transpose:

Solution	Equation
1. B <-- alpha(A-1)B	AX = alphaB
2. B <-- alpha(A-T)B	ATX = alphaB
3. B <-- alphaB(A-1)	XA = alphaB
4. B <-- alphaB(A-T)	XAT = alphaB

CTRSM and ZTRSM perform one of the following solves for a triangular system of equations with multiple right-hand sides, using scalar alpha, rectangular matrix B, and triangular matrix A, its transpose, or its conjugate transpose:

Solution	Equation
1. B <-- alpha(A-1)B	AX = alphaB
2. B <-- alpha(A-T)B	ATX = alphaB
3. B <-- alphaB(A-1)	XA = alphaB
4. B <-- alphaB(A-T)	XAT = alphaB
5. B <-- alpha(A-T)B	AHX = alphaB
6. B <-- alphaB(A-T)	XAH = alphaB

Note:

The term X used in the systems of equations listed above represents the output solution matrix. It is important to note that in these subroutines the solution matrix is actually returned in the input-output argument b.

Table 100. Data Types

A, B, alpha	Subroutine
Short-precision real	STRSM
Long-precision real	DTRSM
Short-precision complex	CTRSM
Long-precision complex	ZTRSM

Syntax

On Entry

side

indicates whether the triangular matrix A is located to the left or right of rectangular matrix B in the system of equations, where:

If side = 'L', A is to the left of B, resulting in solution 1, 2, or 5.

If side = 'R', A is to the right of B, resulting in solution 3, 4, or 6.

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

uplo

indicates whether matrix A is an upper or lower triangular matrix, where:

If uplo = 'U', A is an upper triangular matrix.

If uplo = 'L', A is a lower triangular matrix.

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

transa

indicates the form of matrix A used in the system of equations, where:

If transa = 'N', A is used, resulting in solution 1 or 3.

If transa = 'T', A^T is used, resulting in solution 2 or 4.

If transa = 'C', A^H is used, resulting in solution 5 or 6.

Specified as: a single character. It must be 'N', 'T', or 'C'.

diag

indicates the characteristics of the diagonal of matrix A, where:

If diag = 'U', A is a unit triangular matrix.

If diag = 'N', A is not a unit triangular matrix.

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

m

is the number of rows in rectangular matrix B, and:

If side = 'L', m is the order of triangular matrix A.

Specified as: a fullword integer, where:

If side = 'L', 0 <= m <= lda and m <= ldb.

If side = 'R', 0 <= m <= ldb.

n

is the number of columns in rectangular matrix B, and:

If side = 'R', n is the order of triangular matrix A.

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

If side = 'R', n <= lda.

alpha

is the scalar alpha. Specified as: a number of the data type indicated in Table 100.

a

is the triangular matrix A, of which only the upper or lower triangular portion is used, where:

If side = 'L', A is order m.

If side = 'R', A is order n.

Specified as: a two-dimensional array, containing numbers of the data type indicated in Table 100, where:

If side = 'L', its size must be lda by (at least) m.

If side = 'R', it size must be lda by (at least) n.

lda

is the leading dimension of the array specified for a. Specified as: a fullword integer; lda > 0, and:

If side = 'L', lda >= m.

If side = 'R', lda >= n.

b

is the m by n rectangular matrix B, which contains the right-hand sides of the triangular system to be solved. Specified as: an ldb by (at least) n array, containing numbers of the data type indicated in Table 100.

ldb

is the leading dimension of the array specified for b. Specified as: a fullword integer; ldb > 0 and ldb >= m.

On Return

b

is the m by n matrix B, containing the results of the computation.

Returned as: an ldb by (at least) n array, containing numbers of the data type indicated in Table 100.

Notes

1. These subroutines accept lowercase letters for the transa, side, diag, and uplo arguments.

2. For STRSM and DTRSM, if you specify 'C' for the transa argument, it is interpreted as though you specified 'T'.

3. Matrices A and B must have no common elements or results are unpredictable.

4. If matrix A is upper triangular (uplo = 'U'), these subroutines refer to only the upper triangular portion of the matrix. If matrix A is lower triangular, (uplo = 'L'), these subroutines refer to only the lower triangular portion of the matrix. The unreferenced elements are assumed to be zero.

5. The elements of the diagonal of a unit triangular matrix are always one, so you do not need to set these values. The ESSL subroutines always assume that the values in these positions are 1.0 for STRSM and DTRSM and (1.0, 0.0) for CTRSM and ZTRSM.

6. For a description of triangular matrices and how they are stored, see "Triangular Matrix".

Function

These subroutines solve a triangular system of equations with multiple right-hand sides. The solution B may be any of the following, where A is a triangular matrix and B is a rectangular matrix:

1. B <-- alpha(A^-1)B

2. B <-- alpha(A^-T)B

3. B <-- alphaB(A^-1)

4. B <-- alphaB(A^-T)

5. B <-- alpha(A^-H)B (only for CTRSM and ZTRSM)

6. B <-- alphaB(A^-H) (only for CTRSM and ZTRSM)

where:

alpha is a scalar.

B is an m by n rectangular matrix.

A is an upper or lower triangular matrix, where:

If side = 'L', it has order m, and equation 1, 2, or 5 is performed.

If side = 'R', it has order n, and equation 3, 4, or 6 is performed.

If n or m is 0, no computation is performed. See references [32] and [36].

Error Conditions

Resource Errors

Unable to allocate internal work area.

Computational Errors

None

Note:

If the triangular matrix A is singular, the results returned by this subroutine are unpredictable, and there may be a divide-by-zero program exception message.

Input-Argument Errors

1. m < 0
2. n < 0
3. lda, ldb <= 0
4. side <> 'L' or 'R'
5. uplo <> 'L' or 'U'
6. transa <> 'T', 'N', or 'C'
7. diag <> 'N' or 'U'
8. side = 'L' and m > lda
9. m > ldb
10. side = 'R' and n > lda

Example 1

This example shows the solution B <-- alpha(A^-1)B, where A is a real 5 by 5 upper triangular matrix that is not unit triangular, and B is a real 5 by 3 rectangular matrix.

Call Statement and Input

            SIDE  UPLO TRANSA  DIAG  M   N   ALPHA   A  LDA  B  LDB
             |     |     |      |    |   |     |     |   |   |   |
CALL STRSM( 'L' , 'U' , 'N'  , 'N' , 5 , 3 ,  1.0  , A , 7 , B , 6 )

        *                             *
        | 3.0  -1.0   2.0   2.0   1.0 |
        |  .   -2.0   4.0  -1.0   3.0 |
        |  .     .   -3.0   0.0   2.0 |
A    =  |  .     .     .    4.0  -2.0 |
        |  .     .     .     .    1.0 |
        |  .     .     .     .     .  |
        |  .     .     .     .     .  |
        *                             *

        *                    *
        |   6.0  10.0   -2.0 |
        | -16.0  -1.0    6.0 |
B    =  |  -2.0   1.0   -4.0 |
        |  14.0   0.0  -14.0 |
        |  -1.0   2.0    1.0 |
        |    .     .      .  |
        *                    *

Output

        *                 *
        |  2.0  3.0   1.0 |
        |  5.0  5.0   4.0 |
B    =  |  0.0  1.0   2.0 |
        |  3.0  1.0  -3.0 |
        | -1.0  2.0   1.0 |
        |   .    .     .  |
        *                 *

Example 2

This example shows the solution B <-- alpha(A^-T)B, where A is a real 5 by 5 upper triangular matrix that is not unit triangular, and B is a real 5 by 4 rectangular matrix.

Call Statement and Input

            SIDE  UPLO TRANSA  DIAG  M   N    ALPHA  A  LDA  B  LDB
             |     |     |      |    |   |     |     |   |   |   |
CALL STRSM( 'L' , 'U' , 'T'  , 'N' , 5 , 4 ,  1.0  , A , 7 , B , 6 )

        *                              *
        | -1.0  -4.0  -2.0   2.0   3.0 |
        |   .   -2.0   2.0   2.0   2.0 |
        |   .     .   -3.0  -1.0   4.0 |
A    =  |   .     .     .    1.0   0.0 |
        |   .     .     .     .   -2.0 |
        |   .     .     .     .     .  |
        |   .     .     .     .     .  |
        *                              *

        *                          *
        | -1.0  -2.0   -3.0   -4.0 |
        |  2.0  -2.0  -14.0  -12.0 |
B    =  | 10.0   5.0   -8.0   -7.0 |
        | 14.0  15.0    1.0    8.0 |
        | -3.0   4.0    3.0   16.0 |
        |   .     .      .      .  |
        *                          *

Output

        *                        *
        |  1.0   2.0   3.0   4.0 |
        |  3.0   3.0  -1.0   2.0 |
B    =  | -2.0  -1.0   0.0   1.0 |
        |  4.0   4.0  -3.0  -3.0 |
        |  2.0   2.0   2.0   2.0 |
        |   .     .     .     .  |
        *                        *

Example 3

This example shows the solution B <-- alphaB(A^-1), where A is a real 5 by 5 lower triangular matrix that is not unit triangular, and B is a real 3 by 5 rectangular matrix.

Call Statement and Input

            SIDE  UPLO TRANSA  DIAG  M   N    ALPHA  A  LDA  B  LDB
             |     |     |      |    |   |     |     |   |   |   |
CALL STRSM( 'R' , 'L' , 'N'  , 'N' , 3 , 5 ,  1.0  , A , 7 , B , 4 )

        *                            *
        | 2.0   .     .     .     .  |
        | 2.0  3.0    .     .     .  |
        | 2.0  1.0   1.0    .     .  |
A    =  | 0.0  3.0   0.0  -2.0    .  |
        | 2.0  4.0  -1.0   2.0  -1.0 |
        |  .    .     .     .     .  |
        |  .    .     .     .     .  |
        *                            *

        *                             *
        | 10.0   4.0   0.0  0.0   1.0 |
B    =  | 10.0  14.0  -4.0  6.0  -3.0 |
        | -8.0   2.0  -5.0  4.0  -2.0 |
        |   .     .     .    .     .  |
        *                             *

Output

        *                              *
        |  3.0   4.0  -1.0  -1.0  -1.0 |
B    =  |  2.0   1.0  -1.0   0.0   3.0 |
        | -2.0  -1.0  -3.0   0.0   2.0 |
        |   .     .     .     .     .  |
        *                              *

Example 4

This example shows the solution B <-- alphaB(A^-1), where A is a real 6 by 6 upper triangular matrix that is unit triangular, and B is a real 1 by 6 rectangular matrix.

Note:

Because matrix A is unit triangular, the diagonal elements are not referenced. ESSL assumes a value of 1.0 for the diagonal element.

Call Statement and Input

            SIDE  UPLO TRANSA  DIAG  M   N    ALPHA  A  LDA  B  LDB
             |     |     |      |    |   |     |     |   |   |   |
CALL STRSM( 'R' , 'U' , 'N'  , 'U' , 1 , 6 ,  1.0  , A , 7 , B , 2 )

        *                               *
        | .  2.0  -3.0  1.0   2.0   4.0 |
        | .   .    0.0  1.0   1.0  -2.0 |
        | .   .     .   4.0  -1.0   1.0 |
A    =  | .   .     .    .    0.0  -1.0 |
        | .   .     .    .     .    2.0 |
        | .   .     .    .     .     .  |
        | .   .     .    .     .     .  |
        *                               *

        *                                 *
B    =  | 1.0  4.0  -2.0  10.0  2.0  -6.0 |
        |  .    .     .     .    .     .  |
        *                                 *

Output

        *                                *
B    =  | 1.0  2.0  1.0  3.0  -1.0  -2.0 |
        |  .    .    .    .     .     .  |
        *                                *

Example 5

This example shows the solution B <-- alphaB(A^-1), where A is a complex 5 by 5 lower triangular matrix that is not unit triangular, and B is a complex 3 by 5 rectangular matrix.

Call Statement and Input

            SIDE  UPLO TRANSA  DIAG  M   N    ALPHA   A  LDA  B  LDB
             |     |     |      |    |   |      |     |   |   |   |
CALL CTRSM( 'R' , 'L' , 'N'  , 'N' , 3 , 5 ,  ALPHA , A , 7 , B , 4 )
 
ALPHA    =  (1.0, 0.0)

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

        *                                                                      *
        | (22.0, -41.0)  (7.0, -26.0)  (9.0, 0.0)  (-15.0, -3.0)  (-15.0, 8.0) |
B    =  | (29.0, -18.0) (24.0, -10.0)  (9.0, 6.0) (-12.0, -24.0) (-19.0, -8.0) |
        |  (-15.0, 2.0) (-3.0, -21.0) (-2.0, 4.0)  (-4.0, -12.0) (-10.0, -6.0) |
        |        .           .             .            .            .         |
        *                                                                      *

Output

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

Example 6

This example shows the solution B <-- alpha(A^-T)B, where A is a complex 5 by 5 upper triangular matrix that is not unit triangular, and B is a complex 5 by 1 rectangular matrix.

Call Statement and Input

            SIDE  UPLO TRANSA  DIAG  M   N    ALPHA   A  LDA  B  LDB
             |     |     |      |    |   |      |     |   |   |   |
CALL CTRSM( 'L' , 'U' , 'C'  , 'N' , 5 , 1 ,  ALPHA , A , 6 , B , 6 )
 
ALPHA    =  (1.0, 0.0)

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

        *               *
        | (-8.0, -19.0) |
        |  (8.0,  21.0) |
B    =  | (44.0,  -8.0) |
        | (13.0,  -7.0) |
        | (19.0,   2.0) |
        |      .        |
        *               *

Output

        *             *
        |  (3.0, 4.0) |
        | (-4.0, 2.0) |
B    =  | (-5.0, 0.0) |
        |  (1.0, 3.0) |
        |  (3.0, 1.0) |
        |      .      |
        *             *

STRI, DTRI, STPI, and DTPI--Triangular Matrix Inverse

These subroutines find the inverse of triangular matrix A:

A <-- A^-1

Matrix A can be either upper or lower triangular, where:

• For the _TRI subroutines, it is stored in upper- or lower-triangular storage mode, respectively.

• For the _TPI subroutines, it is stored in upper- or lower-triangular-packed storage mode, respectively.

Table 101. Data Types

A	Subroutine
Short-precision real	STRI and STPI
Long-precision real	DTRI and DTPI

Syntax

On Entry

uplo

indicates whether matrix A is an upper or lower triangular matrix, where:

If uplo = 'U', A is an upper triangular matrix.

If uplo = 'L', A is a lower triangular matrix.

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

diag

indicates the characteristics of the diagonal of matrix A, where:

If diag = 'U', A is a unit triangular matrix.

If diag = 'N', A is not a unit triangular matrix.

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

a

is the upper or lower triangular matrix A of order n, stored in upper- or lower-triangular storage mode, respectively. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 101.

lda

is the leading dimension of the arrays specified for a. Specified as: a fullword integer; lda > 0 and lda >= n.

ap

is the upper or lower triangular matrix A of order n, stored in upper- or lower-triangular-packed storage mode, respectively. Specified as: a one-dimensional array of (at least) length n(n+1)/2, containing numbers of the data type indicated in Table 101.

n

is the order of matrix A. Specified as: a fullword integer; n >= 0, where:

On Return

a

is the inverse of the upper or lower triangular matrix A of order n, stored in upper- or lower-triangular storage mode, respectively. Returned as: an lda by (at least) n array, containing numbers of the data type indicated in Table 101.

ap

is the inverse of the upper or lower triangular matrix A of order n, stored in upper- or lower-triangular-packed storage mode, respectively. Returned as: a one-dimensional array of (at least) length n(n+1)/2, containing numbers of the data type indicated in Table 101.

Notes

1. These subroutines accept lowercase letters for the uplo and diag arguments.

2. If matrix A is upper triangular (uplo = 'U'), these subroutines refer to only the upper triangular portion of the matrix. If matrix A is lower triangular, (uplo = 'L'), these subroutines refer to only the lower triangular portion of the matrix. The unreferenced elements are assumed to be zero.

3. The elements of the diagonal of a unit triangular matrix are always one, so you do not need to set these values.

4. For a description of triangular matrices and how they are stored in upper- and lower-triangular storage mode and in upper- and lower-triangular-packed storage mode, see "Triangular Matrix".

Function

These subroutines find the inverse of triangular matrix A, where A is either upper or lower triangular:

A <-- A^-1

where:

A is the triangular matrix of order n.

A^-1 the inverse of the triangular matrix of order n.

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

Error Conditions

Resource Errors

Unable to allocate internal work area.

Computational Errors

Matrix A is singular.

• One or more of the diagonal elements of matrix A are zero. The first column, i, of matrix A, in which a zero diagonal element is found, 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 2145 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. uplo <> 'U' or 'L'
2. diag <> 'U' or 'N'
3. n < 0
4. lda <= 0
5. lda < n

Example 1

This example shows how the inverse of matrix A is computed, where A is a 5 by 5 upper triangular matrix that is not unit triangular and is stored in upper-triangular storage mode. Matrix A is:

                    *                                *
                    | 1.00  3.00  4.00   5.00   6.00 |
                    | 0.00  2.00  8.00   9.00   1.00 |
                    | 0.00  0.00  4.00   8.00   4.00 |
                    | 0.00  0.00  0.00  -2.00   6.00 |
                    | 0.00  0.00  0.00   0.00  -1.00 |
                    *                                *

and where the following inverse matrix is computed. Matrix A^-1 is:

                  *                                   *
                  | 1.00  -1.50   2.00   3.75   35.00 |
                  | 0.00   0.50  -1.00  -1.75  -14.00 |
                  | 0.00   0.00   0.25   1.00    7.00 |
                  | 0.00   0.00   0.00  -0.50   -3.00 |
                  | 0.00   0.00   0.00   0.00   -1.00 |
                  *                                   *

Call Statement and Input

           UPLO  DIAG  A  LDA  N
            |     |    |   |   |
CALL STRI( 'U' , 'N' , A , 5 , 5)

        *                                *
        | 1.00  3.00  4.00   5.00   6.00 |
        |  .    2.00  8.00   9.00   1.00 |
A    =  |  .     .    4.00   8.00   4.00 |
        |  .     .     .    -2.00   6.00 |
        |  .     .     .      .    -1.00 |
        *                                *

Output

        *                                   *
        | 1.00  -1.50   2.00   3.75   35.00 |
        |  .     0.50  -1.00  -1.75  -14.00 |
A    =  |  .      .     0.25   1.00    7.00 |
        |  .      .      .    -0.50   -3.00 |
        |  .      .      .      .     -1.00 |
        *                                   *

Example 2

This example shows how the inverse of matrix A is computed, where A is a 5 by 5 lower triangular matrix that is unit triangular and is stored in lower-triangular storage mode. Matrix A is:

                       *                         *
                       | 1.0  0.0  0.0  0.0  0.0 |
                       | 3.0  1.0  0.0  0.0  0.0 |
                       | 4.0  8.0  1.0  0.0  0.0 |
                       | 5.0  9.0  8.0  1.0  0.0 |
                       | 6.0  1.0  4.0  6.0  1.0 |
                       *                         *

and where the following inverse matrix is computed. Matrix A^-1 is:

                   *                                 *
                   |    1.0     0.0   0.0   0.0  0.0 |
                   |   -3.0     1.0   0.0   0.0  0.0 |
                   |   20.0    -8.0   1.0   0.0  0.0 |
                   | -138.0    55.0  -8.0   1.0  0.0 |
                   |  745.0  -299.0  44.0  -6.0  1.0 |
                   *                                 *

Call Statement and Input

           UPLO  DIAG  A  LDA  N
            |     |    |   |   |
CALL STRI( 'L' , 'U' , A , 5 , 5)

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

Output

        *                               *
        |     .       .     .     .   . |
        |   -3.0      .     .     .   . |
A    =  |   20.0    -8.0    .     .   . |
        | -138.0    55.0  -8.0    .   . |
        |  745.0  -299.0  44.0  -6.0  . |
        *                               *

Example 3

This example shows how the inverse of matrix A is computed, where A is the same matrix shown in Example 1 and is stored in upper-triangular-packed storage mode. The inverse matrix computed here is the same as the inverse matrix shown in Example 1 and is stored in upper-triangular-packed storage mode.

Call Statement and Input

           UPLO  DIAG  AP   N
            |     |    |    |
CALL STPI( 'U' , 'N' , AP , 5)

AP       =  (1.00, 3.00, 2.00, 4.00, 8.00, 4.00, 5.00, 9.00, 8.00,
             -2.00, 6.00, 1.00, 4.00, 6.00, -1.00)

Output

AP       =  (1.00, -1.50, 0.50, 2.00, -1.00, 0.25, 3.75, -1.75, 1.00,
             -0.50, 35.00, -14.00, 7.00, -3.00, -1.00)

Example 4

This example shows how the inverse of matrix A is computed, where A is the same matrix shown in Example 2 and is stored in lower-triangular-packed storage mode. The inverse matrix computed here is the same as the inverse matrix shown in Example 2 and is stored in lower-triangular-packed storage mode.

Note:

Because matrix A is unit triangular, the diagonal elements are not referenced. ESSL assumes a value of 1.0 for the diagonal elements.

Call Statement and Input

           UPLO  DIAG  AP   N
            |     |    |    |
CALL STPI( 'L' , 'U' , AP , 5)

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

Output

AP       =  ( . , -3.0, 20.0, -138.0, 745.0, . , -8.0, 55.0, -299.0,
             . , -8.0, 44.0, . , -6.0, . )

Banded Linear Algebraic Equation Subroutines

This section contains the banded linear algebraic equation subroutine descriptions.

SGBF and DGBF--General Band Matrix Factorization

These subroutines factor general band matrix A, stored in general-band storage mode, using Gaussian elimination. To solve the system of equations with one or more right-hand sides, follow the call to these subroutines with one or more calls to SGBS or DGBS, respectively.

Table 102. Data Types

A	Subroutine
Short-precision real	SGBF
Long-precision real	DGBF

Note:

The output from these factorization subroutines should be used only as input to the solve subroutines SGBS and DGBS, respectively.

Syntax

On Entry

agb

is the general band matrix A of order n, stored in general-band storage mode, to be factored. It has an upper band width mu and a lower band width ml. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 102, where lda >= 2ml+mu+16.

lda

is the leading dimension of the array specified for agb. Specified as: a fullword integer; lda > 0 and lda >= 2ml+mu+16.

n

is the order of the matrix A. Specified as: a fullword integer; n > ml and n > mu.

ml

is the lower band width ml of the matrix A. Specified as: a fullword integer; 0 <= ml < n.

mu

is the upper band width mu of the matrix A. Specified as: a fullword integer; 0 <= mu < n.

ipvt

See 'On Return'.

On Return

agb

is the transformed matrix A of order n, containing the results of the factorization. See "Function". Returned as: an lda by (at least) n array, containing numbers of the data type indicated in Table 102.

ipvt

is the integer vector ipvt of length n, containing the pivot information necessary to construct matrix L from the information contained in the output array agb. Returned as: a one-dimensional array of (at least) length n, containing fullword integers.

Notes

1. ipvt is not a permutation vector in the strict sense. It is used to record column interchanges in L due to partial pivoting and to improve performance.

2. The entire lda by n array specified for agb must remain unchanged between calls to the factorization and solve subroutines.

3. This subroutine can be used for tridiagonal matrices (ml = mu = 1); however, the tridiagonal subroutines SGTF/DGTF and SGTS/DGTS are faster.

4. For a description of how a general band matrix is stored in general-band storage mode in an array, see "General Band Matrix".

Function

The general band matrix A, stored in general-band storage mode, is factored using Gaussian elimination with partial pivoting to compute the LU factorization of A, where:

ipvt is a vector containing the pivoting information.

L is a unit lower triangular band matrix.

U is an upper triangular band matrix.

The transformed matrix A contains U in packed format, along with the multipliers necessary to construct, with the help of ipvt, a matrix L, such that A = LU. This factorization can then be used by SGBS or DGBS, respectively, to solve the system of equations. See reference [38].

Error Conditions

Resource Errors

Unable to allocate internal work area.

Computational Errors

Matrix A is singular.

• One or more columns of L and the corresponding diagonal of U contain all zeros (all columns of L are checked). The last column, i, of L with a corresponding U = 0 diagonal element 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 2103 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. lda <= 0
2. ml < 0
3. ml >= n
4. mu < 0
5. mu >= n
6. lda < 2ml+mu+16

Example

This example shows a factorization of a general band matrix A of order 9, with a lower band width of 2 and an upper band width of 3. On input matrix A is:

            *                                                *
            | 1.0  1.0  1.0  1.0  0.0  0.0   0.0   0.0   0.0 |
            | 1.0  1.0  1.0  1.0  1.0  0.0   0.0   0.0   0.0 |
            | 4.0  1.0  1.0  1.0  1.0  1.0   0.0   0.0   0.0 |
            | 0.0  5.0  1.0  1.0  1.0  1.0   1.0   0.0   0.0 |
            | 0.0  0.0  6.0  1.0  1.0  1.0   1.0   1.0   0.0 |
            | 0.0  0.0  0.0  7.0  1.0  1.0   1.0   1.0   1.0 |
            | 0.0  0.0  0.0  0.0  8.0  1.0   1.0   1.0   1.0 |
            | 0.0  0.0  0.0  0.0  0.0  9.0   1.0   1.0   1.0 |
            | 0.0  0.0  0.0  0.0  0.0  0.0  10.0  11.0  12.0 |
            *                                                *

Matrix A is stored in general-band storage mode in the two-dimensional array AGB of size LDA by N, where LDA = 2ml+mu+16 = 23. The array AGB is declared as AGB(1:23,1:9).

Note:

Matrix A is the same matrix used in the examples in subroutines SGEF and DGEF (see "Example 1") and SGEFCD and DGEFCD (see "Example").

Call Statement and Input

           AGB  LDA   N   ML  MU  IPVT )
            |    |    |   |   |    |
CALL SGBF( AGB , 23 , 9 , 2 , 3 , IPVT )

        *                                                                           *
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  1.0000  1.0000  1.0000   1.0000   1.0000   1.0000 |
        | 0.0000  0.0000  1.0000  1.0000  1.0000  1.0000   1.0000   1.0000   1.0000 |
        | 0.0000  1.0000  1.0000  1.0000  1.0000  1.0000   1.0000   1.0000   1.0000 |
        | 1.0000  1.0000  1.0000  1.0000  1.0000  1.0000   1.0000   1.0000  12.0000 |
        | 1.0000  1.0000  1.0000  1.0000  1.0000  1.0000   1.0000  11.0000   0.0000 |
        | 4.0000  5.0000  6.0000  7.0000  8.0000  9.0000  10.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
AGB  =  | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000 |
        *                                                                           *

Output

        *                                                                             *
        | 0.0000  0.0000  0.0000  0.0000   0.0000   1.0000   1.0000   1.0000   1.0000 |
        | 0.0000  0.0000  0.0000  0.0000   1.0000   1.0000   1.0000   1.0000   1.0000 |
        | 0.0000  0.0000  0.0000  1.0000   1.0000   1.0000   1.0000   1.0000   1.0000 |
        | 0.0000  0.0000  1.0000  1.0000   1.0000   1.0000   1.0000   1.0000  12.0000 |
        | 0.0000  1.0000  1.0000  1.0000   1.0000   1.0000   1.0000  11.0000   0.3111 |
        | 0.2500  0.2000  0.1600  0.1400   0.1250   0.1100   0.1000   5.5380  -325.00 |
        | 0.0000  0.1500  0.0000  0.0714   0.0000  -0.0556  -0.0306   0.9385   0.0000 |
        | 0.2500  0.1500  0.1000  0.0714  -0.0714  -0.0694  -0.0194   0.0000   0.0000 |
        | 0.2500  0.0000  0.1000  0.0000   0.0536   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
AGB  =  | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        | 0.0000  0.0000  0.0000  0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 |
        *                                                                             *

IPVT     =  (2, -65534, -131070, -196606, -262142, -327678, -327678,
             -327680, -327680)

SGBS and DGBS--General Band Matrix Solve

These subroutines solve the system Ax = b for x, where A is a general band matrix, and x and b are vectors. They use the results of the factorization of matrix A, produced by a preceding call to SGBF or DGBF, respectively.

Table 103. Data Types

A, b, x	Subroutine
Short-precision real	SGBS
Long-precision real	DGBS

Note:

The input to these solve subroutines must be the output from the factorization subroutines SGBF and DGBF, respectively.

Syntax

On Entry

agb

is the factorization of general band matrix A, produced by a preceding call to SGBF or DGBF. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 103, where lda >= 2ml+mu+16.

lda

is the leading dimension of the array specified for agb. Specified as: a fullword integer; lda > 0 and lda >= 2ml+mu+16.

n

is the order of the matrix A. Specified as: a fullword integer; n > ml and n > mu.

ml

is the lower band width ml of the matrix A. Specified as: a fullword integer; 0 <= ml < n.

mu

is the upper band width mu of the matrix A. Specified as: a fullword integer; 0 <= mu < n.

ipvt

is the integer vector ipvt of length n, produced by a preceding call to SGBF or DGBF. It contains the pivot information necessary to construct matrix L from the information contained in the array specified for agb.

Specified as: a one-dimensional array of (at least) length n, containing fullword integers.

bx

is the vector b of length n, containing the right-hand side of the system. Specified as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 103.

On Return

bx

is the solution vector x of length n, containing the results of the computation. Returned as: a one-dimensional array, containing numbers of the data type indicated in Table 103.

Notes

1. The scalar data specified for input arguments lda, n, ml, and mu for these subroutines must be the same as that specified for SGBF and DGBF, respectively.

2. The array data specified for input arguments agb and ipvt for these subroutines must be the same as the corresponding output arguments for SGBF and DGBF, respectively.

3. The entire lda by n array specified for agb must remain unchanged between calls to the factorization and solve subroutines.

4. The vectors and matrices used in this computation must have no common elements; otherwise, results are unpredictable. See "Concepts".

5. This subroutine can be used for tridiagonal matrices (ml = mu = 1); however, the tridiagonal subroutines, SGTF/DGTF and SGTS/DGTS, are faster.

6. For a description of how a general band matrix is stored in general-band storage mode in an array, see "General Band Matrix".

Function

The real system Ax = b is solved for x, where A is a real general band matrix, stored in general-band storage mode, and x and b are vectors. These subroutines use the results of the factorization of matrix A, produced by a preceding call to SGBF or DGBF, respectively. The transformed matrix A, used by this computation, consists of the upper triangular matrix U and the multipliers necessary to construct L using ipvt, as defined in "Function". See reference [38].

Error Conditions

Computational Errors

Input-Argument Errors

1. lda <= 0
2. ml < 0
3. ml >= n
4. mu < 0
5. mu >= n
6. lda < 2ml+mu+16

Example

This example shows how to solve the system Ax = b, where general band matrix A is the same matrix factored in "Example" for SGBF and DGBF. The input for AGB and IPVT in this example is the same as the output for that example.

Call Statement and Input

           AGB  LDA   N   ML  MU  IPVT   BX
            |    |    |   |   |    |     |
CALL SGBS( AGB , 23 , 9 , 2 , 3 , IPVT , BX )

IPVT     =  (2, -65534, -131070, -196606, -262142, -327678, -327678,
             -327680, -327680)
BX       =  (4.0000, 5.0000, 9.0000, 10.0000, 11.0000, 12.0000,
             12.0000, 12.0000, 33.0000)
AGB      =(same as output AGB in
"Example")

Output

BX       =  (1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
             0.9999, 1.0001)

SPBF, DPBF, SPBCHF, and DPBCHF--Positive Definite Symmetric Band Matrix Factorization

These subroutines factor positive definite symmetric band matrix A, stored in lower-band-packed storage mode, using:

• Gaussian elimination for SPBF and DPBF
• Cholesky factorization for SPBCHF and DPBCHF

To solve the system of equations with one or more right-hand sides, follow the call to these subroutines with one or more calls to SPBS, DPBS, SPBCHS, or DPBCHS, respectively.

Table 104. Data Types

A	Subroutine
Short-precision real	SPBF and SPBCHF
Long-precision real	DPBF and DPBCHF

Notes:

1. The output from these factorization subroutines should be used only as input to the solve subroutines SPBS, DPBS, SPBCHS, and DPBCHS, respectively.

2. For optimal performance:
   • For wide band widths, use _PBCHF.
   • For narrow band widths, use either _PBF or _PBCHF.
   • For very narrow band widths:
     • Use either SPBF or SPBCHF.
     • Use DPBF.

Syntax

On Entry

apb

is the positive definite symmetric band matrix A of order n, stored in lower-band-packed storage mode, to be factored. It has a half band width of m. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 104. See "Notes".

lda

is the leading dimension of the array specified for apb. Specified as: a fullword integer; lda > 0 and lda > m.

n

is the order n of matrix A. Specified as: a fullword integer; n > m.

m

is the half band width of the matrix A. Specified as: a fullword integer; 0 <= m < n.