C     ******
C
C     SUBDIV
C
C     **** VERSION OF AUGUST 15, 1983
C
C     SUBDIV IS A COLLECTION OF FOUR SUBROUTINES DEALING WITH
C     THE VARIOUS PIECES OF THE SUBDIVISION
C
C     FSTPCE SETS UP THE INITIAL PIECE
C
C     NEWPCE UPDATES THE PIECE GIVEN A LEAVING VERTEX OR SET OF
C     VERTICES
C
C     LSTPCE CALCULATES THE APPROXIMATE ZERO OF THE FUNCTION
C
C     PRTPCE PRINTS THE CURRENT PIECE
C
C     FIVE SUBDIVISIONS ARE POSSIBLE
C     IF VDLT.EQ.0 THE SUBDIVISION IS MODIFIED J1
C     IF VDLT.GE.1 THE SUBDIVISION IS VAN DER LAAN AND TALMAN'S
C                  FOR VDLT.EQ.1 THE TRIANGULATION ON TOP IS J1
C                  FOR VDLT.EQ.2 IT IS K'
C                  FOR VDLT.GE.3 IT IS K1
C     IF VDLT.LT.0 THE SUBDIVISION IS J-PRIME
C
C     EACH PIECE IS DETERMINED BY THE FOLLOWING:
C
C     VERT(.,1)     INTEGER          STARTING VERTEX
C
C     ITOP          INTEGER          NUMBER OF VERTICES ON TOP
C
C     PI(.)         INTEGER          MAPPING FROM 1,2,...,ITOP-1
C                                    INTO 1,2,...,N: COORDINATES
C                                    USED ON TOP IN ORDER
C
C     SIGN(.)       INTEGER          SIGN VECTOR: SIGN(PI(I)) IS
C                                    SENSE (+ OR -1) USED IN
C                                    COORDINATE PI(I)
C
C     FOR I=1,2,...,ITOP, THE VERTICES ON TOP ARE GIVEN BY (VERT(.,I),1)
C     WHERE
C
C     VERT(.,I+1) = VERT(.,I) + SIGN(PI(I))*E(PI(I))
C
C     WHERE E(K) IS THE KTH UNIT VECTOR IN R**N
C     IF VDLT.LT.0 THE INCREMENT IS DOUBLED FOR I.EQ.1 AND
C     REPLACED BY SIGN(PI(2))*E(PI(2))-SIGN(PI(1))*E(PI(1)) FOR I.EQ.2
C     IF VDLT.EQ.0 THE VERTICES ON THE BOTTOM ARE
C
C     (VERT(.,ITOP) + SUM(+ OR - E(K)), 0),
C
C     WHERE K RUNS OVER THE COORDINATES NOT USED ON TOP
C     IF VDLT.LT.0 THE SAME IS TRUE EXCEPT THAT IF ITOP.EQ.2
C     VERT(.,ITOP) SHOULD BE REPLACED WITH (VERT(.,1)+VERT(.,2))/2
C     IF VDLT.GE.1 VERT(K,I) IS ZERO FOR ALL SUCH K AND ALL I
C     IN THIS CASE THE VERTICES ON THE BOTTOM ARE AS ABOVE EXCEPT
C     THAT VERT(.,ITOP) IS REPLACED WITH S(VERT(.,ITOP)) WHERE
C     S(X(.)) IS THE (0,+1,-1) SIGN VECTOR OF VECTOR X
C
C     IF VDLT.LE.1 SIGN IS ARBITRARY AND ALL VERT(K,1) ARE EVEN INTEGERS
C     EXCEPT THAT IF VDLT.LT.0 SIGN(PI(1)) MUST BE POSITIVE
C     IF VDLT.EQ.2 SIGN=S(VERT(.,ITOP)) AND ALL VERT(K,1) ARE INTEGERS
C     IF VDLT.GE.3 SIGN(K)=1 ALL K AND ALL VERT(K,1) ARE INTEGERS
C
C
C     EACH POINT Y IN THE SUBDIVIDED R**N CORRESPONDS TO A
C     POINT X = XORIG + DELTA*Y IN THE DOMAIN OF THE FUNCTION F
C
C     XORIG(.)     DOUBLE PRECISION     BASE POINT
C
C     DELTA(.)     DOUBLE PRECISION     ELEMENTS OF DIAGONAL SCALING
C                                       MATRIX DELTA
C
C     ******
C
      SUBROUTINE FSTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,DELTA,
     *                  XORIG,IOUT,KOUT,KOUTSG,VTWT,
     *                  XNEW,IIN,KIN,KINSG)
      INTEGER N,NP1,PRTLEV,VDLT,VERT(N,NP1),ITOP,PI(N),SIGN(N)
      INTEGER IOUT,KOUT,KOUTSG,IIN,KIN,KINSG
      DOUBLE PRECISION DELTA(N),XORIG(N),VTWT(NP1),XNEW(N)
