C      ALGORITHM 685, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 16, NO. 4, PP. 325-351.
C-----TOMS SERRG2 (MAIN PORTION OF ALGORITHM)
      SUBROUTINE SERRG2(K, M, X, N, Y, P1, R, Q1, P2, Q2, F, V,
     1   IA, BETAA, GAMA, IB, BETAB, GAMB, IC, BETAC, GAMC, ID, BETAD,
     1   GAMD)
C
C     SERRG2  -  SEPARABLE ELLIPTIC RAYLEIGH-RITZ-GALERKIN TWO
C                 DIMENSIONAL PARTIAL DIFFERENTIAL EQUATION SOLVER.
C
C THIS SUBROUTINE IS DESIGNED TO SOLVE THE SEPARABLE ELLIPTICAL
C P.D.E.PROBLEM
C      (L(X) + M(Y))U(X,Y)=F(X,Y)     (1)
C WHERE
C L(X)= -D (P1(X)D ) + R(X)D + Q1(X)
C         X       X         X
C
C AND M(Y) = -D (P2(Y)D )+Q2(Y)
C              Y       Y
C
C ON THE RECTANGLE A<=X<=B, C<=Y<=D.
C
C THE SUBROUTINE RETURNS THE COEFFICIENTS V TO A RAYLEIGH-RITZ-
C GALERKIN APPROXIMATION TO U OF THE FORM
C
C         M      N
C U(X,Y)=SUM ( SUM V(I,J) B (X) B(Y))      (2)
C        I=1   J=1         I   J
C
C WHERE THE B(X)'S ARE B-SPLINES OF ORDER K ASSOCIATED WITH THE MESH
C A=X =.=X < X   <=.....<=X <X   =...=X   =B  (3)
C    1     K  K+1          M  M+1      M+K
C
C AND THE B(Y)'S ARE B-SPLINES OF ORDER K ASSOCIATED WITH THE MESH
C C=Y =...=Y <Y   <=......<=Y <Y   =...=Y   =D     (4)
C    1      K  K+1           N  N+1      N+K
C
C THIS CODE CAN HANDLE DIRICHLET, NEUMANN, MIXED, AND PERIODIC
C BOUNDARY CONDITIONS. FOR PERIODIC BOUNDARY CONDITIONS, THE
C B SPLINES ARE PERIODIC B 'SPLINES AND THE KNOT SEQUENCE, SAY FOR
C PERIODIC IN X, WOULD BE
C A=X <X <=.....<=X <X   =B            (5)
C    1  2          M  M+1
C
C DIRICHLET BOUNDARY CONDITIONS MAY BE SPECIFIED IN 2 WAYS:
C ONE MAY EITHER SPECIFY THE COEFFICIENTS THEMSELVES ON THE
C BOUNDARY(AN OPTION USUALLY USED FOR HOMOGENEOUS CONDITIONS
C WHEN ONE WISHES THE BOUNDARY TO BE ZERO) OR A FUNCTION.
C E.G. AT X=A ONE MAY SPECIFY   V(1,J), J=1,..N OR A FUNCTION
C GAMA(Y) WHICH MAY BE EVALUATED AT ANY Y BETWEEN C AND D.
C
C MIXED BOUNDARY CONDITIONS FOR X=A ARE ASSUMED TO HAVE THE
C FORM
C  BETAA U(A,Y)-P1(A) D U(A,Y)=GAMA(Y)      (6)
C                      X
C
C WHERE BETAA IS A SCALAR AND GAMA(Y) IS A FUNCTION.
C MIXED CONDITIONS AT Y=C ARE ASSUMED TO HAVE THE FORM
C  BETAC U(X,C)-P2(C) D U(X,C)=GAMC(Y)      (7)
C                      Y
C
C SIMILAR CONDITIONS MAY BE SPECIFIED AT X=B AND Y=D.
C IF THE BETA'S ARE 0, ONE OBTAINS NEUMANN CONDITIONS
C
C INPUT PARAMETERS:
C
C K      ORDER OF THE B SPLINE APPROXIMATION
C M      IN THE X-DIRECTION, M+K KNOTS ARE SPECIFIED FOR A
C        NON-PERIODIC PROBLEM AND M+1 FOR A PERIODIC ONE; ALSO
C        THE NUMBER OF UNKNOWNS V IS M X N.
C X      THE X MESH DEFINED AS IN (3) FOR NON PERIODIC OR (5) FOR
C        PERIODIC
C N      IN THE Y-DIRECTION, N+K KNOTS ARE SPECIFIED FOR A  NONPERIODIC
C        PROBLEM AND N+1 FOR A PERIODIC PROBLEM; THE NUMBER OF UNKNOWNS
C        V IS M X N.
C Y      THE Y MESH DEFINED AS IN (4) FOR NON-PERIODIC AND LIKE (5)
C        FOR PERIODIC
C P1     A FUNCTION WRITTEN BY THE USER, DECLARED EXTERNAL IN THE
C        CALLING PROGRAM, WHICH IS DEFINED FOR A<=X<=B THAT IS USED
C        TO DEFINE THE DIFFERENTIAL EQUATION AS GIVEN IN (1).
C        IT IS USED AS P=P1(X).
C R      A FUNCTION WRITTEN BY THE USER, DECLARED EXTERNAL IN THE
C        CALLING PROGRAM, WHICH WOULD EVALUATE R(X) AS GIVEN IN (1)
C        AT A POINT X WHERE A<=X<=B. IT IS CALLED AS A RR=R(X).
C Q1     A FUNCTION WRITTEN BY THE USER, DECLARED EXTERNAL IN THE
C        CALLING RROGRAM, WHICH WOULD EVALUATE Q1(X) AS GIVEN IN (1) AT
C        A POINT X WHERE A<=X<=B. IT IS CALLED AS Q=Q1(X).
C P2     A FUNCTION WRITTEN BY THE USER, DECLARED EXTERNAL IN THE
C        CALLING PROGRAM, WHICH WOULD EVALUATE P2(Y) AS GIVEN IN (1) AT
C        A POINT Y WHERE C<=Y<=D. IT IS CALLED AS P=P2(Y).
C Q2     A FUNCTION WRITTEN BY THE USER, DECLARED EXTERNAL IN THE
C        CALLLING PROGRAM, WHICH WOULD EVALUATE Q2(Y) AS GIVEN IN (1) AT
C        A POINT Y WHERE C<=Y<=D. IT IS CALLED AS Q=Q2(Y).
C F      A SUBROUTINE WRITTEN BY THE USER, DECLARED EXTERNAL IN THE
C        CALLING PROGRAM WHICH EVALUATES F(X,Y) AS GIVEN IN (1)
C        THE CALLING SEQUENCE IS F(X,Y,FXY,NY)
C        WHERE X IS A SCALAR INPUT PARAMETER AND Y IS A VECTOR INPUT
C        PARAMETER OF LENGTH NY. THE SUBROUTINE F SHOULD RETURN
C        THE VECTOR FXY OF LENGTH NY SUCH THAT
C          FXY(J)=F(X,Y(J)), J=1,2,....NY
C        ON SUBSEQUENT CALLS X WILL CHANGE BUT THE VECTOR Y WILL
C        NOT. HENCE IF PART OF THE COMPUTATION IS DEPENDENT ONLY
C        ON Y, IT NEED BE DONE ONLY ONCE AND STORED FOR SUBSEQUENT
C        CALLS. IN GENERAL NY=(K-1)*(N-K).
C IA     AN INTEGER GIVING THE TYPE OF BOUNDARY CONDITIONS AT X=A
C             VALUE         TYPE
C              1         DIRICHLET DEFINED BY THE FUNCTION GAMA(Y)
C              2         DIRICHLET WITH V(1,J), J=1,..N SPECIFIED
C              3         MIXED AS IN (6) WITH BETAA AND GAMA GIVEN
C              4         PERIODIC IN X U(A,Y)=U(B,Y)
C BETAA  A DOUBLE PRECISION SCALAR DEFINED WHEN A MIXED CONDITION IS
C        SPECIFIED AT X=A AS IN (6). IF NEUMANN IS WANTED, SET BETAA TO
C        0.
C GAMA   A FUNCTION WRITTEN BY THE USER, DECLARED EXTERNAL IN THE
C        CALLING PROGRAM WHICH IF IA =1, SPECIFIES THE FUNCTION
C        AT X=A FOR A GIVEN Y AND IF IA=3, IS USED AS IN (6)
C IB     AN INTEGER GIVING THE TYPE OF BOUNDARY CONDITION AT X=B.
C             VALUE        TYPE
C              1         DIRICHLET DEFINED BY THE FUNCTION GAMB(Y)
C              2         DIRICHLET WITH  V(M,J), J=1,...N SPECIFIED
C              3         MIXED WITH BETAB AND GAMB SPECIFIED. SEE (6)
C              4         PERIODIC IN X U(A,Y)=U(B,Y)
C BETAB  A DOUBLE PRECISION SCALAR DEFINED WHEN A MIXED CONDITION IS
C        SPECIFIED AT X=B. IF NEUMANN, SET BETAB=0.0
C GAMB   A FUNCTION WRITTEN BY THE USER HAVING THE SAME MEANING
C        AS GAMA ABOVE BUT AT X=B.
C IC     AN INTEGER GIVING THE TYPE OF BOUNDARY CONDITIONS AT Y=C
C             VALUE       TYPE
C              1         DIRICHLET DEFINED BY THE FUNCTION GAMC(X)
C              2         DIRICHLET WITH V(I,1), I=1,2,...M SPECIFIED
C              3         MIXED WITH BETAC AND GAMC SPECIFIED
C              4         PERIODIC IN Y U(X,C)=X(X,D)
C BETAC  A DOUBLE PRECISIONL SCALAR DEFINED WHEN A MIXED CONDITION IS
C        GIVEN AT Y=C
C GAMC   A FUNCTION WRITTEN BY THE USER, DECLARED EXTERNAL IN THE
C        CALLING PROGRAM WHICH IF IC=1, SPECIFIES THE FUNCTION AT
C        Y=C FOR A GIVEN X AND IF IC=3, IS USED AS IN (6)
C ID     AN INTEGER GIVING THE TYPE OF BOUNDARY CONDITION AT Y=D,
C        DEFINED LIKE IC ABOVE
C BETAD  A DOUBLE PRECISION SCALAR DEFINED WHEN  ID=3
C GAMD   A FUNCTION WRITTEN BY THE USER FOR EITHER THE DIRICHLET
C        OR MIXED CONDITION ACCORDING TO THE SPECIFICATION OF ID
C
C OUTPUT PARAMETER
C V      THE M TIMES N MATRIX OF B-SPLINE COEFFICIENTS. (SEE (2))
C        V SHOULD BE DECLARED V(M,N) IN THE CALLING PROGRAM
C THIS PROGRAM WAS WRITTEN BY LINDA KAUFMAN AT BELL LABS WITH
C ASSISTANCE FROM DAN WARNER AT CLEMSON UNIVERSITY
      INTEGER M, N
      EXTERNAL P1, R, Q1, P2, Q2, F
      DOUBLE PRECISION GAMA, GAMB, GAMC, GAMD
      EXTERNAL GAMA, GAMB, GAMC, GAMD
      INTEGER K, IA, IB, IC, ID
      DOUBLE PRECISION X(1), Y(1), P1, R, Q1, P2
      INTEGER TEMP
      DOUBLE PRECISION Q2, V(M, N), BETAA, BETAB, BETAC, BETAD
      LOGICAL B2DMSH
      DOUBLE PRECISION DS(500)
      COMMON /CSTAK/ DS
      INTEGER IS(1000)
      EQUIVALENCE (IS(1),DS(1))
C
C CALLS FROM PORT: SETERR, ISTKGT, ENTER, LEAVE
C CALLS FROM BLAS: DCOPY
C ALSO CALLS N2MESH,B2SERG
C
      CALL ENTER(0)
C
C CHECK THE INPUT PARAMETERS
C
        IF(K.LT.2.OR.K.GT.20)
     1   CALL SETERR(17HSERRG-INCORRECT K,17,1,2)
        IF (N.LT.1)
     1   CALL SETERR(17HSERRG-INCORRECT N,17,2,2)
        IF (M.LT.1)
     1   CALL SETERR(17HSERRG-INCORRECT M,17,3,2)
        IF (IA.LT.1.OR.IA.GT.4)
     1   CALL SETERR(18HSERRG-INCORRECT IA,18,4,2)
        IF (IA.NE.4.AND.(IB.LT.1.OR.IB.GT.3))
     1   CALL SETERR(18HSERRG-INCORRECT IB,18,5,2)
        IF (IC.LT.1.OR.IC.GT.4)
     1   CALL SETERR(18HSERRG-INCORRECT IC,18,6,2)
        IF (IC.NE.4.AND.(ID.LT.1.OR.ID.GT.3))
     1   CALL SETERR(18HSERRG-INCORRECT ID,18,7,2)
        IF (B2DMSH(M,K,X,IA))
     1   CALL SETERR(17HSERRG-INCORRECT X,17,8,2)
        IF (B2DMSH(N,K,Y,IC))
     1   CALL SETERR(17HSERRG-INCORRECT Y,17,9,2)
         IF (IA.EQ.4.AND.M.LT.K-1)
     1   CALL SETERR(17HSERRG-M TOO SMALL,17,11,1)
         IF(IC.EQ.4.AND.N.LT.K-1)
     1    CALL SETERR(17HSERRG-N TOO SMALL,17,12,1)
C
C DETERMINE THE AMOUNT OF SPACE NEEDED FOR THE MASS AND STIFFNESS
C MATRICES. FOR PERIODIC BOUNDARY CONDITIONS THE BANDWIDTH IS
C ESSENTIALLY DOUBLED BECAUSE VARIABLES ARE REORDERED ACCORDING TO
C 1,N,2,N-1,3,N-2,....N/2 SO THAT ONE STILL HAS BAND MATRICES
C
      K1=K-1
      KSY = K
      IF (IC .EQ. 4) KSY = MIN0(N,2*K-1)
      KSX = K
      IF (IA .EQ. 4) KSX = MIN0(N,2*K-1)
      K3X = 3*(KSX-1)+1
C
C GET SPACE FOR MASS AND STIFFNESS ARRAYS
C
      IMX = ISTKGT(K3X*M, 4)
      ISX= ISTKGT(K3X*M, 4)
      ISY = ISTKGT(KSY*N, 4)
      IMY = ISTKGT(KSY*N, 4)
C
      NINTY = N-K1
      IF (IC .EQ. 4) NINTY = N
      NINTX = M-K1
      IF (IA .EQ. 4) NINTX = M
      NYK1 = NINTY+K1
      MXK1 = NINTX+K1
C
C THIS SECTION OF CODE WAS DONE SO THAT THE SAME CODE CAN HANDLE
C REGULAR B-SPLINES AND THE EXPANDED MESH FOR A
C PERIODIC BSPLINE
C IPT. IS JUST A POINTER ARRAY WHICH FOR NORMAL
C B- SPLINES CORRESPONDS TO NO SHUFFLING, BUT
C FOR PERIODIC WILL TELL HOW TO PERMUTE THE ROWS
C AND COLUMNS OF THE MATRICES SO THAT BAND STRUCTURE
C IS PRESERVED
      IXX = ISTKGT(MXK1+K, 4)
      IYY = ISTKGT(NYK1+K, 4)
      IPTY = ISTKGT(NYK1, 2)
      IPTX = ISTKGT(MXK1, 2)
      IPTY1 = IPTY-1
      IPTX1 = IPTX-1
      DO  1 I = 1, NYK1
         TEMP = IPTY1+I
         IS(TEMP) = I
   1     CONTINUE
      DO  2 I = 1, MXK1
         TEMP = IPTX1+I
         IS(TEMP) = I
   2     CONTINUE
      CALL DCOPY(MXK1+K, X, 1, DS(IXX), 1)
      CALL DCOPY(NYK1+K, Y, 1, DS(IYY), 1)
C
C CHANGE THE MESH AND SET UP THE POINTER ARRAYS FOR PERIODIC
C MESHES
C
      IF(IC.EQ.4)  CALL N2MESH(N, K, Y, DS(IYY), IS(IPTY),.TRUE.)
      IF (IA.EQ.4) CALL N2MESH(M, K, X, DS(IXX), IS(IPTX),.TRUE.)
      IBB=IB
      IF (IA.EQ.4)IBB=4
      IDD=ID
      IF (IC.EQ.4)IDD=4
      CALL B2SERG(K, M, DS(IXX), N, DS(IYY), P1, R, Q1, P2, Q2, F, V,
     1   IA, BETAA, GAMA,IBB, BETAB, GAMB, IC, BETAC, GAMC,IDD, BETAD,
     1   GAMD,IS(IPTX),IS(IPTY),DS(ISX),DS(IMX),K3X,KSX,
     3   DS(ISY),DS(IMY),KSY,NINTY,NINTX, ISING,IERX, IERY)
       IF (ISING.NE.0)
     1CALL SETERR(22H SERRG-SINGULAR MATRIX,22,12,1)
       IF (IERX.NE.0)
     1CALL SETERR(22H SERRG-P1 NOT POSITIVE,22,13,1)
       IF (IERY.NE.0)
     1CALL SETERR(22H SERRG-P2 NOT POSITIVE,22,14,1)
       CALL LEAVE
      RETURN
      END
      SUBROUTINE B2SERG(K, M, X, N, Y, P1, R, Q1, P2, Q2, F, V,
     1   IA, BETAA, GAMA, IB, BETAB, GAMB, IC, BETAC, GAMC, ID, BETAD,
     1   GAMD,IPTX,IPTY,SX,MX,K3X,KSX,SY,MY,KSY,NINTY,NINTX, ISING,
     2   IERX,IERY)
C
C THIS IS THE REAL WORK HORSE OF THE PACKAGE THAT CALLS ALL THE
C UNDERLYING SUBROUTINES.
C
C    INPUT PARAMETERS
C     K THROUGH F AND IA THROUGH GAMD - HAVE THE SAME MEANING AS THEY
C          HAD ON ENTRY TO SERRG2
C     IPTX  ORDER OF THE X VARIABLES. FOR NON-PERIOD PROBLEMS IPTX(I)=I
C           FOR PERIODIC PROBLEMS, IPTX REFLECTS THE FACT THAT THE
C           UNKNOWNS HAVE TO BE REORDERED AS 1,N,2,N-1,....N/2
C     IPTY  ORDER OF THE Y VARIABLES. FOR NON-PERIOD PROBLEMS IPTY(I)=I
C           FOR PERIODIC PROBLEMS, IPTY REFLECTS THE FACT THAT THE
C           UNKNOWNS HAVE TO BE REORDERED AS 1,N,2,N-1,....N/2
C     K3X   LEADING DIMENSION OF THE SCRATCH SX AND MX ARRAYS
C     KSY   LEADING DIMENSION OF THE SCRATCH SY AND MY ARRAYS
C     NINTX NUMBER OF POINTS IN X DIRECTIONS WITHOUT MULTIPLICITY
C           OF END POINTS
C
C     OUTPUT PARAMETERS
C
C     V       THE M TIMES N MATRIX OF B-SPLINE COEFFICIENTS
C     ISING   IF NONZERO, SINGULAR PROBLEM
C     IERX    IF NONZERO, P1 IS NOT POSITIVE
C     IERY    IF NONZERO, P2 IS NOT POSITIVE
C
C     SCRATCH ARRAYS
C     SX      WILL HOLD STIFFNESS MATRIX IN X
C     MX      WILL HOLD MASS MATRIX IN X DIRECTION
C     SY      WILL HOLD STIFFNESS MATRIX IN Y DIRECTION
C     MY      WILL HOLD MASS MATRIX IN Y DIRECTION
C
C  THIS CODE IS DIVIDED INTO 3 MAIN SECTIONS
C
C  SECTION 1: OBTAIN THE MATRICES DEFINED BY THE PDE
C      CALL R2G1Y TO OBTAIN THE MASS AND STIFFNESS MATRICES IN THE
C          Y DIRECTION. EACH ELEMENT OF THESE MATRICES IS AN INTEGRAL.
C      CALL R2G1X TO OBTAIN THE MASS AND STIFFNESS MATRICES IN THE
C          X DIRECTION AND THE RIGHT HAND SIDE
C      THE LINEAR SYSTEM TO BE SOLVED AS THE FORM
C          H= SX TENSOR MY + SY TENSOR MX = F
C
C  SECTION 2: MODIFY H AND F , I.E. SX, MY, SY, MX AND F TO
C    HANDLE BOUNDARY CONDITIONS
C    FOR EACH SIDE
C       IF THERE ARE DIRICHLET BOUNDARY CONDITIONS DEFINED
C           BY A FUNCTION CALL  D2IRFN TO CHANGE THESE TO B-SPLINE
C           COEFFICIENT SPECIFICATIONS
C       IF THER ARE DIRICHLET BOUNDARY CONDTIONS, CALL D2IRX
C       IF THERE ARE MIXED BOUNDARY CONDITIONS CALL
C          B2SYGA IF THE BOUNDARY IS EITHER AT X=A OR X=B
C          OR CALL B2SXGA IF IT IS A Y BOUNDARY
C
C  SECTION 3: SOLVE THE TENSOR PRODUCT SYSTEM
C       CALL S2EDM TO SOLVE THE BANDED GENERALIZED EIGENVALUE
C          PROBLEM SY Z = MY  Z D
C       CALL  S2ESLV TO USE Z AND D TO SOLVE THE TENSOR SYSTEM
C         BY PRE AND POST TMULTIPLYING BY Z AND SOLVING A BUNCH
C         OF BANDED LINEAR SYSTEMS WITH THE SX AND MX MATRICES
C
C  THIS CODE IS COMPLICATED BECAUSE OF THE DECISION TO PERMIT
C  PERIODIC BOUNDARY CONDITIONS WHICH IS TRANSPARENT TO THE SUBROUTINES
C  R2G1Y, R2G1X, S2EDM AND S2ESLV.IT HAS ALSO BEEN MADE MORE COMPLICATED
C  BECAUSE IT WAS DECIDED TO PERMIT EXTRA GENERALITY IN THE X DIRECTION
C  WHICH LEADS TO AN UNSYMMETRIC STIFFNESS MATRIX AND TO TAKE ADVANTAGE
C  OF THE SYMMETRIC NATURE OF THE STIFFNESS MATRIX IN Y.
C
      INTEGER M, N
      EXTERNAL P1, R, Q1, P2, Q2, F
      EXTERNAL GAMA, GAMB, GAMC, GAMD
      DOUBLE PRECISION GAMA, GAMB, GAMC, GAMD
      INTEGER K, IA, IB, IC, ID
      INTEGER IERX, IERY
      DOUBLE PRECISION X(1), Y(1), P1, R, Q1, P2
      DOUBLE PRECISION Q2, V(M, N), BETAA, BETAB, BETAC, BETAD
      INTEGER KSX,K3X,KSY,IPTX(M),IPTY(N)
      DOUBLE PRECISION SX(K3X,M),SY(KSY,N),MX(K3X,M),MY(KSY,M)
      DOUBLE PRECISION DS(500)
      COMMON /CSTAK/ DS
      INTEGER ICA, ICB, ICC, ICD, IBS, IVB
      INTEGER JVB,  IBY
      INTEGER IZY, IWZ
      INTEGER IZX,  ISTKGT,  K3XB
      INTEGER  MXK1, NYK1
      INTEGER I, J, IMYB, IMXB, IFTR, ISXB
      INTEGER ISYB, IFZY,   I1, K1
      INTEGER   MM, IS(1000), NN
      INTEGER IV, IW, IZ, NINTY, NINTX
      DOUBLE PRECISION  CORN
      INTEGER TEMP,  TEMP2
      LOGICAL RESTAR
      DOUBLE PRECISION ABSI(20,2),WTS(20,2)
      COMMON /GLQNW/ ABSI,WTS
      EQUIVALENCE (DS(1), IS(1))
      DATA IORDER /0/
