./ ADD NAME=README TIME=492783148
FILES ON THIS DISKETTE INCLUDE THE FORTRAN SUBROUTINES

                DIVDIF, XDERIV, YDERIV AND EVAL

THESE ARE REASONABLY THOROUGHLY DOCUMENTED VIA COMMENT CARDS.

ALSO ON THE DISKETTE ARE TWO EXAMPLE TEST PROGRAMS TEST1 AND TEST2.

THE PROGRAM AS CURRENTLY WRITTEN ONLY HANDLES DATA MONOTONE INCREASING

IN BOTH X AND Y VARIABLES.  ENHANCEMENTS FOR DATA OF VARYING MONOTONICITY

ARE IN PROGRESS.  THERE ALSO EXISTS AN EVALUATOR WHICH IS MORE MEMORY

INTENSIVE BUT HOPEFULLY FASTER.  THIS STORES THE COEFFICIENTS

OF THE QUADRATIC ON EACH SUB-TRIANGLE IN A HUGE ARRAY.

I HAVE NOT INCLUDED THIS ON THE DISKETTE BUT IT WILL BE AVAILABLE

SOON.

     PLEASE CONTACT ME ABOUT ANY ERRORS IN THE PAPER, OR CODE. I WOULD

BE VERY INTERESTED IN ANY ENHANCEMENTS YOU MAKE TO IMPROVE SPEED ETC.

                                          YOURS,

                                          RICK BEATSON.


./ ADD NAME=divdif.f TIME=492785136
      SUBROUTINE DIVDIF(N,M,X,Y,F,HX,HY,DDX,DDY)
