# poistg.f

          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