C
C CALLS FROM PORT: SETERR, ISTKGT, ENTER, LEAVE, ISTKRL
C CALLS FROM BLAS: DCOPY
C ALSO CALLS G2QM11, R2G1Y, R2G1X, P2SHM, G2ETBZ, P2SHF, D2IRFN,
C      D2IRX, B2SYGA, B2SXGA, S2EDM, S2ESLV
C
      CALL ENTER(0)
      K1 = K-1
      K2X=2*(KSX-1)+1
C
C  STORE GIVEN COEFFICIENTS DIRICHLET CONDITIONS IN SAFE PLACE
C
      IF (IA .LE. 2) ICA = ISTKGT(N, 4)
      IF (IB .EQ. 1 .OR. IB .EQ. 2) ICB = ISTKGT(N, 4)
      IF (IC .EQ. 1 .OR. IC .EQ. 2) ICC = ISTKGT(M, 4)
      IF (ID .EQ. 1 .OR. ID .EQ. 2) ICD = ISTKGT(M, 4)
      IF (IA .EQ. 2) CALL DCOPY(N, V, M, DS(ICA), 1)
      IF (IB .EQ. 2) CALL DCOPY(N, V(M, 1), M, DS(ICB), 1)
      IF (IC .EQ. 2) CALL DCOPY(M, V, 1, DS(ICC), 1)
      IF (ID .EQ. 2) CALL DCOPY(M, V(1, N), 1, DS(ICD), 1)
C
C
C SECTION 1 : CONSTRUCTING THE MATRICES THAT ARE DEFINED BY THE PDE
C
      NQY=K1
      IF ((IC.EQ.3.AND.ID.EQ.3.OR.IC.EQ.4).AND.K.EQ.2)NQY=K
C
C GET GAUSS QUADRATURE POINTS IF NECESSARY
C
      NCOLY=NQY-K1+1
      IF (IORDER.EQ.K1) GO TO 1
        CALL G2QM11(K1,ABSI,WTS)
        CALL G2QM11(K,ABSI(1,2),WTS(1,2))
        IORDER=K1
  1   CONTINUE
      IBS = ISTKGT(K*K, 4)
      IW = ISTKGT(K*K, 4)
      IZY = ISTKGT(NQY*NINTY, 4)
      IWZ = ISTKGT(NQY, 4)
      IBY = ISTKGT(NQY*NINTY*K, 4)
      NYK1=NINTY+K1
      MXK1=NINTX+K1
C
C GET MASS AND STIFFNESS MATRICES IN Y
C
      IF (IC .GE. 4) GOTO 2
         CALL R2G1Y(Y, N, K, P2, Q2, SY, MY, NQY, NINTY, DS(
     1      IZY), DS(IWZ), DS(IBY), DS(IBS), DS(IW),NCOLY,IERY)
         GOTO  3
C
C WE ARE WORKING WITH A PERIODIC B-SPLINE IN Y
C SOLVE PROBLEM WITH BIGGER MESH
C
   2     CONTINUE
         ISYB = ISTKGT(K*NYK1, 4)
         IMYB = ISTKGT(K*NYK1, 4)
         CALL R2G1Y(Y, NYK1, K, P2, Q2, DS(ISYB), DS(IMYB), NQY, N
     1      , DS(IZY), DS(IWZ), DS(IBY), DS(IBS), DS(IW),NCOLY, IERY)
C
C SHUFFLE INTO MATRICES OF LARGER BANDWIDTH
C BECAUSE FOR PERIODIC WE REORDER THE VARIABLES TO BE 1,N,2,N-1,ETC.
C WHICH DOUBLES THE BANDWIDTH
C
         CALL P2SHM(N, K, DS(ISYB), K, SY, KSY, IPTY, 1, K, K
     1      , KSY, .FALSE.)
         CALL P2SHM(N, K, DS(IMYB), K, MY, KSY, IPTY, 1, K, K
     1      , KSY, .FALSE.)
         CALL ISTKRL(2)
C
C
   3        CONTINUE
         IF (IERY.NE.0) GO TO 501
          IF (NQY.EQ.K1)GO TO 4
C
C FIX UP B'S ETC WHEN WORKING WITH NO. QUADRATURE POINTS N.E.K-1
C
            CALL ISTKRL(3)
            IZY=ISTKGT(K1*NINTY, 4)
            IWZ=ISTKGT(K1, 4)
            IBY=ISTKGT(K1*NINTY*K, 4)
            CALL G2ETBZ(Y,NYK1,K,K1,NINTY,DS(IZY),DS(IWZ),DS(IBY),
     1      DS(IBS),DS(IW))
  4        CONTINUE
      IZX = ISTKGT(K1, 4)
      IFZY = ISTKGT(K1*NINTY, 4)
      IFTR = ISTKGT(NYK1, 4)
      IF (IA .GE. 4) GOTO 7
         IF (IC .GE. 4) GOTO 5
C
C GET MASS AND STIFFNESS MATRICES IN X AND DETERMINE
C RIGHT HAND SIDE
C
            CALL R2G1X(X,M,K,P1,R,Q1,F,K1, NINTY,DS(IZY),DS(
     1         IBY), K3X, SX, MX, N, V, DS(IZX), DS(IWZ), DS(
     1         IBS), DS(IW), DS(IFZY), DS(IFTR), IERX)
            GOTO  6
   5        IV = ISTKGT(M*NYK1, 4)
C
C WE ARE WORKING WITH MESH THAT IS NOT PERIODIC IN X BUT PERIODIC IN Y
C
            CALL R2G1X(X,M,K,P1,R,Q1,F,K1,NINTY,DS(IZY),DS(
     1         IBY), K3X, SX, MX, NYK1, DS(IV), DS(IZX), DS(
     1         IWZ), DS(IBS), DS(IW), DS(IFZY), DS(IFTR), IERX)
C
C RESHUFFLE RIGHT HAND SIDE TO REFLECT ORDERING OF UNKNOWNS AS
C 1,N,2,N-1.... FOR PERIODIC PROBLEM
C
            CALL P2SHF(DS(IV), M, V, M, N, IPTX, IPTY, NYK1, M)
   6     CONTINUE
         GOTO  8
   7     CONTINUE
C WE ARE WORKING WITH A PERIODIC B-SPLINE IN X
C SET UP PROBLEM FOR NEW MESH  AND GET MATRICES AND RIGHT HAND SIDE
         K3XB = 3*K1+1
         IV = ISTKGT(MXK1*NYK1, 4)
         ISXB = ISTKGT(K3XB*MXK1, 4)
         IMXB = ISTKGT(K3XB*MXK1, 4)
         CALL R2G1X(X,MXK1,K,P1,R,Q1,F,K1,NINTY,DS(IZY),
     1      DS(IBY), K3XB, DS(ISXB), DS(IMXB), NYK1, DS(IV), DS(IZX),
     1      DS(IWZ), DS(IBS), DS(IW), DS(IFZY), DS(IFTR),IERX)
C
C RESHUFFLE MATRICES AND RIGHT HAND SIDE TO GET BIGGER BAND WIDTH
C FOR PERIODICITY IN X
C
         CALL P2SHM(M, K, DS(IMXB), K3XB, MX, K3X, IPTX, K,
     1      KSX, KSX, 2*(KSX-1)+1, .FALSE.)
         CALL P2SHM(M, K, DS(ISXB), K3XB, SX, K3X, IPTX, K, 3*(
     1      K-1)+1, KSX, 2*(KSX-1)+1, .TRUE.)
         CALL ISTKRL(2)
         CALL P2SHF(DS(IV), MXK1, V, M, N, IPTX, IPTY, NYK1,
     1      MXK1)
   8  CONTINUE
      IF(IERX.NE.0) GO TO 501
C
      DO  10 I = 2, KSX
         I1 = I-1
         K2XPI1=K2X+I1
         IM1=I-1
         K2XMI1=K2X-I1
         DO  9 J = I, M
             JMIM1=J-IM1
            MX(K2XPI1,JMIM1)=MX(K2XMI1,J)
   9        CONTINUE
  10     CONTINUE
C
C
C SECTION 2:    BOUNDARY CONDITIONS
C
C MM  AND NN ARE SIZES OF THE 1 D PROBLEM TO BE SOLVED
C WITH DIRICHLET BOUNDARY CONDITIONS, THE SIZE
C IS REDUCED
C IVB AND JVB INDICATES WHERE TO BEGIN WITH THE RIGHT HAND SIDE
C AND MATRICES. AGAIN DIRICHLET CONDITIONS CAN CHANGE THEM
C
      MM = M
      NN = N
      IVB=1
      JVB=1
C
C HANDLE DIRICHLET CONDITIONS SPECIFIED BY A FUNCTION
C GET ALL THE DIRICHLET FUNCTIONS FIRST SO THAT ONE CAN MODIFY
C THE CORNERS LATER
C
      IF (IA.GT.1.AND.IB.GT.1) GO TO 80
      IWS=ISTKGT(K*NINTY*(K+1)+K+N*KSY, 4)
      IF(IA.EQ.1)
     1CALL D2IRFN(K,KSY,N,NINTY,Y,IPTY,GAMA,DS(ICA),IC,
     1.TRUE.,DS(IWS))
       RESTAR=IA.NE.IB
       IF (IB.EQ.1)
     1 CALL D2IRFN(K,KSY,N,NINTY,Y,IPTY,GAMB,DS(ICB),IC,RESTAR,DS(IWS))
       CALL ISTKRL(1)
80     CONTINUE
       IF(IC.GT.1.AND.ID.GT.1) GO TO 90
       IWS=ISTKGT(K*NINTX*(K+1)+K+M*KSX, 4)
       IF(IC.EQ.1)
     1 CALL D2IRFN(K,KSX,M,NINTX,X,IPTX,GAMC,DS(ICC),IA,.TRUE.,DS(IWS))
       RESTAR=IC.NE.ID
       IF (ID.EQ.1)
     1 CALL D2IRFN(K,KSX,M,NINTX,X,IPTX,GAMD,DS(ICD),IA,RESTAR,DS(IWS))
       CALL ISTKRL(1)
90     CONTINUE
C
C BOUNDARY CONDITIONS AT X=A
C
       IF (IA-3)100,103,205
C
C HANDLE DIRICHLET CONDITIONS- FIRST MAKE SURE
C CORNERS MATCH
C
100    IF(IC.GT.2)GO TO 101
       CORN=(DS(ICA)+DS(ICC))/2.D0
       DS(ICA)=CORN
       DS(ICC)=0.0D0
101    IF(ID.GT.2) GO TO 102
         ICAPN=ICA+N
         CORN=(DS(ICAPN-1)+DS(ICD))/2.D0
         DS(ICAPN-1)=CORN
         DS(ICD)=0.0D0
102   CONTINUE
C
C FIX UP THE RIGHT HAND SIDE TO REFLECT THE DIRICHLET CONDITIONS
C
         IPC=1
         IF (IA.EQ.2)IPC=IC
         CALL D2IRX(M, N, SY, MY, KSY, K,IPC, IPTY, K, DS(
     1      ICA), V, 1, SX(K2X,1), MX(K2X,1),1, KSX, MM, 1)
         IVB=2
         GO TO 105
C
C HANDLE MIXED CONDITION AT X=A
C
103    CALL B2SYGA(Y, NYK1, K, GAMA, V, M, DS(IZY)
     1   , DS(IBY), K1, NINTY,  IC, IPTY,.TRUE.)
C
C BOUNDARY CONDITIONS AT X=B
C
105      IF (IB.EQ.3) GO TO 203
C
C HANDLE DIRICHLET CONDITIONS
C
         IF (IC.GT.2) GO TO 201
          ICCPM=ICC+M-1
          CORN=(DS(ICB)+DS(ICCPM))/2.0D0
          DS(ICB)=CORN
          DS(ICCPM)=0.0D0
201      IF(ID.GT.2) GO TO 202
          ICBPN=ICB+N
          ICDPM=ICD+M
          CORN=(DS(ICBPN-1)+DS(ICDPM-1))/2.0D0
           DS(ICBPN-1)=CORN
           DS(ICDPM-1)=0.0D0
202     IPC=1
        TEMP = M-KSX
         IF (IB.EQ.2)IPC=IC
         CALL D2IRX(M, N, SY, MY, KSY, K,IPC, IPTY, K, DS(
     1      ICB), V(TEMP, 1), 1, SX(KSX-1,M), MX(KSX-1,M), 1, KSX, MM,
     1      1)
         GO TO 205
C
C MIXED CONDITION AT X=B
C
203   CALL B2SYGA(Y, NYK1, K, GAMB, V(M, 1), M, DS(IZY),
     1    DS(IBY), K1, NINTY,
     1    IC, IPTY,.FALSE.)
C
C
C BOUNDARY CONDITIONS AT Y=C
C
205     CONTINUE
        IF(IC-3)301,303,405
C
C DIRICHLET CONDITONS
C
301    CONTINUE
         IPA=1
         IF(IC.EQ.2)IPA=IA
         CALL D2IRX(1, M, SX, MX, K3X, KSX,IPA, IPTX, K,
     1      DS(ICC), V, M, SY(KSY,1), MY(KSY,1), KSY-1, KSY, NN, 0)
         JVB=2
         GO TO 305
C
C MIXED CONDITION
C
303   CALL B2SXGA(X, MXK1, K, GAMC, V, 1, DS(IZX)
     1   , DS(IW), DS(IBS), K1,  IA, IPTX,.TRUE.)
C
C BOUNDARY CONDITIONS AT Y=D
C
305     IF(ID.EQ.3) GO TO 403
C
C DIRICHLET CONDITION
C
         IPA=1
         TEMP2 = N-KSY
          IF (ID.EQ.2)IPA=IA
         CALL D2IRX(1, M, SX, MX, K3X, KSX,IPA, IPTX, K,
     1      DS(ICD),V(1,TEMP2),M,SY(KSY,N-1),MY(KSY,N-1), 1,KSY,NN
     1      , 0)
         GO TO 405
C
C MIXED CONDITION AT Y=D
C
403   CALL B2SXGA(X, MXK1, K, GAMD, V(1, N), 1, DS(
     1   IZX), DS(IW), DS(IBS), K1,  IA, IPTX,.FALSE.)
405      CONTINUE
         IF(IA.EQ.3)SX(K2X,1)=SX(K2X,1)+BETAA
         IF (IB.EQ.3)SX(K2X,M)=SX(K2X,M)-BETAB
         IF(IC.EQ.3)SY(KSY,1)=SY(KSY,1)+BETAC
         IF(ID.EQ.3)SY(KSY,N)=SY(KSY,N)-BETAD
C
C SECTION 3:    FACTOR AND SOLVE THE TENSOR PRODUCT SYSTEM
C
      IDD = ISTKGT(N, 4)
      IZ = ISTKGT(N*N, 4)
      CALL S2EDM(NN, KSY, SY(1,JVB),MY(1,JVB),DS(IDD), DS(IZ))
      CALL S2ESLV(MM, NN, DS(IZ), V(IVB, JVB), M, SX(1,IVB),MX(1,IVB),
     1   K3X, DS(IDD), KSX,ISING)
         IF(IC.LT.3)CALL DCOPY(M,DS(ICC),1,V,1)
          IF(ID.LT.3)CALL DCOPY(M,DS(ICD),1,V(1,N),1)
         IF (IA.LT.3)CALL DCOPY(N,DS(ICA),1,V,M)
         IF (IB.LT.3)CALL DCOPY(N,DS(ICB),1,V(M,1),M)
C
C FOR PERIODIC BOUNDARY CONDITION THE COEFFICIENTS HAVE BEEN
C DETERMINED IN A PERMUTED ORDER. THEY MUST NOW BE PERMUTED BACK INTO
C NATURAL ORDER
C
      IF (IC.EQ.4.OR.IA.EQ.4)
     1 CALL S2HFB(K,M,N,V,DS(IV),IPTX,IPTY,IC,IA)
C
501   CALL LEAVE
      RETURN
C LAST CARD OF B2SERG
      END
      SUBROUTINE N2MESH(N, K, OLDX, NEWX, IPT, POINT)
C
C THIS SUBROUTINE SETS UP BIGGER MESH FOR PERIODIC PROBLEMS
C AND SETS UP SHUFFLE POINTER ARRAY
C
C      INPUT PARAMETERS
C
C     N    N+1 KNOTS ARE SPECIFIED
C     K    ORDER OF APPROXIMATION
C     OLDX ORIGINAL MESH
C     POINT LOGICAL INDICATING WHETHER POINTER ARRAY SHOULD BE FIXED
C           ALSO
C
C      OUTPUT PARAMETERS
C
C      NEWX EXPANDED MESH WHICH ONE EACH END HAS KM1 NEW POINTS
C      IPT  IF POINT IS TRUE, PERMUTATION ARRAY GIVING NEW ORDER
C           OF UNKNOWNS FOR PERIODIC PROBLEM I.E. 1,N,2,N-1...
C
      INTEGER N
      LOGICAL POINT
      INTEGER K, IPT(N)
      DOUBLE PRECISION OLDX(N), NEWX(N)
      INTEGER I, I1, I2, IPKM1, NPKPI, KM1
      INTEGER NP1, NMK1PI
      INTEGER TEMP
      KM1 = K-1
      NP1 = N+1
C  SHIFT OVER MESH TO MAKE ROOM FOR NEW BEGINNING POINTS
      DO  1 I = 1, NP1
         IPKM1 = I+KM1
         NEWX(IPKM1) = OLDX(I)
   1     CONTINUE
C
C EXPAND MESH AT BOTH ENDPOINTS SO IT APPEARS PERIODIC
C
      DO  2 I = 1, KM1
         NMK1PI = N-KM1+I
         NEWX(I) = OLDX(1)-(OLDX(NP1)-OLDX(NMK1PI))
         NPKPI = N+K+I
         NEWX(NPKPI) = OLDX(N+1)+OLDX(I+1)-OLDX(1)
   2     CONTINUE
       IF (.NOT.POINT)RETURN
C
C FIX UP PERMUTATION ARRAY FOR PERIODIC MESH
C
      I1 = K
      I2 = N+KM1
      DO  3 I = 1, N, 2
         IPT(I1) = I
         IPT(I2) = I+1
         I1 = I1+1
         I2 = I2-1
   3     CONTINUE
      IF (I1 .EQ. I2) IPT(I1) = N
      DO  4 I = 1, KM1
         TEMP = N+I
         IPT(I) = IPT(TEMP)
   4     CONTINUE
      RETURN
      END
      SUBROUTINE P2SHM(N, K, MXB, IMXB, MX, IMX, IPT, IBEG, JEND
     1   , KDO, KDN, U)
      INTEGER IMX, N, IMXB
      INTEGER K, IPT(N), IBEG, JEND, KDO, KDN
      LOGICAL U
      DOUBLE PRECISION MXB(IMXB, N), MX(IMX, N)
      INTEGER MIN0, MAX0, I, J, ICOL, J2
      INTEGER TEMP, TEMP1, TEMP2, TEMP3, TEMP4
C
C THIS SUBROUTINE SHUFFLES THE 1-D MATRICES
C WHEN THERE ARE PERIODIC BOUNDARY CONDITIONS
C
C    INPUT PARAMETERS
C    N       NUMBER OF COLUMNS IN THE MATRIX IS N+K-1
C    K       ORDER OF THE APPROXIMATION
C    MXB     ORIGINAL BAND MATRIX
C    IMX     LEADING DIMENSION OF MX
C    IMXB    LEADING DIMENSION OF MXB
C    IPT     INTEGER ARRAY GIVING PERMUTATION DEFINED BY N2MESH
C    IBEG    FIRST NONZERO ROW IN MXB BY LINPACK CONVENTION
C    JEND    LAST NONZERO ROW IN MXB BY LINPACK CONVENTION
C    KD0     MAIN DIAGONAL OF OLD MATRIX
C    KDN     MAIN DIAGONAL OF NEW MATRIX
C    U       LOGICAL WHICH IF .TRUE. MEANS THE MATRIX IS UNSYMMETRIC
C
C    OUTPUT PARAMETERS
C
C    MX      PERMUTED MATRIX ACCORDING TO IPT
C
C IF U IS .TRUE. WE ARE WORKING WITH AN UNSYMMETRIC MATRIX
      CALL S2ETD(IMX*N, 0.0D0, MX)
      TEMP = N+K-1
      DO  4 I = 1, TEMP
         ICOL = IPT(I)
         TEMP2 = MAX0(KDO+1-I, IBEG)
         TEMP1 = MIN0(JEND, KDO+(N-I+K))
         DO  3 J = TEMP2, TEMP1
            TEMP3 = I+J-KDO
            J2 = IPT(TEMP3)
            IF (J2 .GT. ICOL .AND. (.NOT. U)) GOTO 1
               TEMP3 = KDN+J2-ICOL
               MX(TEMP3, ICOL) = MX(TEMP3, ICOL)+MXB(J, I)
               GOTO  2
   1           TEMP4 = KDN+ICOL-J2
               MX(TEMP4, J2) = MX(TEMP4, J2)+MXB(J, I)
   2        CONTINUE
   3        CONTINUE
   4     CONTINUE
      RETURN
      END
      SUBROUTINE P2SHF(OLDF, IOF, NEWF, M, N, IPTX, IPTY, NCOLO,
     1   NROWO)
      INTEGER IOF, M, N
      INTEGER IPTX(M), IPTY(N), NCOLO, NROWO
      DOUBLE PRECISION OLDF(IOF, N), NEWF(M, N)
      INTEGER I, J, I2, J2
