relibc_openlibm/slatec/blktri.f

*DECK BLKTRI
      SUBROUTINE BLKTRI (IFLG, NP, N, AN, BN, CN, MP, M, AM, BM, CM,
     +   IDIMY, Y, IERROR, W)
C***BEGIN PROLOGUE  BLKTRI
C***PURPOSE  Solve a block tridiagonal system of linear equations
C            (usually resulting from the discretization of separable
C            two-dimensional elliptic equations).
C***LIBRARY   SLATEC (FISHPACK)
C***CATEGORY  I2B4B
C***TYPE      SINGLE PRECISION (BLKTRI-S, CBLKTR-C)
C***KEYWORDS  ELLIPTIC PDE, FISHPACK, TRIDIAGONAL LINEAR SYSTEM
C***AUTHOR  Adams, J., (NCAR)
C           Swarztrauber, P. N., (NCAR)
C           Sweet, R., (NCAR)
C***DESCRIPTION
C
C     Subroutine BLKTRI Solves a System of Linear Equations of the Form
C
C          AN(J)*X(I,J-1) + AM(I)*X(I-1,J) + (BN(J)+BM(I))*X(I,J)
C
C          + CN(J)*X(I,J+1) + CM(I)*X(I+1,J) = Y(I,J)
C
C               for I = 1,2,...,M  and  J = 1,2,...,N.
C
C     I+1 and I-1 are evaluated modulo M and J+1 and J-1 modulo N, i.e.,
C
C          X(I,0) = X(I,N),  X(I,N+1) = X(I,1),
C          X(0,J) = X(M,J),  X(M+1,J) = X(1,J).
C
C     These equations usually result from the discretization of
C     separable elliptic equations.  Boundary conditions may be
C     Dirichlet, Neumann, or Periodic.
C
C
C     * * * * * * * * * *     ON INPUT     * * * * * * * * * *
C
C     IFLG
C       = 0  Initialization only.  Certain quantities that depend on NP,
C            N, AN, BN, and CN are computed and stored in the work
C            array  W.
C       = 1  The quantities that were computed in the initialization are
C            used to obtain the solution X(I,J).
C
C       NOTE   A call with IFLG=0 takes approximately one half the time
C              as a call with IFLG = 1  .  However, the
C              initialization does not have to be repeated unless NP, N,
C              AN, BN, or CN change.
C
C     NP
C       = 0  If AN(1) and CN(N) are not zero, which corresponds to
C            periodic boundary conditions.
C       = 1  If AN(1) and CN(N) are zero.
C
C     N
C       The number of unknowns in the J-direction. N must be greater
C       than 4. The operation count is proportional to MNlog2(N), hence
C       N should be selected less than or equal to M.
C
C     AN,BN,CN
C       One-dimensional arrays of length N that specify the coefficients
C       in the linear equations given above.
C
C     MP
C       = 0  If AM(1) and CM(M) are not zero, which corresponds to
C            periodic boundary conditions.
C       = 1  If AM(1) = CM(M) = 0  .
C
C     M
C       The number of unknowns in the I-direction. M must be greater
C       than 4.
C
C     AM,BM,CM
C       One-dimensional arrays of length M that specify the coefficients
C       in the linear equations given above.
C
C     IDIMY
C       The row (or first) dimension of the two-dimensional array Y as
C       it appears in the program calling BLKTRI.  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 right
C       side of the linear system of equations given above.  Y must be
C       dimensioned at least M*N.
C
C     W
C       A one-dimensional array that must be provided by the user for
C       work space.
C             If NP=1 define K=INT(log2(N))+1 and set L=2**(K+1) then
C                     W must have dimension (K-2)*L+K+5+MAX(2N,6M)
C
C             If NP=0 define K=INT(log2(N-1))+1 and set L=2**(K+1) then
C                     W must have dimension (K-2)*L+K+5+2N+MAX(2N,6M)
C
C       **IMPORTANT** For purposes of checking, the required dimension
C                     of W is computed by BLKTRI and stored in W(1)
C                     in floating point format.
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
C       = 0  No error.
C       = 1  M is less than 5.
C       = 2  N is less than 5.
C       = 3  IDIMY is less than M.
C       = 4  BLKTRI failed while computing results that depend on the
C            coefficient arrays AN, BN, CN.  Check these arrays.
C       = 5  AN(J)*CN(J-1) is less than 0 for some J. Possible reasons
C            for this condition are
C            1. The arrays AN and CN are not correct.
C            2. Too large a grid spacing was used in the discretization
C               of the elliptic equation.
C            3. The linear equations resulted from a partial
C               differential equation which was not elliptic.
C
C     W
C       Contains intermediate values that must not be destroyed if
C       BLKTRI will be called again with IFLG=1.  W(1) contains the
C       number of locations required by W in floating point format.
C
C *Long Description:
C
C     * * * * * * *   Program Specifications    * * * * * * * * * * * *
C
C     Dimension of   AN(N),BN(N),CN(N),AM(M),BM(M),CM(M),Y(IDIMY,N)
C     Arguments      W(See argument list)
C
C     Latest         June 1979
C     Revision
C
C     Required       BLKTRI,BLKTRI,PROD,PRODP,CPROD,CPRODP,COMPB,INDXA,
C     Subprograms    INDXB,INDXC,PPADD,PSGF,PPSGF,PPSPF,BSRH,TEVLS,
C                    R1MACH
C
C     Special        The Algorithm may fail if ABS(BM(I)+BN(J)) is less
C     Conditions     than ABS(AM(I))+ABS(AN(J))+ABS(CM(I))+ABS(CN(J))
C                    for some I and J. The Algorithm will also fail if
C                    AN(J)*CN(J-1) is less than zero for some J.
C                    See the description of the output parameter IERROR.
C
C     Common         CBLKT
C     Blocks
C
C     I/O            None
C
C     Precision      Single
C
C     Specialist     Paul Swarztrauber
C
C     Language       FORTRAN
C
C     History        Version 1 September 1973
C                    Version 2 April     1976
C                    Version 3 June      1979
C
C     Algorithm      Generalized Cyclic Reduction (See Reference below)
C
C     Space
C     Required       Control Data 7600
C
C     Portability    American National Standards Institute Fortran.
C                    The machine accuracy is set using function R1MACH.
C
C     Required       None
C     Resident
C     Routines
C
C     References     Swarztrauber,P. and R. Sweet, 'Efficient FORTRAN
C                    Subprograms For The Solution Of Elliptic Equations'
C                    NCAR TN/IA-109, July, 1975, 138 PP.
C
C                    Swarztrauber P. ,'A Direct Method For The Discrete
C                    Solution Of Separable Elliptic Equations', S.I.A.M.
C                    J. Numer. Anal.,11(1974) PP. 1136-1150.
C
C***REFERENCES  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               P. N. Swarztrauber, A direct method for the discrete
C                 solution of separable elliptic equations, SIAM Journal
C                 on Numerical Analysis 11, (1974), pp. 1136-1150.
C***ROUTINES CALLED  BLKTR1, COMPB, CPROD, CPRODP, PROD, PRODP
C***COMMON BLOCKS    CBLKT
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  BLKTRI
C
      DIMENSION       AN(*)      ,BN(*)      ,CN(*)      ,AM(*)      ,
     1                BM(*)      ,CM(*)      ,Y(IDIMY,*) ,W(*)
      EXTERNAL        PROD       ,PRODP      ,CPROD      ,CPRODP
      COMMON /CBLKT/  NPP        ,K          ,EPS        ,CNV        ,
     1                NM         ,NCMPLX     ,IK