C
C SUBROUTINE TO INITIALIZE ARRAYS OF SUBINTERVAL LENGTHS (HX, HY) AND
C OF FIRST DIVIDED DIFFERENCES IN THE TWO COORDINATE DIRECTIONS (DDX,DDY
C
C
C INPUT      N     NUMBER OF DATA POINTS ON A X GRID LINE.
C            M     NUMBER OF DATA POINTS ON A Y GRID LINE.
C            X     ARRAY OF X VALUES. SUBSCRIPTS 1-N. DIMENSION N.
C            Y     ARRAY OF Y VALUES. SUBSCRIPTS 1-M. DIMENSION M.
C            F     ARRAY OF INPUT Z VALUES TO WHICH THE SPLINE WILL BE
C                  FITTED. DIMENSION N BY M.
C
C
C OUTPUT     HX    ARRAY OF X SUBINTERVAL LENGTHS. DIMENSION N.
C            HY    ARRAY OF Y SUBINTERVAL LENGTHS. DIMENSION M.
C            DDX   ARRAY OF FIRST DIVIDED DIFFERENCES IN THE X DIRECTION
C                  DIMENSION N BY M.
C            DDY   ARRAY OF FIRST DIVIDED DIFFERENCES IN THE Y DIRECTION
C                  DIMENSION N BY M.
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(N),Y(M),F(N,M),HX(N),HY(M),DDX(N,M),DDY(N,M)
      DO 1 I=1,N-1
         IP1=I+1
         H=X(IP1)-X(I)
         HX(I)=H
         DO 1 J=1,M
            DDX(I,J)=(F(IP1,J)-F(I,J))/H
1     CONTINUE
      DO 2 J=1,M-1
         JP1=J+1
         H=Y(JP1)-Y(J)
         HY(J)=H
         DO 2 I=1,N
            DDY(I,J)=(F(I,JP1)-F(I,J))/H
2     CONTINUE
      RETURN
      END

./ ADD NAME=eval.f TIME=492785137
       SUBROUTINE BINSOR(A,N,X,L)
C ARRAY A ASSUMED ORDERED IN ASCENDING ORDER
C RETURNS WITH L THE GREATEST L SUCH THAT A(L) IS LESS THAN OR EQUAL TO
C X, SUBJECT TO 1 < OR = L < OR = N-1.
       IMPLICIT REAL*8 (A-H,O-Z)
       DIMENSION A(N)
       L=1
       J=N
C
1      CONTINUE
       IF (J-L.LE.1) RETURN
C
           I=(L+J)/2
           IF(A(I).LE.X) THEN
                  L=I
           ELSE
                  J=I
           END IF
       GO TO 1
       END
C
C
       FUNCTION EVAL(X,Y,F,XPART,YPART,N,M,UIN,VIN,IFLAG)
C
C     FUNCTION TO EVALUATE A PIECEWISE QUADRATIC SPLINE IN TWO VARIABLES
C     OF THE TYPE SPECIFIED IN THE BEATSON ZIEGLER ARTICLE TO APPEAR
C     IN SIAM J. NUMERICAL ANALYSIS 1985.
C
C     SUBROUTINES DIVDIF, XDERIV, YDERIV ARE TO BE CALLED BEFORE
C     THIS ROUTINE.
C
C     THIS ROUTINE CALLS THE BINARY SORT ROUTINE BINSOR.
C
C
C     PARAMETERS X,Y   ARRAYS SPECIFYING THE DATA MESH
C                F     ARRAY OF Z VALUES TO INTERPOLATE
C                N     NUMBER OF POINTS IN X ARRAY SUBSCRIPTS 1-N.
C                M     NUMBER OF POINTS IN Y ARRAY SUBSCRIPTS 1-M.
C                XPART ARRAY OF XPARTIALS FITTED BY PREVIOUS CALLS TO
C                      DIVDIF AND XDERIV.  DIMENSION N BY M. SUBSCRIPTS
C                      START AT 1.
C                YPART ARRAY OF Y PARTIALS FITTED BY PREVIOUS CALLS TO
C                      DIVDIF AND YDERIV.  DIMENSION N BY M. SUBSCRIPTS
C                      START AT 1.
C                UIN   INPUT X COORDINATE AT WHICH THE SPLINE IS TO BE
C                      EVALUATED.
C                VIN   INPUT Y COORDINATE AT WHICH THE SPLINE IS TO BE
C                      EVALUATED.
C
C
       IMPLICIT REAL*8 (A-H,O-Z)
       DIMENSION X(N),Y(M),F(N,M),XPART(N,M),YPART(N,M)
       INTEGER SUBSQU,TRIANG,OLDTRI,OLDSUB,OLDI,OLDJ
       LOGICAL*1 DOIT
       SAVE OLDI,OLDJ,OLDSUB,OLDTRI
       SAVE F0,F1,F2,F3,F0X,F0Y,F1X,F1Y,F2X,F2Y,F3X,F3Y,XM,YM
       SAVE XL,YB,W0,W1,W2,W3,W0P,W1P,W2P,W3P
       SAVE T0,T1,T2,T01,T12,T20
C
C
C IFLAG EQUAL 0 FORCES REINITIALIZATION OF SAVED VARIABLES AND SHOULD
C BE SET TO O WHEN CHANGING TO A DIFFERENT SURFACE IN THE SAME RUN.
C IT IS CHANGED TO 1 BY THE PROGRAM SO THAT THE SUBSEQUENT CALLS ,
C FOR THE SAME SURFACE ARE SOMEWHAT EFFICIENT.
C
C
C A FIRST STEP TO MAKING EVAL FASTER WOULD BE TO MAKE IT ACCEPT ARRAYS
C OF U'S AND V'S AT WHICH TO EVALUATE THE SPLINE SURFACE.  THIS SHOULD
C GREATLY DECREASE THE OVERHEAD INHERENT IN CALLING A SUBROUTINE
C HUNDREDS ( THOUSANDS ) OF TIMES. A SECOND STEP WOULD BE TO MAKE USE
C OF THE FACT THAT WE WILL USUALLY MOVE FROM ONE SUPERSQUARE TO AN
C ADJACENT ONE.
C
        U=UIN
        V=VIN
       IF(IFLAG.EQ.ZERO)THEN
             DOIT=.TRUE.
             IFLAG=1
       ELSE
             DOIT=.FALSE.
       END IF
C
C
C IDENTIFY THE SUPERSQUARE AND AVOID SOME WORK IF THIS IS THE SAME AS
C DURING THE LAST FUNCTION EVALUATION.
C
C
       CALL BINSOR(X,N,U,I)
       CALL BINSOR(Y,M,V,J)
C
C
       IF(.NOT.((OLDI.EQ.I).AND.(OLDJ.EQ.J)))DOIT=.TRUE.
          IF (DOIT) THEN
C
C
       OLDI=I
       OLDJ=J
       IP1=I+1
       JP1=J+1
       F0=F(I,J)
       F1=F(IP1,J)
       F2=F(IP1,JP1)
       F3=F(I,JP1)
       FACX=(X(I+1)-X(I))/2.0
       FACY=(Y(J+1)-Y(J))/2.0
       XM=X(I)+FACX
       YM=Y(J)+FACY
       F0X=XPART(I,J)*FACX
       F0Y=YPART(I,J)*FACY
       F1X=XPART(IP1,J)*FACX
       F1Y=YPART(IP1,J)*FACY
       F2X=XPART(IP1,JP1)*FACX
       F2Y=YPART(IP1,JP1)*FACY
       F3X=XPART(I,JP1)*FACX
       F3Y=YPART(I,JP1)*FACY
       END IF
C
C
C IDENTIFY THE SUBSQUARE WITHIN THE SUPERSQUARE AND AVOID SOME WORK IF
C THE SAME SUBSQUARE WAS USED IN THE PREVIOUS FUNCTION EVALUATION.
C
C
       IF(U.LE.XM)THEN
           IF(V.LE.YM) THEN
               SUBSQU=1
           ELSE
               SUBSQU=4
           END IF
       ELSE
           IF(V.LE.YM) THEN
               SUBSQU=2
           ELSE
               SUBSQU=3
           END IF
       END IF
C
C
       IF (.NOT.(OLDSUB.EQ.SUBSQU)) DOIT=.TRUE.
          IF(DOIT) THEN
C
C
C
C MAKE THE APPROPRIATE ASSIGNMENTS TO W0,W1,W2,W3,W0P,W1P,W2P,W3P
C THE DUAL FUNCTIONALS FOR THE SUBSQUARE WITHIN THE SUPERSQUARE.
C
C NOTE THAT THESE DUAL FUNCTIONALS ARE INVARIANT UNDER AFFINE
C TRANSFORMATIONS OF THE (X,Y) PLANE.
C
C IN WHAT FOLLOWS THE DOMAIN IS CONVERTED BY AN AFFINE TRANSFORMATION
C TO THE STANDARD SUBSQUARE [0,1]X[0,1].
C
           OLDSUB=SUBSQU
        IF(SUBSQU.EQ.1) THEN
C BOTTOM LEFT SUBSQUARE
            XL=X(I)
            YB=Y(J)
            W0=F0
            W1=0.5*(F0+F1) +0.25*(F0X-F1X)
            W2=0.25*(F0 +F1 +F2 +F3) +0.125*(F0X+F0Y +F1Y + F3X)
     1                              -0.125*(F1X+F2X+ F2Y +F3Y)
            W3=0.5*(F0+F3) +0.25*(F0Y-F3Y)
            W0P=F0X
            W1P=0.5*(F0Y +F1Y)
            W2P=-0.5*((F1-F0) -0.5*(F0X+F1X) +0.5*(F1Y-F0Y))
     1         -0.5*((F2-F3) -0.5*(F2X+F3X) +0.5*(F3Y-F2Y))
            W3P=F0-F3+0.5*(F0Y+F3Y)
         ELSE IF (SUBSQU.EQ.4) THEN
C TOP LEFT SUBSQUARE
            XL=X(I)
            YB=YM
            W0=0.5*(F0+F3) + 0.25*(F0Y-F3Y)
            W1= 0.25*(F0 +F1 +F2 +F3) +0.125*(F0X +F0Y + F1Y +F3X)
     1                               -0.125*(F1X +F2X + F2Y +F3Y)
            W2= 0.5*(F2+F3)  +0.25*(F3X-F2X)
            W3=F3
            W0P=0.5*(F0X+F3X)
            W1P=0.5*((F2-F1)-0.5*(F1Y+F2Y)+0.5*(F1X-F2X))
     1        +0.5*((F3-F0)-0.5*(F3Y+F0Y)+0.5*(F3X-F0X))
            W2P=F3 -F2 +0.5*(F3X+ F2X)
            W3P=-F3Y
       ELSE IF (SUBSQU.EQ.2) THEN
C BOTTOM RIGHT SUBSQUARE
            XL=XM
            YB=Y(J)
            W0=0.5*(F0+F1) +0.25*(F0X-F1X)
            W1=F1
            W2=0.5*(F1+F2)+ 0.25*(F1Y-F2Y)
            W3=0.25*(F0+F1+F2+F3)+ 0.125*(F0X +F0Y+F1Y +F3X)
     1                          - 0.125*(F1X +F2X+F2Y +F3Y)
            W0P= F1-F0-0.5*(F0X +F1X)
            W1P=F1Y
            W2P=-0.5*(F1X+F2X)
            W3P=-0.5*((F2-F1) -0.5*(F1Y+F2Y) +0.5*(F1X- F2X))
     1         -0.5*((F3-F0) -0.5*(F3Y+F0Y)+ 0.5*(F3X- F0X))
         ELSE
C TOP RIGHT SUBSQUARE
            XL=XM
            YB=YM
            W0=0.25*(F0+ F1 +F2 +F3) +0.125*(F0X +F0Y +F1Y +F3X)
     1                              -0.125*(F1X+ F2X +F2Y +F3Y)
            W1=0.5*(F1+F2) + 0.25*(F1Y-F2Y)
            W2=F2
            W3=0.5*(F2+F3)+0.25*(F3X-F2X)
            W0P=+0.5*((F1-F0)-0.5*(F0X+F1X)+0.5*(F1Y-F0Y))
     1         +0.5*((F2-F3)-0.5*(F2X+F3X)+0.5*(F3Y-F2Y))
            W1P=F2-F1 -0.5*(F1Y +F2Y)
            W2P=-F2X
            W3P=-0.5*(F2Y +F3Y)
         END IF
C
       END IF
C
C
C
C
C IDENTIFY TRIANGLE WITHIN SUBSQUARE AND AVOID SOME WORK IF THIS IS THE
C SAME AS FOR THE LAST EVALUATION.
C
C
        U=(U-XL)/FACX
        V=(V-YB)/FACY
        IF(U+V.LE.1.0) THEN
           IF(V-U.LE.0.0) THEN
               TRIANG=1
           ELSE
               TRIANG=4
           END IF
       ELSE
           IF(V-U.LE.0.0) THEN
               TRIANG=2
           ELSE
               TRIANG=3
           END IF
       END IF
C
C
       IF(.NOT.(OLDTRI.EQ.TRIANG)) DOIT=.TRUE.
          IF(DOIT) THEN
C
C
C EVALUATE THE FUNCTIONALS FOR THE SUBTRIANGLE
C
C NOTE THAT IN EACH SUBTRIANGLE
C ALPHA 0 IS TAKEN AS THE LEFTMOST (AND IN THE CASE OF A TIE
C ALSO THE TOPMOST) VERTEX. ALSO NOTE THAT EACH EDGE FUNCTIONAL
C (FOR AN EDGE OF THE SUBTRIANGLE WITH ALPHA 4 AS AN ENDPOINT)
C T(J,J+1) = 2*ZP(X(I)) +DP(X(I))(X(I+1) -X(I))
C          = 2*ZP(X(I+1)) +DP(X(I+1))(X(I)-X(I+1))
C CAN BE EXPRESSED IN TERMS OF THE SUBSQUARE FUNCTIONALS AS
C             W(I-1)+W(I) +0.5*(W(I-1)' +W(I)')
C WHERE ALPHA I IS THE SUBSQUARE VERTEX WHICH THE TRIANGLE
C EDGE JOINS TO ALPHA 4.
C
       OLDTRI=TRIANG
       ZM= 0.25*(W0+W1+W2+W3) +0.125*(W0P+W1P+W2P+W3P)
       IF(TRIANG.EQ.1) THEN
C BOTTOM TRIANGLE
           T0=W0
           T1=W1
           T2=ZM
           T01=2*W0+W0P
           T12=W0+W1+0.5*(W0P+W1P)
           T20=W3+ W0+0.5*(W3P+W0P)
        ELSE IF (TRIANG.EQ.4) THEN
C LEFTMOST TRIANGLE
            T0=W3
            T1=W0
            T2=ZM
            T01=2*W3+W3P
            T12=W3+W0+0.5*(W3P+W0P)
            T20=W2+W3+0.5*(W2P+W3P)
       ELSE IF(TRIANG.EQ.2) THEN
C RIGHTMOST TRIANGLE
            T0=ZM
            T1=W1
            T2=W2
            T01=W0+W1 +0.5*(W0P+W1P)
            T12=2*W1+W1P
            T20=W1+W2 +0.5*(W1P+W2P)
         ELSE
C TOP TRIANGLE
            T0=W3
            T1=ZM
            T2=W2
            T01=W2+W3+0.5*(W2P+W3P)
            T12=W1+W2+0.5*(W1P+W2P)
            T20=2*W2+W2P
       END IF
C
C
       END IF
C MAKE THE APPROPRIATE ASSIGNMENTS TO XLAM0,XLAM1,XLAM2
C FOR THE SUBTRIANGLE WITHIN THE SUBSQUARE. THESE ARE
C THE BARYCENTRIC COORDINATES.
       IF(TRIANG.EQ.1) THEN
C BOTTOM TRIANGLE
            XLAM0 =1.0 -U -V
            XLAM1 =U -V
            XLAM2 =2.0*V
        ELSE IF (TRIANG.EQ.4) THEN
C LEFTMOST TRIANGLE
            XLAM0=V-U
            XLAM1=1.0-U-V
            XLAM2=2.0*U
       ELSE IF(TRIANG.EQ.2) THEN
C RIGHTMOST TRIANGLE
            XLAM0=2.0*(1.0-U)
            XLAM1=U -V
            XLAM2= U+V-1.0
         ELSE
C TOP TRIANGLE
            XLAM0=V-U
            XLAM1=2.0*(1.0-V)
            XLAM2=U+V -1.0
       END IF
       EVAL=(XLAM0*(T0*XLAM0+T01*XLAM1)) +(XLAM1*XLAM1*T1)
     1      +  (XLAM2*(T2*XLAM2+T12*XLAM1+T20*XLAM0))
       RETURN
       END

./ ADD NAME=test1.f TIME=492785138
       FUNCTION FUNC(X,Y)
       IMPLICIT REAL*8 (A-H,O-Z)
       FUNC=X*Y
       RETURN
       END
C
C
       PARAMETER (N=5,M=5)
       IMPLICIT REAL*8 (A-H,O-Z)
       DIMENSION X(N),Y(M),F(N,M),HX(N),HY(M),DDX(N,M),DDY(N,M)
       DIMENSION XPART(N,M),YPART(N,M)
       DO 1 I=1,N
             X(I)=FLOAT(I-1)
1      CONTINUE
       DO 2 J=1,M
             Y(J)=FLOAT(2*J-2)
2      CONTINUE
       WRITE(3,*)'ECHO PRINT THE DATA'
       DO 3 I=1,N
       DO 3 J=1,M
             F(I,J)=FUNC(X(I),Y(J))
             WRITE(3,50)X(I),Y(J),F(I,J)
3      CONTINUE
       CALL DIVDIF(N,M,X,Y,F,HX,HY,DDX,DDY)
       WRITE(3,*)'AFTER CALL TO DIVDIF     X(I) AND Y(J) FOLLOW'
       WRITE(3,50)(X(I),I=1,N)
       WRITE(3,50)(Y(J),J=1,M)
       WRITE(3,*)'HX HY FOLLOW'
       WRITE(3,50)(HX(I),I=1,N-1)
       WRITE(3,50)(HY(J),J=1,M-1)
       CALL XDERIV(N,M,HX,HY,DDX,DDY,XPART)
       CALL YDERIV(N,M,HX,HY,DDX,DDY,YPART)
       WRITE(3,*)'DDX AND XPART FOLLOW   SUBSCRIPT I CHANGES FASTEST'
       DO 56 J=1,M
       WRITE(3,50)(DDX(I,J),I=1,N-1)
       WRITE(3,50)(XPART(I,J),I=1,N)
56     CONTINUE
       WRITE(3,*)'DDY AND YPART FOLLOW    SUBSCRIPT J CHANGES FASTEST'
       DO 55 I=1,N
       WRITE(3,50)(DDY(I,J),J=1,M-1)
       WRITE(3,50)(YPART(I,J),J=1,M)
55     CONTINUE
       H=.25
       WRITE(3,*)'X       Y        S(X,Y)    F(X,Y)   ERROR FOLLOW'
 
       DO 4 I=1,17
          DO 4 J=1,17
              U=FLOAT(I-1)*H
              V=FLOAT(J-1)*H
              Z=EVAL(X,Y,F,XPART,YPART,N,M,U,V,IFLAG)
              WRITE(3,50)U,V,Z,FUNC(U,V),Z-FUNC(U,V)
4     CONTINUE
      STOP
50    FORMAT(1X,5(E13.5,2X))
      END
./ ADD NAME=test2.f TIME=492785139
       FUNCTION FUNC(X,Y)
       IMPLICIT REAL*8 (A-H,O-Z)
       R=SQRT(X*X+Y*Y)
       IF(R.GT.2.0D0)THEN
            FUNC=R-2.0D0
       ELSE
            FUNC=0.0D0
       END IF
       RETURN
       END
C
C
       PARAMETER (N=5,M=5)
       IMPLICIT REAL*8 (A-H,O-Z)
       DIMENSION X(N),Y(M),F(N,M),HX(N),HY(M),DDX(N,M),DDY(N,M)
       DIMENSION XPART(N,M),YPART(N,M)
       DO 1 I=1,N
            X(I)=FLOAT(I-1)
1      CONTINUE
       DO 2 J=1,M
            Y(J)=FLOAT(J-1)
2      CONTINUE
       DO 3 I=1,N
          DO 3 J=1,M
             F(I,J)=FUNC(X(I),Y(J))
3      CONTINUE
       CALL DIVDIF(N,M,X,Y,F,HX,HY,DDX,DDY)
       WRITE(3,*)'AFTER CALL TO DIVDIF     X(I) AND Y(J) FOLLOW'
       WRITE(3,50)(X(I),I=1,N)
       WRITE(3,50)(Y(J),J=1,M)
       WRITE(3,*)'HX HY FOLLOW'
       WRITE(3,50)(HX(I),I=1,N-1)
       WRITE(3,50)(HY(J),J=1,M-1)
       CALL XDERIV(N,M,HX,HY,DDX,DDY,XPART)
       CALL YDERIV(N,M,HX,HY,DDX,DDY,YPART)
       WRITE(3,*)'DDX AND XPART FOLLOW   SUBSCRIPT I CHANGES FASTEST'
       DO 56 J=1,M
       WRITE(3,50)(DDX(I,J),I=1,N-1)
       WRITE(3,50)(XPART(I,J),I=1,N)
56     CONTINUE
       WRITE(3,*)'DDY AND YPART FOLLOW    SUBSCRIPT J CHANGES FASTEST'
       DO 55 I=1,N
       WRITE(3,50)(DDY(I,J),J=1,M-1)
       WRITE(3,50)(YPART(I,J),J=1,M)
55     CONTINUE
       H=.25
       WRITE(3,*)'X       Y        S(X,Y)    F(X,Y)   ERROR FOLLOW'
 
       DO 4 I=1,17
          DO 4 J=1,17
              U=FLOAT(I-1)*H
              V=FLOAT(J-1)*H
              Z=EVAL(X,Y,F,XPART,YPART,N,M,U,V,IFLAG)
              WRITE(3,50)U,V,Z,FUNC(U,V),Z-FUNC(U,V)
4     CONTINUE
50    FORMAT(1X,5(E13.5,2X))
      END

./ ADD NAME=xderiv.f TIME=492785140
      SUBROUTINE XDERIV(N,M,HX,HY,DDX,DDY,XPART)
C THIS SUBROUTINE OBTAINS GOOD X PARTIALS SATISFYING THE SUFFICIENT
C CONDITIONS FOR MONOTONICITY AND APPROXIMATING THE X PARTIALS OF THE
C UNDERLYING FUNCTION.
C
C SUBROUTINE DIVDIF IS TO BE CALLED BEFORE THIS SUBROUTINE.
C
C
C INPUT      N     NUMBER OF DATA POINTS ON A X GRID LINE.
C            M     NUMBER OF DATA POINTS ON A Y GRID LINE.
C            HX    ARRAY OF X SUBINTERVAL LENGTHS. DIMENSION N.
C            HY    ARRAY OF Y SUBINTERVAL LENGTHS. DIMENSION M.
C            DDX   ARRAY OF FIRST DIVIDED DIFFERENCES IN THE X DIRECTION
C                  DIMENSION N BY M.
C            DDY   ARRAY OF FIRST DIVIDED DIFFERENCES IN THE Y DIRECTION
C                  DIMENSION N BY M.
C OUTPUT     XPART OUTPUT ARRAY OF APPROXIMATIONS TO THE FIRST PARTIALS
C                  IN THE X DIRECTION AT THE GRID POINTS. XPART IS
C                  FITTED BY THE ALGORITHM IN THE B.Z. PAPER. DIMENSION
C                  N BY M.
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
      INTRINSIC MAX,MIN
      DIMENSION HX(N),HY(M),DDX(N,M),DDY(N,M),XPART(N,M)
      ZERO=0.0D0
C
C
      N1=N-1
      M1=M-1
C
C
C SETTING UP THE INITIAL X PARTIALS ON THE LEFT AND RIGHT HAND EDGES
C OF THE ORIGINAL RECTANGLE
C
C
      N2=N-2
      H1=HX(1)
      H2=H1+HX(2)
      H3=HX(N1)
      H4=H3+HX(N2)
C
C
      DO 3 J=1,M
         U=DDX(1,J)
         V=(DDX(2,J)-U)/H2
         U=U-H1*V
         XPART(1,J)=MAX(ZERO,U)
         U=DDX(N1,J)
         V=(U-DDX(N2,J))/H4
         U=U+H3*V
         XPART(N,J)=MAX(ZERO,U)
    3 CONTINUE
C
C
C SETTING UP THE INITIAL X PARTIALS AT POINTS OTHER THAN THE LEFT AND
C RIGHT HAND EDGES.
C
C
      DO 4 I=2,N1
         I1=I-1
         H=HX(I1)
         HH=H+HX(I)
C
C NOTE BELOW THAT IF THE DATA IS MONOTONE INCREASING (=NON DECREASING)
C THEN THE INTERPOLATING QUADRATIC AT 3 POINTS HAS A NON-NEGATIVE
C DERIVATIVE AT THE SECOND POINT (OF THE THREE POINTS IN ASCENDING ORDER
C
         DO 4 J=1,M
              U=DDX(I1,J)
              V=(DDX(I,J)-U)/HH
              XPART(I,J)=U+H*V
    4 CONTINUE
C
C
C
C THIS PART OF THE CODE IMPLEMENTS STEP 2 OF THE ALGORITHM.
C RADIAL PROJECTION.
C
      DO 5 I=1,N1
         IP1=I+1
         DO 5 J=1,M
              TEMP1=3.0*DDX(I,J)
              TEMP2=XPART(I,J)+XPART(IP1,J)
              IF(TEMP2.GT.TEMP1)THEN
                   TEMP2=TEMP1/TEMP2
                   XPART(I,J)=TEMP2*XPART(I,J)
                   XPART(IP1,J)=TEMP2*XPART(IP1,J)
               END IF
    5 CONTINUE
C
C
C THIS PART OF THE CODE DOES THE DOWNWARD SWEEP - THAT IS STEP 3 OF
C THE ALGORITHM. SEE EQUATION (5D) OF B.Z.
C
C
      DO 6 I=1,N1
         DO 6 K=1,M1
              J=M-K
              JP1=J+1
              XPART(I,J)=MIN(XPART(I,J),XPART(I,JP1)+
     1                                       (HY(J)/HX(I))*DDY(I,J))
    6 CONTINUE
C
C
C THIS PART OF THE CODE DOES THE UPWARD SWEEP - THAT IS STEP 4 OF
C THE ALGORITHM. SEE EQUATION (5B) OF B.Z.
C
C
        DO 7 I=2,N
           DO 7 J=1,M1
              JP1=J+1
              XPART(I,JP1)=MIN(XPART(I,JP1),XPART(I,J)+
     1                                       (HY(J)/HX(I-1))*DDY(I,J))
    7 CONTINUE
      RETURN
      END

./ ADD NAME=yderiv.f TIME=492785141
      SUBROUTINE YDERIV(N,M,HX,HY,DDX,DDY,YPART)
C THIS SUBROUTINE OBTAINS GOOD Y PARTIALS SATISFYING THE SUFFICIENT
C CONDITIONS FOR MONOTONICITY AND APPROXIMATING THE Y PARTIALS OF THE
C UNDERLYING FUNCTION.
C
C
C SUBROUTINE DIVDIF IS TO BE CALLED BEFORE THIS SUBROUTINE.
C
C
C INPUT      N     NUMBER OF DATA POINTS ON A X GRID LINE.
C            M     NUMBER OF DATA POINTS ON A Y GRID LINE.
C            HX    ARRAY OF X SUBINTERVAL LENGTHS. DIMENSION N.
C            HY    ARRAY OF Y SUBINTERVAL LENGTHS. DIMENSION M.
C            DDX   ARRAY OF FIRST DIVIDED DIFFERENCES IN THE X DIRECTION
C                  DIMENSION N BY M.
C            DDY   ARRAY OF FIRST DIVIDED DIFFERENCES IN THE Y DIRECTION
C                  DIMENSION N BY M.
C OUTPUT    YPART  OUTPUT ARRAY OF APPROXIMATIONS TO THE FIRST PARTIALS
C                  IN THE Y DIRECTION AT THE GRID POINTS. YPART IS
C                  FITTED BY THE ALGORITHM IN THE B.Z. PAPER. DIMENSION
C                  N BY M.
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
      INTRINSIC MAX,MIN
      DIMENSION HX(N),HY(M),DDX(N,M),DDY(N,M),YPART(N,M)
C
C
      ZERO=0.0D0
      N1=N-1
      M1=M-1
C
C
C SETTING UP THE INITIAL Y PARTIALS ON THE TOP AND BOTTOM EDGES
C OF THE ORIGINAL RECTANGLE
C
C
      M2=M-2
      H1=HY(1)
      H2=H1+HY(2)
      H3=HY(M1)
      H4=H3+HY(M2)
C
C
      DO 3 I=1,N
         U=DDY(I,1)
         V=(DDY(I,2)-U)/H2
         U=U-H1*V
         YPART(I,1)=MAX(ZERO,U)
         U=DDY(I,M1)
         V=(U-DDY(I,M2))/H4
         U=U+H3*V
         YPART(I,M)=MAX(ZERO,U)
    3 CONTINUE
C
C
C SETTING UP THE INITIAL YPARTIALS AT POINTS OTHER THAN THOSE ON THE
C BOTTOM AND TOP EDGES.
C
C
      DO 4 J=2,M1
         J1=J-1
         H=HY(J1)
         HH=H+HY(J)
C
C NOTE BELOW THAT IF THE DATA IS MONOTONE INCREASING (=NON DECREASING)
C THEN THE INTERPOLATING QUADRATIC AT 3 POINTS HAS A NON-NEGATIVE
C DERIVATIVE AT THE SECOND POINT (OF THE THREE POINTS IN ASCENDING ORDER
C
         DO 4 I=1,N
              U=DDY(I,J1)
              V=(DDY(I,J)-U)/HH
              YPART(I,J)=U+H*V
    4 CONTINUE
C
C
C
C THIS PART OF THE CODE IMPLEMENTS STEP 2 OF THE ALGORITHM.
C RADIAL PROJECTION.
C
      DO 5 J=1,M1
         JP1=J+1
         DO 5 I=1,N
              TEMP1=3.0*DDY(I,J)
              TEMP2=YPART(I,J)+YPART(I,JP1)
              IF(TEMP2.GT.TEMP1)THEN
                  TEMP2=TEMP1/TEMP2
                  YPART(I,J)=TEMP2*YPART(I,J)
                  YPART(I,JP1)=TEMP2*YPART(I,JP1)
              END IF
    5 CONTINUE
C
C
C THIS PART OF THE CODE DOES THE RIGHT TO LEFT SWEEP. SEE EQUATION
C (5A) IN B.Z.
C
C
      DO 6 J=1,M1
         DO 6 K=1,N1
              I=N-K
              IP1=I+1
              YPART(I,J)=MIN(YPART(I,J),YPART(IP1,J)+
     1                                    (HX(I)/HY(J))*DDX(I,J))
    6 CONTINUE
C
C
C THIS PART OF THE CODE DOES THE LEFT TO RIGHT SWEEP. SEE EQUATION (5C)
C IN B.Z.
C
C
       DO 7 J=2,M
         DO 7 I=1,N1
              IP1=I+1
              YPART(IP1,J)=MIN(YPART(IP1,J),YPART(I,J)+
     1                                        (HX(I)/HY(J-1))*DDX(I,J))
    7 CONTINUE
      RETURN
      END

./ ENDUP

.