*DECK SEPELI
      SUBROUTINE SEPELI (INTL, IORDER, A, B, M, MBDCND, BDA, ALPHA, BDB,
     +   BETA, C, D, N, NBDCND, BDC, GAMA, BDD, XNU, COFX, COFY, GRHS,
     +   USOL, IDMN, W, PERTRB, IERROR)
C***BEGIN PROLOGUE  SEPELI
C***PURPOSE  Discretize and solve a second and, optionally, a fourth
C            order finite difference approximation on a uniform grid to
C            the general separable elliptic partial differential
C            equation on a rectangle with any combination of periodic or
C            mixed boundary conditions.
C***LIBRARY   SLATEC (FISHPACK)
C***CATEGORY  I2B1A2
C***TYPE      SINGLE PRECISION (SEPELI-S)
C***KEYWORDS  ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, SEPARABLE
C***AUTHOR  Adams, J., (NCAR)
C           Swarztrauber, P. N., (NCAR)
C           Sweet, R., (NCAR)
C***DESCRIPTION
C
C Dimension of           BDA(N+1), BDB(N+1), BDC(M+1), BDD(M+1),
C Arguments              USOL(IDMN,N+1), GRHS(IDMN,N+1),
C                        W (see argument list)
C
C Latest Revision        March 1977
C
C Purpose                SEPELI solves for either the second-order
C                        finite difference approximation or a
C                        fourth-order approximation to a separable
C                        elliptic equation.
C
C                                    2    2
C                             AF(X)*d U/dX + BF(X)*dU/dX  + CF(X)*U +
C                                    2    2
C                             DF(Y)*d U/dY  + EF(Y)*dU/dY + FF(Y)*U
C
C                             = G(X,Y)
C
C                        on a rectangle (X greater than or equal to A
C                        and less than or equal to B; Y greater than
C                        or equal to C and less than or equal to D).
C                        Any combination of periodic or mixed boundary
C                        conditions is allowed.
C
C Purpose                The possible boundary conditions are:
C                        in the X-direction:
C                         (0) Periodic, U(X+B-A,Y)=U(X,Y) for all Y,X
C                         (1) U(A,Y), U(B,Y) are specified for all Y
C                         (2) U(A,Y), dU(B,Y)/dX+BETA*U(B,Y) are
C                             specified for all Y
C                         (3) dU(A,Y)/dX+ALPHA*U(A,Y),dU(B,Y)/dX+
C                             BETA*U(B,Y) are specified for all Y
C                         (4) dU(A,Y)/dX+ALPHA*U(A,Y),U(B,Y) are
C                             specified for all Y
C
C                        in the Y-direction:
C                         (0) Periodic, U(X,Y+D-C)=U(X,Y) for all X,Y
C                         (1) U(X,C),U(X,D) are specified for all X
C                         (2) U(X,C),dU(X,D)/dY+XNU*U(X,D) are specified
C                             for all X
C                         (3) dU(X,C)/dY+GAMA*U(X,C),dU(X,D)/dY+
C                             XNU*U(X,D) are specified for all X
C                         (4) dU(X,C)/dY+GAMA*U(X,C),U(X,D) are
C                             specified for all X
C
C Arguments
C
C On Input               INTL
C                          = 0 On initial entry to SEPELI or if any of
C                              the arguments C, D, N, NBDCND, COFY are
C                              changed from a previous call
C                          = 1 If C, D, N, NBDCND, COFY are unchanged
C                              from the previous call.
C
C                        IORDER
C                          = 2 If a second-order approximation is sought
C                          = 4 If a fourth-order approximation is sought
C
C                        A,B
C                          The range of the X-independent variable;
C                          i.e., X is greater than or equal to A and
C                          less than or equal to B.  A must be less than
C                          B.
C
C                        M
C                          The number of panels into which the interval
C                          [A,B] is subdivided.  Hence, there will be
C                          M+1 grid points in the X-direction given by
C                          XI=A+(I-1)*DLX for I=1,2,...,M+1 where
C                          DLX=(B-A)/M is the panel width.  M must be
C                          less than IDMN and greater than 5.
C
C                        MBDCND
C                          Indicates the type of boundary condition at
C                          X=A and X=B
C                          = 0 If the solution is periodic in X; i.e.,
C                              U(X+B-A,Y)=U(X,Y) for all Y,X
C                          = 1 If the solution is specified at X=A and
C                              X=B; i.e., U(A,Y) and U(B,Y) are
C                              specified for all Y
C                          = 2 If the solution is specified at X=A and
C                              the boundary condition is mixed at X=B;
C                              i.e., U(A,Y) and dU(B,Y)/dX+BETA*U(B,Y)
C                              are specified for all Y
C                          = 3 If the boundary conditions at X=A and X=B
C                              are mixed; i.e., dU(A,Y)/dX+ALPHA*U(A,Y)
C                              and dU(B,Y)/dX+BETA*U(B,Y) are specified
C                              for all Y
C                          = 4 If the boundary condition at X=A is mixed
C                              and the solution is specified at X=B;
C                              i.e., dU(A,Y)/dX+ALPHA*U(A,Y) and U(B,Y)
C                              are specified for all Y
C
C                        BDA
C                          A one-dimensional array of length N+1 that
C                          specifies the values of dU(A,Y)/dX+
C                          ALPHA*U(A,Y) at X=A, when MBDCND=3 or 4.
C                               BDA(J) = dU(A,YJ)/dX+ALPHA*U(A,YJ);
C                               J=1,2,...,N+1
C                          when MBDCND has any other value, BDA is a
C                          dummy parameter.
C
C On Input               ALPHA
C                          The scalar multiplying the solution in case
C                          of a mixed boundary condition at X=A (see
C                          argument BDA).  If MBDCND = 3,4 then ALPHA is
C                          a dummy parameter.
C
C                        BDB
C                          A one-dimensional array of length N+1 that
C                          specifies the values of dU(B,Y)/dX+
C                          BETA*U(B,Y) at X=B.  When MBDCND=2 or 3
C                               BDB(J) = dU(B,YJ)/dX+BETA*U(B,YJ);
C                               J=1,2,...,N+1
C                          When MBDCND has any other value, BDB is a
C                          dummy parameter.
C
C                        BETA
C                          The scalar multiplying the solution in case
C                          of a mixed boundary condition at X=B (see
C                          argument BDB).  If MBDCND=2,3 then BETA is a
C                          dummy parameter.
C
C                        C,D
C                          The range of the Y-independent variable;
C                          i.e., Y is greater than or equal to C and
C                          less than or equal to D.  C must be less than
C                          D.
C
C                        N
C                          The number of panels into which the interval
C                          [C,D] is subdivided.  Hence, there will be
C                          N+1 grid points in the Y-direction given by
C                          YJ=C+(J-1)*DLY for J=1,2,...,N+1 where
C                          DLY=(D-C)/N is the panel width.  In addition,
C                          N must be greater than 4.
C
C                        NBDCND
C                          Indicates the types of boundary conditions at
C                          Y=C and Y=D
C                          = 0 If the solution is periodic in Y; i.e.,
C                              U(X,Y+D-C)=U(X,Y) for all X,Y
C                          = 1 If the solution is specified at Y=C and
C                              Y = D, i.e., U(X,C) and U(X,D) are
C                              specified for all X
C                          = 2 If the solution is specified at Y=C and
C                              the boundary condition is mixed at Y=D;
C                              i.e., U(X,C) and dU(X,D)/dY+XNU*U(X,D)
C                              are specified for all X
C                          = 3 If the boundary conditions are mixed at
C                              Y=C and Y=D; i.e., dU(X,D)/dY+GAMA*U(X,C)
C                              and dU(X,D)/dY+XNU*U(X,D) are specified
C                              for all X
C                          = 4 If the boundary condition is mixed at Y=C
C                              and the solution is specified at Y=D;
C                              i.e. dU(X,C)/dY+GAMA*U(X,C) and U(X,D)
C                              are specified for all X
C
C                        BDC
C                          A one-dimensional array of length M+1 that
C                          specifies the value of dU(X,C)/dY+GAMA*U(X,C)
C                          at Y=C.  When NBDCND=3 or 4
C                             BDC(I) = dU(XI,C)/dY + GAMA*U(XI,C);
C                             I=1,2,...,M+1.
C                          When NBDCND has any other value, BDC is a
C                          dummy parameter.
C
C                        GAMA
C                          The scalar multiplying the solution in case
C                          of a mixed boundary condition at Y=C (see
C                          argument BDC).  If NBDCND=3,4 then GAMA is a
C                          dummy parameter.
C
C                        BDD
C                          A one-dimensional array of length M+1 that
C                          specifies the value of dU(X,D)/dY +
C                          XNU*U(X,D) at Y=C.  When NBDCND=2 or 3
C                            BDD(I) = dU(XI,D)/dY + XNU*U(XI,D);
C                            I=1,2,...,M+1.
C                          When NBDCND has any other value, BDD is a
C                          dummy parameter.
C
C                        XNU
C                          The scalar multiplying the solution in case
C                          of a mixed boundary condition at Y=D (see
C                          argument BDD).  If NBDCND=2 or 3 then XNU is
C                          a dummy parameter.
C
C                        COFX
C                          A user-supplied subprogram with
C                          parameters X, AFUN, BFUN, CFUN which
C                          returns the values of the X-dependent
C                          coefficients AF(X), BF(X), CF(X) in
C                          the elliptic equation at X.
C
C                        COFY
C                          A user-supplied subprogram with
C                          parameters Y, DFUN, EFUN, FFUN which
C                          returns the values of the Y-dependent
C                          coefficients DF(Y), EF(Y), FF(Y) in
C                          the elliptic equation at Y.
C
C                        NOTE:  COFX and COFY must be declared external
C                        in the calling routine.  The values returned in
C                        AFUN and DFUN must satisfy AFUN*DFUN greater
C                        than 0 for A less than X less than B,
C                        C less than Y less than D (see IERROR=10).
C                        The coefficients provided may lead to a matrix
C                        equation which is not diagonally dominant in
C                        which case solution may fail (see IERROR=4).
C
C                        GRHS
C                          A two-dimensional array that specifies the
C                          values of the right-hand side of the elliptic
C                          equation; i.e., GRHS(I,J)=G(XI,YI), for
C                          I=2,...,M; J=2,...,N.  At the boundaries,
C                          GRHS is defined by
C
C                          MBDCND   GRHS(1,J)   GRHS(M+1,J)
C                          ------   ---------   -----------
C                            0      G(A,YJ)     G(B,YJ)
C                            1         *           *
C                            2         *        G(B,YJ)  J=1,2,...,N+1
C                            3      G(A,YJ)     G(B,YJ)
C                            4      G(A,YJ)        *
C
C                          NBDCND   GRHS(I,1)   GRHS(I,N+1)
C                          ------   ---------   -----------
C                            0      G(XI,C)     G(XI,D)
C                            1         *           *
C                            2         *        G(XI,D)  I=1,2,...,M+1
C                            3      G(XI,C)     G(XI,D)
C                            4      G(XI,C)        *
C
C                          where * means these quantities are not used.
C                          GRHS should be dimensioned IDMN by at least
C                          N+1 in the calling routine.
C
C                        USOL
C                          A two-dimensional array that specifies the
C                          values of the solution along the boundaries.
C                          At the boundaries, USOL is defined by
C
C                          MBDCND   USOL(1,J)   USOL(M+1,J)
C                          ------   ---------   -----------
C                            0         *           *
C                            1      U(A,YJ)     U(B,YJ)
C                            2      U(A,YJ)        *     J=1,2,...,N+1
C                            3         *           *
C                            4         *        U(B,YJ)
C
C                          NBDCND   USOL(I,1)   USOL(I,N+1)
C                          ------   ---------   -----------
C                            0         *           *
C                            1      U(XI,C)     U(XI,D)
C                            2      U(XI,C)        *     I=1,2,...,M+1
C                            3         *           *
C                            4         *        U(XI,D)
C
C                          where * means the quantities are not used in
C                          the solution.
C
C                          If IORDER=2, the user may equivalence GRHS
C                          and USOL to save space.  Note that in this
C                          case the tables specifying the boundaries of
C                          the GRHS and USOL arrays determine the
C                          boundaries uniquely except at the corners.
C                          If the tables call for both G(X,Y) and
C                          U(X,Y) at a corner then the solution must be
C                          chosen.  For example, if MBDCND=2 and
C                          NBDCND=4, then U(A,C), U(A,D), U(B,D) must be
C                          chosen at the corners in addition to G(B,C).
C
C                          If IORDER=4, then the two arrays, USOL and
C                          GRHS, must be distinct.
C
C                          USOL should be dimensioned IDMN by at least
C                          N+1 in the calling routine.
C
C                        IDMN
C                          The row (or first) dimension of the arrays
C                          GRHS and USOL as it appears in the program
C                          calling SEPELI.  This parameter is used to
C                          specify the variable dimension of GRHS and
C                          USOL.  IDMN must be at least 7 and greater
C                          than or equal to M+1.
C
C                        W
C                          A one-dimensional array that must be provided
C                          by the user for work space.  Let
C                          K=INT(log2(N+1))+1 and set  L=2**(K+1).
C                          then (K-2)*L+K+10*N+12*M+27 will suffice
C                          as a length of W.  THE actual length of W in
C                          the calling routine must be set in W(1) (see
C                          IERROR=11).
C
C On Output              USOL
C                          Contains the approximate solution to the
C                          elliptic equation.  USOL(I,J) is the
C                          approximation to U(XI,YJ) for I=1,2...,M+1
C                          and J=1,2,...,N+1.  The approximation has
C                          error O(DLX**2+DLY**2) if called with
C                          IORDER=2 and O(DLX**4+DLY**4) if called with
C                          IORDER=4.
C
C                        W
C                          Contains intermediate values that must not be
C                          destroyed if SEPELI is called again with
C                          INTL=1.  In addition W(1) contains the exact
C                          minimal length (in floating point) required
C                          for the work space (see IERROR=11).
C
C                        PERTRB
C                          If a combination of periodic or derivative
C                          boundary conditions (i.e., ALPHA=BETA=0 if
C                          MBDCND=3; GAMA=XNU=0 if NBDCND=3) is
C                          specified and if the coefficients of U(X,Y)
C                          in the separable elliptic equation are zero
C                          (i.e., CF(X)=0 for X greater than or equal to
C                          A and less than or equal to B; FF(Y)=0 for
C                          Y greater than or equal to C and less than
C                          or equal to D) then a solution may not exist.
C                          PERTRB is a constant calculated and
C                          subtracted from the right-hand side of the
C                          matrix equations generated by SEPELI which
C                          insures that a solution exists.  SEPELI then
C                          computes this solution which is a weighted
C                          minimal least squares solution to the
C                          original problem.
C
C                        IERROR
C                          An error flag that indicates invalid input
C                          parameters or failure to find a solution
C                          = 0 No error
C                          = 1 If A greater than B or C greater than D
C                          = 2 If MBDCND less than 0 or MBDCND greater
C                              than 4
C                          = 3 If NBDCND less than 0 or NBDCND greater
C                              than 4
C                          = 4 If attempt to find a solution fails.
C                              (the linear system generated is not
C                              diagonally dominant.)
C                          = 5 If IDMN is too small (see discussion of
C                              IDMN)
C                          = 6 If M is too small or too large (see
C                              discussion of M)
C                          = 7 If N is too small (see discussion of N)
C                          = 8 If IORDER is not 2 or 4
C                          = 9 If INTL is not 0 or 1
C                          = 10 If AFUN*DFUN less than or equal to 0 for
C                               some interior mesh point (XI,YJ)
C                          = 11 If the work space length input in W(1)
C                               is less than the exact minimal work
C                               space length required output in W(1).
C
C                          NOTE (concerning IERROR=4):  for the
C                          coefficients input through COFX, COFY, the
C                          discretization may lead to a block
C                          tridiagonal linear system which is not
C                          diagonally dominant (for example, this
C                          happens if CFUN=0 and BFUN/(2.*DLX) greater
C                          than AFUN/DLX**2).  In this case solution may
C                          fail.  This cannot happen in the limit as
C                          DLX, DLY approach zero.  Hence, the condition
C                          may be remedied by taking larger values for M
C                          or N.
C
C Entry Points           SEPELI, SPELIP, CHKPRM, CHKSNG, ORTHOG, MINSOL,
C                        TRISP, DEFER, DX, DY, BLKTRI, BLKTR1, INDXB,
C                        INDXA, INDXC, PROD, PRODP, CPROD, CPRODP,
C                        PPADD, PSGF, BSRH, PPSGF, PPSPF, COMPB,
C                        TRUN1, STOR1, TQLRAT
C
C Special Conditions     NONE
C
C Common Blocks          SPLP, CBLKT
C
C I/O                    NONE
C
C Precision              Single
C
C Specialist             John C. Adams, NCAR, Boulder, Colorado  80307
C
C Language               FORTRAN
C
C History                Developed at NCAR during 1975-76.
C
C Algorithm              SEPELI automatically discretizes the separable
C                        elliptic equation which is then solved by a
C                        generalized cyclic reduction algorithm in the
C                        subroutine, BLKTRI.  The fourth-order solution
C                        is obtained using 'Deferred Corrections' which
C                        is described and referenced in sections,
C                        references and method.
C
C Space Required         14654 (octal) = 6572 (decimal)
C
C Accuracy and Timing    The following computational results were
C                        obtained by solving the sample problem at the
C                        end of this write-up on the Control Data 7600.
C                        The op count is proportional to M*N*log2(N).
C                        In contrast to the other routines in this
C                        chapter, accuracy is tested by computing and
C                        tabulating second- and fourth-order
C                        discretization errors.  Below is a table
C                        containing computational results.  The times
C                        given do not include initialization (i.e.,
C                        times are for INTL=1).  Note that the
C                        fourth-order accuracy is not realized until the
C                        mesh is sufficiently refined.
C
C              Second-order    Fourth-order   Second-order  Fourth-order
C    M    N   Execution Time  Execution Time    Error         Error
C               (M SEC)         (M SEC)
C     6    6         6              14          6.8E-1        1.2E0
C    14   14        23              58          1.4E-1        1.8E-1
C    30   30       100             247          3.2E-2        9.7E-3
C    62   62       445           1,091          7.5E-3        3.0E-4
C   126  126     2,002           4,772          1.8E-3        3.5E-6
C
C Portability            There are no machine-dependent constants.
C
C Required Resident      SQRT, ABS, LOG
C Routines
C
C References             Keller, H.B., 'Numerical Methods for Two-point
C                          Boundary-value Problems', Blaisdel (1968),
C                          Waltham, Mass.
C
C                        Swarztrauber, P., and R. Sweet (1975):
C                          'Efficient FORTRAN Subprograms for The
C                          Solution of Elliptic Partial Differential
C                          Equations'.  NCAR Technical Note
C                          NCAR-TN/IA-109, pp. 135-137.
C
C***REFERENCES  H. B. Keller, Numerical Methods for Two-point
C                 Boundary-value Problems, Blaisdel, Waltham, Mass.,
C                 1968.
C               P. N. Swarztrauber and R. Sweet, Efficient Fortran
C                 subprograms for the solution of elliptic equations,
C                 NCAR TN/IA-109, July 1975, 138 pp.
C***ROUTINES CALLED  CHKPRM, SPELIP
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   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  SEPELI
C
      DIMENSION       GRHS(IDMN,*)           ,USOL(IDMN,*)
      DIMENSION       BDA(*)     ,BDB(*)     ,BDC(*)     ,BDD(*)     ,
     1                W(*)
      EXTERNAL        COFX       ,COFY