C***FIRST EXECUTABLE STATEMENT  BLKTRI
      NM = N
      IERROR = 0
      IF (M-5) 101,102,102
  101 IERROR = 1
      GO TO 119
  102 IF (NM-3) 103,104,104
  103 IERROR = 2
      GO TO 119
  104 IF (IDIMY-M) 105,106,106
  105 IERROR = 3
      GO TO 119
  106 NH = N
      NPP = NP
      IF (NPP) 107,108,107
  107 NH = NH+1
  108 IK = 2
      K = 1
  109 IK = IK+IK
      K = K+1
      IF (NH-IK) 110,110,109
  110 NL = IK
      IK = IK+IK
      NL = NL-1
      IWAH = (K-2)*IK+K+6
      IF (NPP) 111,112,111
C
C     DIVIDE W INTO WORKING SUB ARRAYS
C
  111 IW1 = IWAH
      IWBH = IW1+NM
      W(1) = IW1-1+MAX(2*NM,6*M)
      GO TO 113
  112 IWBH = IWAH+NM+NM
      IW1 = IWBH
      W(1) = IW1-1+MAX(2*NM,6*M)
      NM = NM-1
C
C SUBROUTINE COMP B COMPUTES THE ROOTS OF THE B POLYNOMIALS
C
  113 IF (IERROR) 119,114,119
  114 IW2 = IW1+M
      IW3 = IW2+M
      IWD = IW3+M
      IWW = IWD+M
      IWU = IWW+M
      IF (IFLG) 116,115,116
  115 CALL COMPB (NL,IERROR,AN,BN,CN,W(2),W(IWAH),W(IWBH))
      GO TO 119
  116 IF (MP) 117,118,117
C
C SUBROUTINE BLKTR1 SOLVES THE LINEAR SYSTEM
C
  117 CALL BLKTR1 (NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),W(IW1),W(IW2),
     1             W(IW3),W(IWD),W(IWW),W(IWU),PROD,CPROD)
      GO TO 119
  118 CALL BLKTR1 (NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),W(IW1),W(IW2),
     1             W(IW3),W(IWD),W(IWW),W(IWU),PRODP,CPRODP)
  119 CONTINUE
      RETURN
      END