C
C     FSTPCE INITIALIZES THE PIECE OF THE SUBDIVISION AND
C     RETURNS IN XNEW THE FIRST POINT TO BE EVALUATED
C
C     SUBROUTINE USED -- PRTPCE
C
      INTEGER J
C
      IF ((N.LT.2).AND.(VDLT.LT.0)) VDLT=0
C
C     THERE IS ONE VERTEX ON TOP
C
      ITOP=1
C
C     THIS VERTEX IS THE ORIGIN. THE NEW POINT XNEW IS THE BASE POINT
C
      DO 10 J=1,N
         VERT(J,1)=0
         XNEW(J)=XORIG(J)
   10    CONTINUE
C
C     THE NEW POINT IS ON TOP AND THE FIRST VERTEX
C
      KIN=0
      IIN=1
C
C     IF DESIRED PRINT THE FIRST PIECE
C
      IF (PRTLEV.GE.5) CALL PRTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,
     *                             DELTA,XORIG,IOUT,KOUT,KOUTSG,VTWT,
     *                             XNEW,IIN,KIN,KINSG)
      RETURN
C
C     LAST CARD OF FSTPCE FOLLOWS
C
      END
      SUBROUTINE NEWPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,IPOINT,
     *                  DELTA,XORIG,IOUT,KOUT,KOUTSG,VTWT,
     *                  XNEW,IIN,KIN,KINSG,MOSTV,FLAG)
      INTEGER N,NP1,PRTLEV,VDLT,VERT(N,NP1),ITOP,PI(N),SIGN(N)
      INTEGER IPOINT(NP1),MOSTV(NP1),FLAG
      INTEGER IOUT,KOUT,KOUTSG,IIN,KIN,KINSG
      DOUBLE PRECISION DELTA(N),XORIG(N),VTWT(NP1),XNEW(N),DBLE  
C
C     NEWPCE UPDATES THE PIECE OF THE SUBDIVISION
C
C     ON ENTRY,
C
C     VERT,ITOP,PI,SIGN      DETERMINE THE CURRENT PIECE
C
C     IOUT                   IF NONZERO INDICATES THE INDEX OF
C                            THE VERTEX ON TOP TO BE DROPPED
C                            IF ZERO VERTICES ON BOTTOM DROPPED
C
C     KOUT                   IF ZERO VERTEX ON TOP TO BE DROPPED
C                            OTHERWISE DETERMINES WITH KOUTSG
C                            THE SET OF VERTICES TO BE DROPPED
C
C     KOUTSG                 THE VERTICES TO BE DROPPED ARE
C
C                            (VERT(.,ITOP)+SUM(+ OR - E(KOUT)), 0 )
C
C                            WHERE THE SIGN OF E(KOUT) IS KOUTSG
C                            IF VDLT.GE.1 REPLACE VERT BY S(VERT)
C
C     ON RETURN,
C
C     VERT,ITOP,PI,SIGN      DETERMINE THE NEW PIECE
C
C     IIN                    IF NONZERO INDICATES INDEX OF NEW VERTEX
C                            ON TOP
C
C     XNEW                   IF IIN NONZERO XNEW IS CORRESPONDING POINT
C                            IN THE DOMAIN OF F
C
C     KIN                    IF IIN NONZERO KIN IS ZERO. OTHERWISE
C                            DETERMINES WITH KINSG THE NEW SET OF
C                            VERTICES
C
C     KINSG                  THE NEW VERTICES ARE
C
C                            (VERT(.,ITOP)+SUM(+ OR - E(KIN)), 0 )
C
C                            WHERE E(KIN) HAS SIGN KINSG
C                            IF VDLT.GE.1 REPLACE VERT BY S(VERT)
C
C     SUBROUTINES AND FUNCTIONS USED -- PRTPCE, DBLE   
C
C
      INTEGER I,MO,MOM,MOP,IM1,II,IOUTM1,ITOPM1,J,J1,J2,K,SIGN1,SIGN2
