LCHRG

Computes the Cholesky decomposition of a symmetric positive definite matrix with optional column pivoting.

Required Arguments

A — N by N symmetric positive definite matrix to be decomposed. (Input)

Only the upper triangle of A is referenced.

Only the upper triangle of A is referenced.

FACT — N by N matrix containing the Cholesky factor of the permuted matrix in its upper triangle. (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.

Optional Arguments

N — Order of the matrix A. (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).

PIVOT — Logical variable. (Input)

PIVOT = .TRUE. means column pivoting is done. PIVOT = .FALSE. means no pivoting is done.

Default: PIVOT = .TRUE.

PIVOT = .TRUE. means column pivoting is done. PIVOT = .FALSE. means no pivoting is done.

Default: PIVOT = .TRUE.

IPVT — Integer vector of length N containing information that controls the selection of the pivot columns. (Input/Output)

On input, if IPVT(K) > 0, then the K-th column of A is an initial column; if IPVT(K) = 0, then the K-th column of A is a free column; if IPVT(K) < 0, then the K-th column of A is a final column. See Comments. On output, IPVT(K) contains the index of the diagonal element of A that was moved into the K-th position. IPVT is only referenced when PIVOT is equal to .TRUE..

On input, if IPVT(K) > 0, then the K-th column of A is an initial column; if IPVT(K) = 0, then the K-th column of A is a free column; if IPVT(K) < 0, then the K-th column of A is a final column. See Comments. On output, IPVT(K) contains the index of the diagonal element of A that was moved into the K-th position. IPVT is only referenced when PIVOT is equal to .TRUE..

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 LCHRG (A, FACT [, …])

Specific: The specific interface names are S_LCHRG and D_LCHRG.

FORTRAN 77 Interface

Single: CALL LCHRG (N, A, LDA, PIVOT, IPVT, FACT, LDFACT)

Double: The double precision name is DLCHRG.

Description

Routine LCHRG is based on the LINPACK routine SCHDC; see Dongarra et al. (1979).

Before the decomposition is computed, initial elements are moved to the leading part of A and final elements to the trailing part of A. During the decomposition only rows and columns corresponding to the free elements are moved. The result of the decomposition is an upper triangular matrix R and a permutation matrix P that satisfy PT AP = RT R, where P is represented by IPVT.

Comments

1. Informational error

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

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

2. Before the decomposition is computed, initial elements are moved to the leading part of A and final elements to the trailing part of A. During the decomposition only rows and columns corresponding to the free elements are moved. The result of the decomposition is an upper triangular matrix R and a permutation matrix P that satisfy PT AP = RT R, where P is represented by IPVT.

3. LCHRG can be used together with subroutines PERMU and LSLDS to solve the positive definite linear system AX = B with the solution X overwriting the right-hand side B as follows:

CALL ISET (N, 0, IPVT, 1)

CALL LCHRG (A, FACT, N, LDA,.TRUE, IPVT, LDFACT)

CALL PERMU (B, IPVT, B, N, 1)

CALL LSLDS (FACT, B, B, N, LDFACT)

CALL PERMU (B, IPVT, B, N, 2)

CALL LCHRG (A, FACT, N, LDA,.TRUE, IPVT, LDFACT)

CALL PERMU (B, IPVT, B, N, 1)

CALL LSLDS (FACT, B, B, N, LDFACT)

CALL PERMU (B, IPVT, B, N, 2)

Example

Routine LCHRG can be used together with the IMSL routines PERMU (see Chapter 11) and LFSDS to solve a positive definite linear system Ax = b. Since A = PRT RP, the system Ax = b is equivalent to RT R(Px) = Pb. LFSDS is used to solve RT Ry = Pb for y. The routine PERMU is used to compute both Pb and x = Py.

USE IMSL_LIBRARIES

! Declare variables

PARAMETER (N=3, LDA=N, LDFACT=N)

INTEGER IPVT(N)

REAL A(LDA,N), FACT(LDFACT,N), B(N), X(N)

LOGICAL PIVOT

!

! Set values for A and B

!

! A = ( 1 -3 2 )

! ( -3 10 -5 )

! ( 2 -5 6 )

!

! B = ( 27 -78 64 )

!

DATA A/1.,-3.,2.,-3.,10.,-5.,2.,-5.,6./

DATA B/27.,-78.,64./

! Pivot using all columns

PIVOT = .TRUE.

IPVT = 0

! Compute Cholesky factorization

CALL LCHRG (A, FACT, PIVOT=PIVOT, IPVT=IPVT)

! Permute B and store in X

CALL PERMU (B, IPVT, X, IPATH=1)

! Solve for X

CALL LFSDS (FACT, X, X)

! Inverse permutation

CALL PERMU (X, IPVT, X, IPATH=2)

! Print X

CALL WRRRN (’X’, X, 1, N, 1)

!

END

Output

X

1 2 3

1.000 -4.000 7.000