./ 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
ARE IN PROGRESS. THERE ALSO EXISTS AN EVALUATOR WHICH IS MORE MEMORY
OF THE QUADRATIC ON EACH SUB-TRIANGLE IN A HUGE ARRAY.
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
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
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
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)
DO 2 J=1,M
Y(J)=FLOAT(2*J-2)
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)
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)
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)
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)
STOP
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)
DO 2 J=1,M
Y(J)=FLOAT(J-1)
DO 3 I=1,N
DO 3 J=1,M
F(I,J)=FUNC(X(I),Y(J))
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)
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)
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)
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
.
