Computes the LDU factorization of a complex tridiagonal matrix A using a cyclic reduction algorithm.

Required Arguments

C — Complex 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 — Complex array of size 2N containing the diagonal of the N by N tridiagonal matrix in the entries A(1), …, A(N).   (Input/Output)

B — Complex 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 — Complex array of size 2N containing the right-hand side of the system Ax = y in the order Y(1),…,Y(N).   (Input/Output)
The vector x overwrites Y in storage.

U — Real 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 these arrays 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              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 LSLCQ (C, A, B, Y, U, IR, IS [,…])

Specific:                             The specific interface names are S_LSLCQ and D_LSLCQ.

FORTRAN 77 Interface

Single:            CALL LSLCQ (N, C, A, B, IJOB, Y, U, IR, IS)

Double:                              The double precision name is DLSLCQ.

Description

Routine LSLCQ factors and solves the complex tridiagonal linear system Ax = y. The matrix is decomposed in the form A = LDU, 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).

LSLCQ is intended just for tridiagonal systems. The coefficient matrix does not have to be Hermitian. 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 LSLCQ will outperform either LSLTQ or LSLQB 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 real skew-symmetric tridiagonal matrix, A, of dimension n = 1000 is given by ck = −k, ak = 0, and bk = k, k = 1,…, n − 1, an = 0. This matrix will have eigenvalues that are purely imaginary. The eigenvalue closest to the imaginary unit is required. This number is obtained by using inverse iteration to approximate a complex eigenvector y. The eigenvalue is approximated  by 
λ = yH Ay/yH y. (This example is contrived in the sense that the given tridiagonal skew-symmetric matrix eigenvalue problem is essentially equivalent to the tridiagonal symmetic eigenvalue problem where the ck = k and the other data are unchanged.)

 

      USE LSLCQ_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), K, NOUT

      REAL       AIMAG, U(N2)

      COMPLEX    A(N2), B(N2), C(N2), CMPLX, CONJG, S, T, Y(N2)

      INTRINSIC  AIMAG, CMPLX, CONJG

!                                 Define entries of skew-symmetric

!                                 matrix, A:

      DO 10  I=1, N - 1

         C(I) = -I

!                                 This amounts to subtracting the

!                                 positive imaginary unit from the

!                                 diagonal.  (The eigenvalue closest

!                                 to this value is desired.)

         A(I) = CMPLX(0.E0,-1.0E0)

         B(I) = I

!                                 This initializes the approximate

!                                 eigenvector.

         Y(I) = 1.E0

   10 CONTINUE

      A(N) = CMPLX(0.E0,-1.0E0)

      Y(N) = 1.E0

!                                 First step of inverse iteration

!                                 follows. Obtain decomposition of

!                                 matrix and solve the first system:

      IJOB = 1

      CALL LSLCQ (C, A, B, Y, U, IR, IS, N=N, IJOB=IJOB)

!

!                                 Next steps of inverse iteration

!                                 follow. Solve the system again with

!                                 the decomposition ready:

      IJOB = 2

      DO 20  K=1, 3

         CALL LSLCQ (C, A, B, Y, U, IR, IS, N=N, IJOB=IJOB)

   20 CONTINUE

!

!                                 Compute the Raleigh quotient to

!                                 estimate the eigenvalue closest to

!                                 the positive imaginary unit. After

!                                 the approximate eigenvector, y, is

!                                 computed, the estimate of the

!                                 eigenvalue is ctrans(y)*A*y/t,

!                                 where t = ctrans(y)*y.

      S = -CONJG(Y(1))*Y(2)

      T = CONJG(Y(1))*Y(1)

      DO 30  I=2, N - 1

         S = S + CONJG(Y(I))*((I-1)*Y(I-1)-I*Y(I+1))

         T = T + CONJG(Y(I))*Y(I)

   30 CONTINUE

      S = S + CONJG(Y(N))*(N-1)*Y(N-1)

      T = T + CONJG(Y(N))*Y(N)

      S = S/T

      CALL UMACH (2, NOUT)

      WRITE (NOUT,*) ' The value of n is:  ', N

      WRITE (NOUT,*) ' Value of approximate imaginary eigenvalue:',&

                   AIMAG(S)

      STOP

      END

Output

 

The value of n is:    1000
Value of approximate imaginary eigenvalue:    1.03811

Visual Numerics, Inc. Visual Numerics - Developers of IMSL and PV-WAVE http://www.vni.com/ PHONE: 713.784.3131 FAX:713.781.9260