C
C     BEGIN PROGRAM
C
C     DETERMINE WHERE LEAVING VERTEX (VERTICES) ARE
C
      IOUTM1=IOUT-1
      IF (IOUT.EQ.0) GO TO 300
         IF (IOUT.EQ.ITOP) GO TO 200
            IF (IOUT.EQ.1) GO TO 100
C
C     INTERNAL VERTEX LEAVES
C
               J1=PI(IOUTM1)
               J2=PI(IOUT)
               SIGN1=SIGN(J1)
               SIGN2=SIGN(J2)
               IF ((IOUT.EQ.2).AND.(VDLT.LT.0)) GO TO 50
               VERT(J1,IOUT)=VERT(J1,IOUT)-SIGN1
               VERT(J2,IOUT)=VERT(J2,IOUT)+SIGN2
               PI(IOUTM1)=J2
               PI(IOUT)=J1
               IIN=IOUT
               GO TO 500
   50          CONTINUE
               IF (SIGN2.LT.0) GO TO 80
               VERT(J1,2)=VERT(J1,2)-2
               VERT(J2,2)=VERT(J2,2)+2
               PI(1)=J2
               PI(2)=J1
               IIN=2
               GO TO 500
   80          CONTINUE
               VERT(J1,2)=VERT(J1,2)-2
               VERT(J2,1)=VERT(J2,1)-2
               PI(1)=J2
               PI(2)=J1
               SIGN(J2)=1
               DO 90 J=1,NP1
                  IF (IPOINT(J).EQ.1) IPOINT(J)=2
   90             CONTINUE
               MOSTV(2)=MOSTV(1)
               IIN=1
               GO TO 500
C
C     FIRST VERTEX LEAVES
C
  100       J1=PI(1)
            IF (VDLT.GE.2) GO TO 120
            IF (VDLT.GE.0) GO TO 115
            VERT(J1,1)=VERT(J1,1)+2
            SIGN(J1)=-1
            SIGN2=1
            IF (ITOP.GT.2) J2=PI(2)
            IF (ITOP.GT.2) SIGN2=SIGN(J2)
            IF (SIGN2.LT.0) GO TO 110
            DO 105 J=1,NP1
               IF (IPOINT(J).EQ.2) IPOINT(J)=1
  105          CONTINUE
            MOSTV(1)=MOSTV(2)
            ITOP=ITOP-1
            IIN=0
            KIN=J1
            KINSG=1
            IF (ITOP.EQ.1) GO TO 600
            ITOP=ITOP+1
            VERT(J2,2)=VERT(J2,2)+2
            IIN=2
            PI(1)=J2
            PI(2)=J1
            GO TO 500
  110       CONTINUE
            VERT(J2,1)=VERT(J2,1)-2
            IIN=1
            PI(1)=J2
            PI(2)=J1
            SIGN(J2)=1
            GO TO 500
  115       CONTINUE
            SIGN1=SIGN(J1)
            VERT(J1,1)=VERT(J1,1)+2*SIGN1
            SIGN(J1)=-SIGN1
            IIN=1
            GO TO 500
  120       J=J1
            SIGN1=SIGN(J)
            VERT(J,ITOP)=VERT(J,ITOP)+SIGN1
            MO=MOSTV(ITOP)
            DO 130 I=2,ITOP
               II=ITOP+1-I
               K=PI(II)
               PI(II)=J
               J=K
               SIGN1=SIGN(J)
               VERT(J,II)=VERT(J,II)+SIGN1
               MOM=MOSTV(II)
               MOSTV(II)=MO
               MO=MOM
  130          CONTINUE
            DO 140 J=1,NP1
               IF (IPOINT(J).GT.0) IPOINT(J)=IPOINT(J)-1
  140          CONTINUE
            IF ((VDLT.GE.3).AND.(VERT(J1,1).EQ.0)) GO TO 150
            IIN=ITOP
            GO TO 500
  150       ITOP=ITOP-1
            KIN=J1
            IIN=0
            KINSG=1
            GO TO 600
