SUBROUTINE POIS3D (LPEROD,L,C1,MPEROD,M,C2,NPEROD,N,A,B,C,LDIMF, 1 MDIMF,F,IERROR,W) C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * F I S H P A K * C * * C * * C * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * C * * C * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * C * * C * (VERSION 3.1 , OCTOBER 1980) * 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 * * * * * * * * * PURPOSE * * * * * * * * * * * * * * * * * * C C SUBROUTINE POIS3D SOLVES THE LINEAR SYSTEM OF EQUATIONS C C C1*(X(I-1,J,K)-2.*X(I,J,K)+X(I+1,J,K)) C + C2*(X(I,J-1,K)-2.*X(I,J,K)+X(I,J+1,K)) C + A(K)*X(I,J,K-1)+B(K)*X(I,J,K)+C(K)*X(I,J,K+1) = F(I,J,K) C C FOR I=1,2,...,L , J=1,2,...,M , AND K=1,2,...,N . C C THE INDICES K-1 AND K+1 ARE EVALUATED MODULO N, I.E. C X(I,J,0) = X(I,J,N) AND X(I,J,N+1) = X(I,J,1). THE UNKNOWNS C X(0,J,K), X(L+1,J,K), X(I,0,K), AND X(I,M+1,K) ARE ASSUMED TO TAKE C ON CERTAIN PRESCRIBED VALUES DESCRIBED BELOW. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C C * * * * * * * * PARAMETER DESCRIPTION * * * * * * * * * * C C C * * * * * * ON INPUT * * * * * * C C LPEROD INDICATES THE VALUES THAT X(0,J,K) AND X(L+1,J,K) ARE C ASSUMED TO HAVE. C C = 0 IF X(0,J,K) = X(L,J,K) AND X(L+1,J,K) = X(1,J,K). C = 1 IF X(0,J,K) = X(L+1,J,K) = 0. C = 2 IF X(0,J,K) = 0 AND X(L+1,J,K) = X(L-1,J,K). C = 3 IF X(0,J,K) = X(2,J,K) AND X(L+1,J,K) = X(L-1,J,K). C = 4 IF X(0,J,K) = X(2,J,K) AND X(L+1,J,K) = 0. C C L THE NUMBER OF UNKNOWNS IN THE I-DIRECTION. L MUST BE AT C LEAST 3. C C C1 THE REAL CONSTANT THAT APPEARS IN THE ABOVE EQUATION. C C MPEROD INDICATES THE VALUES THAT X(I,0,K) AND X(I,M+1,K) ARE C ASSUMED TO HAVE. C C = 0 IF X(I,0,K) = X(I,M,K) AND X(I,M+1,K) = X(I,1,K). C = 1 IF X(I,0,K) = X(I,M+1,K) = 0. C = 2 IF X(I,0,K) = 0 AND X(I,M+1,K) = X(I,M-1,K). C = 3 IF X(I,0,K) = X(I,2,K) AND X(I,M+1,K) = X(I,M-1,K). C = 4 IF X(I,0,K) = X(I,2,K) AND X(I,M+1,K) = 0. C C M THE NUMBER OF UNKNOWNS IN THE J-DIRECTION. M MUST BE AT C LEAST 3. C C C2 THE REAL CONSTANT WHICH APPEARS IN THE ABOVE EQUATION. C C NPEROD = 0 IF A(1) AND C(N) ARE NOT ZERO. C = 1 IF A(1) = C(N) = 0. C C N THE NUMBER OF UNKNOWNS IN THE K-DIRECTION. N MUST BE AT C LEAST 3. C C C A,B,C ONE-DIMENSIONAL ARRAYS OF LENGTH N THAT SPECIFY THE C COEFFICIENTS IN THE LINEAR EQUATIONS GIVEN ABOVE. C C IF NPEROD = 0 THE ARRAY ELEMENTS MUST NOT DEPEND UPON THE C INDEX K, BUT MUST BE CONSTANT. SPECIFICALLY,THE C SUBROUTINE CHECKS THE FOLLOWING CONDITION C C A(K) = C(1) C C(K) = C(1) C B(K) = B(1) C C FOR K=1,2,...,N. C C LDIMF THE ROW (OR FIRST) DIMENSION OF THE THREE-DIMENSIONAL C ARRAY F AS IT APPEARS IN THE PROGRAM CALLING POIS3D. C THIS PARAMETER IS USED TO SPECIFY THE VARIABLE DIMENSION C OF F. LDIMF MUST BE AT LEAST L. C C MDIMF THE COLUMN (OR SECOND) DIMENSION OF THE THREE-DIMENSIONAL C ARRAY F AS IT APPEARS IN THE PROGRAM CALLING POIS3D. C THIS PARAMETER IS USED TO SPECIFY THE VARIABLE DIMENSION C OF F. MDIMF MUST BE AT LEAST M. C C F A THREE-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUES OF C THE RIGHT SIDE OF THE LINEAR SYSTEM OF EQUATIONS GIVEN C ABOVE. F MUST BE DIMENSIONED AT LEAST L X M X N. C C W A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE C USER FOR WORK SPACE. THE LENGTH OF W MUST BE AT LEAST C 30 + L + M + 2*N + MAX(L,M,N) + C 7*(INT((L+1)/2) + INT((M+1)/2)). C C C * * * * * * ON OUTPUT * * * * * * C C F CONTAINS THE SOLUTION X. C C IERROR AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. C EXCEPT FOR NUMBER ZERO, A SOLUTION IS NOT ATTEMPTED. C = 0 NO ERROR C = 1 IF LPEROD .LT. 0 OR .GT. 4 C = 2 IF L .LT. 3 C = 3 IF MPEROD .LT. 0 OR .GT. 4 C = 4 IF M .LT. 3 C = 5 IF NPEROD .LT. 0 OR .GT. 1 C = 6 IF N .LT. 3 C = 7 IF LDIMF .LT. L C = 8 IF MDIMF .LT. M C = 9 IF A(K) .NE. C(1) OR C(K) .NE. C(1) OR B(I) .NE.B(1) C FOR SOME K=1,2,...,N. C = 10 IF NPEROD = 1 AND A(1) .NE. 0 OR C(N) .NE. 0 C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY C INCORRECT CALL TO POIS3D, THE USER SHOULD TEST IERROR C AFTER THE CALL. C C C * * * * * * * PROGRAM SPECIFICATIONS * * * * * * * * * * * * C C DIMENSION OF A(N),B(N),C(N),F(LDIMF,MDIMF,N), C ARGUMENTS W(SEE ARGUMENT LIST) C C LATEST DECEMBER 1, 1978 C REVISION C C SUBPROGRAMS POIS3D,POS3D1,TRID,RFFTI,RFFTF,RFFTF1,RFFTB, C REQUIRED RFFTB1,COSTI,COST,SINTI,SINT,COSQI,COSQF,COSQF1 C COSQB,COSQB1,SINQI,SINQF,SINQB,CFFTI,CFFTI1, C CFFTB,CFFTB1,PASSB2,PASSB3,PASSB4,PASSB,CFFTF, C CFFTF1,PASSF1,PASSF2,PASSF3,PASSF4,PASSF,PIMACH, C C SPECIAL NONE C CONDITIONS C C COMMON VALUE 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 JULY,1977 C C ALGORITHM THIS SUBROUTINE SOLVES THREE-DIMENSIONAL BLOCK C TRIDIAGONAL LINEAR SYSTEMS ARISING FROM FINITE C DIFFERENCE APPROXIMATIONS TO THREE-DIMENSIONAL C POISSON EQUATIONS USING THE FOURIER TRANSFORM C PACKAGE SCLRFFTPAK WRITTEN BY PAUL SWARZTRAUBER. C C SPACE 6561(DECIMAL) = 14641(OCTAL) LOCATIONS ON THE C REQUIRED NCAR CONTROL DATA 7600 C C TIMING AND THE EXECUTION TIME T ON THE NCAR CONTROL DATA C ACCURACY 7600 FOR SUBROUTINE POIS3D IS ROUGHLY PROPORTIONAL C TO L*M*N*(LOG2(L)+LOG2(M)+5), BUT ALSO DEPENDS ON C INPUT PARAMETERS LPEROD AND MPEROD. SOME TYPICAL C VALUES ARE LISTED IN THE TABLE BELOW WHEN NPEROD=0. 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(K) = C(K) = -0.5*B(K) = 1, K=1,2,...,N C C AND, WHEN NPEROD = 1 C C A(1) = C(N) = 0 C A(N) = C(1) = 2. C C THE SOLUTION X WAS SUBSTITUTED INTO THE GIVEN SYS- C TEM AND, USING DOUBLE PRECISION, A RIGHT SIDE Y WAS C COMPUTED. USING THIS ARRAY Y SUBROUTINE POIS WAS C CALLED TO PRODUCE AN APPROXIMATE SOLUTION Z. THEN C THE RELATIVE ERROR, DEFINED AS C C E = MAX(ABS(Z(I,J,K)-X(I,J,K)))/MAX(ABS(X(I,J,K))) C C WHERE THE TWO MAXIMA ARE TAKEN OVER I=1,2,...,L, C J=1,2,...,M AND K=1,2,...,N, WAS COMPUTED. THE C VALUE OF E IS GIVEN IN THE TABLE BELOW FOR SOME C TYPICAL VALUES OF L,M AND N. C C C L(=M=N) LPEROD MPEROD T(MSECS) E C ------ ------ ------ -------- ------ C C 16 0 0 272 1.E-13 C 15 1 1 287 4.E-13 C 17 3 3 338 2.E-13 C 32 0 0 1755 2.E-13 C 31 1 1 1894 2.E-12 C 33 3 3 2042 7.E-13 C C C PORTABILITY AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN. C THE MACHINE DEPENDENT CONSTANT PI IS DEFINED IN C FUNCTION PIMACH. C C REQUIRED COS,SIN,ATAN C RESIDENT C ROUTINES C C REFERENCE NONE C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DIMENSION A(1) ,B(1) ,C(1) , 1 F(LDIMF,MDIMF,1) ,W(1) ,SAVE(6) LP = LPEROD+1 MP = MPEROD+1 NP = NPEROD+1 C C CHECK FOR INVALID INPUT. C IERROR = 0 IF (LP.LT.1 .OR. LP.GT.5) IERROR = 1 IF (L .LT. 3) IERROR = 2 IF (MP.LT.1 .OR. MP.GT.5) IERROR = 3 IF (M .LT. 3) IERROR = 4 IF (NP.LT.1 .OR. NP.GT.2) IERROR = 5 IF (N .LT. 3) IERROR = 6 IF (LDIMF .LT. L) IERROR = 7 IF (MDIMF .LT. M) IERROR = 8 IF (NP .NE. 1) GO TO 103 DO 101 K=1,N IF (A(K) .NE. C(1)) GO TO 102 IF (C(K) .NE. C(1)) GO TO 102 IF (B(K) .NE. B(1)) GO TO 102 101 CONTINUE GO TO 104 102 IERROR = 9 103 IF (NPEROD.EQ.1 .AND. (A(1).NE.0. .OR. C(N).NE.0.)) IERROR = 10 104 IF (IERROR .NE. 0) GO TO 122 IWYRT = L+1 IWT = IWYRT+M IWD = IWT+MAX0(L,M,N)+1 IWBB = IWD+N IWX = IWBB+N IWY = IWX+7*((L+1)/2)+15 GO TO (105,114),NP C C REORDER UNKNOWNS WHEN NPEROD = 0. C 105 NH = (N+1)/2 NHM1 = NH-1 NODD = 1 IF (2*NH .EQ. N) NODD = 2 DO 111 I=1,L DO 110 J=1,M DO 106 K=1,NHM1 NHPK = NH+K NHMK = NH-K W(K) = F(I,J,NHMK)-F(I,J,NHPK) W(NHPK) = F(I,J,NHMK)+F(I,J,NHPK) 106 CONTINUE W(NH) = 2.*F(I,J,NH) GO TO (108,107),NODD 107 W(N) = 2.*F(I,J,N) 108 DO 109 K=1,N F(I,J,K) = W(K) 109 CONTINUE 110 CONTINUE 111 CONTINUE SAVE(1) = C(NHM1) SAVE(2) = A(NH) SAVE(3) = C(NH) SAVE(4) = B(NHM1) SAVE(5) = B(N) SAVE(6) = A(N) C(NHM1) = 0. A(NH) = 0. C(NH) = 2.*C(NH) GO TO (112,113),NODD 112 B(NHM1) = B(NHM1)-A(NH-1) B(N) = B(N)+A(N) GO TO 114 113 A(N) = C(NH) 114 CONTINUE CALL POS3D1 (LP,L,MP,M,N,A,B,C,LDIMF,MDIMF,F,W,W(IWYRT),W(IWT), 1 W(IWD),W(IWX),W(IWY),C1,C2,W(IWBB)) GO TO (115,122),NP 115 DO 121 I=1,L DO 120 J=1,M DO 116 K=1,NHM1 NHMK = NH-K NHPK = NH+K W(NHMK) = .5*(F(I,J,NHPK)+F(I,J,K)) W(NHPK) = .5*(F(I,J,NHPK)-F(I,J,K)) 116 CONTINUE W(NH) = .5*F(I,J,NH) GO TO (118,117),NODD 117 W(N) = .5*F(I,J,N) 118 DO 119 K=1,N F(I,J,K) = W(K) 119 CONTINUE 120 CONTINUE 121 CONTINUE C(NHM1) = SAVE(1) A(NH) = SAVE(2) C(NH) = SAVE(3) B(NHM1) = SAVE(4) B(N) = SAVE(5) A(N) = SAVE(6) 122 CONTINUE RETURN END .