C***FIRST EXECUTABLE STATEMENT  SEPELI
      CALL CHKPRM (INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,COFX,COFY,
     1             IDMN,IERROR)
      IF (IERROR .NE. 0) RETURN
C
C     COMPUTE MINIMUM WORK SPACE AND CHECK WORK SPACE LENGTH INPUT
C
      L = N+1
      IF (NBDCND .EQ. 0) L = N
      LOGB2N = INT(LOG(L+0.5)/LOG(2.0))+1
      LL = 2**(LOGB2N+1)
      K = M+1
      L = N+1
      LENGTH = (LOGB2N-2)*LL+LOGB2N+MAX(2*L,6*K)+5
      IF (NBDCND .EQ. 0) LENGTH = LENGTH+2*L
      IERROR = 11
      LINPUT = INT(W(1)+0.5)
      LOUTPT = LENGTH+6*(K+L)+1
      W(1) = LOUTPT
      IF (LOUTPT .GT. LINPUT) RETURN
      IERROR = 0
C
C     SET WORK SPACE INDICES
C
      I1 = LENGTH+2
      I2 = I1+L
      I3 = I2+L
      I4 = I3+L
      I5 = I4+L
      I6 = I5+L
      I7 = I6+L
      I8 = I7+K
      I9 = I8+K
      I10 = I9+K
      I11 = I10+K
      I12 = I11+K
      I13 = 2
      CALL SPELIP (INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,N,
     1             NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,W(I1),W(I2),W(I3),
     2             W(I4),W(I5),W(I6),W(I7),W(I8),W(I9),W(I10),W(I11),
     3             W(I12),GRHS,USOL,IDMN,W(I13),PERTRB,IERROR)
      RETURN
      END