LSLCR
Computes the L DU factorization of a real tridiagonal matrix A using a cyclic reduction algorithm.
Required Arguments
C — Array of size 2N containing the upper codiagonal of the N by N tridiagonal matrix in the entries C(1), …, C(N − 1). (Input/Output)
A — Array of size 2N containing the diagonal of the N by N tridiagonal matrix in the entries A(1), …, A(N). (Input/Output)
B — Array of size 2N containing the lower codiagonal of the N by N tridiagonal matrix in the entries B(1), …, B(N − 1). (Input/Output)
Y — Array of size 2N containing the right hand side for the system Ax = y in the order Y(1), …, Y(N). (Input/Output)
The vector x overwrites Y in storage.
U — Array of size 2N of flags that indicate any singularities of A. (Output)
A value U(I) = 1. means that a divide by zero would have occurred during the factoring. Otherwise U(I) = 0.
IR — Array of integers that determine the sizes of loops performed in the cyclic reduction algorithm. (Output)
IS — Array of integers that determine the sizes of loops performed in the cyclic reduction algorithm. (Output)
The sizes of IR and IS must be at least log2 (N) + 3.
Optional Arguments
N — Order of the matrix. (Input)
N must be greater than zero
Default: N = size (C,1).
IJOB — Flag to direct the desired factoring or solving step. (Input)
Default: IJOB = 1.
|
IJOB |
Action |
|
1 |
Factor the matrix A and solve the system Ax = y, where y is stored in array Y. |
|
2 |
Do the solve step only. Use y from array Y. (The factoring step has already been done.) |
|
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 the validity of any input value. |
FORTRAN 90 Interface
Generic: CALL LSLCR (C, A, B, Y, U, IR, IS [, …])
Specific: The specific interface names are S_LSLCR and D_LSLCR.
FORTRAN 77 Interface
Single: CALL LSLCR (N, C, A, B, IJOB, Y, U, IR, IS)
Double: The double precision name is DLSLCR.
Description
Routine LSLCR factors and solves the real tridiagonal linear system Ax = y. The matrix is decomposed in the form A = L DU, where L is unit lower triangular, U is unit upper triangular, and D is diagonal. The algorithm used for the factorization is effectively that described in Kershaw (1982). More details, tests and experiments are reported in Hanson (1990).
LSLCR is intended just for tridiagonal systems. The coefficient matrix does not have to be symmetric. The algorithm amounts to Gaussian elimination, with no pivoting for numerical stability, on the matrix whose rows and columns are permuted to a new order. See Hanson (1990) for details. The expectation is that LSLCR will outperform either LSLTR or LSLPB on vector or parallel computers. Its performance may be inferior for small values of n, on scalar computers, or high-performance computers with non-optimizing compilers.
Example
A system of n = 1000 linear equations is solved. The coefficient matrix is the symmetric matrix of the second difference operation, and the right-hand-side vector y is the first column of the identity matrix. Note that an,n= 1. The solution vector will be the first column of the inverse matrix of A. Then a new system is solved where y is now the last column of the identity matrix. The solution vector for this system will be the last column of the inverse matrix.
USE LSLCR_INT
USE UMACH_INT
! Declare variables
INTEGER LP, N, N2
PARAMETER (LP=12, N=1000, N2=2*N)
!
INTEGER I, IJOB, IR(LP), IS(LP), NOUT
REAL A(N2), B(N2), C(N2), U(N2), Y1(N2), Y2(N2)
!
! Define matrix entries:
DO 10 I=1, N - 1
C(I) = -1.E0
A(I) = 2.E0
B(I) = -1.E0
Y1(I+1) = 0.E0
Y2(I) = 0.E0
10 CONTINUE
A(N) = 1.E0
Y1(1) = 1.E0
Y2(N) = 1.E0
!
! Obtain decomposition of matrix and
! solve the first system:
IJOB = 1
CALL LSLCR (C, A, B, Y1, U, IR, IS, IJOB=IJOB)
!
! Solve the second system with the
! decomposition ready:
IJOB = 2
CALL LSLCR (C, A, B, Y2, U, IR, IS, IJOB=IJOB)
CALL UMACH (2, NOUT)
WRITE (NOUT,*) ’ The value of n is: ’, N
WRITE (NOUT,*) ’ Elements 1, n of inverse matrix columns 1 ’//&
’and n:’, Y1(1), Y2(N)
END
The value of n is: 1000
Elements 1, n of inverse matrix columns 1 and n: 1.00000 1000.000
