SUBROUTINE CBLKTR (IFLG,NP,N,AN,BN,CN,MP,M,AM,BM,CM,IDIMY,Y, 1 IERROR,W) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * F I S H P A C K * C * * C * * C * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * C * * C * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * C * * C * (VERSION 4.0 , JUNE 1999) * C * * C * BY * C * * C * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * C * * C * OF * C * * C * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * C * * C * BOULDER, COLORADO (80307) U.S.A. * C * * C * WHICH IS SPONSORED BY * C * * C * THE NATIONAL SCIENCE FOUNDATION * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C 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 REVISION NOVEMBER 1988 C C PURPOSE CBLKTR SOLVES A SYSTEM OF LINEAR EQUATIONS C OF THE FORM C C AN(J)*X(I,J-1) + AM(I)*X(I-1,J) + C (BN(J)+BM(I))*X(I,J) + CN(J)*X(I,J+1) + C 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 C 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 C DISCRETIZATION OF SEPARABLE ELLIPTIC C EQUATIONS. BOUNDARY CONDITIONS MAY BE C DIRICHLET, NEUMANN, OR PERIODIC. C C CBLKTRI IS A COMPLEX VERSION OF PACKAGE C BLKTRI ON ULIB. C C USAGE CALL CBLKTR (IFLG,NP,N,AN,BN,CN,MP,M,AM,BM, C CM,IDIMY,Y,IERROR,W) C C ARGUMENTS C C ON INPUT IFLG C C = 0 INITIALIZATION ONLY. C CERTAIN QUANTITIES THAT DEPEND ON NP, C N, AN, BN, AND CN ARE COMPUTED AND C STORED IN THE WORK ARRAY W. C C = 1 THE QUANTITIES THAT WERE COMPUTED C IN THE INITIALIZATION ARE USED C TO OBTAIN THE SOLUTION X(I,J). C C NOTE: C A CALL WITH IFLG=0 TAKES C APPROXIMATELY ONE HALF THE TIME C AS A CALL WITH IFLG = 1. C HOWEVER, THE INITIALIZATION DOES C NOT HAVE TO BE REPEATED UNLESS NP, C N, AN, BN, OR CN CHANGE. C C NP C = 0 IF AN(1) AND CN(N) ARE NOT ZERO, C WHICH CORRESPONDS TO PERIODIC C BOUNARY CONDITIONS. C C = 1 IF AN(1) AND CN(N) ARE ZERO. C C N C THE NUMBER OF UNKNOWNS IN THE J-DIRECTION. C N MUST BE GREATER THAN 4. C THE OPERATION COUNT IS PROPORTIONAL TO C MNLOG2(N), HENCE N SHOULD BE SELECTED C LESS THAN OR EQUAL TO M. C C AN,BN,CN C ONE-DIMENSIONAL ARRAYS OF LENGTH N C THAT SPECIFY THE COEFFICIENTS IN THE C LINEAR EQUATIONS GIVEN ABOVE. C C MP C = 0 IF AM(1) AND CM(M) ARE NOT ZERO, C WHICH CORRESPONDS TO PERIODIC C BOUNDARY CONDITIONS. C C = 1 IF AM(1) = CM(M) = 0 . C C M C THE NUMBER OF UNKNOWNS IN THE I-DIRECTION. C M MUST BE GREATER THAN 4. C C AM,BM,CM C COMPLEX ONE-DIMENSIONAL ARRAYS OF LENGTH M C THAT SPECIFY THE COEFFICIENTS IN THE LINEAR C EQUATIONS GIVEN ABOVE. C C IDIMY C THE ROW (OR FIRST) DIMENSION OF THE C TWO-DIMENSIONAL ARRAY Y AS IT APPEARS C IN THE PROGRAM CALLING CBLKTR. C THIS PARAMETER IS USED TO SPECIFY THE C VARIABLE DIMENSION OF Y. C IDIMY MUST BE AT LEAST M. C C Y C A COMPLEX TWO-DIMENSIONAL ARRAY THAT C SPECIFIES THE VALUES OF THE RIGHT SIDE OF C THE LINEAR SYSTEM OF EQUATIONS GIVEN ABOVE. C Y MUST BE DIMENSIONED Y(IDIMY,N) WITH C IDIMY .GE. M. C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE C PROVIDED BY THE USER FOR WORK SPACE. C C IF NP=1 DEFINE K=INT(LOG2(N))+1 AND C SET L=2**(K+1) THEN W MUST HAVE DIMENSION C (K-2)*L+K+5+MAX(2N,12M) C C IF NP=0 DEFINE K=INT(LOG2(N-1))+1 AND C SET L=2**(K+1) THEN W MUST HAVE DIMENSION C (K-2)*L+K+5+2N+MAX(2N,12M) C C **IMPORTANT** C FOR PURPOSES OF CHECKING, THE REQUIRED C DIMENSION OF W IS COMPUTED BY CBLKTR AND C STORED IN W(1) IN FLOATING POINT FORMAT. C C ARGUMENTS C C ON OUTPUT Y C CONTAINS THE SOLUTION X. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID C INPUT PARAMETERS. EXCEPT FOR NUMBER ZER0, C 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 CBLKTR FAILED WHILE COMPUTING RESULTS C THAT DEPEND ON THE COEFFICIENT ARRAYS C AN, BN, CN. CHECK THESE ARRAYS. C = 5 AN(J)*CN(J-1) IS LESS THAN 0 FOR SOME J. C C POSSIBLE REASONS FOR THIS CONDITION ARE C 1. THE ARRAYS AN AND CN ARE NOT CORRECT C 2. TOO LARGE A GRID SPACING WAS USED C IN THE DISCRETIZATION OF THE ELLIPTIC C EQUATION. C 3. THE LINEAR EQUATIONS RESULTED FROM A C PARTIAL DIFFERENTIAL EQUATION WHICH C WAS NOT ELLIPTIC. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST C NOT BE DESTROYED IF CBLKTR WILL BE CALLED C AGAIN WITH IFLG=1. W(1) CONTAINS THE C NUMBER OF LOCATIONS REQUIRED BY W IN C FLOATING POINT FORMAT. C C C SPECIAL CONDITIONS THE ALGORITHM MAY FAIL IF ABS(BM(I)+BN(J)) C IS LESS THAN ABS(AM(I))+ABS(AN(J))+ C ABS(CM(I))+ABS(CN(J)) C FOR SOME I AND J. THE ALGORITHM WILL ALSO C FAIL IF AN(J)*CN(J-1) IS LESS THAN ZERO FOR C SOME J. C SEE THE DESCRIPTION OF THE OUTPUT PARAMETER C IERROR. C C C I/O NONE C C PRECISION SINGLE C C REQUIRED LIBRARY COMF FROM FISHPAK C FILES C C LANGUAGE FORTRAN C C HISTORY WRITTEN BY PAUL SWARZTRAUBER AT NCAR IN C THE EARLY 1970'S. REWRITTEN AN RELEASED C ON NCAR'S PUBLIC SOFTWARE LIBRARIES IN C JANUARY, 1980. C C ALGORITHM GENERALIZED CYCLIC REDUCTION C (SEE REFERENCE BELOW) C C PORTABILITY FORTRAN 77 C THE APPROXIMATE MACHINE ACCURACY IS COMPUTED C IN FUNCTION EPMACH C C REFERENCES SWARZTRAUBER,P. AND R. SWEET, 'EFFICIENT C FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC EQUATIONS' C NCAR TN/IA-109, JULY, 1975, 138 PP. C C SWARZTRAUBER P. N.,A DIRECT METHOD FOR C THE DISCRETE SOLUTION OF SEPARABLE C ELLIPTIC EQUATIONS, S.I.A.M. C J. NUMER. ANAL.,11(1974) PP. 1136-1150. C C***********************************************************************