*DECK BLKTRI SUBROUTINE BLKTRI (IFLG, NP, N, AN, BN, CN, MP, M, AM, BM, CM, + IDIMY, Y, IERROR, W) C***BEGIN PROLOGUE BLKTRI C***PURPOSE Solve a block tridiagonal system of linear equations C (usually resulting from the discretization of separable C two-dimensional elliptic equations). C***LIBRARY SLATEC (FISHPACK) C***CATEGORY I2B4B C***TYPE SINGLE PRECISION (BLKTRI-S, CBLKTR-C) C***KEYWORDS ELLIPTIC PDE, FISHPACK, TRIDIAGONAL LINEAR SYSTEM C***AUTHOR Adams, J., (NCAR) C Swarztrauber, P. N., (NCAR) C Sweet, R., (NCAR) C***DESCRIPTION C C Subroutine BLKTRI Solves a System of Linear Equations of the Form C C AN(J)*X(I,J-1) + AM(I)*X(I-1,J) + (BN(J)+BM(I))*X(I,J) C C + CN(J)*X(I,J+1) + CM(I)*X(I+1,J) = Y(I,J) C C for I = 1,2,...,M and J = 1,2,...,N. C C I+1 and I-1 are evaluated modulo M and J+1 and J-1 modulo N, i.e., C C X(I,0) = X(I,N), X(I,N+1) = X(I,1), C X(0,J) = X(M,J), X(M+1,J) = X(1,J). C C These equations usually result from the discretization of C separable elliptic equations. Boundary conditions may be C Dirichlet, Neumann, or Periodic. C C C * * * * * * * * * * ON INPUT * * * * * * * * * * C C IFLG C = 0 Initialization only. Certain quantities that depend on NP, C N, AN, BN, and CN are computed and stored in the work C array W. C = 1 The quantities that were computed in the initialization are C used to obtain the solution X(I,J). C C NOTE A call with IFLG=0 takes approximately one half the time C as a call with IFLG = 1 . However, the C initialization does not have to be repeated unless NP, N, C AN, BN, or CN change. C C NP C = 0 If AN(1) and CN(N) are not zero, which corresponds to C periodic boundary conditions. C = 1 If AN(1) and CN(N) are zero. C C N C The number of unknowns in the J-direction. N must be greater C than 4. The operation count is proportional to MNlog2(N), hence C N should be selected less than or equal to M. C C AN,BN,CN C One-dimensional arrays of length N that specify the coefficients C in the linear equations given above. C C MP C = 0 If AM(1) and CM(M) are not zero, which corresponds to C periodic boundary conditions. C = 1 If AM(1) = CM(M) = 0 . C C M C The number of unknowns in the I-direction. M must be greater C than 4. C C AM,BM,CM C One-dimensional arrays of length M that specify the coefficients C in the linear equations given above. C C IDIMY C The row (or first) dimension of the two-dimensional array Y as C it appears in the program calling BLKTRI. 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 array that specifies the values of the right C side of the linear system of equations given above. Y must be C dimensioned at least M*N. C C W C A one-dimensional array that must be provided by the user for C work space. C If NP=1 define K=INT(log2(N))+1 and set L=2**(K+1) then C W must have dimension (K-2)*L+K+5+MAX(2N,6M) C C If NP=0 define K=INT(log2(N-1))+1 and set L=2**(K+1) then C W must have dimension (K-2)*L+K+5+2N+MAX(2N,6M) C C **IMPORTANT** For purposes of checking, the required dimension C of W is computed by BLKTRI and stored in W(1) C in floating point format. C C * * * * * * * * * * On Output * * * * * * * * * * C C Y C Contains the solution X. C C IERROR C An error flag that indicates invalid input parameters. Except C for number zero, a solution is not attempted. C C = 0 No error. C = 1 M is less than 5. C = 2 N is less than 5. C = 3 IDIMY is less than M. C = 4 BLKTRI failed while computing results that depend on the C coefficient arrays AN, BN, CN. Check these arrays. C = 5 AN(J)*CN(J-1) is less than 0 for some J. Possible reasons C for this condition are C 1. The arrays AN and CN are not correct. C 2. Too large a grid spacing was used in the discretization C of the elliptic equation. C 3. The linear equations resulted from a partial C differential equation which was not elliptic. C C W C Contains intermediate values that must not be destroyed if C BLKTRI will be called again with IFLG=1. W(1) contains the C number of locations required by W in floating point format. C C *Long Description: C C * * * * * * * Program Specifications * * * * * * * * * * * * C C Dimension of AN(N),BN(N),CN(N),AM(M),BM(M),CM(M),Y(IDIMY,N) C Arguments W(See argument list) C C Latest June 1979 C Revision C C Required BLKTRI,BLKTRI,PROD,PRODP,CPROD,CPRODP,COMPB,INDXA, C Subprograms INDXB,INDXC,PPADD,PSGF,PPSGF,PPSPF,BSRH,TEVLS, C R1MACH C C Special The Algorithm may fail if ABS(BM(I)+BN(J)) is less C Conditions than ABS(AM(I))+ABS(AN(J))+ABS(CM(I))+ABS(CN(J)) C for some I and J. The Algorithm will also fail if C AN(J)*CN(J-1) is less than zero for some J. C See the description of the output parameter IERROR. C C Common CBLKT C Blocks C C I/O None C C Precision Single C C Specialist Paul Swarztrauber C C Language FORTRAN C C History Version 1 September 1973 C Version 2 April 1976 C Version 3 June 1979 C C Algorithm Generalized Cyclic Reduction (See Reference below) C C Space C Required Control Data 7600 C C Portability American National Standards Institute Fortran. C The machine accuracy is set using function R1MACH. C C Required None C Resident C Routines C C References Swarztrauber,P. and R. Sweet, 'Efficient FORTRAN C Subprograms For The Solution Of Elliptic Equations' C NCAR TN/IA-109, July, 1975, 138 PP. C C Swarztrauber P. ,'A Direct Method For The Discrete C Solution Of Separable Elliptic Equations', S.I.A.M. C J. Numer. Anal.,11(1974) PP. 1136-1150. C C***REFERENCES P. N. Swarztrauber and R. Sweet, Efficient Fortran C subprograms for the solution of elliptic equations, C NCAR TN/IA-109, July 1975, 138 pp. C P. N. Swarztrauber, A direct method for the discrete C solution of separable elliptic equations, SIAM Journal C on Numerical Analysis 11, (1974), pp. 1136-1150. C***ROUTINES CALLED BLKTR1, COMPB, CPROD, CPRODP, PROD, PRODP C***COMMON BLOCKS CBLKT 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 BLKTRI C DIMENSION AN(*) ,BN(*) ,CN(*) ,AM(*) , 1 BM(*) ,CM(*) ,Y(IDIMY,*) ,W(*) EXTERNAL PROD ,PRODP ,CPROD ,CPRODP COMMON /CBLKT/ NPP ,K ,EPS ,CNV , 1 NM ,NCMPLX ,IK C***FIRST EXECUTABLE STATEMENT BLKTRI NM = N IERROR = 0 IF (M-5) 101,102,102 101 IERROR = 1 GO TO 119 102 IF (NM-3) 103,104,104 103 IERROR = 2 GO TO 119 104 IF (IDIMY-M) 105,106,106 105 IERROR = 3 GO TO 119 106 NH = N NPP = NP IF (NPP) 107,108,107 107 NH = NH+1 108 IK = 2 K = 1 109 IK = IK+IK K = K+1 IF (NH-IK) 110,110,109 110 NL = IK IK = IK+IK NL = NL-1 IWAH = (K-2)*IK+K+6 IF (NPP) 111,112,111 C C DIVIDE W INTO WORKING SUB ARRAYS C 111 IW1 = IWAH IWBH = IW1+NM W(1) = IW1-1+MAX(2*NM,6*M) GO TO 113 112 IWBH = IWAH+NM+NM IW1 = IWBH W(1) = IW1-1+MAX(2*NM,6*M) NM = NM-1 C C SUBROUTINE COMP B COMPUTES THE ROOTS OF THE B POLYNOMIALS C 113 IF (IERROR) 119,114,119 114 IW2 = IW1+M IW3 = IW2+M IWD = IW3+M IWW = IWD+M IWU = IWW+M IF (IFLG) 116,115,116 115 CALL COMPB (NL,IERROR,AN,BN,CN,W(2),W(IWAH),W(IWBH)) GO TO 119 116 IF (MP) 117,118,117 C C SUBROUTINE BLKTR1 SOLVES THE LINEAR SYSTEM C 117 CALL BLKTR1 (NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),W(IW1),W(IW2), 1 W(IW3),W(IWD),W(IWW),W(IWU),PROD,CPROD) GO TO 119 118 CALL BLKTR1 (NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),W(IW1),W(IW2), 1 W(IW3),W(IWD),W(IWW),W(IWU),PRODP,CPRODP) 119 CONTINUE RETURN END