C
C THIS SUBROUTINE SHUFFLES THE RIGHT HAND SIDE
C IF THERE ARE PERIODIC BOUNDARY CONDITIONS
C
C     INPUT PARAMETERS
C
C     OLDF  ORIGINAL RIGHT HAND SIDE STORED AS A 2 DIMENSIONAL ARRAY
C     IOF   LEADING DIMENSION OF OLDF
C     M     THE NEW F IS PUT IN M X N ARRAY
C     N     THE NEW F IS PUT IN AN M X N ARRAY
C     IPTX  PERMUTATION ARRAY IN X DIRECTION
C     IPTY  PERMUTATION ARRAY IN Y DIRECTION
C     NCOLO IN 2 D ARRAY OLDF GIVES NUMBER OF COLUMNS TO LOOK AT
C     NROWO IN 2 D ARRAY OLDF GIVES NUMBER OF ROWS TO LOOK AT
C
C     OUTPUT PARAMETER
C     NEWF  THE PERMUTED RIGHT HAND SIDE WHICH REFLECTS THE SYMMETRIC
C           REORDERING OF THE BAND MATRICES. NOTE THAT MORE THAN
C           ELEMENT IN THE OLDF CAN CONTRIBUTE TO AN ELEMENT IN THE
C           NEW RIGHT HAND SIDE
C
      CALL S2ETD(M*N, 0.0D0, NEWF)
      DO  2 J = 1, NCOLO
         J2 = IPTY(J)
         DO  1 I = 1, NROWO
            I2 = IPTX(I)
            NEWF(I2, J2) = NEWF(I2, J2)+OLDF(I, J)
   1        CONTINUE
   2     CONTINUE
      RETURN
      END
      SUBROUTINE D2IRX(M, N, SY, MY, KSY, KY, IC, IPTY, K, C, V,
     1   IV, SXV, MXV, ISX, KSX, MM, IX)
C
C THIS SUBROUTINE FIXES UP THE RIGHT HAND SIDE FOR
C DIRICHLET BOUNDARY B-SPLINE CONDITIONS ONCE THE COEFFICIENTS
C HAVE BEEN SPECIFIED
C THIS ENTAILS MULTIPLYING THE KNOWN SOLUTION POINTS BY THE
C COEFFICIENT MATRIX WHICH IS IN TENSOR PRODUCT FORM AND
C SUBTRACTING THEM FROM THE CURRENT RIGHT HAND SIDE
C
C      INPUT PARAMETERS
C      M     INCREMENT IN RIGHT HAND SIDE WHICH IS 1 IN Y DIRECTION
C            BUT M IN X DIRECTION
C      N     BLOCK SIZE FOR BANDED MATRIX
C      SY    STIFFNESS MATRIX WHICH WOULD BE SX IF CONDITION ALONG
C            X BOUNDARY AND OTHERWISE SY OF B2SERG
C      MY    MASS MATRIX WHICH MAY BE MX OF B2SERG IF CONDITION IS
C            ALONG X BOUNDARY
C      KSY   LEADING DIMENSION OF MASS AND STIFFNESS MATRICES
C      KY    IF MATRICES ARE IN LINPACK'S UNSYMMETRIC FORM, THE
C            FIRST NONZERO ROW
C      IC    INTEGER WHICH IF 4 MEANS THE KNOWN COEFFICIENTS MUST BE
C            FIRST PERMUTED TO REFLECT A PERIODIC BOUNDARY CONDITION
C            IN THE OTHER DIRECTION
C      IPTY  PERMUTATION VECTOR IF THERE ARE PERIODIC CONDITIONS IN
C            THE OTHER DIRECTION
C      K     THE ORDER OF THE APPROXIMATION
C      C     THE VECTOR OF KNOWN B-SPLINE COEFFICIENTS
C      V     THE ORIGINAL RIGHT HAND SIDE BEFORE IT IS MANGLED
C      IV    THE LEADING DIMENSION OF V WHICH IS STORED AS A 2 D ARRAY
C      SXV   THE OTHER STIFFNESS MATRIX THAN SY
C      MXV   THE MASS MATRIX IN THE OTHER DIRECTION
C      ISX   LEADING DIMENSION OF THE SXV AND MXV
C      KSX   NUMBER OF AFFECTED BLOCKS IN THE MASS AND STIFFNESS ARRAYS
C      MM    THE SIZE OF THE PROBLEM IN THE CURRENT DIRECTION WHICH IS
C            ABOUT TO BE DECREASED
C      IX    IF IX=1 SY AND MY ARE SYMMETRIC BANDED MATRICES, I. E.
C            REALLY MATRICES IN THE Y DIRECTION. IF IX.NE.1. SY AND MY
C            ARE UNSYMMETRIC BANDED MATRICES AND THEY ARE REALLY
C            THE MATRICES IN THE X DIRECTION
C
C      OUTPUT PARAMETERS
C        V   THE CHANGED RIGHT HAND SIDE
C        MM   THE NEW SIZE OF THE PROBLEM IN THE RELEVANT DIRECTION
C             I.E. THE SIZE HAS BEEN DECREASED BY 1 SINCE THE
C             COEFFICIENTS ARE KNOWN ALONG ONE BOUNDARY
C
C
      INTEGER ISX, KSY, M, N, IV
      INTEGER KY, IC, IPTY(N), K, KSX, MM
      INTEGER IX
      DOUBLE PRECISION SY(KSY, N), MY(KSY, M), C(N), V(IV, 1), SXV(ISX
     1   , 1), MXV(ISX, 1)
      DOUBLE PRECISION A
      DOUBLE PRECISION DS(500)
      COMMON /CSTAK/ DS
      INTEGER IWW, ISTKGT, ITM1, I, I2, IT
      INTEGER MU
      INTEGER TEMP, TEMP1
      CALL ENTER(0)
C
C CALLS FROM PORT: ENTER, ISTKGT, LEAVE
C CALLS FROM BLAS: DCOPY, DAXPY
C OTHER CALLS: G2X, G2Y
C
      MM = MM-1
      IWW = ISTKGT(2*N, 4)
      IF (IC .NE. 4) GOTO 2
C
C IF C'S HAVE NOT BEEN PERMUTED FOR PERIODIC BOUNDARY
C CONDITIONS IN THE OTHER VARIABLE, PERMUTE THEM NOW
C
         IT = ISTKGT(N, 4)
         CALL DCOPY(N, C, 1, DS(IT), 1)
         ITM1 = IT-K
         TEMP = N+K-1
         DO  1 I = K, TEMP
            I2 = IPTY(I)
            TEMP1 = ITM1+I
            C(I2) = DS(TEMP1)
   1        CONTINUE
C
C DO A MATRIX VECTOR MULTIPLICATION WITH THE KNOWN
C B-SPLINE COEFFICIENTS WITH BOTH THE STIFFNESS AND
C MASS MATRICES. NOTE WE ARE EXPLICITLY USING THE TENSOR
C PRODUCT FORM OF THE MAIN MATRIX TO ELIMINATE THE KNOWN
C BOUNDARY COMPONENTS
C
   2  IF (IX .NE. 1) GOTO 3
         CALL G2X(N, C, SY, MY, KSY, DS(IWW))
         GOTO  4
   3     CALL G2Y(N, C, SY, MY, KY, KSY, DS(IWW))
C
C MULTIPLY THE RESULTS OF THE BANDED MULTIPLICATIONS BY THE
C APPROPRIATE SCALAR TO CHANGE THE RIGHT HAND SIDE
C
   4  DO  5 MU = 2, KSX
         A=-SXV(1,MU)
         CALL DAXPY(N,A,DS(IWW), 1, V(1, MU), M)
         TEMP = IWW+N
         A=-MXV(1,MU)
         CALL DAXPY(N, A, DS(TEMP), 1, V(1,MU), M)
   5     CONTINUE
      CALL LEAVE
      RETURN
      END
      SUBROUTINE G2X(N, V, SY, MY, K, W)
      INTEGER K, N
      DOUBLE PRECISION V(N), SY(K, N), MY(K, N), W(N, 2)
      INTEGER I2P1, I, J,  I2, KM1
      INTEGER TEMP
C THIS SUBROUTINE IS USED WHEN DIRICHLET BOUNDARY CONDITIONS
C ARE NEEDED .IT COMPUTES W(1)=MY*V AND W(2)=SY*V WHERE SY AND MY ARE
C SYMMETRIC BANDED MATRICES STORED IN LINPACK'S PACKED FORM
      DO  1 J = 1, N
         W(J, 1) = MY(K, J)*V(J)
         W(J, 2) = SY(K, J)*V(J)
   1     CONTINUE
      KM1 = K-1
      DO  4 I = 1, KM1
         I2 = K-I
         I2P1 = I2+1
         IF (I2P1.GT.N) GO TO 4
         DO  3 J = I2P1, N
            TEMP = J-I2
            W(TEMP, 1) = W(TEMP, 1)+MY(I, J)*V(J)
            W(TEMP, 2) = W(TEMP, 2)+SY(I, J)*V(J)
            W(J, 1) = W(J, 1)+MY(I, J)*V(TEMP)
            W(J, 2) = W(J, 2)+SY(I, J)*V(TEMP)
   3        CONTINUE
   4     CONTINUE
      RETURN
      END
      SUBROUTINE G2Y(M, V, SX, MX, K, K3, W)
      INTEGER M, K3
      INTEGER K
      DOUBLE PRECISION V(M), SX(K3, M), MX(K3, M), W(M, 2)
      INTEGER MIN0, MAX0, I, J, JBEG, IDIF
      INTEGER JEND
      INTEGER TEMP
C THIS SUBROUTINE IS USED WHEN DIRICHLET BOUNDARY CONDITIONS
C ARE NEEDED .IT COMPUTES W(1)=MY*V AND W(2)=SY*V WHERE SY AND MY ARE
C UNSYMMETRIC BANDED MATRICES STORED IN LINPACK'S PACKED FORM
C
C CALLS S2ETD
C
      CALL S2ETD(2*M, 0.0D0, W)
      DO  2 I = K, K3
         IDIF = I-(2*K-1)
         JBEG = MAX0(1, 1-IDIF)
         JEND = MIN0(M, M-IDIF)
         IF (JBEG.GT.JEND) GO TO 2
         DO  1 J = JBEG, JEND
            TEMP = J+IDIF
            W(TEMP, 1) = W(TEMP, 1)+MX(I, J)*V(J)
            W(TEMP, 2) = W(TEMP, 2)+SX(I, J)*V(J)
   1        CONTINUE
   2     CONTINUE
      RETURN
      END
      SUBROUTINE B2SXGA(X, M, K, GAM, V, IV, ZX, WX, BX, KM1,
     1    IA, IPTX,ADD)
      INTEGER  M, IV, KM1
      INTEGER K,  IA, IPTX(M)
      DOUBLE PRECISION X(M), GAM, V(IV, 1), ZX(KM1), WX(KM1), BX(KM1)
      EXTERNAL GAM
      INTEGER J, LEFT, IC, JJ
      DOUBLE PRECISION GW, ZZ
      LOGICAL ADD
      DOUBLE PRECISION ABSI(20,2),WTS(20,2)
      COMMON /GLQNW/ ABSI,WTS
C
C THIS SUBROUTINE FIXES UP THE RIGHT HAND SIDE WHEN A MIXED CONDITION
C IS ENCOUNTERED IN X I.E. WHEN THE B-SPLINE COEFFICIENTS HAVE
C NOT BEEN SAVED
C SEE KAUFMAN AND WARNER, SINUM 1984.
C THE FIX UP INVOLVES ADDING THE INTEGRAL OF THE FUNCTION GAM
C MULTIPLIED BY THE B-SPLINES IN THE X DIRECTION
C
C CALLS BSPLVB FROM DE BOORS'S B SPLINE PACKAGE
C CALLS G2QAB
C
      DO  3 LEFT = K, M
C        DO NOT PERFORM QUADRATURE ON EMPTY INTERVALS.  EMPTY INTERVALS
C        WILL ARISE WHEN THERE ARE INTERIOR MULTIPLE KNOTS.
         IF (X(LEFT) .EQ. X(LEFT+1))  GO TO 3
C        GET GAUSS QUADRATURE POINTS FOR THIS INTERVAL
         CALL G2QAB(X(LEFT), X(LEFT+1), KM1, ZX, WX,1)
         DO  2 J = 1, KM1
            ZZ = ZX(J)
            GW = GAM(ZZ)*WX(J)
            IF(.NOT.ADD)GW=-GW
            CALL BSPLVB(X, K, 1, ZZ, LEFT, BX)
            DO  1 JJ = 1, K
               IC = LEFT-K+JJ
               IF (IA .EQ. 4) IC = IPTX(IC)
               V(1, IC) = V(1, IC)+GW*BX(JJ)
   1           CONTINUE
   2        CONTINUE
   3     CONTINUE
      RETURN
      END
      SUBROUTINE B2SYGA(Y,N, K, GAM, V, IV, ZY, BY, KM1, NINT,
     1    IC, IPTY,ADD)
      INTEGER K, N, NINT, IV, KM1
      INTEGER  IC, IPTY(N)
      DOUBLE PRECISION GAM, V(IV, 1), ZY(KM1, NINT),BY(
     1   KM1, NINT, K)
      DOUBLE PRECISION Y(N)
      EXTERNAL GAM
      INTEGER ICC, J, LEFT, JJ, IY
      DOUBLE PRECISION GW, ZZ
      LOGICAL ADD
C
C THIS SUBROUTINE FIXES UP THE RIGHT HAND SIDE WHEN A MIXED CONDITION
C IS ENCOUNTERED IN Y WHEN THE B SPLINE COEFFICIENTS HAVE
C BEEN SAVED IN BY.
C
      DO  3 LEFT = K, N
         IY = LEFT-K+1
         IF (Y(LEFT).EQ.Y(LEFT+1)) GO TO 3
         DO  2 J = 1, KM1
            ZZ = ZY(J, IY)
            GW = GAM(ZZ)
            IF(.NOT.ADD)GW=-GW
            DO  1 JJ = 1, K
               ICC = LEFT-K+JJ
               IF (IC .EQ. 4) ICC = IPTY(ICC)
               V(1, ICC) = V(1, ICC)+GW*BY(J, IY, JJ)
   1           CONTINUE
   2        CONTINUE
   3     CONTINUE
      RETURN
      END
        SUBROUTINE S2HFB(K,M,N,F,SCRF,IPTX,IPTY,IC,IA)
C
C THIS SUBROUTINE SHUFFLES THE SOLUTION BACK TO THE CORRECT ORDER
C WHEN A PERIODIC CONDITION IS FOUND
C
         INTEGER M,N,IC,IA,K
         INTEGER IPTX(M),IPTY(N)
         DOUBLE PRECISION F(M,N),SCRF(M,N)
         INTEGER I,J,I2,J2,IDISP,JDISP
C
C CALLS FROM BLAS: DCOPY
C
         CALL DCOPY(M*N,F,1,SCRF,1)
         IDISP=0
         IF (IA.EQ.4)IDISP=K-1
         JDISP=0
         IF (IC.EQ.4)JDISP=K-1
         DO 20 I=1,M
           IDI=I+IDISP
           I2=IPTX(IDI)
           DO 10 J=1,N
             JDJ=J+JDISP
             J2=IPTY(JDJ)
             F(I,J)=SCRF(I2,J2)
10          CONTINUE
20       CONTINUE
         RETURN
         END
        SUBROUTINE D2IRFN(K,KSY,N,NINTY,Y,IPTY,GAM,C,IC,RESTAR,WS)
C
C THIS SUBROUTINE GETS THE COEFFICIENTS WHEN DIRICHLET BOUNDARY
C CONDITIONS ARE DEFINED BY THE FUNCTION GAM
C TO DETERMINE THE B-SPLINE COEFFICENTS A LEAST SQUARES
C PROBLEMS IS SOLVED WHERE THE COEFFICIENT MATRIX TURNS
C OUT TO BE THE MASS MATRIX WHICH MAY BE COMPUTED BEFORE
C
C    INPUT PARAMETERS
C
C    K     ORDER OF THE APPROXIMATION
C    KSY   LEADING DEMSION OF THE MASS MATRIX
C    N     NUMBER OF B-SPLINE COEFFICIENTS TO BE DETERMINED
C    NINTY NUMBER OF Y INTERVALS FOR THE QUADRATURE
C    Y     THE MESH TO BE USED FOR THE QUADRATURE
C    IPTY  PERMUTATION MATRIX TO BE USED IF THE OTHER BOUNDARY
C          IS PERIODIC
C    GAM   FUNCITION DEFINING THE MIXED CONDITION
C    IC    IF IC.EQ.4 THE OTHER DIRECTION IS PERIODIC
C    RESTAR LOGICAL VARIABLE TELLING WHETHER THIS IS SECOND
C          SIDE WHICH NEEDS THE SAME MASS MATRIX
C    WS    SCRATCH SPACE
C
C    OUTPUT PARAMETERS
C
C    C     THE B-SPLINE COEFFICIENTS WHICH HAVE BEEN DETERMED FROM
C          LEAST SQUARES COMPUTATION
C
        EXTERNAL GAM
        DOUBLE PRECISION GAM
        LOGICAL RESTAR
        INTEGER K,KSY,N,NINTY,IC
        INTEGER IPTY(NINTY)
        DOUBLE PRECISION Y(N),C(N)
        INTEGER IMY, IMYB,IWZ,IBY,IBS,IW
        DOUBLE PRECISION WS(1)
        INTEGER NYK1
        DOUBLE PRECISION DS(500)
        INTEGER IS (1000)
        DOUBLE PRECISION ABSI(20,2),WTS(20,2)
        COMMON /GLQNW/ ABSI,WTS
        COMMON /CSTAK/DS
        EQUIVALENCE (DS(1), IS(1))
C
C CALLS FROM PORT: ENTER, ISTKGT, SETERR, ISTKRL, LEAVE
C CALLS FROM LINPACK: DPBSL, DPBFA
C OTHER CALLS : R6GGDY, S2ETD, B2SYGA, P2SHM
C
        CALL ENTER(0)
        NYK1=NINTY+K-1
        K1=K-1
        IMY =1
      IZY = IMY+N*KSY
      IWZ = IZY+K*NINTY
      IBY = IWZ+K
C
C CHECK IF COEFFICENT MATRIX HAS ALREADY BEEN COMPUTED
C AND FACTORED
C
      IF (.NOT.RESTAR)GO TO 5
      IBS = ISTKGT(K*K, 4)
      IW = ISTKGT(K*K, 4)
      IF (IC .GE. 4) GOTO 3
C
C COMPUTE COEFFICIENT MATRIX IN MY
C FOR NON- PERIODIC PROBLEM
C
         CALL R6GGDY(Y, N, K, WS(IMY), K, NINTY, WS(
     1      IZY), WS(IWZ), WS(IBY), DS(IBS), DS(IW))
         GOTO  4
C
C WE ARE WORKING WITH A PERIODIC B-SPLINE IN Y
C GET COEFFICIENT MATRIX WITH BIGGER MESH
C
    3     CONTINUE
         IMYB = ISTKGT(K*NYK1, 4)
         CALL R6GGDY(Y, NYK1, K,  DS(IMYB), K, N
     1      , WS(IZY), WS(IWZ), WS(IBY), DS(IBS), DS(IW))
C
C SHUFFLE INTO MATRIX OF LARGER BANDWIDTH
C
         CALL P2SHM(N, K, DS(IMYB), K, WS(IMY), KSY, IPTY, 1, K, K
     1      , KSY, .FALSE.)
         CALL ISTKRL(1)
C
4       CONTINUE
C
C FACTOR THE COEFFICIENT MATRIX
C
         CALL DPBFA(WS(IMY),KSY,N,KSY-1,INFO)
        IF (INFO.EQ.0) GO TO 5
          CALL SETERR(26HSERRG-CANNOT GET DIRICHLET,26,17,1)
          CALL LEAVE
          RETURN
C
C GENERATE RIGHT HAND SIDE
C
5       CONTINUE
        CALL S2ETD(N,0.0D0,C)
        CALL B2SYGA(Y, NYK1,K,GAM,C,1,WS(IZY),WS(IBY),
     1   K,NINTY,IC,IPTY,.TRUE.)
C
C SOLVE LEAST SQUARES PROBLEM
C
         IF (INFO.EQ.0)CALL DPBSL(WS(IMY),KSY,N,KSY-1,C)
         CALL LEAVE
         RETURN
         END
      SUBROUTINE R6GGDY(Y, N, K,  MY, KM1, NINT, ZY, WY,
     1                  BY, BS, WORK)
      INTEGER           N, K, KM1, NINT
      DOUBLE PRECISION  Y(1), MY(K,N),
     1                  ZY(KM1,NINT), WY(KM1), BY(KM1,NINT,K),
     1                  BS(K,2), WORK(K,K)
C
      INTEGER           I, J, LEFT, IY, II, JJ, IR, IC, LMK, KMJJ
      DOUBLE PRECISION   ZZ, WW, BJW
      DOUBLE PRECISION ABSI(20,2),WTS(20,2)
      COMMON /GLQNW/ ABSI,WTS
C
C THIS SUBROUTINE COMPUTES THE MASS MATRIX  ONLY
C IT IS USED WHEN A DIRICHLET BOUNDARY CONDITION IS ENCOUNTERED
C DEFINED BY A FUNCTION
C
      DO 10 I=1,K
         DO 10 J=1,N
   10       MY(I,J) = 0.0D0
      DO 20 I=1,K
         DO 20 J=1,KM1
            DO 20 IY=1,NINT
   20          BY(J,IY,I) = 0.0D0
C
C
C FOR EACH INTERVAL GET THE QUADRATURE POINTS AND WEIGHTS
C AND ADD THE RESULT OF MULTIPLYING THE B SPLINES TOGETHER
C BECAUSE OF LOCAL SUPPORT THE MASS MATRIX IS BANDED
C
      DO 50 LEFT=K,N
         LMK = LEFT - K
         IY  = LMK + 1
C        GET THE GAUSSIAN QUADRATURE POINTS
         CALL G2QAB(Y(LEFT), Y(LEFT+1), KM1, ZY(1,IY), WY,2)