C
C     LAST VERTEX LEAVES
C
  200    CONTINUE
         IF (ITOP.GT.1) GO TO 210
            FLAG=2
            WRITE (6,205)
  205       FORMAT (45H0**** ATTEMPTING TO REMOVE ONLY VERTEX ON TOP )
            GO TO 700
  210       CONTINUE
         ITOPM1=ITOP-1
         J1=PI(ITOPM1)
         SIGN1=SIGN(J1)
C
C     BRANCH IF VDL&T SUBDIVISION USED AND
C     DIMENSION OF SIMPLEX ON TOP UNCHANGED
C
         IF ((VDLT.GE.1).AND.(VERT(J1,1).NE.0)) GO TO 215
            ITOP=ITOP-1
            KIN=J1
            IIN=0
            KINSG=-SIGN1
            GO TO 600
  215    IF (VDLT.EQ.1) GO TO 250
            J=J1
            VERT(J,1)=VERT(J,1)-SIGN1
            MO=MOSTV(1)
            DO 220 I=2,ITOP
               IM1=I-1
               K=PI(IM1)
               PI(IM1)=J
               J=K
               SIGN1=SIGN(J)
               VERT(J,I)=VERT(J,I)-SIGN1
               MOP=MOSTV(I)
               MOSTV(I)=MO
               MO=MOP
  220          CONTINUE
            DO 230 J=1,NP1
               IF (IPOINT(J).GT.0) IPOINT(J)=IPOINT(J)+1
  230          CONTINUE
            IIN=1
            GO TO 500
  250    VERT(J1,ITOP)=VERT(J1,ITOP)-2*SIGN1
         SIGN(J1)=-SIGN1
         IIN=ITOP
         GO TO 500
C
C     VERTICES ON BOTTOM LEAVE
C
  300 PI(ITOP)=KOUT
      IF ((VDLT.GE.3).AND.(KOUTSG.EQ.1)) GO TO 360
      SIGN(KOUT)=-KOUTSG
      ITOPM1=ITOP
      ITOP=ITOP+1
      DO 310 J=1,N
         VERT(J,ITOP)=VERT(J,ITOPM1)
  310    CONTINUE
      VERT(KOUT,ITOP)=VERT(KOUT,ITOPM1)-KOUTSG
      IIN=ITOP
      IF ((VDLT.GE.0).OR.(ITOP.GE.4)) GO TO 500
      IF (ITOP.EQ.3) GO TO 340
         IF (KOUTSG.EQ.1) GO TO 320
            VERT(KOUT,2)=VERT(KOUT,2)+1
            GO TO 500
  320       CONTINUE
         DO 330 J=1,NP1
            IF (IPOINT(J).EQ.1) IPOINT(J)=2
  330       CONTINUE
         MOSTV(2)=MOSTV(1)
         VERT(KOUT,2)=VERT(KOUT,2)+1
         VERT(KOUT,1)=VERT(KOUT,1)-2
         SIGN(KOUT)=1
         IIN=1
         GO TO 500
  340    CONTINUE
      J1=PI(1)
      VERT(J1,3)=VERT(J1,3)-1
      GO TO 500
  360 SIGN(KOUT)=1
      J=KOUT
      SIGN1=1
      VERT(J,1)=VERT(J,1)-SIGN1
      MO=MOSTV(1)
      IF (ITOP.EQ.1) GO TO 375
      DO 370 I=2,ITOP
         IM1=I-1
         K=PI(IM1)
         PI(IM1)=J
         J=K
         SIGN1=SIGN(J)
         VERT(J,I)=VERT(J,I)-SIGN1
         MOP=MOSTV(I)
         MOSTV(I)=MO
         MO=MOP
  370    CONTINUE
  375 CONTINUE
      PI(ITOP)=J
      ITOPM1=ITOP
      ITOP=ITOP+1
      MOSTV(ITOP)=MO
      DO 380 K=1,N
         VERT(K,ITOP)=VERT(K,ITOPM1)
  380    CONTINUE
      VERT(J,ITOP)=VERT(J,ITOP)+SIGN1
      DO 390 J=1,NP1
         IF (IPOINT(J).GT.0) IPOINT(J)=IPOINT(J)+1
  390    CONTINUE
      IIN=1
  500 KIN=0
      DO 550 J=1,N
         XNEW(J)=XORIG(J)+DELTA(J)*  DBLE(VERT(J,IIN))
  550    CONTINUE
  600 IF (PRTLEV.GE.5) CALL PRTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,
     *                             DELTA,XORIG,IOUT,KOUT,KOUTSG,VTWT,
     *                             XNEW,IIN,KIN,KINSG)
  700 RETURN
