LFCQS

Computes the RT R Cholesky factorization of a real symmetric positive definite matrix in band symmetric storage mode and estimate its L1condition number.

Required Arguments

A — NCODA + 1 by N array containing the N by N positive definite band coefficient matrix in band symmetric storage mode to be factored. (Input)

NCODA — Number of upper codiagonals of A. (Input)

FACT — NCODA + 1 by N array containing the RTR factorization of the matrix A in band symmetric form. (Output)

If A is not needed, A and FACT can share the same storage locations.

If A is not needed, A and FACT can share the same storage locations.

RCOND — Scalar containing an estimate of the reciprocal of the L1condition number of A. (Output)

Optional Arguments

N — Order of the matrix. (Input)

Default: N = size (A,2).

Default: N = size (A,2).

LDA — Leading dimension of A exactly as specified in the dimension statement of the calling program. (Input)

Default: LDA = size (A,1).

Default: LDA = size (A,1).

LDFACT — Leading dimension of FACT exactly as specified in the dimension statement of the calling program. (Input)

Default: LDFACT = size (FACT,1).

Default: LDFACT = size (FACT,1).

FORTRAN 90 Interface

Generic: CALL LFCQS (A, NCODA, FACT, RCOND [, …])

Specific: The specific interface names are S_LFCQS and D_LFCQS.

FORTRAN 77 Interface

Single: CALL LFCQS (N, A, LDA, NCODA, FACT, LDFACT, RCOND)

Double: The double precision name is DLFCQS.

Description

Routine LFCQS computes an RTR Cholesky factorization and estimates the condition number of a real symmetric positive definite band coefficient matrix. R is an upper triangular band matrix.

The L1condition number of the matrix A is defined to be κ(A) = ∥A∥1∥A-1∥1. Since it is expensive to compute ∥A-1∥1, the condition number is only estimated. The estimation algorithm is the same as used by LINPACK and is described by Cline et al. (1979).

If the estimated condition number is greater than 1/ɛ (where ɛ is machine precision), a warning error is issued. This indicates that very small changes in A can cause very large changes in the solution x. Iterative refinement can sometimes find the solution to such a system.

LFCQS fails if any submatrix of R is not positive definite or if R has a zero diagonal element. These errors occur only if A is very close to a singular matrix or to a matrix which is not positive definite.

The RTR factors are returned in a form that is compatible with routines LFIQS, LFSQS and LFDQS. To solve systems of equations with multiple right-hand-side vectors, use LFCQS followed by either LFIQS or LFSQS called once for each right-hand side. The routine LFDQS can be called to compute the determinant of the coefficient matrix after LFCQS has performed the factorization.

LFCQS is based on the LINPACK routine SPBCO; see Dongarra et al. (1979).

Comments

1. Workspace may be explicitly provided, if desired, by use of L2CQS/DL2CQS. The reference is:

CALL L2CQS (N, A, LDA, NCODA, FACT, LDFACT, RCOND, WK)

The additional argument is:

WK — Work vector of length N.

2. Informational errors

Type | Code | Description |
---|---|---|

3 | 3 | The input matrix is algorithmically singular. |

4 | 2 | The input matrix is not positive definite. |

Example

The inverse of a 4 × 4 symmetric positive definite band matrix with one codiagonal is computed. LFCQS is called to factor the matrix and to check for nonpositive definiteness or ill-conditioning. LFIQS is called to determine the columns of the inverse.

USE LFCQS_INT

USE LFIQS_INT

USE UMACH_INT

USE WRRRN_INT

! Declare variables

INTEGER LDA, LDFACT, N, NCODA, NOUT

PARAMETER (LDA=2, LDFACT=2, N=4, NCODA=1)

REAL A(LDA,N), AINV(N,N), RCOND, FACT(LDFACT,N),&

RES(N), RJ(N)

!

! Set values for A in band symmetric form

!

! A = ( 0.0 1.0 1.0 1.0 )

! ( 2.0 2.5 2.5 2.0 )

!

DATA A/0.0, 2.0, 1.0, 2.5, 1.0, 2.5, 1.0, 2.0/

! Factor the matrix A

CALL LFCQS (A, NCODA, FACT, RCOND)

! Set up the columns of the identity

! matrix one at a time in RJ

RJ = 0.0E0

DO 10 J=1, N

RJ(J) = 1.0E0

! RJ is the J-th column of the identity

! matrix so the following LFIQS

! reference places the J-th column of

! the inverse of A in the J-th column

! of AINV

CALL LFIQS (A, NCODA, FACT, RJ, AINV(:,J), RES)

RJ(J) = 0.0E0

10 CONTINUE

! Print the results

CALL UMACH (2, NOUT)

WRITE (NOUT,99999) RCOND, 1.0E0/RCOND

CALL WRRRN (’AINV’, AINV)

99999 FORMAT (’ RCOND = ’,F5.3,/,’ L1 Condition number = ’,F6.3)

END

Output

RCOND = 0.160

L1 Condition number = 6.239

AINV

1 2 3 4

1 0.6667 -0.3333 0.1667 -0.0833

2 -0.3333 0.6667 -0.3333 0.1667

3 0.1667 -0.3333 0.6667 -0.3333

4 -0.0833 0.1667 -0.3333 0.6667