Computes the RH DR Cholesky factorization of a complex Hermitian positive-definite matrix A in codiagonal band Hermitian storage mode. Solve a system Ax = b.
A Array
containing the N
by N
positive-definite band coefficient matrix and the right hand side in codiagonal
band Hermitian storage mode. (Input/Output)
The number of array
columns must be at least 2 * NCODA + 3. The number
of columns is not an input to this subprogram.
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.
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 the real part of b is stored in column 2 * NCODA + 2 and the imaginary part of b is stored in column 2 * NCODA + 3 of array A. The real and imaginary parts of b are overwritten by the real and imaginary parts of x.
2 solve step only. Use the real part of b as column 2 * NCODA + 2 and the imaginary part of b as column 2 * NCODA + 3 of A. (The factorization step has already been done.) The real and imaginary parts of b are overwritten by the real and imaginary parts of x.
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(*).
Generic: CALL LSLQB (A, NCODA, U [, ])
Specific: The specific interface names are S_LSLQB and D_LSLQB.
Single: CALL LSLQB (N, A, LDA, NCODA, IJOB, U)
Double: The double precision name is DLSLQB.
Routine LSLQB factors and solves the Hermitian positive definite banded linear system Ax = b. The matrix is factored so that A = RH DR, where R is unit upper triangular and D is diagonal and real. 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 nonpositive diagonal element. Such events occur only if A is very close to a singular matrix or is not positive definite.
LSLQB is efficient for problems with a small band width. The particular cases NCODA = 0, 1 are done with special loops within the code. These cases will give good performance. See Hanson (1989) for more on the algorithm. When solving tridiagonal systems, NCODA = 1, the cyclic reduction code LSLCQ should be considered as an alternative. The expectation is that LSLCQ will outperform LSLQB on vector or parallel computers. It may be inferior on scalar computers or even parallel computers with non-optimizing compilers.
1. Workspace may be explicitly provided, if desired, by use of L2LQB/DL2LQB The reference is:
CALL L2LQB (N, A, LDA, NCODA, IJOB, U, WK1, WK2)
The additional arguments are as follows:
WK1 Work vector of length NCODA.
WK2 Work vector of length NCODA.
2. Informational error
Type Code
4 2 The input matrix is not positive definite.
A system of five linear equations is solved. The coefficient matrix has real positive definite codiagonal Hermitian band form and the right-hand-side vector b has five elements.
USE
LSLQB_INT
USE WRRRN_INT
INTEGER LDA, N, NCODA
PARAMETER (N=5, NCODA=1, LDA=N+NCODA)
!
INTEGER I, IJOB, J
REAL A(LDA,2*NCODA+3), U(N)
!
! Set values for A and right hand side
! in codiagonal band Hermitian form:
!
! ( * * * * * )
! ( 2.0 * * 1.0 5.0)
! A = ( 4.0 -1.0 1.0 12.0 -6.0)
! (10.0 1.0 2.0 1.0 -16.0)
! ( 6.0 0.0 4.0 -3.0 -3.0)
! ( 9.0 1.0 1.0 25.0 16.0)
!
DATA ((A(I+NCODA,J),I=1,N),J=1,2*NCODA+3)/2.0, 4.0, 10.0, 6.0,&
9.0, 0.0, -1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 2.0, 4.0, 1.0,&
1.0, 12.0, 1.0, -3.0, 25.0, 5.0, -6.0, -16.0, -3.0, 16.0/
!
! Factor and solve A*x = b.
!
IJOB = 1
CALL LSLQB (A, NCODA, U)
!
! Print results
!
CALL WRRRN (REAL(X), A((NCODA+1):,(2*NCODA+2):), 1, N, 1)
CALL WRRRN (IMAG(X), A((NCODA+1):,(2*NCODA+3):), 1, N, 1)
END
REAL(X)
1 2
3 4
5
2.000 3.000 -1.000 0.000
3.000
IMAG(X)
1 2 3 4 5
1.000 0.000 -1.000 -2.000 2.000
PHONE: 713.784.3131 FAX:713.781.9260 |