*DECK POIS3D SUBROUTINE POIS3D (LPEROD, L, C1, MPEROD, M, C2, NPEROD, N, A, B, + C, LDIMF, MDIMF, F, IERROR, W) C***BEGIN PROLOGUE POIS3D C***PURPOSE Solve a three-dimensional block tridiagonal linear system C which arises from a finite difference approximation to a C three-dimensional Poisson equation using the Fourier C transform package FFTPAK written by Paul Swarztrauber. C***LIBRARY SLATEC (FISHPACK) C***CATEGORY I2B4B C***TYPE SINGLE PRECISION (POIS3D-S) C***KEYWORDS ELLIPTIC PDE, FISHPACK, HELMHOLTZ, POISSON C***AUTHOR Adams, J., (NCAR) C Swarztrauber, P. N., (NCAR) C Sweet, R., (NCAR) C***DESCRIPTION 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 *Long Description: 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,TRIDQ,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 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 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 FFTPAK 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 POIS3D 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 C***REFERENCES (NONE) C***ROUTINES CALLED POS3D1 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***END PROLOGUE POIS3D