C        DO NOT PERFORM QUADRATURE ON EMPTY INTERVALS
C        EMPTY INTERVALS WILL ARISE WHEN THERE ARE MULTIPLE KNOTS
         IF (Y(LEFT) .EQ. Y(LEFT+1)) GO TO 50
         DO 40 J=1,KM1
            ZZ = ZY(J,IY)
            WW = WY(J)
            CALL B2SVD1(Y, K, ZZ, LEFT,  BS)
            DO 40 JJ=1,K
               BJW  = BS(JJ,1)*WW
               BY(J,IY,JJ) = BJW
               IC   = LMK + JJ
               KMJJ = K - JJ
               DO 40 II=1,JJ
                  IR = KMJJ + II
   40             MY(IR,IC) = MY(IR,IC) + BS(II,1)*BJW
   50    CONTINUE
C
      RETURN
      END
      SUBROUTINE G2ETBZ(Y, N, K,   KM1, NINT, ZY, WY,
     1                  BY, BS, WORK)
      INTEGER           N, K, KM1, NINT
      DOUBLE PRECISION  Y(1),
     1                  ZY(KM1,NINT), WY(KM1), BY(KM1,NINT,K),
     1                  BS(K,2), WORK(K,K)
C
      INTEGER           I, J, LEFT, IY,  JJ,   LMK
      DOUBLE PRECISION   ZZ, WW, BJW
      DOUBLE PRECISION ABSI(20,2),WTS(20,2)
      COMMON /GLQNW/ ABSI,WTS
C
C THIS SUBROUTINE COMPUTES GQUSSIAN QUADRATURE POINTS
C EVALUATES THE B-SPLINES AT THEM BUT DOES NOT GO AHEAD AND
C FORM THE M MATRIX. IT IS USED WHEN THE M MATRIX IN Y HAS
C ALREADY BEEN FORMED BUT WITH MORE GAUSSIAN POINTS
C
C
C CALLS  G2QAB, B2SVD1
C
      DO 20 I=1,K
         DO 20 J=1,KM1
            DO 20 IY=1,NINT
   20          BY(J,IY,I) = 0.0D0
C
      DO 50 LEFT=K,N
         LMK = LEFT - K
         IY  = LMK + 1
         CALL G2QAB(Y(LEFT), Y(LEFT+1), KM1, ZY(1,IY), WY,1)
         DO 40 J=1,KM1
            ZZ = ZY(J,IY)
            WW = WY(J)
            CALL B2SVD1(Y, K, ZZ, LEFT,  BS)
            DO 40 JJ=1,K
               BJW  = BS(JJ,1)*WW
               BY(J,IY,JJ) = BJW
  40       CONTINUE
   50    CONTINUE
C
      RETURN
      END
          LOGICAL FUNCTION B2DMSH(N,K,X,IA)
C
C THIS FUNCTION DETERMINES IF A MESH IS SPECIFIED INCORRECTLY
C IA INDICATES WHETHER THE MESH IS NORMAL OR PERIODIC
C
          INTEGER N,K,IA
          DOUBLE PRECISION X(N),XX,XB
          B2DMSH=.TRUE.
          I=1
          IF(IA.EQ.4) GO TO 20
            XX=X(1)
            DO 10 I=2,K
               IF(X(I).NE.XX) RETURN
10          CONTINUE
            I=K
20        IF (X(I).GE.X(I+1)) RETURN
          IB=I+1
          IEND=N-1
          IF (IEND.LT.IB) GO TO 40
          DO 30 I=IB,IEND
             IF (X(I).GT.X(I+1)) RETURN
30        CONTINUE
40        IF (X(N).GE.X(N+1)) RETURN
          IF (IA.EQ.4) GO TO 100
          IBEG=N+2
          IEND=N+K
          XB=X(N+1)
          DO 50 I=IBEG,IEND
             IF (X(I).NE.XB)RETURN
50        CONTINUE
100       B2DMSH=.FALSE.
          RETURN
          END
      SUBROUTINE R2G1Y(Y, N, K, P2, Q2, SY, MY, KM1, NINT, ZY, WY,
     1                  BY, BS, WORK, NCOLY, IERY)
C
C THIS SUBROUTINE COMPUTES THE S AND M MATRICES IN THE Y DIRECTION
C IT SAVES IN BY THE WEIGHTED BASIS SPLINES AT THE QUADRATURE POINTS
C
C     INPUT PARAMETERS
C
C     Y       THE Y MESH
C     N       NUMBER OF ROWS IN THE MASS AND STIFFNESS MATRICES
C     K       ORDER OF THE APPROXIMATION, SY AND KY WILL HAVE 2K-1
C             NONZERO DIAGONALS
C     P2      A FUNCTION WHICH HELPS DEFINE THE PDE
C     Q2      ANOTHER FUNCTION WHICH HELPS DEFINE THE PDE, SEE THE
C             BEGINNING COMMENTS OF SERRG2
C     KM1     THE NUMBER OF QUADRATURE POINTS IN EACH INTERVAL
C     NINT    NUMBER OF INTERVALS
C     NCOLY   THE WEIGHTS AND POINTS FOR GUASSIAN QUADRATURE HAVE BEEN
C             PRECOMPUTED FOR A GENERAL INTERVAL AND LEFT IN THE NCOLY
C             COLUMN IN THE NAMED COMMON
C
C     OUTPUT PARAMETERS
C
C     SY      THE STIFFNESS MATRIX, EACH OF WHOSE ELEMENTS INVOLVES
C             INTEGRATING THE PRODUCT OF 2 B-SPLINES, THEIR DERIVATIVES
C             AND THE P2 AND Q2 FUNCTIONS. BECAUSE EACH OF THE B SPLINES
C             HAS FINITE SUPPORT, THE SY MATRIX IS BANDED. IT IS STORED
C             AS A SYMMETRIC BANDED MATRIX IN THE LINPACK FORMAT
C     MY      THE MASS MATRIX, EACH OF WHOSE ELEMENTS IS THE INTEGRAL OF
C             THE PRODUCT OF 2 B-SPLINES. BECAUSE EACH OF THE B SPLINES
C             HAS FINITE SUPPORT, THE MY MATRIX IS BANDED AND IT IS
C             STORED AS A SYMMETRIC BANDED MATRIX IN THE LINPACK FORMAT
C      ZY     THE GAUSSIAN QUADRATURE POINTS
C      BY     THE BSPLINES EVALUATED AT THE GAUSSIAN QUADRATURE POINTS
C      IERY   IF NONZERO, THEN P2 IS NEGATIVE
C
C      WORK SPACE
C
C      WY     WILL BE FILLED WITH GAUSSIAN QUADRATURE WEIGHTS ON AN
C             INTERVAL
C      BS     WILL BE FILLED WITH THE EVALUTION OF A B-SPLINE AND ITS
C             DERIVATIVE AT A GAUSSIAN QUADRATURE POINT
C      WORK   SPACE NEEDED BY THE UNDERLYING ROUTINE THAT EVALUATES
C             B-SPLINES AT SPECIFIC POINTS.
C
      INTEGER           N, K, KM1, NINT, IERY
      DOUBLE PRECISION  Y(1), P2, Q2, SY(K,N), MY(K,N),
     1                  ZY(KM1,NINT), WY(KM1), BY(KM1,NINT,K),
     1                  BS(K,2), WORK(K,K)
C
      INTEGER           I, J, LEFT, IY, II, JJ, IR, IC, LMK, KMJJ
      DOUBLE PRECISION  PY, QY, ZZ, WW, BJW, PBJW, QBJW
      EXTERNAL P2, Q2
      DOUBLE PRECISION ABSI(20,2),WTS(20,2)
      COMMON /GLQNW/ ABSI,WTS
C
C CALLS G2QAB, B2SVD1
C
      DO 10 I=1,K
         DO 10 J=1,N
            SY(I,J) = 0.0D0
   10       MY(I,J) = 0.0D0
      DO 20 I=1,K
         DO 20 J=1,KM1
            DO 20 IY=1,NINT
   20          BY(J,IY,I) = 0.0D0
C
C
C DO EACH INTERVAL AT A TIME RATHER THAN A PARTICULAR ELEMENT
C OF THE SY AND MY MATRICES. THUS SEVERAL ELEMENTS MAY BE UPDATED
C AT ONCE.
C
      IERY = 0
      DO 50 LEFT=K,N
         LMK = LEFT - K
         IY  = LMK + 1
C
C DO NOT PERFORM QUADRATURE ON EMPTY INTERVALS.  EMPTY INTERVALS
C WILL ARISE WHEN THERE ARE INTERIOR MULTIPLE KNOTS.
         IF (Y(LEFT) .EQ. Y(LEFT+1))  GO TO 50
C GET GAUSS QUADRATURE POINTS FOR THIS INTERVAL
C
         CALL G2QAB(Y(LEFT), Y(LEFT+1), KM1, ZY(1,IY), WY, NCOLY)
         DO 40 J=1,KM1
            ZZ = ZY(J,IY)
            WW = WY(J)
            PY = P2(ZZ)
            IF (PY.LE.0.0D0) GO TO 100
            QY = Q2(ZZ)
C
C EVALUATE THE B-SPLINES AND THEIR DERIVATIVES AT THE GAUSSIAN Q. POINTS
C
            CALL B2SVD1(Y, K, ZZ, LEFT,  BS)
            DO 40 JJ=1,K
               BJW  = BS(JJ,1)*WW
               BY(J,IY,JJ) = BJW
               PBJW = PY*BS(JJ,2)*WW
               QBJW = QY*BJW
               IC   = LMK + JJ
               KMJJ = K - JJ
               DO 40 II=1,JJ
                  IR = KMJJ + II
                  SY(IR,IC) = SY(IR,IC) + PBJW*BS(II,2) + QBJW*BS(II,1)
   40             MY(IR,IC) = MY(IR,IC) + BS(II,1)*BJW
C
   50    CONTINUE
C
      RETURN
100   IERY=1
      RETURN
      END
      SUBROUTINE R2G1X(X,M,K,P1,R,Q1,F,KM1,NY,ZY,BY,KX,
     1                  SX, MX, N, FXY, ZX, WX, BS, WORK, FZY, FT,
     2                  IERX)
C
C THIS SUBROUTINE DETERMINES THE RIGHT HAND SIDE AND THE S AND M
C MATRICES IN THE X DIRECTION. IT USES THE MATRIX BY COMPUTED BY R2G1Y
C
C
C     INPUT PARAMETERS
C
C     X       THE X MESH
C     M       NUMBER OF ROWS IN THE MASS AND STIFFNESS MATRICES
C     K       ORDER OF THE APPROXIMATION, SX AND KX WILL HAVE 2K-1
C             NONZERO DIAGONALS
C     P1      A FUNCTION WHICH HELPS DEFINE THE PDE
C     R       A FUNCTION WHICH HELPS DEFINE THE PDE
C     Q1      ANOTHER FUNCTION WHICH HELPS DEFINE THE PDE, SEE THE
C             BEGINNING COMMENTS OF SERRG2
C     KM1     THE NUMBER OF QUADRATURE POINTS IN EACH INTERVAL
C     NY      NUMBER OF GAUSSIAN QUADRATURE POINTS IN THE Y DIRECTION
C     ZY      GAUSSIAN QUADRATURE POINTS IN THE Y DIRECTION
C     BY      B-SPLINES IN Y EVALUATED AT THE ZYS
C     F       THE USER DEFINED FUNCTION WHICH EVALUATES THE RIGHT HAND
C             SIDE AT ONE VALUE X AND ALL THE Y GAUSSIAN QUADRATURE
C             POINTS
C     KX      ROW DIMENSION OF THE MASS AND STIFFNESS MATRICES
C
C     OUTPUT PARAMETERS
C
C     SX      THE STIFFNESS MATRIX, EACH OF WHOSE ELEMENTS INVOLVES
C             INTEGRATING THE PRODUCT OF 2 B-SPLINES, THEIR DERIVATIVES
C             AND THE P1, R AND Q1 FUNCTIONS. BECAUSE EACH OF THE B
C             SPLINES HAS FINITE SUPPORT, THE SY MATRIX IS BANDED.IT IS
C             STORED AS AN UNSYMMETRIC BANDED MATRIX IN THE LINPACK
C             FORMAT THE R FUNCTION PRODUCES AN UNSYMMETRIC MATRIX
C     MX      THE MASS MATRIX, EACH OF WHOSE ELEMENTS IS THE INTEGRAL OF
C             THE PRODUCT OF 2 B-SPLINES. BECAUSE EACH OF THE B SPLINES
C             HAS FINITE SUPPORT, THE MY MATRIX IS BANDED AND IT IS
C             STORED AS AN UNSYMMETRIC BANDED MATRIX IN THE LINPACK
C             FORMAT
C     FXY     THE RIGHT HAND SIDE EACH OF WHOSE ELEMENTS IS THE
C             INTEGRAL OF THE PRODUCT OF F AND THE B SPLINES IN
C             BOTH X AND Y
C     IERX    IF NONZERO, P1 IS NEGATIVE SOMEWHERE
C
C      WORK SPACE
C
C      WX     WILL BE FILLED WITH GAUSSIAN QUADRATURE WEIGHTS ON AN
C             INTERVAL
C      ZX     WILL BE FILLED WITH GAUSSIAN QUADRATURE POINTS ON AN
C             INTERVAL
C      BS     WILL BE FILLED WITH THE EVALUTION OF A B-SPLINE AND ITS
C             DERIVATIVE AT A GAUSSIAN QUADRATURE POINT IN X
C      WORK   SPACE NEEDED BY THE UNDERLYING ROUTINE THAT EVALUATES
C             B-SPLINES AT SPECIFIC POINTS.
C      FZY    WILL CONTAIN THE FUNCTIONF AT ONE X AND ALL THE ZY'S
C      FT     WILL CONTAIN THE PRODUCT OF FZY AND THE Y BSPLINES FOR
C             EACH X
C
C   FOR HIGHER ORDER APPROXIMATIONS THIS IS MOST TIME CONSUMING PORTION
C   OF THE WHOLE PACKAGE. IT HAS BEEN VECTORIZED IN THE SENSE THAT
C   INNER LOOPS NOW RUN OVER THE NUMBER OF GAUSSIAN QUADRATURE POINTS
C   IN Y RATHER THAN OVER THE ORDER OF APPROXIMATION.
      INTEGER           M, K, KM1, NY, KX, N
      DOUBLE PRECISION X(1), P1, R, Q1, ZY(KM1,NY), BY(KM1,NY,K),
     1                  SX(KX,M), MX(KX,M), FXY(M,N), ZX(KM1), WX(KM1),
     1                  BS(K,2), WORK(K,K), FZY(KM1,NY), FT(N)
      DOUBLE PRECISION BJW,PX,RX,QX,WW,ZZ,QBJW,BJW2
      EXTERNAL P1,Q1,R,F
      DOUBLE PRECISION ABSI(20,2),WTS(20,2)
      COMMON /GLQNW/ ABSI,WTS
C
C CALLS G2QAB, B2SVD1
C
      DO 10 J = 1,M
         DO 10 I = 1,KX
            SX(I,J) = 0.0D0
  10        MX(I,J) = 0.0D0
      DO 20 J = 1,N
         DO 20 I = 1,M
  20        FXY(I,J) = 0.0D0
C
C
C PROCESS ONE INTERVAL AT A TIME RATHER THAN ONE ELEMENT OF F
C THIS INCREASE THE LOOP LENGTH OF THE INNERMOST LOOPS
C
      IERX=0
      DO 100 LEFT = K, M
         LMK = LEFT-K
C
C GET THE GAUSSIAN QUADRATURE POINTS FOR THIS X INTERVAL
C
         CALL G2QAB(X(LEFT), X(LEFT+1), KM1, ZX, WX, 1)
C
C DO NOT PERFORM QUADRATURE ON EMPTY INTERVALS.  EMPTY INTERVALS
C WILL ARISE WHEN THERE ARE INTERIOR MULTIPLE KNOTS.
C
         IF (X(LEFT) .EQ. X(LEFT+1))  GO TO 100
         DO 90 J = 1, KM1
            ZZ = ZX(J)
            WW = WX(J)
            PX = P1(ZZ)
            IF (PX.LE.0.D0) GO TO 110
            RX = R(ZZ)
            QX = Q1(ZZ)
C
C EVALUATE THE B SPLINES IN X AT THE GAUSSIAN QUAD. POINTS
C
            CALL B2SVD1(X, K, ZZ, LEFT,  BS )
C
C THIS IS A CALL TO F WITH ONE X AND ALL THE YS
C
            CALL F(ZZ, ZY, FZY, KM1*NY)
            DO 30 JJ= 1,N
C
C FOR EACH Y INTERVAL ADD ALL THE INFORMATION ONE HAS ABOUT THIS X
C
  30           FT(JJ) = 0.0D0
            DO 40 I3 = 1,K
               I3M1 = I3-1
               DO 40 I1 = 1,KM1
                  DO 40 I2 = 1,NY
                     I2I3M1=I2+I3M1
  40                 FT(I2I3M1) = FT(I2I3M1) + FZY(I1,I2)*BY(I1,I2,I3)
C
C NOW WORK ON S AND M WHERE S IS NONSYMMETRIC AND M IS SYMMETRIC
C
            DO 80 JJ = 1,K
               BJW  = BS(JJ,1)*WW
               BJW2  = BS(JJ,2)*WW
               QBJW= BJW*QX
               IC   = LMK + JJ
               IR0  = (2*K-1) - JJ
               DO 50 II = 1, K
                  IR = II + IR0
  50              SX(IR,IC) = SX(IR,IC) + QBJW*BS(II,1) +
     1                        BJW2*(RX*BS(II,1) + PX*BS(II,2))
               DO 60 II = 1,JJ
                  IR = II + IR0
  60              MX(IR,IC) = MX(IR,IC) + BS(II,1)*BJW
               DO 70 I2=1,N
  70              FXY(IC,I2) = FXY(IC,I2) + FT(I2)*BJW
  80           CONTINUE
  90        CONTINUE
 100     CONTINUE
C
      RETURN
110   IERX=1
      RETURN
      END
      SUBROUTINE G2QAB (A, B, N, NODE, WEIGHT, NCOL)
      INTEGER          N
      DOUBLE PRECISION A, B, NODE(N), WEIGHT(N)
C
C     THIS ROUTINE RETURNS THE NODES AND WEIGHTS FOR THE N-POINT
C     GAUSS-LEGENDRE QUADRATURE RULE SCALED TO THE INTERVAL (A,B)
C
C
      DOUBLE PRECISION  Z(20,2),W(20,2)
      COMMON  /GLQNW/   Z,W
      DOUBLE PRECISION  BA
C
      IF (N .LE. 0  .OR.  11 .LT. N) CALL SETERR(
     1                 22HG2QAB - N OUT OF RANGE, 22,10, 2)
C
      BA = (B-A)/2.
      DO 10 I=1,N
         NODE(I)     = BA*(1.0D0 + Z(I,NCOL)) + A
         WEIGHT(I)   = BA*W(I,NCOL)
10    CONTINUE
C
      RETURN
      END
           SUBROUTINE B2SVD1(T,K,X,LEFT,DBIATX)
C
C   THIS SUBROUTINE COMPUTES THE VALUE AND 1 DERIVATIVE OF THE BSPLINES
C   WHICH DO NOT VANISH AT X
C   IT IS A FASTER VERSION OF BSPLVD TAILORED TO ONE DERIVATIVE
C   INPUT PARAMETERS
C
C   T        MESH IN 1 D
C   K        ORDER OF THE APPROXIMATION
C   X        POINT AT WHICH B-SPLINES AND FIRST DERIVATIVE ARE TO BE
C            COMPUTED
C   LEFT     PARAMETER NEEDED BY B2SLB, DE BOORS B-SPLINE EVALUATOR
C   OUTPUT PARAMETERS
C   DBIATX   ALL THE B SPLINES THAT ARE NONZERO AT X ARE CONTAINED IN
C            ITS FIRST COLUMNS, AND THEIR DERIVATIVES AT X IN THE
C            SECOND COLUMN
        DOUBLE PRECISION DBIATX(K,2),T(1),X,FKP,F2,F3
C
C CALLS BSPLVB FROM DE BOOR'S B SPLINE PACKAGE
C
        KM1=K-1
        CALL BSPLVB(T,KM1,1,X,LEFT,DBIATX)
        DO 5 I=1,KM1
            DBIATX(I+1,2)=DBIATX(I,1)
5        CONTINUE
        CALL BSPLVB(T,K,2,X,LEFT,DBIATX)
        FKP=DBLE(FLOAT(KM1))
        IL=LEFT-K+2
        ILPKM1=IL+KM1
        F2=FKP/(T(ILPKM1)-T(IL))
        DBIATX(1,2)=-DBIATX(2,2)*F2
        IF (KM1.LT.2) GO TO 20
        DO 10 I=2,KM1
           IL=IL+1
            ILPKM1=IL+KM1
            F3=FKP/(T(ILPKM1)-T(IL))
            DBIATX(I,2)=DBIATX(I,2)*F2-F3*DBIATX(I+1,2)
            F2=F3
10      CONTINUE
20      CONTINUE
        LEKM1=LEFT+KM1
        DBIATX(K,2)=FKP*DBIATX(K,2)/(T(LEKM1)-T(LEFT))
        RETURN
        END
      SUBROUTINE G2QM11(N,X,W)
C
C  TO COMPUTE THE ABCISSAE (X) AND WEIGHTS (W) FOR AN N POINT
C  GAUSS-LEGENDRE QUADRATURE RULE ON (-1,+1).
C
C  SCRATCH SPACE ALLOCATED - 17*N DOUBLE PRECISION WORDS.
C
C   THIS SUBROUTINE CHECKS PARAMETERS AND OBTAINS SCRATCH SPACE
C   AND CALLS A SUBROUTINE WHICH IMPLEMENTS SACK AND DONOVANS 1972
C   ALGORITHM. IT AND ITS CALLEES WERE LIFTED FROM THE PORT LIBRARY
C   LOWER DOWN AND EIGENSYSTEM MUST BE SOLVED FOR ALL THE
C   EIGENVALUES AND ONE VECTOR. THE LINPACK SUBROUTINE TQL2 IS
C   MODIFIED TO SOLVE THIS SYSTEM.
C
C  ERROR STATES -
C
C    1 - N.LT.1.
C
      DOUBLE PRECISION X(N),W(N)
C
      COMMON /CSTAK/S
      DOUBLE PRECISION S(500)
C
C CALLS FROM PORT: SETERR, ISTKGT, ISTKRL
C OTHER CALLS: D2GAUS
C
      IF (N.LT.1)
     1   CALL SETERR(15HG2QM11 - N.LT.1,15,1,2)
C
      N2=2*N
