*DECK POISTG SUBROUTINE POISTG (NPEROD, N, MPEROD, M, A, B, C, IDIMY, Y, + IERROR, W) C***BEGIN PROLOGUE POISTG C***PURPOSE Solve a block tridiagonal system of linear equations C that results from a staggered grid finite difference C approximation to 2-D elliptic PDE's. C***LIBRARY SLATEC (FISHPACK) C***CATEGORY I2B4B C***TYPE SINGLE PRECISION (POISTG-S) C***KEYWORDS ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, TRIDIAGONAL C***AUTHOR Adams, J., (NCAR) C Swarztrauber, P. N., (NCAR) C Sweet, R., (NCAR) C***DESCRIPTION C C Subroutine POISTG solves the 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 The indices I+1 and I-1 are evaluated modulo M, i.e. C X(0,J) = X(M,J) and X(M+1,J) = X(1,J), and X(I,0) may be equal to C X(I,1) or -X(I,1) and X(I,N+1) may be equal to X(I,N) or -X(I,N) C depending on an input parameter. C C C * * * * * * * * Parameter Description * * * * * * * * * * C C * * * * * * On Input * * * * * * C C NPEROD C Indicates the values which X(I,0) and X(I,N+1) are assumed C to have. C = 1 If X(I,0) = -X(I,1) and X(I,N+1) = -X(I,N) C = 2 If X(I,0) = -X(I,1) and X(I,N+1) = X(I,N) C = 3 If X(I,0) = X(I,1) and X(I,N+1) = X(I,N) C = 4 If X(I,0) = X(I,1) and X(I,N+1) = -X(I,N) C C N C The number of unknowns in the J-direction. N must C 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. M must C be greater than 2. C C A,B,C C One-dimensional arrays of length M that specify the coefficients C in the linear equations given above. If MPEROD = 0 the array C elements must not depend on the index I, but must be constant. C Specifically, the subroutine checks the following condition C C A(I) = C(1) C B(I) = B(1) C C(I) = C(1) C C for I = 1, 2, ..., M. C C IDIMY C The row (or first) dimension of the two-dimensional array Y as C it appears in the program calling POISTG. 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 C right side of the linear system of equations given above. C Y must be dimensioned at least M X N. C C W C A one-dimensional work array that must be provided by the user C for work space. W may require up to 9M + 4N + M(INT(log2(N))) C locations. The actual number of locations used is computed by C POISTG and returned in location W(1). C 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 = 0 No error C = 1 If M .LE. 2 C = 2 If N .LE. 2 C = 3 IDIMY .LT. M C = 4 If NPEROD .LT. 1 or NPEROD .GT. 4 C = 5 If MPEROD .LT. 0 or MPEROD .GT. 1 C = 6 If MPEROD = 0 and C A(I) .NE. C(1) or B(I) .NE. B(1) or C(I) .NE. C(1) C for some I = 1, 2, ..., M. C = 7 If MPEROD .EQ. 1 .AND. (A(1).NE.0 .OR. C(M).NE.0) C C W C W(1) contains the required length of W. C C *Long Description: C C * * * * * * * Program Specifications * * * * * * * * * * * * C C Dimension of A(M),B(M),C(M),Y(IDIMY,N), C Arguments W(see argument list) C C Latest June 1, 1977 C Revision C C Subprograms POISTG,POSTG2,COSGEN,MERGE,TRIX,TRI3,PIMACH C Required 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 in 1973 C Revised by Roland Sweet in 1977 C C C Space 3297(decimal) = 6341(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 POISTG is roughly proportional C to M*N*log2(N). Some typical values are listed C in the table below. More comprehensive timing C charts may be found in the reference. 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(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 B(1) = B(M) =-1. 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 POISTG was C called to produce an approximate solution Z. Then C the relative error, defined as C C E = MAX(ABS(Z(I,J)-X(I,J)))/MAX(ABS(X(I,J))) C C where the two maxima are taken over all I=1,2,...,M C and J=1,2,...,N, was computed. The value of E is C given in the table below for some typical values of C M and N. C C C M (=N) MPEROD NPEROD T(MSECS) E C ------ ------ ------ -------- ------ C C 31 0-1 1-4 45 9.E-13 C 31 1 1 21 4.E-13 C 31 1 3 41 3.E-13 C 32 0-1 1-4 51 3.E-12 C 32 1 1 32 3.E-13 C 32 1 3 48 1.E-13 C 33 0-1 1-4 42 1.E-12 C 33 1 1 30 4.E-13 C 33 1 3 34 1.E-13 C 63 0-1 1-4 186 3.E-12 C 63 1 1 91 1.E-12 C 63 1 3 173 2.E-13 C 64 0-1 1-4 209 4.E-12 C 64 1 1 128 1.E-12 C 64 1 3 199 6.E-13 C 65 0-1 1-4 143 2.E-13 C 65 1 1 160 1.E-11 C 65 1 3 138 4.E-13 C C Portability American National Standards Institute FORTRAN. C The machine dependent constant PI is defined in C function PIMACH. C C Required COS C Resident C Routines C C Reference Schumann, U. and R. Sweet,'A Direct Method for C the Solution of Poisson's Equation With Neumann C Boundary Conditions on a Staggered Grid of C Arbitrary Size,' J. Comp. Phys. 20(1976), C pp. 171-182. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C***REFERENCES U. Schumann and R. Sweet, A direct method for the C solution of Poisson's equation with Neumann boundary C conditions on a staggered grid of arbitrary size, C Journal of Computational Physics 20, (1976), C pp. 171-182. C***ROUTINES CALLED POSTG2 C***REVISION HISTORY (YYMMDD) C 801001 DATE WRITTEN C 861211 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 POISTG