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

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

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 indices. 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.
4. On both input and output, matrix A conforms to LAPACK format.

Function

Matrix A is factored using Gaussian elimination with partial pivoting (ipvt) to compute the LU factorization of A, where (A=PLU):

L is a unit lower triangular matrix.

U is an upper triangular matrix.

P is the permutation matrix.

On output, the transformed matrix A contains U in the upper triangle and L in the strict lower triangle where ipvt contains the pivots representing permutation P, such that A = PLU.

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 [69] and [70].

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  )

Output