C
C ... ALLOCATE SCRATCH SPACE FOR A,B,C AND NU.
C
      IA=ISTKGT(8*N,4)-1
      IB=IA+N2
      IC=IB+N2
      INU=IC+N2
C
C ... COMPUTE A,B,C AND NU FOR CALL TO D2GAUS.
C
      DO 10 I=1,N2
         IDXA=IA+I
         IDXB=IB+I
         IDXC=IC+I
         IDXNU=INU+I
         S(IDXA)=DBLE(FLOAT(I))/DBLE(FLOAT(2*I-1))
         S(IDXB)=0.0D0
         S(IDXC)=DBLE(FLOAT(I-1))/DBLE(FLOAT(2*I-1))
 10      S(IDXNU)=0.0D0
      S(INU+1)=2.0D0
C
C ... GET X AND W.
C
      CALL D2GAUS(N,S(IA+1),S(IB+1),S(IC+1),S(INU+1),X,W)
C
      CALL ISTKRL(1)
C
      RETURN
C
      END
      SUBROUTINE D2GAUS(N,A,B,C,NU,X,W)
C
C  LET THE RECURRENCE RELATION
C
C  X*P(X,L)=A(L)*P(X,L+1)+B(L)*P(X,L)+C(L)*P(X,L-1), L=1,...,2*N
C
C  DEFINE POLYNOMIALS OF DEGREE L-1, WITH P(X,0)=0 AND P(X,1)=1.
C
C  LET NU(L)=INTEGRAL(A TO B)(W(X)*P(X,L))DX, L=1,...,2*N,
C
C  WHERE W(X) IS SOME NON-NEGATIVE WEIGHT FUNCTION ON (A,B).
C
C  THIS ROUTINE THEN RETURNS (X(I),W(I)), I=1,...,N, SUCH THAT
C  THE QUADRATURE RULE GIVEN BY
C
C      INTEGRAL(A TO B)(W(X)*F(X))DX = SUM(I=1,...,N)(W(I)*F(X(I)))
C
C  IS EXACT FOR ALL POLYNOMIALS OF DEGREE LESS THAN 2*N.
C
C  SCRATCH SPACE ALLOCATED - 9*N DOUBLE PRECISION WORDS.
C
C  ERROR STATES -
C
C    1 - N.LT.1.
C    2 - A(I)=0.
C    3 - ALL P(X,L) ARE IDENTICALLY ZERO.
C    4 - CANNOT OBTAIN X AND W.
C    5 - CANNOT OBTAIN X AND W.
C
C  SACK AND DONOVAN, NUM. MATH. 18, 465-478(1972).
C
      DOUBLE PRECISION A(1),B(1),C(1),NU(1),X(N),W(N)
C     DOUBLE PRECISION (A,B,C,NU)(2*N)
C
      COMMON /CSTAK/S
      DOUBLE PRECISION S(500)
      DOUBLE PRECISION ABCMAX,MEPS
C
C CALLS TO PORT: I1MACH, SETERR, ISTKGT, ISTKRL
C OTHER CALLS D6AUSQ, N2TQL
C
C
C ... MEPS IS THE ROUNDING ERROR LEVEL OF THE MACHINE.
C
      MEPS=DBLE(FLOAT(I1MACH(10)))**(1-I1MACH(14))
C
      IF (N.LT.1) CALL SETERR(15HD2GAUS - N.LT.1,15,1,2)
C
      NM1=N-1
      N2=2*N
C
C ... SCALE TO PREVENT OVER OR UNDER FLOW.
C
      ABCMAX=0.0D0
      DO 10 I=1,N2
         IF (A(I).EQ.0.0D0) CALL SETERR(15HD2GAUS - A(I)=0,15,2,2)
 10      ABCMAX=DMAX1(ABCMAX,DABS(A(I)),DABS(B(I)),DABS(C(I)))
      IF (ABCMAX.EQ.0.0D0)
     1   CALL SETERR(40HD2GAUS - ALL P(X,L) ARE IDENTICALLY ZERO,
     2               40,3,2)
      DO 20 I=1,N2
         A(I)=A(I)/ABCMAX
         B(I)=B(I)/ABCMAX
         C(I)=C(I)/ABCMAX
 20      NU(I)=NU(I)/ABCMAX
C
C ... ALLOCATE SCRATCH SPACE.
C
      IAC=ISTKGT(9*N,4)-1
      IB=IAC+N2
      IA2=IB+N
      IE=IA2+N
      IR1=IE+N
      IR2=IR1+N2
C
      DO 30 L=2,N2
         IDXAC=IAC+L
 30      S(IDXAC)=A(L-1)*C(L)
C
C ... COMPUTE MATRIX.
C
      CALL D6AUSQ(A,B,S(IAC+1),NU,S(IB+1),S(IA2+1),N,S(IR1+1),S(IR2+1))
C
      IF (1.GT.NM1) GO TO 50
      DO 40 L=1,NM1
         IDXE=IE+L+1
         IDXA2=IA2+L
         IDXB=IB+L
         IF (S(IDXA2).LT.0.0D0) CALL SETERR
     1      (30HD2GAUS - CANNOT OBTAIN X AND W,30,4,2)
         S(IDXE)=DSQRT(S(IDXA2))
 40      X(L)=S(IDXB)
 50   X(N)=S(IA2)
C
C ... SOLVE THE EIGENSYSTEM.
C
      CALL N2TQL(1,N,X,S(IE+1),W,IERR,1)
C
      IF (IERR.NE.0)
     1   CALL SETERR(30HD2GAUS - CANNOT OBTAIN X AND W,30,5,2)
C
      DO 60 L=1,N
         X(L)=X(L)*ABCMAX
 60      W(L)=NU(1)*(W(L)**2)*ABCMAX
C
      CALL ISTKRL(1)
C
      RETURN
C
      END
      SUBROUTINE D6AUSQ(A,B,AC,NU,BETA,ALPHA2,N,R1,R2)
C
C THIS IS A TRANSLATION OF THE ALGOL PROGRAM IN
C SACK AND DONOVAN, NUM. MATH. 18, 465-478(1972).
C IT HAS BEEN LIFTED EN MASSE FROM THE PORT LIBRARY
C
      DOUBLE PRECISION A(1),B(1),AC(1),NU(1),BETA(1),ALPHA2(1),R1(2)
      DOUBLE PRECISION R2(1),TERM,SIGMA
      LOGICAL EVEN
      IMAX=N
      JMAX=IMAX+N
      TERM=1.0D0/NU(1)
      DO 10 J=1,JMAX
         R1(J)=TERM*NU(J)
         R2(J)=0.0D0
 10      TERM=TERM*A(J)
      SIGMA=R1(2)+B(1)
      BETA(1)=SIGMA
      EVEN=.TRUE.
      IF (2.GT.IMAX) GO TO 80
      DO 70 I=2,IMAX
         IP1=I+1
         IF (EVEN) GO TO 40
         JMAX=JMAX-1
         EVEN=.NOT.EVEN
         IF (I.GT.JMAX) GO TO 25
         DO 20 J=I,JMAX
 20         R1(J)=(B(J)-SIGMA)*R2(J)+R2(J+1)+AC(J)*R2(J-1)-R1(J)
 25      TERM=R1(I)
         ALPHA2(I-1)=TERM
         TERM=1.0D0/TERM
         R1(I)=1.0D0
         IF (IP1.GT.JMAX) GO  TO 35
         DO 30 J=IP1,JMAX
 30         R1(J)=R1(J)*TERM
 35      SIGMA=R1(IP1)-R2(I)+B(I)
         GO TO 70
 40      JMAX=JMAX-1
         EVEN=.NOT.EVEN
         IF (I.GT.JMAX) GO TO 55
         DO 50 J=I,JMAX
 50         R2(J)=(B(J)-SIGMA)*R1(J)+R1(J+1)+AC(J)*R1(J-1)-R2(J)
 55      TERM=R2(I)
         ALPHA2(I-1)=TERM
         TERM=1.0D0/TERM
         R2(I)=1.0D0
         IF (IP1.GT.JMAX) GO TO 65
         DO 60 J=IP1,JMAX
 60         R2(J)=R2(J)*TERM
 65      SIGMA=R2(IP1)-R1(I)+B(I)
 70      BETA(I)=SIGMA
 80   ALPHA2(IMAX)=0.0D0
      RETURN
      END
      SUBROUTINE S2ETD(N,V,B)
C
C     S2ETD SETS THE N DOUBLE PRECISION ITEMS IN B TO V
C
      DOUBLE PRECISION B(N),V
C
      IF(N .LE. 0) RETURN
C
      DO 10 I = 1, N
 10     B(I) = V
C
      RETURN
C
      END
      SUBROUTINE N2TQL(NM,N,D,E,Z,IERR,NZ)
C
      INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR
      DOUBLE PRECISION D(N),E(N),Z(NM,NZ)
      DOUBLE PRECISION C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,P2THAG
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
C     WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
C ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C     IT HAS BEEN MODIFIED SO THAT ONE CAN GET ONE EIGENVECTOR WHICH
C     IS NEEDED BY THE GAUSS-QUADRATURE RULE DEFINING ROUTINE. THUS
C     NZ IS THE NUMBER OF EIGENVECTORS REQUESTED.
C :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C
C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
C     FULL MATRIX TO TRIDIAGONAL FORM.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
C          THE IDENTITY MATRIX.
C
C       NZ THE NUMBER OF EIGENVECTORS REQUESTED
C
C      ON OUTPUT
C
C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
C          UNORDERED FOR INDICES 1,2,...,IERR-1.
C
C        E HAS BEEN DESTROYED.
C
C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
C          EIGENVALUES.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     CALLS P2THAG FOR  DSQRT(A*A + B*B) .
C
C
C     THIS VERSION DATED APRIL 1985.
C
C     ------------------------------------------------------------------
C
      IERR = 0
      IF(NZ.EQ.1)Z(1,1)=1.0D0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      IF (NZ.GT.1) GO TO 103
          DO 102 I=2,N
            Z(1,I)=0.0D0
102       CONTINUE
103   CONTINUE
      F = 0.0D0
      TST1 = 0.0D0
      E(N) = 0.0D0
C
      DO 240 L = 1, N
         J = 0
         H = DABS(D(L)) + DABS(E(L))
         IF (TST1 .LT. H) TST1 = H
C     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
         DO 110 M = L, N
            TST2 = TST1 + DABS(E(M))
            IF (TST2 .EQ. TST1) GO TO 120
C     .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP ..........
  110    CONTINUE
C
  120    IF (M .EQ. L) GO TO 220
  130    IF (J .EQ. 30) GO TO 1000
         J = J + 1
C     .......... FORM SHIFT ..........
         L1 = L + 1
         L2 = L1 + 1
         G = D(L)
         P = (D(L1) - G) / (2.0D0 * E(L))
         R = P2THAG(P,1.0D0)
         D(L) = E(L) / (P + DSIGN(R,P))
         D(L1) = E(L) * (P + DSIGN(R,P))
         DL1 = D(L1)
         H = G - D(L)
         IF (L2 .GT. N) GO TO 145
C
         DO 140 I = L2, N
  140    D(I) = D(I) - H
C
  145    F = F + H
C     .......... QL TRANSFORMATION ..........
         P = D(M)
         C = 1.0D0
         C2 = C
         EL1 = E(L1)
         S = 0.0D0
         MML = M - L
C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
         DO 200 II = 1, MML
            C3 = C2
            C2 = C
            S2 = S
            I = M - II
            G = C * E(I)
            H = C * P
            R = P2THAG(P,E(I))
            E(I+1) = S * R
            S = E(I) / R
            C = P / R
            P = C * D(I) - S * G
            D(I+1) = H + S * (C * G + S * D(I))
C     .......... FORM VECTOR ..........
            DO 180 K = 1, NZ
               H = Z(K,I+1)
               Z(K,I+1) = S * Z(K,I) + C * H
               Z(K,I) = C * Z(K,I) - S * H
  180       CONTINUE
C
  200    CONTINUE
C
         P = -S * S2 * C3 * EL1 * E(L) / DL1
         E(L) = S * P
         D(L) = C * P
         TST2 = TST1 + DABS(E(L))
         IF (TST2 .GT. TST1) GO TO 130
  220    D(L) = D(L) + F
  240 CONTINUE
C     .......... ORDER EIGENVALUES AND EIGENVECTORS ..........
      DO 300 II = 2, N
         I = II - 1
         K = I
         P = D(I)
C
         DO 260 J = II, N
            IF (D(J) .GE. P) GO TO 260
            K = J
            P = D(J)
  260    CONTINUE
C
         IF (K .EQ. I) GO TO 300
         D(K) = D(I)
         D(I) = P
C
         DO 280 J = 1, NZ
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P
  280    CONTINUE
C
  300 CONTINUE
C
      GO TO 1001
C     .......... SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ..........
 1000 IERR = L
 1001 RETURN
      END
      DOUBLE PRECISION FUNCTION P2THAG(A,B)
      DOUBLE PRECISION A,B
C
C     FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
C     THIS IS NOT THE VERSION THAT TRIES TO AVOID DIVISION THAT
C     IS IN EISPACK
C
      DOUBLE PRECISION P,Q,R
      P = DMAX1(DABS(A),DABS(B))
      Q = DMIN1(DABS(A),DABS(B))
      IF (Q .EQ. 0.0D0) GO TO 20
         R = (Q/P)**2 +1.D0
         P2THAG=P*DSQRT(R)
         RETURN
   20 P2THAG = P
      RETURN
      END
      SUBROUTINE S2EDM(N, K, SY, MY, D, Z)
      INTEGER N, K
      DOUBLE PRECISION  SY(K, N), MY(K, N), D(N), Z(N, N)
C
C THIS SUBROUTINE SOLVES SY Z = MY D Z, THE GENERALIZED EIGENVALUE
C PROBLEM WHERE SY AND MY ARE N X N MATRICES, SYMMETRIC AND BANDED
C WITH 2*K-1 DIAGONALS STORED AS SYMMETRIC MATRICES ACCORDING TO THE
C LINPAC SCHEME. MY IS POSITIVE DEFINITE AND SY IS POSITIVE SEMIDEFINITE
C      MODIFIED BY D.D.WARNER  3/29/82
C WE DID NOT USE A CODE FROM EISPACK BECAUSE THEY HAVE NO BANDED
C GENERALIZED EIGENVALUE SOLVERS.
C
C INPUT PARAMETERS
C N      NUMBER OF COLUMNS IN MY AND SY
C K      SY AND MY HAVE 2K-1 DIAGONALS
C SY      ARRAY WITH ROW DIMENSION K CONTAINING SY MATRIX
C MY      ARRAY WITH ROW DIMENSION K CONTAINING MY MATRIX
C OUTPUT PARAMETERS
C D       THE N-VECTOR OF EIGENVALUES
C Z       THE NXN MATRIX OF EIGENVECTORS
C
      INTEGER ISTKGT, IFV1, IFV2,  IERR
      INTEGER  IA, IB
C
      DOUBLE PRECISION R(500)
      COMMON  /CSTAK / R
      INTEGER          ISTAK(1000)
      EQUIVALENCE (ISTAK(1), R(1))
C
C CALLS TO PORT: ENTER, SETERR, ISTKGT, LEAVE
C CALLS TO EISPACK: TRED2
C OTHER CALLS: S2ETD, S2ETIN, R2BDUC, N2TQL, R2BBAK
C
C
C CHECK INPUT
      IF (N .LT. 1) CALL SETERR(
     1   16H ERROR IN SEDCMP, 16, 1, 2)
C
C GET SPACE FOR GENERAL MATRICES A AND B, SET IT ALL TO ZERO
C AND THEN PUT SY AND MY INTO IT
C
      CALL ENTER(0)
      IA = ISTKGT(N*N, 4)
      IB = ISTKGT(N*N, 4)
      CALL S2ETD(N*N, 0.0D0, R(IA))
      CALL S2ETD(N*N, 0.0D0, R(IB))
      CALL S2ETIN(N, R(IA), N, SY, K)
      CALL S2ETIN(N, R(IB), N, MY, K)
C
C NOW SOLVE THE GENERALIZED EIGENVALUE PROBLEM
C
C     CALL  R2BDUC  A VERSION OF THE EISPAC SUBROUTINE REDUC
C           WHICH TAKES INTO CONSIDERATION THAT THE ORIGINAL
C           MATRICES WERE BANDED BUT RETURNS A GENERAL MATRIX
C     CALL TRED2 FROM EISPACK
C     CALL N2TQL, A VERSION OF TQL2 FROM EISPACK WHICH PERMITS
C            A VARIABLE NUMBER OF Z VECTORS
C     CALL R2BBAK, A BANDED VERSION OF EISPACK'S REBAK
C
      IFV1 = ISTKGT(N, 4)
      IFV2 = ISTKGT(N, 4)
      ML1 = K-1
      CALL R2BDUC(N, N, R(IA), R(IB), R(IFV2), IERR, ML1, R(IFV1))
      IF (IERR.EQ.0) GO TO 10
          CALL SETERR(26HSERRG-SINGULAR MASS MATRIX,26,13,1)
          CALL LEAVE
          RETURN
10    CONTINUE
      CALL TRED2(N, N, R(IA), D, R(IFV1), Z)
      CALL N2TQL(N, N, D, R(IFV1), Z, IERR, N)
      CALL R2BBAK(N, N, R(IB), R(IFV2), N, Z, ML1, R(IFV1))
      IF (IERR .NE. 0) CALL SETERR(20H  SERRG- EIGEN ERROR, 20,14, 1)
      CALL LEAVE
      RETURN
      END
      SUBROUTINE S2ETIN(N, A, IA, B, K)
      INTEGER N, IA, K
      DOUBLE PRECISION  A(IA, N), B(K, N)
C
C     THIS PROCEDURE TAKES A BANDED SYMMETRIC MATRIX AND MOVES IT TO
C     AN UNBANDED SYMMETRIC MATRIX WITH ONLY THE UPPER TRIANGLE STORED
C
      DO 1 I=1,K
         IF (I .GT. N)  GOTO 2
         KMI=K-I
         DO 1 J=I,N
            JMI=J-I
   1        A(JMI+1,J) = B(KMI+1,J)
C
   2  RETURN
      END
      SUBROUTINE S2ESLV(NX, NY, Z, RHS, IRHS, SX, MX, ISX, D, KX, INFO)
      INTEGER NX, NY, IRHS, ISX, KX
      DOUBLE PRECISION Z(NY,NY),RHS(IRHS,NY),SX(ISX,NX),MX(ISX,NX)
      DOUBLE PRECISION D(NY)
C
C THIS PROCEDURE DOES A SEPARABLE SOLVE USING THE RESULTS OF S2EDM
C
C IT HAS 3 MAIN SECTIONS
C
C  1: MULTIPLY Z(TRANSPOSE) BY RHS (TRANSPOSE) AND CALL RESULT G
C  2: SOLVE THE N BANDED SYSTEMS
C        (SX TENSOR I + MX TENSOR D) H = G
C  3: MULTIPLY Z BY H  AND PUT THE TRANSPOSE OF THE RESULT BACK INTO RHS
C
C INPUT PARAMETERS
C NX    NUMBER OF X POINTS THAT ARE RELEVANT
C NY    NUMBER OF Y POINTS THAT ARE RELEVANT
C Z     EIGENVECTOR MATRIX FROM S2EDM
C RHS   RIGHT HAND SIDE MATRIX
C IRHS  ROW DIMENSION OF RIGHT HAND SIDE MATRIX
C       MUST BE GREATER THAN NX
C SX    STIFFNESS MATRIX IN X STORED IN BANDED FORM
C MX    MASS MATRIX IN X STORED IN BANDED FORM
C ISX   ROW DIMENSION OF MASS AND STIFFNESS MATRICES
C KX    MASS AND STIFFNESS MATRICES HAVE 2KX-1 DIAGONALS
C D     EIGENVALUE MATRIX COMPUTED BY S2EDM
C
C OUTPUT PARAMETERS
C
C RHS   SOLUTION
C INFO  IF 0, NO PROBLEMS, ELSE SINGULAR PROBLEM
C
      INTEGER JOB, ISTKGT,  M, INFO
      INTEGER IPVT, IA, ML, IW
      INTEGER IAB
C
      DOUBLE PRECISION DSTAK(500)
      COMMON  /CSTAK/  DSTAK
      INTEGER          ISTAK(1000)
      EQUIVALENCE (DSTAK(1), ISTAK(1))
C
C CALLS TO PORT: ENTER , ISTKGT, LEAVE
C OTHER CALLS: S2ET3, G2BFAM, G2BSLM, M2M
C
      CALL ENTER(0)
      IW = ISTKGT(NX*NY, 4)
C
C DO MATRIX MULTIPLY Z(T)*RHS(T)
C
      CALL M2M(NY,NX,NY,Z,NY,RHS,IRHS,DSTAK(IW),NY)
C
C
C SOLVE BANDED SYSTEMS
C
      ML   = KX-1
      M    = 2*KX-1
      JOB  = 0
      IAB  = 3*KX-2
      IA   = ISTKGT(IAB*NX*NY, 4)
      IPVT = ISTKGT(NX*NY, 2)
C
C SET UP NY BANDED SYSTEMS IN LINPACK BANDED FORMAT
C
      CALL S2ET3(NX, M, ML, DSTAK(IA), IAB, SX, MX, ISX, D, NY)
      ILL  = ISTKGT(NY, 2)
      ITM  = ISTKGT(NY, 4)
C
C SOLVE ALL THE SYSTEMS SIMULTANESOULY
C
      CALL G2BFAM(DSTAK(IA), IAB, NX, ML, ML, ISTAK(IPVT), INFO,
     1            NY, NY, DSTAK(ITM), ISTAK(ILL))
      IF (INFO.EQ.0) GO TO 10
      CALL LEAVE
      RETURN
10    CONTINUE
      CALL G2BSLM(DSTAK(IA), IAB, NX, ML, ML, ISTAK(IPVT), DSTAK(IW),
     1            JOB, NY, NY, DSTAK(ITM))
C
C DO MATRIX-MATRIX MULTIPLICATIONS Z W =RHS (TRANSPOSE)
C
      CALL M2M(NX, NY, NY, DSTAK(IW), NY,Z,NY, RHS, IRHS)
C
      CALL LEAVE
      RETURN
      END
         SUBROUTINE M2M(M,N,K,A,LDA,B,LDB,C,LDC)
         INTEGER M,N,K,LDA,LDB,LDC
         DOUBLE PRECISION C(LDC,N),A(LDA,M),B(LDB,K)
C
C THIS SUBROUTINE COMPUTES C=A(TRANSPOSE)*B(TRANSPOSE)
C
C   INPUT PARAMATERS
C
C M      NUMBER OF ROWS IN C
C N      NO. OF COLUMNS IN C
C K      NO.OF ROWS IN A
C A      FIRST MATRIX, K X M
C LDA    LEADING DIMENSION OF A, MUST BE GREATER THAN K
C B      SECOND MATRIX, N X K
C LDB    LEADING DIMENSION OF B, MUST BE GREATER THAN N
C LDC    LEADING DIMENSION OF C, MUST BE GREATER THAN M
C
C OUTPUT PARAMETER
C
C C      C = A(TRANSPOSE) * B( TRANSPOSE)
C
          DOUBLE PRECISION BB,BB1
          INTEGER I,J,J2,JMIN
C
C SET THE C MATRIX TO 0.
C
          DO 20 J=1,N
             DO 10 I=1,M
                C(I,J)=0.0D0
10           CONTINUE
20        CONTINUE
C
C WE WILL DO LOOP UNROLLING ONCE
C
           J2=MOD(K,2)
           JMIN=J2+2
           DO 100 I1=1,N
              IF (J2.EQ.0) GO TO 40
              BB=B(I1,1)
              DO 30 I=1,M
                 C(I,I1)=C(I,I1)+A(1,I)*BB
30            CONTINUE
40            CONTINUE
              DO 60 J=JMIN,K,2
                 BB1=B(I1,J-1)
                 BB=B(I1,J)
                 DO 50 I=1,M
                    C(I,I1)=C(I,I1)+A(J-1,I)*BB1+A(J,I)*BB
50               CONTINUE
60            CONTINUE
100        CONTINUE
           RETURN
           END
      SUBROUTINE S2ET3(N, M, ML, A, IA, SY, MY, ISY, D, NMA)
      INTEGER N, M, ML, IA, ISY, NMA
      DOUBLE PRECISION  A(NMA,IA, N), SY(ISY, N), MY(ISY, N), D(NMA)
      DOUBLE PRECISION YMJI,SYJI
C
C THIS SUBROUTINE CREATES NMA BANDED SYSTEMS IN THE ARRAY A
C WHERE THE KTH SYSTEM IS SY + D(K) TIMES MY
C IN LINPACK UNSYMMETRIC BAND FORM
C
C     INPUT PARAMETERS
C
C N    NUMBER OF COLUMNS IN THE SY AND MY ARRAYS
C M    NUMBER OF NONZERO BANDS IN SY AND MY
C ML   NUMBER OF NONZERO DIAGONALS BELOW THE MAIN DIAGONAL
C IA   SECOND DIMENSION OF THE A MATRIY. MUST BE GREATER THAN ML+M
C SY   BANDED MATRIY WITH BANDWIDTH M
C MY   BANDED MATRIY WITH BANDWIDTH M
C D    VECTOR OF LENGTH NMA
C NMA  LENGTH OF D AND LEADING DIMENSION OF A MATRIX
C
C     OUTPUT PARAMETERS
C A    A(K,.,.) CONTAINS THE BANDED SYSTEM SY +D(K) TIMES MY
C      IN LINPACK UNSYMMETRIC BAND FORM
C
      INTEGER MLP1, I, J, MLPM
C
      MLP1 = ML+1
      MLPM = ML+M
      DO 3 I=1,N
         DO 2 J=MLP1,MLPM
            SYJI = SY(J,I)
            YMJI = MY(J,I)
            DO 1 K=1,NMA
   1           A(K,J,I) = SYJI + D(K)*YMJI
   2        CONTINUE
   3     CONTINUE
      RETURN
      END
      SUBROUTINE R2BDUC(NM, N, A, B, DL, IERR, ML, T)
C
      INTEGER I,J,K,N,I1,NM,NN,IERR
      DOUBLE PRECISION  A(NM,N),B(NM,N),DL(N),T(N)
      DOUBLE PRECISION  X,Y,BKI,DLI
      INTEGER IABS
C
C
C     THIS SUBROUTINE REDUCES THE GENERALIZED SYMMETRIC EIGENPROBLEM
C     AX=(LAMBDA)BX, WHERE B IS POSITIVE DEFINITE, TO THE STANDARD
C     SYMMETRIC EIGENPROBLEM USING THE CHOLESKY FACTORIZATION OF B.
C     AND A AND B ARE BANDED WITH BANDWIDTH 2ML+1
C
C   THIS SUBROUTINE IS BASED ON REDUC IN EISPACK BUT ADVANTAGE
C   IS TAKEN THAT THE MATRICES A AND B ARE INITIALLY ASSUMED TO BE
C   BANDED
C
C     ON INPUT-
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRICES A AND B.  IF THE CHOLESKY
C          FACTOR L OF B IS ALREADY AVAILABLE, N SHOULD BE PREFIXED
C          WITH A MINUS SIGN,
C
C        A AND B CONTAIN THE  SYMMETRIC INPUT MATRICES.  ONLY THE
C          FULL UPPER TRIANGLES OF THE MATRICES NEED BE SUPPLIED.  IF
C          N IS NEGATIVE, THE STRICT LOWER TRIANGLE OF B CONTAINS,
C          INSTEAD, THE STRICT LOWER TRIANGLE OF ITS CHOLESKY FACTOR L,
C
C        DL CONTAINS, IF N IS NEGATIVE, THE DIAGONAL ELEMENTS OF L.
C
C        ML THE BANDWIDTHS OF A AND B IS 2*ML -1
C
C     ON OUTPUT-
C
C        A CONTAINS IN ITS FULL LOWER TRIANGLE THE FULL LOWER TRIANGLE
C          OF THE SYMMETRIC MATRIX DERIVED FROM THE REDUCTION TO THE
C          STANDARD FORM.  THE STRICT UPPER TRIANGLE OF A IS UNALTERED,
C
C        B CONTAINS IN ITS STRICT LOWER TRIANGLE THE STRICT LOWER
C          TRIANGLE OF ITS CHOLESKY FACTOR L.  THE FULL UPPER
C          TRIANGLE OF B IS UNALTERED,
C
C        DL CONTAINS THE DIAGONAL ELEMENTS OF L,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          7*N+1      IF B IS NOT POSITIVE DEFINITE.
C    SCRATCH SPACE
C
C    T        LENGTH N
C
C CALLS FROM BLAS: DCOPY
C
C     ------------------------------------------------------------------
C
      IERR = 0
      NN = IABS(N)
      IF (N .LT. 0) GO TO 100
C     1********* FORM L IN THE ARRAYS B AND DL **********
      DO 80 I = 1, N
         I1 = I - 1
C
         JEND= MIN0(N,I+ML)
         DO 80 J = I, JEND
            X = B(I,J)
            IBEG=MAX0(1,J-ML)
            IF(I1.LT.IBEG)GO TO 40
            DO 20 K=IBEG,I1
C
   20       X = X - B(I,K) * B(J,K)
C
   40       IF (J .NE. I) GO TO 60
            IF (X .LE. 0.0D0) GO TO 1000
            Y = DSQRT(X)
            DL(I) = Y
            GO TO 80
   60       B(J,I) = X / Y
   80 CONTINUE
C     1********* FORM THE TRANSPOSE OF THE UPPER TRIANGLE OF INV(L)*A
C                IN THE LOWER TRIANGLE OF THE ARRAY A **********
  100 DO 200 I = 1, NN
             IPML=MIN0(NN,I+ML)
              IMML=MAX0(1,I-ML)
         I1 = I - 1
         Y = DL(I)
C
         DO 200 J = I, IPML
            X = A(I,J)
            IF (I1.LT.IMML)GO TO 180
C
            DO 160 K = IMML, I1
  160       X = X - B(I,K) * A(J,K)
C
  180       A(J,I) = X / Y
  200 CONTINUE
C     1********* PRE-MULTIPLY BY INV(L) AND OVERWRITE **********
C
         DO 300 I = 1, NN
            X=A(I,I)
            IF (I.EQ.1) GO TO 297
           IIM=MAX0(1,I-ML)
           CALL DCOPY(I,A(I,1),NM,T,1)
            I1 = I - 1
C
            IF (I1.LT.IIM) GO TO 280
            DO 220 K = IIM, I1
               BKI=-B(I,K)
                 CALL DAXPY(K,BKI,A(K,1),NM,T,1)
                KP1=K+1
                IF (I1.LT.KP1) GO TO 220
                DO 210 J=KP1,I1
                   T(J)=T(J)+BKI*A(J,K)
210             CONTINUE
220          CONTINUE
C
C
 280       DLI=1.0D0/DL(I)
           DO 290 J=1,I1
              A(I,J)=T(J)*DLI
290        CONTINUE
           IF (I1.LT.IIM) GO TO 297
           DO 295 K=IIM, I1
              X=X-B(I,K)*A(I,K)
295        CONTINUE
297        A(I,I)=X/DL(I)
  300 CONTINUE
C
      GO TO 1001
C     1********* SET ERROR -- B IS NOT POSITIVE DEFINITE **********
 1000 IERR = 7 * N + 1
 1001 RETURN
C     1********* LAST CARD OF R2BDUC **********
      END
      SUBROUTINE R2BBAK(NM, N, B, DL, M, Z, ML, T)
C
      INTEGER I,J,K,M,N,I1,II,NM
      DOUBLE PRECISION B(NM,N),DL(N),Z(NM,M)
      DOUBLE PRECISION T(M),BKI,DLI
C
C
C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A GENERALIZED
C     SYMMETRIC EIGENSYSTEM BY BACK TRANSFORMING THOSE OF THE
C     DERIVED SYMMETRIC MATRIX DETERMINED BY  R2BDUC.
C
C     THIS SUBROUTINE HAS THE SAME FUNCTIOLN AS REBAK IN EISPAC
C     BUT CONSIDERS THE FACT THAT B WAS ORIGINALLY BANDED
C
C     ON INPUT-
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX SYSTEM,
C
C        B CONTAINS INFORMATION ABOUT THE SIMILARITY TRANSFORMATION
C          (CHOLESKY DECOMPOSITION) USED IN THE REDUCTION BY  R2BDUC
C          IN ITS STRICT LOWER TRIANGLE,
C          IT IS BANDED BUT STORED AS A GENERAL MATRIX
C
C        DL CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATION,
C
C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
C
C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
C          IN ITS FIRST M COLUMNS.
C
C        ML BANDWIDTH OF B IS 2*ML-1
C
C     ON OUTPUT-
C
C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
C          IN ITS FIRST M COLUMNS.
C
C
C    SCRATCH SPACE
C
C        T DOUBLE PRECISION VECTOR OF LENGTH M
C
C    THIS SUBROUTINE VECTORIZES WELL. THE INNER LOOP IS OVER M
C    THE NUMBER OF EIGENVECTORS, WHICH FOR SERRG2 IS N AND NOT
C    OVER ML
C
C CALLS FROM BLAS: DCOPY, DAXPY
C     ------------------------------------------------------------------
C
      IF (M .EQ. 0) GO TO 200
C
C     1********* FOR I=N STEP -1 UNTIL 1 DO -- **********
         DO 100 II = 1, N
            I = N + 1 - II
            I1 = I + 1
            I2=MIN0(I+ML,N)
            CALL DCOPY(M,Z(I,1),NM,T,1)
            IF (I2.LT.I1)GO TO 80
C
            DO 60 K = I1, I2
                 BKI=-B(K,I)
                 CALL DAXPY(M,BKI,Z(K,1),NM,T,1)
60          CONTINUE
C
   80       DLI=1.0D0/DL(I)
            DO 90 J=1,M
90             Z(I,J)=DLI*T(J)
  100 CONTINUE
C
  200 RETURN
C     1********* LAST CARD OF R2BBAK **********
      END
      SUBROUTINE G2BFAM(ABD,LDA,N,ML,MU,IPVT,INFO,NMA,LNMA,TM,LL)
      INTEGER LDA,N,ML,MU,IPVT(LNMA,N),INFO
      DOUBLE PRECISION ABD(LNMA,LDA,1)
C
       DOUBLE PRECISION TM(NMA)
       INTEGER LL(NMA)
C
C     G2BFAM FACTORS NMA DOUBLE PRECISION BAND MATRICES SIMULTANEOUSLY
C     WHICH HAVE BEEN STORED ACCORDING TO THE UNSYMMETRIC LINPACK
C     BAND STRUCTURE. WHEN POSSIBLE INNER LOOPS ARE OVER NMA
C     RATHER THAN THE BANDWIDTH.
C
C     ON ENTRY
C
C        ABD     DOUBLE PRECISION(NMA,LDA, N). ABD(I,*.*) CONTAINS
C                THE ITH MATRIX IN BAND STORAGE.  THE COLUMNS
C                OF THE MATRIX ARE STORED IN THE COLUMNS OF  ABD  AND
C                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
C                ML+1 THROUGH 2*ML+MU+1 OF  ABD(I,*,*) .
C                SEE THE COMMENTS BELOW FOR DETAILS.
C
C        LDA     INTEGER
C                THE SECOND  DIMENSION OF THE ARRAY  ABD .
C                LDA MUST BE .GE. 2*ML + MU + 1 .
C
C        N       INTEGER
C                THE ORDER OF THE ORIGINAL MATRIX.
C
C        ML      INTEGER
C                NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
C                0 .LE. ML .LT. N .
C
C        MU      INTEGER
C                NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
C                0 .LE. MU .LT. N .
C                MORE EFFICIENT IF  ML .LE. MU .
C         NMA    NUMBER OF BANDED MATRICES CONTAINED IN ABD
C
C         LNMA   LEADING DIMENSION OF ABD. MUST BE AT LEAST NMA
C
C     ON RETURN
C
C        ABD     ABD(I,.,.) CONTAINS AN UPPER TRIANGULAR MATRIX IN
C                BAND STORAGE FORM FOR THE ITH SYSTEM
C                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
C
C        IPVT    INTEGER(LNMA,N)
C                IPVT(I,.) IS ANN INTEGER VECTOR OF PIVOT INDICES
C                FOR SYSTEM I
C
C        INFO    INTEGER
C                = 0  NORMAL VALUE.
C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
C                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C                     INDICATE THAT G2BSLM WILL DIVIDE BY ZERO IF
C                     CALLED.
C
C   SCRATCH STORAGE
C
C        LL      INTEGER(NMA)
C        TM      DOUBLE PRECISION (NMA)
C
C     BAND STORAGE
C
C           IF  A  IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT
C           WILL SET UP THE INPUT.
C
C                   ML = (BAND WIDTH BELOW THE DIAGONAL)
C                   MU = (BAND WIDTH ABOVE THE DIAGONAL)
C                   M = ML + MU + 1
C                   DO 20 J = 1, N
C                      I1 = MAX0(1, J-MU)
C                      I2 = MIN0(N, J+ML)
C                      DO 10 I = I1, I2
C                         K = I - J + M
C                         ABD(K,J) = A(I,J)
C                10    CONTINUE
C                20 CONTINUE
C
C           THIS USES ROWS  ML+1  THROUGH  2*ML+MU+1  OF  ABD .
C           IN ADDITION, THE FIRST  ML  ROWS IN  ABD  ARE USED FOR
C           ELEMENTS GENERATED DURING THE TRIANGULARIZATION.
C           THE TOTAL NUMBER OF ROWS NEEDED IN  ABD  IS  2*ML+MU+1 .
C           THE  ML+MU BY ML+MU  UPPER LEFT TRIANGLE AND THE
C           ML BY ML  LOWER RIGHT TRIANGLE ARE NOT REFERENCED.
C
C
C     SUBROUTINES AND FUNCTIONS
C
C     FORTRAN MAX0,MIN0
C     CALLS S2XPY
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION T
      INTEGER I,I0,J,JU,JZ,J0,J1,K,KP1,LM,M,MM,NM1
C
C
      M = ML + MU + 1
      MLP1=ML+1
      MP1=M+1
      MM1=M-1
      INFO = 0
C
C     ZERO INITIAL FILL-IN COLUMNS
C
      J0 = MU + 2
      J1 = MIN0(N,M) - 1
      IF (J1 .LT. J0) GO TO 30
      DO 20 JZ = J0, J1
         I0 = M + 1 - JZ
         DO 10 I = I0, ML
            DO 5 MN=1,NMA
            ABD(MN,I,JZ) = 0.0D0
5             CONTINUE
   10    CONTINUE
   20 CONTINUE
   30 CONTINUE
      JZ = J1
      JU = 0
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      NM1 = N - 1
      IF (NM1 .LT. 1) GO TO 130
      DO 120 K = 1, NM1
         KP1 = K + 1
C
C        ZERO NEXT FILL-IN COLUMN
C
         JZ = JZ + 1
         IF (JZ .GT. N) GO TO 50
         IF (ML .LT. 1) GO TO 50
            DO 40 I = 1, ML
               DO 35 MN=1,NMA
               ABD(MN,I,JZ) = 0.0D0
35              CONTINUE
   40       CONTINUE
   50    CONTINUE
C
C        FIND L = PIVOT INDEX
C
         LM = MIN0(ML,N-K)
         LPE=M+LM
          DO 51 MN=1,NMA
             LL(MN)=M
             TM(MN) = DABS(ABD(MN,M,K))
51        CONTINUE
          IPVTK=M
         IF (LPE.LT.MP1) GO TO 54
         DO 53 LP=MP1,LPE
              DO 52 MN=1,NMA
                   IF (DABS(ABD(MN,LP,K)).LT.TM(MN)) GO TO 52
                   TM(MN) = DABS(ABD(MN,LP,K))
                   LL(MN)=LP
                   IPVTK=LP
52            CONTINUE
53       CONTINUE
54       CONTINUE
         KMM= K-M
         IPVTK=IPVTK+KMM
         JU = MIN0(MAX0(JU,MU+IPVTK),N)
         DO 56 MN=1,NMA
            L=LL(MN)
            IPVT(MN,K) = L + KMM
C
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
            IF (ABD(MN,L,K) .EQ. 0.0D0) GO TO 200
C
C           INTERCHANGE IF NECESSARY
C
            IF (L .EQ. M) GO TO 56
               MM=M
                DO 55 KM=K,JU
                   T = ABD(MN,L,KM)
                   ABD(MN,L,KM) = ABD(MN,MM,KM)
                   ABD(MN,MM,KM) = T
                   MM=MM-1
                   L=L-1
55             CONTINUE
56           CONTINUE
C
C           COMPUTE MULTIPLIERS
C           INNER LOOP IS OVER NMA
C
            DO 59 MN=1,NMA
59            TM(MN)=-1.0D0/ABD(MN,M,K)
            MPLM=M+LM
            DO 61 LPML=MP1,MPLM
               DO 62 MN=1,NMA
                  ABD(MN,LPML,K)=ABD(MN,LPML,K)*TM(MN)
62             CONTINUE
61         CONTINUE
C
C           ROW ELIMINATION WITH COLUMN INDEXING
C           S2XPY IS A BLA THAT SHOULD HAVE BEEN CREATED
C           IT IS ESSENTIAL FOR PARALLEL COMPUTATION
C           ESSENTIALLY IT SAYS Y(I)=Y(I)+A(I)*X(I)
C           I.E. IT IS LIKE SAXPY BUT A DOES NOT HAVE TO BE
C           A SCALAR.
C           WITH S2XPY THE INNER LOOP IS OVER THE NUMBER OF BANDED
C           SYSTEMS RATHER THAN THE BANDWIDTHS OF THE SYSTEMS
C
            MM = M
            IF (JU .LT. KP1) GO TO 90
            DO 80 J = KP1, JU
               MM = MM - 1
               MMP=MM
               DO 71 LPML=MP1,MPLM
                  MMP=MMP+1
                 CALL S2XPY(NMA,ABD(1,MMP,J),ABD(1,MM,J),ABD(1,LPML,K))
71             CONTINUE
   80       CONTINUE
   90       CONTINUE
  120 CONTINUE
130    CONTINUE
      DO 145 J=1,NMA
         IF(ABD(J,M,N).EQ.0.0)GO TO 200
145   CONTINUE
      RETURN
200   INFO=K
      RETURN
      END
      SUBROUTINE G2BSLM(ABD,LDA,N,ML,MU,IPVT,B,JOB,LNMA,NMA,TM)
      INTEGER LDA,N,ML,MU,IPVT(LNMA,N),JOB
      DOUBLE PRECISION ABD(LNMA,LDA,1),B(LNMA,1),TM(NMA)
C
C     G2BSLM SOLVES NMA DOUBLE PRECISION BAND SYSTEMS
C     A * X = B
C     USING THE FACTORS COMPUTED BY  G2BFAM.
C
C     THIS SUBROUTINE IS BASED ON LINPACK'S BAND SOLVER
C     BUT ALL THE SYSTEMS ARE HANDLED SIMULTANEOUSLY AND
C     INNER LOOPS ARE OVER THE NUMBER OF SYSTEMS RATHER THAN
C     THEIR BANDWIDTHS
C
C
C     ON ENTRY
C
C        ABD     DOUBLE PRECISION(LNMA,LDA, N)
C                THE OUTPUT FROM  G2BFAM.
C
C        LDA     INTEGER
C                THE SECOND DIMENSION OF THE ARRAY  ABD .
C
C        N       INTEGER
C                THE ORDER OF THE ORIGINAL MATRICES.
C
C        ML      INTEGER
C                NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL IN THE
C                SYSTEMS.
C
C        MU      INTEGER
C                NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
C
C        IPVT    INTEGER(LNMA,N)
C                THE PIVOT MATRIX FROM  G2BFAM.
C
C        B       DOUBLE PRECISION(LNMA,N)
C                THE RIGHT HAND SIDE FOR EACH SYSTEM
C
C        JOB     INTEGER
C                = 0         TO SOLVE  A*X = B ,
C
C        LNMA    THE LEADING DIMENSION OF THE ABD MATRIX
C                MUST BE AT LEAST NMA
C
C        NMA     NUMBER OF SYSTEMS TO BE SOLVED
C     ON RETURN
C
C        B       THE SOLUTION VECTORS  X .
C
C     SCRATCH SPACE
C
C        T       DOUBLE PRECISION (NMA)
C     ERROR CONDITION
C
C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
C        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
C        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
C        SETTING OF LDA.
C
C     SUBROUTINES AND FUNCTIONS
C
C     FORTRAN MIN0
C     CALLS S2XPY
C
C     INTERNAL VARIABLES
C
      INTEGER K,KB,L,LA,LB,LM,M,NM1
C
      M = MU + ML + 1
      MM1=M-1
      MP1 = M+ 1
      NM1 = N - 1
      IF (JOB .NE. 0) GO TO 50
