*DECK CMGNBN
SUBROUTINE CMGNBN (NPEROD, N, MPEROD, M, A, B, C, IDIMY, Y,
+ IERROR, W)
C***BEGIN PROLOGUE CMGNBN
C***PURPOSE Solve a complex block tridiagonal linear system of
C equations by a cyclic reduction algorithm.
C***LIBRARY SLATEC (FISHPACK)
C***CATEGORY I2B4B
C***TYPE COMPLEX (GENBUN-S, CMGNBN-C)
C***KEYWORDS CYCLIC REDUCTION, ELLIPTIC PDE, FISHPACK,
C TRIDIAGONAL LINEAR SYSTEM
C***AUTHOR Adams, J., (NCAR)
C Swarztrauber, P. N., (NCAR)
C Sweet, R., (NCAR)
C***DESCRIPTION
C
C Subroutine CMGNBN solves the complex linear system of equations
C
C A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J)
C
C + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J)
C
C For I = 1,2,...,M and J = 1,2,...,N.
C
C The indices I+1 and I-1 are evaluated modulo M, i.e.,
C X(0,J) = X(M,J) and X(M+1,J) = X(1,J), and X(I,0) may be equal to
C 0, X(I,2), or X(I,N) and X(I,N+1) may be equal to 0, X(I,N-1), or
C X(I,1) depending on an input parameter.
C
C
C * * * * * * * * Parameter Description * * * * * * * * * *
C
C * * * * * * On Input * * * * * *
C
C NPEROD
C Indicates the values that X(I,0) and X(I,N+1) are assumed to
C have.
C
C = 0 If X(I,0) = X(I,N) and X(I,N+1) = X(I,1).
C = 1 If X(I,0) = X(I,N+1) = 0 .
C = 2 If X(I,0) = 0 and X(I,N+1) = X(I,N-1).
C = 3 If X(I,0) = X(I,2) and X(I,N+1) = X(I,N-1).
C = 4 If X(I,0) = X(I,2) and X(I,N+1) = 0.
C
C N
C The number of unknowns in the J-direction. N must be greater
C than 2.
C
C MPEROD
C = 0 If A(1) and C(M) are not zero
C = 1 If A(1) = C(M) = 0
C
C M
C The number of unknowns in the I-direction. N must be greater
C than 2.
C
C A,B,C
C One-dimensional complex arrays of length M that specify the
C coefficients in the linear equations given above. If MPEROD = 0
C the array elements must not depend upon the index I, but must be
C constant. Specifically, the subroutine checks the following
C condition
C
C A(I) = C(1)
C C(I) = C(1)
C B(I) = B(1)
C
C For I=1,2,...,M.
C
C IDIMY
C The row (or first) dimension of the two-dimensional array Y as
C it appears in the program calling CMGNBN. This parameter is
C used to specify the variable dimension of Y. IDIMY must be at
C least M.
C
C Y
C A two-dimensional complex array that specifies the values of the
C right side of the linear system of equations given above. Y
C must be dimensioned at least M*N.
C
C W
C A one-dimensional complex array that must be provided by the
C user for work space. W may require up to 4*N +
C (10 + INT(log2(N)))*M LOCATIONS. The actual number of locations
C used is computed by CMGNBN and is returned in location W(1).
C
C
C * * * * * * On Output * * * * * *
C
C Y
C Contains the solution X.
C
C IERROR
C An error flag which indicates invalid input parameters. Except
C for number zero, a solution is not attempted.
C
C = 0 No error.
C = 1 M .LE. 2
C = 2 N .LE. 2
C = 3 IDIMY .LT. M
C = 4 NPEROD .LT. 0 or NPEROD .GT. 4
C = 5 MPEROD .LT. 0 or MPEROD .GT. 1
C = 6 A(I) .NE. C(1) or C(I) .NE. C(1) or B(I) .NE. B(1) for
C some I=1,2,...,M.
C = 7 A(1) .NE. 0 or C(M) .NE. 0 and MPEROD = 1
C
C W
C W(1) contains the required length of W.
C
C *Long Description:
C
C * * * * * * * Program Specifications * * * * * * * * * * * *
C
C Dimension of A(M),B(M),C(M),Y(IDIMY,N),W(see parameter list)
C Arguments
C
C Latest June 1979
C Revision
C
C Subprograms CMGNBN,CMPOSD,CMPOSN,CMPOSP,CMPCSG,CMPMRG,
C Required CMPTRX,CMPTR3,PIMACH
C
C Special None
C Conditions
C
C Common None
C Blocks
C
C I/O None
C
C Precision Single
C
C Specialist Roland Sweet
C
C Language FORTRAN
C
C History Written by Roland Sweet at NCAR in June, 1977
C
C Algorithm The linear system is solved by a cyclic reduction
C algorithm described in the reference.
C
C Space 4944(DECIMAL) = 11520(octal) locations on the NCAR
C Required Control Data 7600
C
C Timing and The execution time T on the NCAR Control Data
C Accuracy 7600 for subroutine CMGNBN is roughly proportional
C to M*N*log2(N), but also depends on the input
C parameter NPEROD. Some typical values are listed
C in the table below.
C To measure the accuracy of the algorithm a
C uniform random number generator was used to create
C a solution array X for the system given in the
C 'PURPOSE' with
C
C A(I) = C(I) = -0.5*B(I) = 1, I=1,2,...,M
C
C and, when MPEROD = 1
C
C A(1) = C(M) = 0
C A(M) = C(1) = 2.
C
C The solution X was substituted into the given sys-
C tem and a right side Y was computed. Using this
C array Y subroutine CMGNBN was called to produce an
C approximate solution Z. Then the relative error,
C defined as
C
C E = MAX(ABS(Z(I,J)-X(I,J)))/MAX(ABS(X(I,J)))
C
C where the two maxima are taken over all I=1,2,...,M
C and J=1,2,...,N, was computed. The value of E is
C given in the table below for some typical values of
C M and N.
C
C
C M (=N) MPEROD NPEROD T(MSECS) E
C ------ ------ ------ -------- ------
C
C 31 0 0 77 1.E-12
C 31 1 1 45 4.E-13
C 31 1 3 91 2.E-12
C 32 0 0 59 7.E-14
C 32 1 1 65 5.E-13
C 32 1 3 97 2.E-13
C 33 0 0 80 6.E-13
C 33 1 1 67 5.E-13
C 33 1 3 76 3.E-12
C 63 0 0 350 5.E-12
C 63 1 1 215 6.E-13
C 63 1 3 412 1.E-11
C 64 0 0 264 1.E-13
C 64 1 1 287 3.E-12
C 64 1 3 421 3.E-13
C 65 0 0 338 2.E-12
C 65 1 1 292 5.E-13
C 65 1 3 329 1.E-11
C
C Portability American National Standards Institute Fortran.
C The machine dependent constant PI is defined in
C function PIMACH.
C
C Required COS
C Resident
C Routines
C
C Reference Sweet, R., 'A Cyclic Reduction Algorithm for
C Solving Block Tridiagonal Systems Of Arbitrary
C Dimensions,' SIAM J. on Numer. Anal.,
C 14(SEPT., 1977), PP. 706-720.
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C***REFERENCES R. Sweet, A cyclic reduction algorithm for solving
C block tridiagonal systems of arbitrary dimensions,
C SIAM Journal on Numerical Analysis 14, (September
C 1977), pp. 706-720.
C***ROUTINES CALLED CMPOSD, CMPOSN, CMPOSP
C***REVISION HISTORY (YYMMDD)
C 801001 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890531 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE CMGNBN
C
C
COMPLEX A ,B ,C ,Y ,
1 W ,A1
DIMENSION Y(IDIMY,*)
DIMENSION W(*) ,B(*) ,A(*) ,C(*)
C***FIRST EXECUTABLE STATEMENT CMGNBN
IERROR = 0
IF (M .LE. 2) IERROR = 1
IF (N .LE. 2) IERROR = 2
IF (IDIMY .LT. M) IERROR = 3
IF (NPEROD.LT.0 .OR. NPEROD.GT.4) IERROR = 4
IF (MPEROD.LT.0 .OR. MPEROD.GT.1) IERROR = 5
IF (MPEROD .EQ. 1) GO TO 102
DO 101 I=2,M
IF (ABS(A(I)-C(1)) .NE. 0.) GO TO 103
IF (ABS(C(I)-C(1)) .NE. 0.) GO TO 103
IF (ABS(B(I)-B(1)) .NE. 0.) GO TO 103
101 CONTINUE
GO TO 104
102 IF (ABS(A(1)).NE.0. .AND. ABS(C(M)).NE.0.) IERROR = 7
GO TO 104
103 IERROR = 6
104 IF (IERROR .NE. 0) RETURN
IWBA = M+1
IWBB = IWBA+M
IWBC = IWBB+M
IWB2 = IWBC+M
IWB3 = IWB2+M
IWW1 = IWB3+M
IWW2 = IWW1+M
IWW3 = IWW2+M
IWD = IWW3+M
IWTCOS = IWD+M
IWP = IWTCOS+4*N
DO 106 I=1,M
K = IWBA+I-1
W(K) = -A(I)
K = IWBC+I-1
W(K) = -C(I)
K = IWBB+I-1
W(K) = 2.-B(I)
DO 105 J=1,N
Y(I,J) = -Y(I,J)
105 CONTINUE
106 CONTINUE
MP = MPEROD+1
NP = NPEROD+1
GO TO (114,107),MP
107 GO TO (108,109,110,111,123),NP
108 CALL CMPOSP (M,N,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
1 W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
2 W(IWP))
GO TO 112
109 CALL CMPOSD (M,N,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWW1),
1 W(IWD),W(IWTCOS),W(IWP))
GO TO 112
110 CALL CMPOSN (M,N,1,2,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
1 W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
2 W(IWP))
GO TO 112
111 CALL CMPOSN (M,N,1,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
1 W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
2 W(IWP))
112 IPSTOR = REAL(W(IWW1))
IREV = 2
IF (NPEROD .EQ. 4) GO TO 124
113 GO TO (127,133),MP
114 CONTINUE
C
C REORDER UNKNOWNS WHEN MP =0
C
MH = (M+1)/2
MHM1 = MH-1
MODD = 1
IF (MH*2 .EQ. M) MODD = 2
DO 119 J=1,N
DO 115 I=1,MHM1
MHPI = MH+I
MHMI = MH-I
W(I) = Y(MHMI,J)-Y(MHPI,J)
W(MHPI) = Y(MHMI,J)+Y(MHPI,J)
115 CONTINUE
W(MH) = 2.*Y(MH,J)
GO TO (117,116),MODD
116 W(M) = 2.*Y(M,J)
117 CONTINUE
DO 118 I=1,M
Y(I,J) = W(I)
118 CONTINUE
119 CONTINUE
K = IWBC+MHM1-1
I = IWBA+MHM1
W(K) = (0.,0.)
W(I) = (0.,0.)
W(K+1) = 2.*W(K+1)
GO TO (120,121),MODD
120 CONTINUE
K = IWBB+MHM1-1
W(K) = W(K)-W(I-1)
W(IWBC-1) = W(IWBC-1)+W(IWBB-1)
GO TO 122
121 W(IWBB-1) = W(K+1)
122 CONTINUE
GO TO 107
C
C REVERSE COLUMNS WHEN NPEROD = 4
C
123 IREV = 1
NBY2 = N/2
124 DO 126 J=1,NBY2
MSKIP = N+1-J
DO 125 I=1,M
A1 = Y(I,J)
Y(I,J) = Y(I,MSKIP)
Y(I,MSKIP) = A1
125 CONTINUE
126 CONTINUE
GO TO (110,113),IREV
127 CONTINUE
DO 132 J=1,N
DO 128 I=1,MHM1
MHMI = MH-I
MHPI = MH+I
W(MHMI) = .5*(Y(MHPI,J)+Y(I,J))
W(MHPI) = .5*(Y(MHPI,J)-Y(I,J))
128 CONTINUE
W(MH) = .5*Y(MH,J)
GO TO (130,129),MODD
129 W(M) = .5*Y(M,J)
130 CONTINUE
DO 131 I=1,M
Y(I,J) = W(I)
131 CONTINUE
132 CONTINUE
133 CONTINUE
C
C RETURN STORAGE REQUIREMENTS FOR W ARRAY.
C
W(1) = CMPLX(REAL(IPSTOR+IWP-1),0.)
RETURN
END