SUBROUTINE GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,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 A(M),B(M),C(M),Y(IDIMY,N), C W(SEE PARAMETER LIST) C ARGUMENTS C C LATEST REVISION NOVEMBER 1988 C C PURPOSE THE NAME OF THIS PACKAGE IS A MNEMONIC FOR THE C GENERALIZED BUNEMAN ALGORITHM. C C IT SOLVES THE REAL LINEAR SYSTEM OF EQUATIONS C C A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J) 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 INDICES I+1 AND I-1 ARE EVALUATED MODULO M, C I.E., X(0,J) = X(M,J) AND X(M+1,J) = X(1,J), C AND X(I,0) MAY EQUAL 0, X(I,2), OR X(I,N), C AND X(I,N+1) MAY EQUAL 0, X(I,N-1), OR X(I,1) C DEPENDING ON AN INPUT PARAMETER. C C USAGE CALL GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y, C IERROR,W) C C ARGUMENTS C C ON INPUT NPEROD C C INDICATES THE VALUES THAT X(I,0) AND C X(I,N+1) ARE ASSUMED TO HAVE. C C = 0 IF X(I,0) = X(I,N) AND X(I,N+1) = C 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) = C 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. C N MUST BE GREATER 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. C N MUST BE GREATER THAN 2. C C A,B,C C ONE-DIMENSIONAL ARRAYS OF LENGTH M THAT C SPECIFY THE COEFFICIENTS IN THE LINEAR C EQUATIONS GIVEN ABOVE. IF MPEROD = 0 C THE ARRAY ELEMENTS MUST NOT DEPEND UPON C THE INDEX I, BUT MUST BE CONSTANT. C SPECIFICALLY, THE SUBROUTINE CHECKS THE C FOLLOWING 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 C TWO-DIMENSIONAL ARRAY Y AS IT APPEARS C IN THE PROGRAM CALLING GENBUN. 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 TWO-DIMENSIONAL COMPLEX ARRAY THAT C SPECIFIES THE VALUES OF THE RIGHT SIDE C OF THE LINEAR SYSTEM OF EQUATIONS GIVEN C ABOVE. C Y MUST BE DIMENSIONED AT LEAST M*N. C C W C A ONE-DIMENSIONAL ARRAY THAT MUST C BE PROVIDED BY THE USER FOR WORK C SPACE. W MAY REQUIRE UP TO 4*N + C (10 + INT(LOG2(N)))*M LOCATIONS. C THE ACTUAL NUMBER OF LOCATIONS USED IS C COMPUTED BY GENBUN AND IS RETURNED IN C LOCATION W(1). C C C ON OUTPUT Y C C CONTAINS THE SOLUTION X. C C IERROR C AN ERROR FLAG WHICH INDICATES INVALID C INPUT PARAMETERS EXCEPT FOR NUMBER C 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 C B(I) .NE. B(1) FOR C SOME I=1,2,...,M. C = 7 A(1) .NE. 0 OR C(M) .NE. 0 AND C MPEROD = 1 C C W C W(1) CONTAINS THE REQUIRED LENGTH OF W. C C SPECIAL CONDITONS NONE C C I/O NONE C C PRECISION SINGLE C C REQUIRED LIBRARY COMF AND GNBNAUX FROM FISHPAK C FILES C C LANGUAGE FORTRAN C C HISTORY WRITTEN IN 1979 BY ROLAND SWEET OF NCAR'S C SCIENTIFIC COMPUTING DIVISION. MADE AVAILABLE C ON NCAR'S PUBLIC LIBRARIES IN JANUARY, 1980. C C ALGORITHM THE LINEAR SYSTEM IS SOLVED BY A CYCLIC C REDUCTION ALGORITHM DESCRIBED IN THE C REFERENCE. C C PORTABILITY FORTRAN 77 -- C THE MACHINE DEPENDENT CONSTANT PI IS C DEFINED IN FUNCTION PIMACH. C C REFERENCES SWEET, R., "A CYCLIC REDUCTION ALGORITHM FOR C SOLVING BLOCK TRIDIAGONAL SYSTEMS OF ARBITRARY C DIMENSIONS," SIAM J. ON NUMER. ANAL., 14 (1977) C PP. 706-720. C C ACCURACY THIS TEST WAS PERFORMED ON A CDC 7600: C C A UNIFORM RANDOM NUMBER GENERATOR WAS USED C TO CREATE A SOLUTION ARRAY X FOR THE SYSTEM C GIVEN IN THE 'PURPOSE' DESCRIPTION ABOVE C WITH 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 C GIVEN SYSTEM AND, USING DOUBLE PRECISION C A RIGHT SIDE Y WAS COMPUTED. C USING THIS ARRAY Y, SUBROUTINE GENBUN C WAS CALLED TO PRODUCE APPROXIMATE C SOLUTION Z. THEN RELATIVE ERROR C E = MAX(ABS(Z(I,J)-X(I,J)))/ C MAX(ABS(X(I,J))) C WAS COMPUTED, WHERE THE TWO MAXIMA ARE TAKEN C OVER I=1,2,...,M AND J=1,...,N. C C THE VALUE OF E IS GIVEN IN THE TABLE C BELOW FOR SOME TYPICAL VALUES OF M AND N. C C M (=N) MPEROD NPEROD T(MSECS) E C ------ ------ ------ -------- ------ C C 31 0 0 36 6.E-14 C 31 1 1 21 4.E-13 C 31 1 3 41 3.E-13 C 32 0 0 29 9.E-14 C 32 1 1 32 3.E-13 C 32 1 3 48 1.E-13 C 33 0 0 36 9.E-14 C 33 1 1 30 4.E-13 C 33 1 3 34 1.E-13 C 63 0 0 150 1.E-13 C 63 1 1 91 1.E-12 C 63 1 3 173 2.E-13 C 64 0 0 122 1.E-13 C 64 1 1 128 1.E-12 C 64 1 3 199 6.E-13 C 65 0 0 143 2.E-13 C 65 1 1 120 1.E-12 C 65 1 3 138 4.E-13 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *