LSLPB
Computes the RTDR Cholesky factorization of a real symmetric positive definite matrix A in codiagonal band symmetric storage mode. Solve a system Ax = b.
Required Arguments
A — Array containing the N by N positive definite band coefficient matrix and right hand side in codiagonal band symmetric storage mode. (Input/Output)
The number of array columns must be at least NCODA + 2. The number of column is not an input to this subprogram.
On output, A contains the solution and factors. See Comments section for details.
NCODA — Number of upper codiagonals of matrix A. (Input)
Must satisfy NCODA ≥ 0 and NCODA < N.
U — Array of flags that indicate any singularities of A, namely loss of positive-definiteness of a leading minor. (Output)
A value U(I) = 0. means that the leading minor of dimension I is not positive-definite. Otherwise, U(I) = 1.
Optional Arguments
N — Order of the matrix. (Input)
Must satisfy N > 0.
Default: N = size (A,2).
LDA — Leading dimension of A exactly as specified in the dimension statement of the calling program. (Input)
Must satisfy LDA ≥ N + NCODA.
Default: LDA = size (A,1).
IJOB — Flag to direct the desired factorization or solving step. (Input)
Default: IJOB = 1.
IJOB Meaning
1 factor the matrix A and solve the system Ax = b, where b is stored in column NCODA + 2 of array A. The vector x overwrites b in storage.
2 solve step only. Use b as column NCODA + 2 of A. (The factorization step has already been done.) The vector x overwrites b in storage.
3 factor the matrix A but do not solve a system.
4,5,6 same meaning as with the value IJOB - 3. For efficiency, no error checking is done on values LDA, N, NCODA, and U(*).
FORTRAN 90 Interface
Generic: CALL LSLPB (A, NCODA, U [, …])
Specific: The specific interface names are S_LSLPB and D_LSLPB.
FORTRAN 77 Interface
Single: CALL LSLPB (N, A, LDA, NCODA, IJOB, U)
Double: The double precision name is DLSLPB.
Description
Routine LSLPB factors and solves the symmetric positive definite banded linear system Ax = b. The matrix is factored so that A = RTDR, where R is unit upper triangular and D is diagonal. The reciprocals of the diagonal entries of D are computed and saved to make the solving step more efficient. Errors will occur if D has a non-positive diagonal element. Such events occur only if A is very close to a singular matrix or is not positive definite.
LSLPB is efficient for problems with a small band width. The particular cases
NCODA = 0, 1, 2 are done with special loops within the code. These cases will give good performance. See Hanson (1989) for details. When solving tridiagonal systems,
NCODA = 1, the cyclic reduction code
LSLCR should be considered as an alternative. The expectation is that
LSLCR will outperform
LSLPB on vector or parallel computers. It may be inferior on scalar computers or even parallel computers with non-optimizing compilers.
Comments
1. Workspace may be explicitly provided, if desired, by use of L2LPB/DL2LPB. The reference is:
CALL L2LPB (N, A, LDA, NCODA, IJOB, U, WK)
The additional argument is:
WK — Work vector of length NCODA.
2. If IJOB=1, 3, 4, or 6, A contains the factors R and D on output. These are stored in codiagonal band symmetric storage mode. Column 1 of A contains the reciprocal of diagonal matrix D. Columns 2 through NCODA+1 contain the upper diagonal values for upper unit diagonal matrix R. If IJOB=1,2, 4, or 5, the last column of A contains the solution on output, replacing b.
3. Informational error
Type | Code | Description |
|---|
4 | 2 | The input matrix is not positive definite. |
Example
A system of four linear equations is solved. The coefficient matrix has real positive definite codiagonal band form and the right-hand-side vector b has four elements.
USE LSLPB_INT
USE WRRRN_INT
! Declare variables
INTEGER LDA, N, NCODA
PARAMETER (N=4, NCODA=2, LDA=N+NCODA)
!
INTEGER IJOB
REAL A(LDA,NCODA+2), U(N)
REAL R(N,N), RT(N,N), D(N,N), WK(N,N), AA(N,N)!!
!
! Set values for A and right side in
! codiagonal band symmetric form:
!
! A = ( * * * * )
! ( * * * * )
! (2.0 * * 6.0)
! (4.0 0.0 * -11.0)
! (7.0 2.0 -1.0 -11.0)
! (3.0 -1.0 1.0 19.0)
!
DATA ((A(I+NCODA,J),I=1,N),J=1,NCODA+2)/2.0, 4.0, 7.0, 3.0, 0.0,&
0.0, 2.0, -1.0, 0.0, 0.0, -1.0, 1.0, 6.0, -11.0, -11.0,&
19.0/
DATA R/16*0.0/, D/16*0.0/, RT/16*0.0/
! Factor and solve A*x = b.
CALL LSLPB(A, NCODA, U)
! Print results
CALL WRRRN ('X', A((NCODA+1):,(NCODA+2):), NRA=1, NCA=N, LDA=1)
END
Output
X
1 2 3 4
4.000 -6.000 2.000 9.000
