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