C
C     LAST CARD OF NEWPCE FOLLOWS
C
      END
      SUBROUTINE LSTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,DELTA,
     *                  XORIG,IOUT,KOUT,KOUTSG,VTWT,
     *                  XNEW,IIN,KIN,KINSG)
      INTEGER N,NP1,PRTLEV,VDLT,VERT(N,NP1),ITOP,PI(N),SIGN(N)
      INTEGER IOUT,KOUT,KOUTSG,IIN,KIN,KINSG
      DOUBLE PRECISION DELTA(N),XORIG(N),VTWT(NP1),XNEW(N),DBLE  
C
C     LSTPCE CALCULATES THE APPROXIMATE ZERO BY WEIGHTING THE
C     CURRENT VERTICES BY THE VECTOR VTWT
C
C     SUBROUTINES AND FUNCTIONS USED -- PRTPCE, DBLE  
C
      INTEGER I,II,IP1,J
      DOUBLE PRECISION TEMP
C
C
      DO 10 J=1,N
         XNEW(J)=  DBLE(VERT(J,1))
   10    CONTINUE
      IF (ITOP.EQ.1) GO TO 25
      TEMP=0.0D0
      I=ITOP
      DO 20 II=2,ITOP
         TEMP=TEMP+VTWT(I)
         I=I-1
         J=PI(I)
         IF (SIGN(J).GT.0) XNEW(J)=XNEW(J)+TEMP
         IF (SIGN(J).LT.0) XNEW(J)=XNEW(J)-TEMP
   20    CONTINUE
      IF (VDLT.LT.0) XNEW(J)=XNEW(J)+VTWT(2)
   25 CONTINUE
      DO 30 J=1,N
         XNEW(J)=XORIG(J)+DELTA(J)*XNEW(J)
   30    CONTINUE
C
C     IF DESIRED PRINT LAST PIECE
C
      IF (PRTLEV.GE.5) CALL PRTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,
     *                             DELTA,XORIG,IOUT,KOUT,KOUTSG,VTWT,
     *                             XNEW,IIN,KIN,KINSG)
      RETURN
C
C     LAST CARD OF LSTPCE FOLLOWS
C
      END
      SUBROUTINE PRTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,DELTA,
     *                  XORIG,IOUT,KOUT,KOUTSG,VTWT,XNEW,IIN,KIN,KINSG)
      INTEGER N,NP1,PRTLEV,VDLT,VERT(N,NP1),ITOP,PI(N),SIGN(N)
      INTEGER IOUT,KOUT,KOUTSG,IIN,KIN,KINSG
      DOUBLE PRECISION DELTA(N),XORIG(N),VTWT(NP1),XNEW(N)