C
C        JOB = 0 , SOLVE  A * X = B
C        FIRST SOLVE L*Y = B
C
         IF (ML .EQ. 0) GO TO 30
         IF (NM1 .LT. 1) GO TO 30
            DO 20 K = 1, NM1
               LM = MIN0(ML,N-K)
               DO 5 MN=1,NMA
               L = IPVT(MN,K)
               TM(MN) = B(MN,L)
C              IF (L .EQ. K) GO TO 5
                  B(MN,L) = B(MN,K)
                  B(MN,K)=TM(MN)
    5              CONTINUE
               MPLM=M+LM
               KP=K
C
C AGAIN S2XPY IS LIKE THE BLA SAXPY BUT A MAY BE
C A VECTOR RATHER THAN A SCALAR. USING S2XPY MEANS
C THAT THE INNER LOOP IS ON THE NUMBER OF SYSTEMS
C BEING SOLVED SIMULTANEOUSLY RATHER THAN ON THE
C BAND WIDTH
C
               DO 11 MP=MP1,MPLM
                  KP=KP+1
                  CALL S2XPY(NMA,B(1,KP),TM,ABD(1,MP,K))
11            CONTINUE
   20       CONTINUE
   30    CONTINUE
C
C        NOW SOLVE  U*X = Y
C
         DO 40 KB = 1, N
            K = N + 1 - KB
            DO 31 MN=1,NMA
            B(MN,K) = B(MN,K)/ABD(MN,M,K)
            TM(MN)=-B(MN,K)
31           CONTINUE
            LM = MIN0(K,M) - 1
            LA = M - LM
            LB = K - LM
            IF (MM1.LT.LA) GO TO 40
C
C BY USING S2XPY WE RUN THE INNER LOOP OVER THE
C NUMBER OF BANDED SYSTEMS RATHER THAN THE BAND WIDTH
C ON THE CRAY THE SPEEDUP IS SO IMPRESSIVE THAT IT
C IS IMPERATIVE THAT A BLA IS CREATED THAT DOES THE
C FUNCTION OF S2XPY
C
            DO 35 LAP=LA,MM1
               CALL S2XPY(NMA,B(1,LB),TM,ABD(1,LAP,K))
               LB=LB+1
35         CONTINUE
   40    CONTINUE
      GO TO 100
   50 CONTINUE
C
C
  100 CONTINUE
      RETURN
      END
      SUBROUTINE S2XPY(N,Y,A,X)
      INTEGER          N
      DOUBLE PRECISION Y(N), A(N), X(N)
C
C    THIS SUBROUTINE IS A RELATIVE OF THE BLA SAXPY
C    BUT HERE A IS PERMITTED TO BE A VECTOR
C    IT IS ESSENTIAL FOR PARALLEL COMPUTATION
C
C
      DO 10 I=1,N
   10    Y(I) = Y(I) + A(I)*X(I)
C
      RETURN
      END
*
      SUBROUTINE B2EVAL(K, M, X, IX, N, Y, IY, V, MP, XPTS, NP, YPTS,
     1  IDERIV, IS, S, SX, SY)
      INTEGER           K, M, N, MP, NP, IDERIV, IS
      DOUBLE PRECISION  X(1), Y(1), V(M,N), XPTS(MP), YPTS(NP),
     1                  S(IS,NP),SX(IS,NP),SY(IS,NP)
      DOUBLE PRECISION DSTAK(500)
      INTEGER ISTAK(500)
      COMMON /CSTAK/DSTAK
      EQUIVALENCE (DSTAK(1),ISTAK(1))
C
C THIS SUBROUTINE COMPUTES THE VALUE OF THE B SPLINE AND
C PERHAPS ITS DERIVATIVES DEFINED BY THE TENSOR PRODUCT
C MESH X AND Y AND THE B SPLINE COEFICIENTS IN V AT THE
C POINTS XPTS(I),YPTS(J),I=1,...MP,J=1,...NP, I.E
C
C  S(F,G)= SUM SUM V(I,J) B (XPTS(F)) B (YPTS(G))
C           I   J          I           J
C
C   INPUT PARAMETERS
C
C      K      IS THE ORDER OF THE B SPLINE
C      M      IS THE DIMENSION OF BI(X)
C      X      IS THE KNOT SEQUENCE IN THE X DIRECTION
C      IX     IF 0, X IS A NORMAL MESH AND OTHERWISE A PERIODIC MESH
C      N      IS THE DIMENSION OF BJ(X)
C      Y      IS THE KNOT SEQUENCE IN THE Y DIRECTION
C      IY      IF 0, Y IS A NORMAL MESH AND OTHERWISE A PERIODIC ONE
C      V      IS THE M BY N ARRAY OF COEFFICIENTS
C      MP     THE NUMBER OF EVALUATION POINTS IN THE X DIRECTION
C      XPTS   THE MP EVALUATION POINTS IN THE X DIRECTION
C      NP     THE NUMBER OF EVALUATION POINTS IN THE Y DIRECTION
C      YPTS   THE NP EVALUATION POINTS IN THE Y DIRECTION
C      IS     ROW DIMENSION OF S,SX,AND SY
C      IDERIV INTEGER IF 0, FUNCTION VALUES ARE COMPUTED, IF 1, FUNCTION
C             AND DERIVATIVE ARE COMPUTED AND IF 2, ONLY DERIVATIVES
C             ARE COMPUTED
C
C    OUTPUT PARAMETERS
C
C      S      THE COMPUTED VALUES OF THE SPLINE ON THE TENSOR PRODUCT
C             MESH AT THE POINTS SPECIFIED BY XPTS AND YPTS.
C      SX     THE COMPUTED VALUES OF THE DERIVATIVE OF THE B SPLINE
C             WITH RESPECT TO X ON THE TENSOR PRODUCT MESH SPECIFIED
C             BY XPTS AND YPTS
C      SY     THE COMPUTED VALUES OF THE DERIVATIVE OF THE B SPLINE
C             WITH RESPECT TO Y ON THE TENSOR PRODUCT MESH SPECIFIED
C             BY XPTS AND YPTS
C
C CALLS FROM PORT: SETERR, ISTLGT, LEAVE
C OTHER CALLS: B6EVAL, N2MESH, S2HOVV
C
C       CHECK THE INPUT PARAMETERS
C
      IF (K.LT.2)CALL SETERR(15H B2EVAL- K.LT.2,15,1,2)
      IF (M.LT.1)CALL SETERR(15H B2EVAL- M.LT.1,15,2,2)
      IF (N.LT.1)CALL SETERR(15H B2EVAL- N.LT.1,15,3,2)
      IF (MP.LT.1)CALL SETERR(16H B2EVAL- MP.LT.1,16,4,2)
      IF (NP.LT.1)CALL SETERR(15H B2EVAL-NP.LT.1,1,15,5,2)
      IF (IDERIV.LT.0.OR.IDERIV.GT.2) CALL
     1SETERR(24H B2EVAL-INCORRECT IDERIV,24,6,2)
      IF(IS.LT.MP)CALL SETERR(17H B2EVAL- IS.LT.MP,17,7,2)
C
C GET SCRATCH SPACE
C
      CALL ENTER(0)
      ILY=ISTKGT(NP,2)
      IBY=ISTKGT(K*2*NP,4)
      IBX=ISTKGT(K*2,4)
      IF (IX.NE.0.OR.IY.NE.0) GO TO 20
C
C FIND THE VALUES AND/OR DERIVATIVES WHEN THERE IS NOT A PERIODIC
C MESH
C
      CALL B6EVAL(K,M,X,N,Y,V,MP,XPTS,NP,
     1 YPTS,IDERIV,IS,S,SX,SY,ISTAK(ILY),DSTAK(IBY),DSTAK(IBX))
      GO TO 40
20    CONTINUE
C
C THERE IS A PERIOD IC MESH. THUS EXPAND THE MESH SOMEWHAT
C DUPLICATE THE B SPLINE COEFFICIENTS AT THE END POINTS AND
C CALL THE UNDERLYING ROUTINE WHICH DOES NOT KNOW THAT THE
C MESH MIGHT BE PERIODIC
C
      NUMX=M+2*K
      NUMY=N+2*K
      IXM=ISTKGT(NUMX, 4)
      IYN=ISTKGT(NUMY, 4)
      IF (IX.NE.0)CALL N2MESH(M,K,X,DSTAK(IXM),IPT,.FALSE.)
      IF (IY.NE.0) CALL N2MESH(N,K,Y,DSTAK(IYN),IPT,.FALSE.)
      NE=M
      IF (IX.NE.0)NE=M+K-1
      NEN=N
      IF (IY.NE.0)NEN=N+K-1
      IV=ISTKGT(NE*NEN, 4)
      CALL S2HOVV(NE,M,N,K,V,DSTAK(IV),IX,IY)
      IF (IX.EQ.0) CALL DCOPY(NUMX-K,X,1,DSTAK(IXM),1)
      IF (IY.EQ.0)CALL DCOPY(NUMY-K,Y,1,DSTAK(IYN),1)
      CALL B6EVAL(K,NE,DSTAK(IXM),NEN,DSTAK(IYN),DSTAK(IV),MP,XPTS,
     1 NP,YPTS,IDERIV,IS,S,SX,SY,ISTAK(ILY),DSTAK(IBY),DSTAK(IBX))
40    CALL LEAVE
      RETURN
      END
      SUBROUTINE B6EVAL(K, M, X, N, Y, V, MP, XPTS, NP, YPTS,
     1  IDERIV, IS, S, SX, SY, LY,BY,BX)
C
C THIS IS THE WORK HORSE OF THE B SPLINE EVALUATOR. ON ENTRY
C PARAMETERS K THROUGH IS HAVE THE SAME MEANING AS THEY HAVE
C ON ENTRY TO B2EVAL. ON EXIT, S,SX,AND SY HAVE THE
C SAME MEANING AS THEY HAVE ON LEAVING B2EVAL.
C
C   WORK SPACE
C     LY        INTEGER VECTOR OF LENGTH NP POINTING OUT
C               WHICH INTERVAL IN THE Y MESH EACH POINT
C               YPTS(J),J=1,...NP RESIDES
C     BY       DOUBLE PRECISION VECTOR LENGTH K X 2 X NP WHICH
C              WILL CONTAIN THE NONZERO B -SPLINES AND Y
C              DERIVATIVE AT THE YPTS
C     BX       DOUBLE PRECISION VECTOR OF LENGTH K X 2
C              WHICH WILL CONTAIN THE NONZERO X SPLINES
C              AT A POINT .
C
      INTEGER           K, M, N, MP, NP, IDERIV, IS
      INTEGER          I, J, IX, JY, LX, LXK, LYK, I2NTRV
C
      DOUBLE PRECISION  X(1), Y(1), V(M,N), XPTS(MP), YPTS(NP),
     1                  S(IS,NP),SX(IS,NP),SY(IS,NP)
      DOUBLE PRECISION BY(K,2,NP),BX(K,2)
      INTEGER LY(K)
      DOUBLE PRECISION SUM, SUMX,SUMY
C
C CALLS B2SVD1, I2NTRV
C
C
C     COMPUTE THE B-SPLINES ON THE Y MESH
C
      DO 10 JY=1,NP
         LY(JY) = I2NTRV(Y, N, YPTS(JY))
   10    CALL B2SVD1(Y, K, YPTS(JY), LY(JY),  BY(1,1,JY))
C
C     COMPUTE UC.
C
      DO 80 IX=1,MP
C
C FOR EACH POINT IN XPTS FIND ITS INTERVAL AND NONZERO B SPLINES
C
         LX  = I2NTRV(X, M, XPTS(IX))
         LXK = LX - K
         CALL B2SVD1(X, K, XPTS(IX), LX, BX)
         DO 30 JY=1,NP
            LYK = LY(JY) - K
            IF (IDERIV.GT.1) GO TO 25
C
C COMPUTE THE B SPLINE AT XPTS(IX),YPTS(JY)
C
            SUM = 0.0D0
            DO 20 I=1,K
               LXKPI=LXK+I
               DO 20 J=1,K
                  LYKPJ=LYK+J
   20             SUM = SUM + V(LXKPI,LYKPJ)*BX(I,1)*BY(J,1,JY)
            S(IX,JY) = SUM
   25       CONTINUE
            IF (IDERIV.LT.1) GO TO 30
C
C COMPUTE THE B SPLINE DERIVATIVES AT XPTS(IX),JPTS(JY)
C
            SUMX=0.0
            SUMY=0.0
            DO 27 I=1,K
               LXKPI=LXK+I
               DO 27 J=1,K
                  LYKPJ=LYK+J
                  SUMX=SUMX+V(LXKPI,LYKPJ)*BX(I,2)*BY(J,1,JY)
                  SUMY=SUMY+V(LXKPI,LYKPJ)*BX(I,1)*BY(J,2,JY)
  27        CONTINUE
            SX(IX,JY)=SUMX
            SY(IX,JY)=SUMY
  30     CONTINUE
  80  CONTINUE
C
      RETURN
      END
      INTEGER FUNCTION I2NTRV(T, N, X)
      INTEGER N
      DOUBLE PRECISION  T(N), X
C
C     I2NTRV IS DETERMINED SO THAT
C        IF           X .LT. T(1)   THEN I2NTRV = 0
C        IF T(I) .LE. X .LT. T(I+1) THEN I2NTRV = I
C        IF T(N) .LE. X             THEN I2NTRV = N
C
C     INPUT
C        T - A MONOTONE INCREASING MESH
C        N - THE NUMBER OF MESH POINTS
C        X - THE POINT TO BE BRACKETED
C
      INTEGER LO, HI, STEP, MID
      DATA    LO/1/
C
      IF (N .LE. 1)       CALL SETERR(17HI2NTRV - N .LE. 1, 17, 1, 2)
      IF (T(N) .LE. T(1)) CALL SETERR(23HI2NTRV - T(N) .LE. T(1),
     1                                 23, 2, 2)
      STEP = 1
      LO = MIN0(LO, N)
      IF (T(LO) .LE. X) GOTO 200
C
C     (X .LT. T(LO)) THEREFORE DOUBLE DOWN
  100 HI = LO
      LO = MAX0(1, HI-STEP)
      IF (T(LO) .LE. X) GOTO 300
      IF (LO .NE. 1) GOTO 110
         I2NTRV = 0
         RETURN
  110 STEP = 2*STEP
      GOTO 100
C
C     (T(LO) .LE. X) THEREFORE DOUBLE UP
  200 HI = MIN0(N, LO+STEP)
      IF (X .LT. T(HI)) GOTO 300
      IF (HI .NE. N) GOTO 210
         I2NTRV = N
         RETURN
  210 STEP = 2*STEP
      LO   = HI
      GOTO 200
C
C     T(LO) .LE. X .LT. T(HI)  THEREFORE BISECT
  300 MID = (LO + HI)/2
      IF (MID .NE. LO) GOTO 310
         I2NTRV = LO
         RETURN
  310 IF (T(MID) .LE. X) GOTO 320
         HI = MID
      GOTO 300
  320    LO = MID
      GOTO 300
C
      END
          SUBROUTINE S2HOVV(NE,M,N,K,V,NEWV,IX,IY)
C
C FOR A PROBLEM WITH A PERIODIC MESH
C THIS SUBROUTINE ADJUSTS THE MATRIX OF COEFFICENTS V
C FOR THE B SPLINES SO THAT THE EVALUATOR DOES NOT HAVE TO KNOW
C THAT THE MESH HAS BEEN ENLARGED FOR A PERIODIC CASE
C
C INPUT PARAMETERS
C
C NE    ROW DIMENSION OF THE NEW MATRIX OF CEOFFICIENTS NEWV
C M     ORIGINAL MATRIX V WAS N X M
C N     ORIGINAL MATRIX V WAS V X M
C V     ORIGINAL COEFFICIENT MATRIX
C IX    IF IX IS NONZERO, THE PERIODICITY IS IN THE X MESH
C IY    IF IY IS NONZERO, THE PERIODICITY IS IN THE Y MESH
C
C OUTPUT PARAMETERS
C
C NEWV  THE NEW COEFFICIENT MATRIX WHICH IF IX IS NONZERO
C       HAS THE LAST K-1 ROWS OF V PLACED IN ITS FIRST
C       K-1 ROWS AND THE REST OF THE V MATRIX SHOVED OVER
C       K-1 ROWS. IF IY IS NONZERO, THE CORRESPONDING OPERATIONS
C       ARE DONE TO THE COLUMNS OF V.
C
          DOUBLE PRECISION V(M,N),NEWV(NE,N)
          K1=K-1
          KI=0
          KJ=0
          IF(IX.NE.0)KI=K1
C
C PLACE V IN NEWV WITH ROW DISPLACEMENT KI AND
C COLUMN DISPLACEMENT KJ
C
          IF(IY.NE.0)KJ=K1
          DO 20 I=1,M
             I2=I+KI
             DO 10 J=1,N
                J2=J+KJ
                NEWV(I2,J2)=V(I,J)
10           CONTINUE
20        CONTINUE
          IF (IX.EQ.0) GO TO 50
C
C FOR PERIODICITY IN X DUPLICATE THE LAST K-1 ROWS OF V
C IN NEWV
C
          DO 40 I=1,K1
             I2=M-K1+I
             DO 30 J=1,N
               JPKJ=J+KJ
               NEWV(I,JPKJ)=V(I2,J)
30           CONTINUE
40       CONTINUE
         IF (IY.EQ.0)RETURN
C
C FOR PERIODICITY IN Y DUPLICATE THE LAST K-1 COLUMNS OF V
C IN NEWV
C
50       DO 70 I=1,K1
            I2=N+I
            DO 60 J=1,NE
               NEWV(J,I)=NEWV(J,I2)
60          CONTINUE
70       CONTINUE
         RETURN
         END
C-----TOMS SUB ( SUBSIDIARY ROUTINES FROM LINPACK, EISPACK, PPPACK)
      SUBROUTINE DPBFA(ABD,LDA,N,M,INFO)
      INTEGER LDA,N,M,INFO
      DOUBLE PRECISION ABD(LDA,1)
C
C     DPBFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE
C     MATRIX STORED IN BAND FORM.
C
C     DPBFA IS USUALLY CALLED BY DPBCO, BUT IT CAN BE CALLED
C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
C
C     ON ENTRY
C
C        ABD     DOUBLE PRECISION(LDA, N)
C                THE MATRIX TO BE FACTORED.  THE COLUMNS OF THE UPPER
C                TRIANGLE ARE STORED IN THE COLUMNS OF ABD AND THE
C                DIAGONALS OF THE UPPER TRIANGLE ARE STORED IN THE
C                ROWS OF ABD .  SEE THE COMMENTS BELOW FOR DETAILS.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  ABD .
C                LDA MUST BE .GE. M + 1 .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C        M       INTEGER
C                THE NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
C                0 .LE. M .LT. N .
C
C     ON RETURN
C
C        ABD     AN UPPER TRIANGULAR MATRIX  R , STORED IN BAND
C                FORM, SO THAT  A = TRANS(R)*R .
C
C        INFO    INTEGER
C                = 0  FOR NORMAL RETURN.
C                = K  IF THE LEADING MINOR OF ORDER  K  IS NOT
C                     POSITIVE DEFINITE.
C
C     BAND STORAGE
C
C           IF  A  IS A SYMMETRIC POSITIVE DEFINITE BAND MATRIX,
C           THE FOLLOWING PROGRAM SEGMENT WILL SET UP THE INPUT.
C
C                   M = (BAND WIDTH ABOVE DIAGONAL)
C                   DO 20 J = 1, N
C                      I1 = MAX0(1, J-M)
C                      DO 10 I = I1, J
C                         K = I-J+M+1
C                         ABD(K,J) = A(I,J)
C                10    CONTINUE
C                20 CONTINUE
C
C     LINPACK.  THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS DDOT
C     FORTRAN MAX0,DSQRT
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION DDOT,T
      DOUBLE PRECISION S
      INTEGER IK,J,JK,K,MU
C     BEGIN BLOCK WITH ...EXITS TO 40
C
C
         DO 30 J = 1, N
            INFO = J
            S = 0.0D0
            IK = M + 1
            JK = MAX0(J-M,1)
            MU = MAX0(M+2-J,1)
            IF (M .LT. MU) GO TO 20
            DO 10 K = MU, M
               T = ABD(K,J) - DDOT(K-MU,ABD(IK,JK),1,ABD(MU,J),1)
               T = T/ABD(M+1,JK)
               ABD(K,J) = T
               S = S + T*T
               IK = IK - 1
               JK = JK + 1
   10       CONTINUE
   20       CONTINUE
            S = ABD(M+1,J) - S
C     ......EXIT
            IF (S .LE. 0.0D0) GO TO 40
            ABD(M+1,J) = DSQRT(S)
   30    CONTINUE
         INFO = 0
   40 CONTINUE
      RETURN
      END
*
      SUBROUTINE DPBSL(ABD,LDA,N,M,B)
      INTEGER LDA,N,M
      DOUBLE PRECISION ABD(LDA,1),B(1)
C
C     DPBSL SOLVES THE DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE
C     BAND SYSTEM  A*X = B
C     USING THE FACTORS COMPUTED BY DPBCO OR DPBFA.
C
C     ON ENTRY
C
C        ABD     DOUBLE PRECISION(LDA, N)
C                THE OUTPUT FROM DPBCO OR DPBFA.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  ABD .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C        M       INTEGER
C                THE NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
C
C        B       DOUBLE PRECISION(N)
C                THE RIGHT HAND SIDE VECTOR.
C
C     ON RETURN
C
C        B       THE SOLUTION VECTOR  X .
C
C     ERROR CONDITION
C
C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
C        A ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES
C        SINGULARITY BUT IT IS USUALLY CAUSED BY IMPROPER SUBROUTINE
C        ARGUMENTS.  IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED
C        CORRECTLY AND  INFO .EQ. 0 .
C
C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
C     WITH  P  COLUMNS
C           CALL DPBCO(ABD,LDA,N,RCOND,Z,INFO)
C           IF (RCOND IS TOO SMALL .OR. INFO .NE. 0) GO TO ...
C           DO 10 J = 1, P
C              CALL DPBSL(ABD,LDA,N,C(1,J))
C        10 CONTINUE
C
C     LINPACK.  THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DDOT
C     FORTRAN MIN0
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION DDOT,T
      INTEGER K,KB,LA,LB,LM
C
C     SOLVE TRANS(R)*Y = B
C
      DO 10 K = 1, N
         LM = MIN0(K-1,M)
         LA = M + 1 - LM
         LB = K - LM
         T = DDOT(LM,ABD(LA,K),1,B(LB),1)
         B(K) = (B(K) - T)/ABD(M+1,K)
   10 CONTINUE
C
C     SOLVE R*X = Y
C
      DO 20 KB = 1, N
         K = N + 1 - KB
         LM = MIN0(K-1,M)
         LA = M + 1 - LM
         LB = K - LM
         B(K) = B(K)/ABD(M+1,K)
         T = -B(K)
         CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1)
   20 CONTINUE
      RETURN
      END
      SUBROUTINE TRED2(NM,N,A,D,E,Z)
C
      INTEGER I,J,K,L,N,II,NM,JP1
      DOUBLE PRECISION A(NM,N),D(N),E(N),Z(NM,N)
      DOUBLE PRECISION F,G,H,HH,SCALE
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A
C     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C     ON OUTPUT
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
C
C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
C          PRODUCED IN THE REDUCTION.
C
C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED APRIL 1983.
C
C     ------------------------------------------------------------------
C
      DO 100 I = 1, N
C
         DO 80 J = I, N
   80    Z(J,I) = A(J,I)
C
         D(I) = A(N,I)
  100 CONTINUE
C
      IF (N .EQ. 1) GO TO 510
C     .......... FOR I=N STEP -1 UNTIL 2 DO -- ..........
      DO 300 II = 2, N
         I = N + 2 - II
         L = I - 1
         H = 0.0D0
         SCALE = 0.0D0
         IF (L .LT. 2) GO TO 130
C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
         DO 120 K = 1, L
  120    SCALE = SCALE + DABS(D(K))
C
         IF (SCALE .NE. 0.0D0) GO TO 140
  130    E(I) = D(L)
C
         DO 135 J = 1, L
            D(J) = Z(L,J)
            Z(I,J) = 0.0D0
            Z(J,I) = 0.0D0
  135    CONTINUE
C
         GO TO 290
C
  140    DO 150 K = 1, L
            D(K) = D(K) / SCALE
            H = H + D(K) * D(K)
  150    CONTINUE
C
         F = D(L)
         G = -DSIGN(DSQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         D(L) = F - G
C     .......... FORM A*U ..........
         DO 170 J = 1, L
  170    E(J) = 0.0D0
C
         DO 240 J = 1, L
            F = D(J)
            Z(J,I) = F
            G = E(J) + Z(J,J) * F
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
               G = G + Z(K,J) * D(K)
               E(K) = E(K) + Z(K,J) * F
  200       CONTINUE
C
  220       E(J) = G
  240    CONTINUE
C     .......... FORM P ..........
         F = 0.0D0
C
         DO 245 J = 1, L
            E(J) = E(J) / H
            F = F + E(J) * D(J)
  245    CONTINUE
C
         HH = F / (H + H)
C     .......... FORM Q ..........
         DO 250 J = 1, L
  250    E(J) = E(J) - HH * D(J)
C     .......... FORM REDUCED A ..........
         DO 280 J = 1, L
            F = D(J)
            G = E(J)
C
            DO 260 K = J, L
  260       Z(K,J) = Z(K,J) - F * E(K) - G * D(K)
C
            D(J) = Z(L,J)
            Z(I,J) = 0.0D0
  280    CONTINUE
C
  290    D(I) = H
  300 CONTINUE
C     .......... ACCUMULATION OF TRANSFORMATION MATRICES ..........
      DO 500 I = 2, N
         L = I - 1
         Z(N,L) = Z(L,L)
         Z(L,L) = 1.0D0
         H = D(I)
         IF (H .EQ. 0.0D0) GO TO 380
C
         DO 330 K = 1, L
  330    D(K) = Z(K,I) / H
C
         DO 360 J = 1, L
            G = 0.0D0
C
            DO 340 K = 1, L
  340       G = G + Z(K,I) * Z(K,J)
C
            DO 360 K = 1, L
               Z(K,J) = Z(K,J) - G * D(K)
  360    CONTINUE
C
  380    DO 400 K = 1, L
  400    Z(K,I) = 0.0D0
C
  500 CONTINUE
C
  510 DO 520 I = 1, N
         D(I) = Z(N,I)
         Z(N,I) = 0.0D0
  520 CONTINUE
C
      Z(N,N) = 1.0D0
      E(1) = 0.0D0
      RETURN
      END
      SUBROUTINE BSPLVB ( T, JHIGH, INDEX, X, LEFT, BIATX )
CALCULATES THE VALUE OF ALL POSSIBLY NONZERO B-SPLINES AT  X  OF ORDER
C
C               JOUT  =  MAX( JHIGH , (J+1)*(INDEX-1) )
C
C  WITH KNOT SEQUENCE  T .
C
C******  I N P U T  ******
C  T.....KNOT SEQUENCE, OF LENGTH  LEFT + JOUT  , ASSUMED TO BE NONDE-
C        CREASING.  A S S U M P T I O N . . . .
C                       T(LEFT)  .LT.  T(LEFT + 1)   .
C   D I VI S I O N  B Y  Z E R O   WILL RESULT IF  T(LEFT) = T(LEFT+1)
C  JHIGH,
C  INDEX.....INTEGERS WHICH DETERMINE THE ORDER  JOUT = MAX(JHIGH,
C        (J+1)*(INDEX-1))   OF THE B-SPLINES WHOSE VALUES AT  X  ARE TO
C        BE RETURNED.  INDEX  IS USED TO AVOID RECALCULATIONS WHEN SEVE-
C        RAL COLUMNS OF THE TRIANGULAR ARRAY OF B-SPLINE VALUES ARE NEE-
C        DED (E.G., IN  BVALUE  OR IN  BSPLVD ). PRECISELY,
C                     IF  INDEX = 1 ,
C        THE CALCULATION STARTS FROM SCRATCH AND THE ENTIRE TRIANGULAR
C        ARRAY OF B-SPLINE VALUES OF ORDERS 1,2,...,JHIGH  IS GENERATED
C        ORDER BY ORDER, I.E., COLUMN BY COLUMN.
C                     IF  INDEX = 2 ,
C        ONLY THE B-SPLINE VALUES OF ORDER  J+1, J+2,..., JOUT  ARE GE-
C        NERATED, THE ASSUMPTION BEING THAT  BIATX, J, DELTAL, DELTAR
C        ARE, ON ENTRY, AS THEY WERE ON EXIT AT THE PREVIOUS CALL.
C           IN PARTICULAR, IF  JHIGH = 0, THEN  JOUT = J+1, I.E., JUST
C        THE NEXT COLUMN OF B-SPLINE VALUES IS GENERATED.
C
C  W A R N I N G . . .  THE RESTRICTION  JOUT .LE. JMAX (= 20) IS IM-
C        POSED ARBITRARILY BY THE DIMENSION STATEMENT FOR  DELTAL  AND
C        DELTAR  BELOW, BUT IS  N O W H E R E   C H E C K E D   FOR .
C
C  X.....THE POINT AT WHICH THE B-SPLINES ARE TO BE EVALUATED.
C  LEFT.....AN INTEGER CHOSEN (USUALLY) SO THAT
C                  T(LEFT) .LE. X .LE. T(LEFT+1)  .
C
C******  O U T P U T  ******
C  BIATX.....ARRAY OF LENGTH  JOUT , WITH BIATX(I) CONTAINING THE VA-
C        LUE AT  X  OF THE POLYNOMIAL OF ORDER  JOUT  WHICH AGEES WITH
C        THE B-SPLINE  B(LEFT-JOUT+I,JOUT,T)  ON THE INTERVAL (T(LEFT),
C        T(LEFT+1))  .
C
C******  M E T H O D  *******
C  THE RECURRENCE RELATION
C
C                       X - T(I)              T(I+J+1) - X
C     B(I,J+1)(X)  =  -----------B(I,J)(X) + ---------------B(I+1,J)(X)
C                     T(I+J)-T(I)            T(I+J+1)-T(I+1)
C
C  IS USED (REPEATEDLY) TO GENERATE THE (J+1)-VECTOR  B(LEFT-J,J+1)(X),
C  ...,B(LEFT,J+1)(X)  FROM THE J-VECTOR  B(LEFT-J+1,J)(X),...,
C  B(LEFT,J)(X), STORING THE NEW VALUES IN BIATX  OVER THE OLD. THE
C  FACTS THAT
C            B(I,1) = 1  IF  T(I) .LE. X .LT. T(I+1)
C  AND THAT
C            B(I,J)(X) = 0  UNLESS  T(I) .LE. X .LT. T(I+J)
C  ARE USED. THE PARTICULAR ORGANIZATION OF THE CALCULATIONS FOLLOWS AL-
C  GORITHM  (8)  IN CHAPTER X OF THE TEXT.
C
C     PARAMETER JMAX = 20
      INTEGER INDEX,JHIGH,LEFT,   I,J,JP1
      DOUBLE PRECISION BIATX(JHIGH),T(1),X,   DELTAL(20),DELTAR(20),
     1                 SAVED,TERM
C     DIMENSION BIATX(JOUT), T(LEFT+JOUT)
CURRENT FORTRAN STANDARD MAKES IT IMPOSSIBLE TO SPECIFY THE LENGTH OF
C  T  AND OF  BIATX  PRECISELY WITHOUT THE INTRODUCTION OF OTHERWISE
C  SUPERFLUOUS ADDITIONAL ARGUMENTS.
      DATA J/1/
C     SAVE J,DELTAL,DELTAR (VALID IN FORTRAN 77)
C
                                        GO TO (10,20), INDEX
   10 J = 1
      BIATX(1) = 1.0D0
      IF (J .GE. JHIGH)                 GO TO 99
C
   20    JP1 = J + 1
         DELTAR(J) = T(LEFT+J) - X
         DELTAL(J) = X - T(LEFT+1-J)
         SAVED = 0.0D0
         DO 26 I=1,J
            TERM = BIATX(I)/(DELTAR(I) + DELTAL(JP1-I))
            BIATX(I) = SAVED + DELTAR(I)*TERM
  26        SAVED = DELTAL(JP1-I)*TERM
        BIATX(JP1) = SAVED
        J = JP1
        IF (J .LT. JHIGH)               GO TO 20
C
  99                                    RETURN
      END
      SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
C
C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DA
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF (DA .EQ. 0.0D0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DY(IY) + DA*DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,4)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DY(I) + DA*DX(I)
   30 CONTINUE
      IF( N .LT. 4 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
        DY(I) = DY(I) + DA*DX(I)
        DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
        DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
        DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
   50 CONTINUE
      RETURN
      END
      SUBROUTINE  DCOPY(N,DX,INCX,DY,INCY)
C
C     COPIES A VECTOR, X, TO A VECTOR, Y.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1)
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,7)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DX(I)
   30 CONTINUE
      IF( N .LT. 7 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        DY(I) = DX(I)
        DY(I + 1) = DX(I + 1)
        DY(I + 2) = DX(I + 2)
        DY(I + 3) = DX(I + 3)
        DY(I + 4) = DX(I + 4)
        DY(I + 5) = DX(I + 5)
        DY(I + 6) = DX(I + 6)
   50 CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DTEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      DDOT = 0.0D0
      DTEMP = 0.0D0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DTEMP = DTEMP + DX(IX)*DY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      DDOT = DTEMP
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DTEMP + DX(I)*DY(I)
   30 CONTINUE
      IF( N .LT. 5 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
     *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
   50 CONTINUE
   60 DDOT = DTEMP
      RETURN
      END
C-----TOMS PORT (SUBSIDIARY ROUTINES FROM PUBLIC PART OF PORT)
      SUBROUTINE ENTER(IRNEW)
C
C  THIS ROUTINE SAVES
C
C    1) THE CURRENT NUMBER OF OUTSTANDING STORAGE ALLOCATIONS, LOUT, AND
C    2) THE CURRENT RECOVERY LEVEL, LRECOV,
C
C  IN AN ENTER-BLOCK IN THE STACK.
C
C  IT ALSO SETS LRECOV = IRNEW IF IRNEW = 1 OR 2.
C  IF IRNEW = 0, THEN THE RECOVERY LEVEL IS NOT ALTERED.
C
C  SCRATCH SPACE ALLOCATED - 3 INTEGER WORDS ARE LEFT ON THE STACK.
C
C  ERROR STATES -
C
C    1 - MUST HAVE IRNEW = 0, 1 OR 2.
C
      COMMON /CSTAK/DSTACK
      DOUBLE PRECISION DSTACK(500)
      INTEGER ISTACK(1000)
      EQUIVALENCE (DSTACK(1),ISTACK(1))
      EQUIVALENCE (ISTACK(1),LOUT)
C
C/6S
C     IF (0.GT.IRNEW .OR. IRNEW.GT.2)
C    1  CALL SETERR(35HENTER - MUST HAVE IRNEW = 0, 1 OR 2,35,1,2)
C/7S
      IF (0.GT.IRNEW .OR. IRNEW.GT.2)
     1  CALL SETERR('ENTER - MUST HAVE IRNEW = 0, 1 OR 2',35,1,2)
C/
C
C  ALLOCATE SPACE FOR SAVING THE ABOVE 2 ITEMS
C  AND A BACK-POINTER FOR CHAINING THE ENTER-BLOCKS TOGETHER.
C
      INOW=ISTKGT(3,2)
C
C  SAVE THE CURRENT NUMBER OF OUTSTANDING ALLOCATIONS.
C
      ISTACK(INOW)=LOUT
C
C  SAVE THE CURRENT RECOVERY LEVEL.
C
      CALL ENTSRC(ISTACK(INOW+1),IRNEW)
C
C  SAVE A BACK-POINTER TO THE START OF THE PREVIOUS ENTER-BLOCK.
C
      ISTACK(INOW+2)=I8TSEL(INOW)
C
      RETURN
C
      END
      SUBROUTINE ENTSRC(IROLD,IRNEW)
C
C  THIS ROUTINE RETURNS IROLD = LRECOV AND SETS LRECOV = IRNEW.
C
C  IF THERE IS AN ACTIVE ERROR STATE, THE MESSAGE IS PRINTED
C  AND EXECUTION STOPS.
C
C  IRNEW = 0 LEAVES LRECOV UNCHANGED, WHILE
C  IRNEW = 1 GIVES RECOVERY AND
C  IRNEW = 2 TURNS RECOVERY OFF.
C
C  ERROR STATES -
C
C    1 - ILLEGAL VALUE OF IRNEW.
C    2 - CALLED WHILE IN AN ERROR STATE.
C
C/6S
C     IF (IRNEW.LT.0 .OR. IRNEW.GT.2)
C    1   CALL SETERR(31HENTSRC - ILLEGAL VALUE OF IRNEW,31,1,2)
C/7S
      IF (IRNEW.LT.0 .OR. IRNEW.GT.2)
     1   CALL SETERR('ENTSRC - ILLEGAL VALUE OF IRNEW',31,1,2)
C/
C
      IROLD=I8SAVE(2,IRNEW,IRNEW.NE.0)
C
C  IF HAVE AN ERROR STATE, STOP EXECUTION.
C
C/6S
C     IF (I8SAVE(1,0,.FALSE.) .NE. 0) CALL SETERR
C    1   (39HENTSRC - CALLED WHILE IN AN ERROR STATE,39,2,2)
C/7S
      IF (I8SAVE(1,0,.FALSE.) .NE. 0) CALL SETERR
     1   ('ENTSRC - CALLED WHILE IN AN ERROR STATE',39,2,2)
C/
C
      RETURN
C
      END
      INTEGER FUNCTION I1MACH(I)
C
C  I/O UNIT NUMBERS.
C
C    I1MACH( 1) = THE STANDARD INPUT UNIT.
C
C    I1MACH( 2) = THE STANDARD OUTPUT UNIT.
C
C    I1MACH( 3) = THE STANDARD PUNCH UNIT.
C
C    I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT.
C
C  WORDS.
C
C    I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT.
C
C    I1MACH( 6) = THE NUMBER OF CHARACTERS PER CHARACTER STORAGE UNIT.
C                 FOR FORTRAN 77, THIS IS ALWAYS 1.  FOR FORTRAN 66,
C                 CHARACTER STORAGE UNIT = INTEGER STORAGE UNIT.
C
C  INTEGERS.
C
C    ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM
C
C               SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
C
C               WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1.
C
C    I1MACH( 7) = A, THE BASE.
C
C    I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS.
C
C    I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE.
C
C  FLOATING-POINT NUMBERS.
C
C    ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT,
C    BASE-B FORM
C
C               SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C
C               WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T,
C               0 .LT. X(1), AND EMIN .LE. E .LE. EMAX.
C
C    I1MACH(10) = B, THE BASE.
C
C  SINGLE-PRECISION
C
C    I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
C
C    I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
C
C    I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
C
C  DOUBLE-PRECISION
C
C    I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
C
C    I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
C
C    I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
C
C  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
C  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
C  REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
C  I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
C  WITH THE LOCAL OPERATING SYSTEM.  FOR FORTRAN 77, YOU MAY WISH
C  TO ADJUST THE DATA STATEMENT SO IMACH(6) IS SET TO 1, AND
C  THEN TO COMMENT OUT THE EXECUTABLE TEST ON I .EQ. 6 BELOW.
C
C  FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), THE FIRST
C  SET OF CONSTANTS BELOW SHOULD BE APPROPRIATE, EXCEPT PERHAPS
C  FOR IMACH(1) - IMACH(4).
C
      INTEGER IMACH(16),OUTPUT,SANITY
C
      EQUIVALENCE (IMACH(4),OUTPUT)
C
C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
C     3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
C     PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300).
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -125 /
C      DATA IMACH(13) /  128 /
C      DATA IMACH(14) /   53 /
C      DATA IMACH(15) / -1021 /
C      DATA IMACH(16) /  1024 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  32 /
C      DATA IMACH( 6) /   4 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /  16 /
C      DATA IMACH(11) /   6 /
C      DATA IMACH(12) / -64 /
C      DATA IMACH(13) /  63 /
C      DATA IMACH(14) /  14 /
C      DATA IMACH(15) / -64 /
C      DATA IMACH(16) /  63 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C      DATA IMACH( 1) /    7 /
C      DATA IMACH( 2) /    2 /
C      DATA IMACH( 3) /    2 /
C      DATA IMACH( 4) /    2 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   33 /
C      DATA IMACH( 9) / Z1FFFFFFFF /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -256 /
C      DATA IMACH(13) /  255 /
C      DATA IMACH(14) /   60 /
C      DATA IMACH(15) / -256 /
C      DATA IMACH(16) /  255 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  48 /
C      DATA IMACH( 6) /   6 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  39 /
C      DATA IMACH( 9) / O0007777777777777 /
C      DATA IMACH(10) /   8 /
C      DATA IMACH(11) /  13 /
C      DATA IMACH(12) / -50 /
C      DATA IMACH(13) /  76 /
C      DATA IMACH(14) /  26 /
C      DATA IMACH(15) / -50 /
C      DATA IMACH(16) /  76 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  48 /
C      DATA IMACH( 6) /   6 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  39 /
C      DATA IMACH( 9) / O0007777777777777 /
C      DATA IMACH(10) /   8 /
C      DATA IMACH(11) /  13 /
C      DATA IMACH(12) / -50 /
C      DATA IMACH(13) /  76 /
C      DATA IMACH(14) /  26 /
C      DATA IMACH(15) / -32754 /
C      DATA IMACH(16) /  32780 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   60 /
C      DATA IMACH( 6) /   10 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   48 /
C      DATA IMACH( 9) / 00007777777777777777B /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   47 /
C      DATA IMACH(12) / -929 /
C      DATA IMACH(13) / 1070 /
C      DATA IMACH(14) /   94 /
C      DATA IMACH(15) / -929 /
C      DATA IMACH(16) / 1069 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   60 /
C      DATA IMACH( 6) /   10 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   48 /
C      DATA IMACH( 9) / O"00007777777777777777" /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   47 /
C      DATA IMACH(12) / -929 /
C      DATA IMACH(13) / 1070 /
C      DATA IMACH(14) /   94 /
C      DATA IMACH(15) / -929 /
C      DATA IMACH(16) / 1069 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR CONVEX C-1.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -128 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   53 /
C      DATA IMACH(15) /-1024 /
C      DATA IMACH(16) / 1023 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /   102 /
C      DATA IMACH( 4) /     6 /
C      DATA IMACH( 5) /    64 /
C      DATA IMACH( 6) /     8 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    63 /
C      DATA IMACH( 9) /  777777777777777777777B /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    47 /
C      DATA IMACH(12) / -8189 /
C      DATA IMACH(13) /  8190 /
C      DATA IMACH(14) /    94 /
C      DATA IMACH(15) / -8099 /
C      DATA IMACH(16) /  8190 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200.
C
C      DATA IMACH( 1) /   11 /
C      DATA IMACH( 2) /   12 /
C      DATA IMACH( 3) /    8 /
C      DATA IMACH( 4) /   10 /
C      DATA IMACH( 5) /   16 /
C      DATA IMACH( 6) /    2 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   15 /
C      DATA IMACH( 9) /32767 /
C      DATA IMACH(10) /   16 /
C      DATA IMACH(11) /    6 /
C      DATA IMACH(12) /  -64 /
C      DATA IMACH(13) /   63 /
C      DATA IMACH(14) /   14 /
C      DATA IMACH(15) /  -64 /
C      DATA IMACH(16) /   63 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7.
C
C      DATA IMACH( 1) /       5 /
C      DATA IMACH( 2) /       6 /
C      DATA IMACH( 3) /       0 /
C      DATA IMACH( 4) /       6 /
C      DATA IMACH( 5) /      24 /
C      DATA IMACH( 6) /       3 /
C      DATA IMACH( 7) /       2 /
C      DATA IMACH( 8) /      23 /
C      DATA IMACH( 9) / 8388607 /
C      DATA IMACH(10) /       2 /
C      DATA IMACH(11) /      23 /
C      DATA IMACH(12) /    -127 /
C      DATA IMACH(13) /     127 /
C      DATA IMACH(14) /      38 /
C      DATA IMACH(15) /    -127 /
C      DATA IMACH(16) /     127 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /   43 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / O377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   63 /
C      DATA IMACH(15) / -127 /
C      DATA IMACH(16) /  127 /, SANITY/987/
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  32 /
C      DATA IMACH( 6) /   4 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  31 /
C      DATA IMACH( 9) / Z7FFFFFFF /
C      DATA IMACH(10) /  16 /
C      DATA IMACH(11) /   6 /
C      DATA IMACH(12) / -64 /
C      DATA IMACH(13) /  63 /
C      DATA IMACH(14) /  14 /
C      DATA IMACH(15) / -64 /
C      DATA IMACH(16) /