C
C     PRTPCE PRINTS THE CURRENT PIECE, WHAT VERTEX (VERTICES) WAS
C     LAST DROPPED, AND WHAT VERTEX (VERTICES) HAVE JUST ENTERED
C
      INTEGER I,INT0,INT1,ITOPM1,J,K,SVERT,PII
      CHARACTER CHAR1,CHAR2,CHAR3,CHAR4
      CHARACTER IPLUS,ISLASH,IMINUS,IONE,IBLANK
      DATA IBLANK,IPLUS,ISLASH,IMINUS,IONE/' ','+','/','-','1'/
      DATA INT0,INT1/0,1/
      WRITE (6,100) ITOP
      WRITE (6,200)
      ITOPM1=ITOP-1
      IF (ITOPM1.LE.0) GO TO 5
      DO 2 I=1,ITOPM1
         PII=PI(I)
         WRITE (6,300) PI(I),SIGN(PII)
    2    CONTINUE
    5 IF (VDLT.GE.1) GO TO 30
         IF (VDLT.EQ.0) WRITE (6,400)
         IF (VDLT.LT.0) WRITE (6,450)
         DO 20 J=1,N
            CHAR1=IBLANK
            CHAR2=IBLANK
            CHAR3=IBLANK
            CHAR4=IBLANK
            IF (VERT(J,1).NE.VERT(J,ITOP)) GO TO 10
               CHAR1=IPLUS
               CHAR2=ISLASH
               CHAR3=IMINUS
               CHAR4=IONE
   10       CONTINUE
            WRITE (6,600) (VERT(J,I),(IBLANK,K=1,4),I=1,ITOP),
     *       VERT(J,ITOP),CHAR1,CHAR2,CHAR3,CHAR4
   20       CONTINUE
         WRITE (6,600)(INT1,(IBLANK,K=1,4),I=1,ITOP),INT0,(IBLANK,K=1,4)
         GO TO 90
   30 IF (VDLT.EQ.1) WRITE (6,500)
      IF (VDLT.EQ.2) WRITE (6,510)
      IF (VDLT.GE.3) WRITE (6,520)
      DO 80 J=1,N
         CHAR1=IBLANK
         CHAR2=IBLANK
         CHAR3=IBLANK
         CHAR4=IBLANK
         IF (VERT(J,ITOP)) 40,50,60
   40       SVERT=-1
            GO TO 70
   50       SVERT=0
            CHAR1=IPLUS
            CHAR2=ISLASH
            CHAR3=IMINUS
            CHAR4=IONE
            GO TO 70
   60       SVERT=1
   70    CONTINUE
         WRITE (6,600) (VERT(J,I),(IBLANK,K=1,4),I=1,ITOP),
     *    SVERT,CHAR1,CHAR2,CHAR3,CHAR4
   80    CONTINUE
      WRITE (6,600) (INT1,(IBLANK,K=1,4),I=1,ITOP),INT0,(IBLANK,K=1,4)
   90 CONTINUE
      WRITE (6,700) IOUT,KOUT,KOUTSG,IIN,KIN,KINSG
      WRITE (6,800) (XNEW(J),J=1,N)
C
C     FORMATS
C
  100 FORMAT (40H0     THE NUMBER OF VERTICES ON TOP IS        ,I5)
  200 FORMAT (53H0     THE COORDINATES USED ON TOP AND THEIR SIGNS ARE/)
  300 FORMAT (6X,14I5)
  400 FORMAT (35H0     USING J1 THE VERTICES ARE        /)
  450 FORMAT (45H0     USING J-PRIME THE VERTICES ARE                 /)
  500 FORMAT (53H0     USING VDLTS SUBDIVISION WITH J1 THE VERTICES    ,
     * 4H ARE/)
  510 FORMAT(55H0     USING VDLTS SUBDIVISION WITH K-PRIME THE VERTICES,
     * 4H ARE/)
  520 FORMAT (53H0     USING VDLTS SUBDIVISION WITH K1 THE VERTICES    ,
     * 4H ARE/)
  600 FORMAT (8(1X,I5,4A1))
  700 FORMAT (11H0     IOUT= ,I3,6H KOUT= ,I3,8H KOUTSG= ,I3,
     * 5H IIN= ,I3,5H KIN= ,I3,7H KINSG= ,I3)
  800 FORMAT (15H0     XNEW=     ,7D10.3)
      RETURN
C
C     LAST CARD OF PRTPCE FOLLOWS
C
      END