C ALGORITHM 791, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 25,NO. 1, March, 1999, P. 74-- 77. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc/ # Doc/PackingList # Fortran77/ # Fortran77/Drivers/ # Fortran77/Drivers/Sp/ # Fortran77/Drivers/Sp/data1 # Fortran77/Drivers/Sp/driver1.f # Fortran77/Drivers/Sp/driver2.f # Fortran77/Drivers/Sp/res1 # Fortran77/Drivers/Sp/res1_2 # Fortran77/Drivers/Sp/res2 # Fortran77/Drivers/Sp/src_792.f # Fortran77/Src/ # Fortran77/Src/Sp/ # Fortran77/Src/Sp/src.f # This archive created: Tue Jul 27 14:09:45 1999 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'PackingList' then echo shar: will not over-write existing file "'PackingList'" else cat << SHAR_EOF > 'PackingList' Algorithm:791 Language:Fortran77 Precision:Sp Src:src.f Driver_0:driver1.f src_792.f @Src DriverLib_0:port Data_0:stdin=data1 Res_0:stdout=res1 Driver_1:driver2.f @Src DriverLib_1:port SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << SHAR_EOF > 'data1' 1 1 10 17 30 6 2 1 10 17 30 6 3 1 10 17 24 3 4 1 10 17 30 6 5 1 10 17 30 5 SHAR_EOF fi # end of overwriting check if test -f 'driver1.f' then echo shar: will not over-write existing file "'driver1.f'" else cat << SHAR_EOF > 'driver1.f' C C TS2TST C 11/20/98 C C This program computes interpolation errors using the C scattered data package TSHEP2D for each of ten test C functions and a 33 by 33 uniform grid of interpolation C points in the unit square. C C This program uses Subroutines TESTDT and TSTFN2 from C ACM Algorithm SURVEY to generate a node set and and the C test function values. C INTEGER NMAX, NRMAX, NI PARAMETER (NMAX=100, NRMAX=10, NI=33) C C Array storage: C DOUBLE PRECISION X(NMAX), Y(NMAX), W(NMAX), RW(NMAX), . A(10,NMAX), P(NI), FT(NI,NI) INTEGER LCELL(NRMAX,NRMAX), LNEXT(NMAX) C DOUBLE PRECISION DEL, DUM, DX, DY, ERMAX, ERMEAN, PW, . RMAX, SSA, SSE, SSM, SUM, SX, SY, . XMIN, YMIN DOUBLE PRECISION TS2VAL INTEGER I, IER, J, K, KF, KFF, KFL, KS, LOUT, . N, NC, NFUN, NP, NR, NSET, NW, NWMAX C C Data: C C LOUT = Logical unit number for output. C NSET = Number of node sets. C NFUN = Number of test functions. C DATA LOUT/2/, NSET/5/, NFUN/10/ OPEN (LOUT,FILE='res1_2') C C Input format: C 100 FORMAT (I2) C C Get a user-specified node set number KS. C 1 WRITE (*,110) NSET 110 FORMAT (///13X,'TS2TST: TSHEP2D Test Program'// . 5X,'Specify a node set number in the range 1', . ' to ',I2,':'/) READ (*,100,ERR=1) KS IF (KS .LT. 1 .OR. KS .GT. NSET) GO TO 1 C C Copy N and the nodal coordinates for node set KS. C CALL TESTDT (KS, N,X,Y) IF (N .LT. 10 .OR. N .GT. NMAX) GO TO 20 C C Allow the user to specify a range of function numbers. C WRITE (*,120) NFUN 120 FORMAT (//5X,'Specify the first test function ', . '(1 to ',I2,'):'/) READ (*,100,ERR=1) KFF IF (KFF .LT. 1 .OR. KFF .GT. NFUN) GO TO 1 C WRITE (*,130) KFF, NFUN 130 FORMAT (//5X,'Specify the last test function (', . I2,' to ',I2,'):'/) READ (*,100,ERR=1) KFL IF (KFL .LT. KFF .OR. KFL .GT. NFUN) GO TO 1 NFUN = KFL-KFF+1 C C Input NC, NW, and NR from the console. C NWMAX = MIN(40,N-1) 2 WRITE (*,140) N 140 FORMAT (//5X,'N =',I4//5X, . 'Specify the number of nodes NC for the ', . 'least squares fits.'/5X,'NC = 18 is ', . 'recommended. NC GE 9.'// . 10X,'NC = '/) READ (*,100,ERR=2) NC IF (NC .LT. 9 .OR. NC .GT. NWMAX) GO TO 2 C 3 WRITE (*,150) 150 FORMAT (///5X,'Specify the number of nodes NW for ', . 'the weights. NW = 32 is'/5X,'recommended. ', . ' 1 LE NW LE MIN(40,N-1).'//10X,'NW = '/) READ (*,100,ERR=2) NW IF (1 .GT. NW .OR. NW .GT. NWMAX) GO TO 3 C 4 WRITE (*,160) NRMAX 160 FORMAT (///5X,'Specify the number of rows and column', . 's NR in the uniform grid'/5X,'of cells used', . ' to locate nearest neighbors. NR = Sqrt(N/', . '3) is'/5X,'recommended. 1 LE NR LE ',I2, . '.'//10X,'NR = '/) READ (*,100,ERR=3) NR IF (NR .LT. 1 .OR. NR .GT. NRMAX) GO TO 4 C C Set up uniform grid points. C DEL = 1./DBLE(NI-1) DO 5 I = 1,NI P(I) = DBLE(I-1)*DEL 5 CONTINUE C C Initialize the average SSE/SSM value to zero. C SSA = 0. C C Print a heading and loop on test functions. C WRITE (LOUT,200) KS, N, NI, NC, NW, NR DO 11 KF = KFF,KFL C C Compute true function values at the nodes. C DO 6 K = 1,N CALL TSTFN1 (KF,X(K),Y(K),0, W(K),DUM,DUM) 6 CONTINUE C C Compute true function values FT on the uniform grid, and C accumulate the sum of values SUM and sum of squared C values SSM. C SUM = 0. SSM = 0. DO 8 I = 1,NI DO 7 J = 1,NI CALL TSTFN1 (KF,P(I),P(J),0, FT(I,J),DUM,DUM) SUM = SUM + FT(I,J) SSM = SSM + FT(I,J)**2 7 CONTINUE 8 CONTINUE C C Compute the sum of squared deviations from the mean SSM. C SSM = SSM - SUM*SUM/DBLE(NI*NI) C C Compute parameters A and RW defining the interpolant. C CALL TSHEP2 (N,X,Y,W,NC,NW,NR, LCELL,LNEXT,XMIN, . YMIN,DX,DY,SX,SY,RMAX,RW,A,IER) IF (IER .NE. 0) GO TO 21 C C Compute interpolation errors. C ERMEAN = 0. ERMAX = 0. SSE = 0. DO 10 I = 1,NI DO 9 J = 1,NI PW = TS2VAL (P(I),P(J),N,X,Y,W,NR,LCELL,LNEXT, . XMIN,YMIN,DX,DY,SX,SY,RMAX,RW,A) - . FT(I,J) ERMEAN = ERMEAN + ABS(PW) ERMAX = MAX(ERMAX,ABS(PW)) SSE = SSE + PW*PW 9 CONTINUE 10 CONTINUE NP = NI*NI ERMEAN = ERMEAN/DBLE(NP) SSE = SSE/SSM SSA = SSA + SSE WRITE (LOUT,210) KF, ERMAX, ERMEAN, SSE 11 CONTINUE C C Print the average SSE/SSM value (averaged over the test C functions). C WRITE (LOUT,220) SSA/DBLE(NFUN) STOP C C N is outside its valid range. C 20 WRITE (LOUT,500) N, NMAX STOP C C Error in TSHEP2. C 21 IF (IER .EQ. 2) WRITE (LOUT,510) IF (IER .EQ. 3) WRITE (LOUT,520) STOP C C Print formats: C 200 FORMAT (30X,'TS2TST Output:'// . 1X,16X,'TSHEP2D Interpolation Errors for ', . 'N nodes and an'/ . 1X,6X,'NI by NI Uniform Grid of Interpolation', . 'ion Points in the Unit Square'//1X, . 6X,'Node set ',I2,4X,'N =',I4,4X,'NI = ',I2, . 4X,'NC = ',I2,4X,'NW = ',I2,4X,'NR = ',I2/// . 1X,16X,'Function',4X,'Max Error',4X, . 'Mean Error',4X,'SSE/SSM'/) 210 FORMAT (1X,19X,I2,9X,F7.4,6X,F8.5,2X,F9.6) 220 FORMAT (//1X,11X,'Average SSE/SSM over the test ', . 'functions = ',F9.6) C C Error message formats: C 500 FORMAT (///1X,10X,'*** Error in data -- N = ',I4, . ', Maximum value =',I4,' ***') 510 FORMAT (///1X,14X,'*** Error in TSHEP2 -- duplicate ', . 'nodes encountered ***') 520 FORMAT (///1X,14X,'*** Error in TSHEP2 -- all nodes ', . 'are collinear ***') END SHAR_EOF fi # end of overwriting check if test -f 'driver2.f' then echo shar: will not over-write existing file "'driver2.f'" else cat << SHAR_EOF > 'driver2.f' C C TS2TEST C 02/20/97 C C C This program tests the scattered data interpolation C package TSHEP2D by printing the maximum errors associated C with interpolated values and gradients on a 20 by 20 C uniform grid in the unit square. The data set consists C of 100 nodes with data values taken from a cosine series C function for which the method is exact. The ratio of C maximum interpolation error relative to the machine C precision is also printed. This should be O(1). The C interpolated values from TS2VAL, TS2GRD, and TS2HES are C compared for agreement. C INTEGER N, NR PARAMETER (N=100, NR=6) C C Array storage: C INTEGER LCELL(NR,NR), LNEXT(N) DOUBLE PRECISION X(N), Y(N), F(N), RW(N), A(10,N), . P(20) C INTEGER I, IER, J, K, LOUT, NC, NW DOUBLE PRECISION C1, C2, C3, CX, CY, CXX, CXY, CYY, . DX, DY, EC, ECX, ECY, ECXX, ECXY, . ECYY, EP1, EPS, PI, PX, PY, RC, RMAX, . SX, SY, XMIN, YMIN, YK, XX, YY DOUBLE PRECISION STORE, TS2VAL, FC, FX, FY, FXX, FXY, . FYY C C TSHEP2 parameters and logical unit for output: C DATA NC/17/, NW/30/, LOUT/1/ C C Cosine series test function and partial derivatives: C FC(XX,YY) = COS(2.*PI*XX)*COS(PI*YY) FX(XX,YY) = -2.*PI*SIN(2.*PI*XX)*COS(PI*YY) FY(XX,YY) = -PI*COS(2.*PI*XX)*SIN(PI*YY) FXX(XX,YY) = -4.*PI*PI*COS(2.*PI*XX)*COS(PI*YY) FXY(XX,YY) = 2.*PI*PI*SIN(2.*PI*XX)*SIN(PI*YY) FYY(XX,YY) = -PI*PI*COS(2.*PI*XX)*COS(PI*YY) C PI = ACOS(-1.D0) C C Output file: C OPEN (LOUT,FILE='res2') C C Generate a 10 by 10 grid of nodes in the unit square with C the natural ordering. C K = 0 DO 2 J = 1,10 YK = DBLE(10-J)/9. DO 1 I = 1,10 K = K + 1 X(K) = DBLE(I-1)/9. Y(K) = YK 1 CONTINUE 2 CONTINUE C C Compute the data values. C DO 3 K = 1,N F(K) = FC(X(K),Y(K)) 3 CONTINUE C C Compute parameters defining the interpolant C. C CALL TSHEP2 (N,X,Y,F,NC,NW,NR, LCELL,LNEXT,XMIN,YMIN, . DX,DY,SX,SY,RMAX,RW,A,IER) IF (IER .NE. 0) GO TO 10 C C Generate a 20 by 20 uniform grid of interpolation points C (P(I),P(J)) in the unit square. The four corners coin- C cide with nodes. C DO 4 I = 1,20 P(I) = DBLE(I-1)/19. 4 CONTINUE C C Compute the machine precision EPS. C EPS = 1. 5 EPS = EPS/2. EP1 = EPS + 1. IF (STORE(EP1) .GT. 1.) GO TO 5 EPS = EPS*2. C C Compute interpolation errors and test for agreement in the C C values returned by TS2VAL, TS2GRD, and TS2HES. C EC = 0. ECX = 0. ECY = 0. ECXX = 0. ECXY = 0. ECYY = 0. DO 7 J = 1,20 PY = P(J) DO 6 I = 1,20 PX = P(I) C1 = TS2VAL (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN, . YMIN,DX,DY,SX,SY,RMAX,RW,A) CALL TS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN, . YMIN,DX,DY,SX,SY,RMAX,RW,A, C2,CX, . CY,IER) CALL TS2HES (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN, . YMIN,DX,DY,SX,SY,RMAX,RW,A, C3,CX,CY, . CXX,CXY,CYY,IER) IF (IER .NE. 0) GO TO 11 IF (ABS(C1-C2) .GT. 5.*ABS(C1)*EPS .OR. . ABS(C1-C3) .GT. 5.*ABS(C1)*EPS) THEN C C Values returned by TS2VAL, TS2GRD, and TS2HES differ by a C relative amount greater than 3*EPS. C WRITE (LOUT,90) C1, C2, C3 90 FORMAT (///1X,'*** Error -- interpolated ', . 'values C1 (TS2VAL), C2 (TS2GRD), and'/ . 5X,'C3 (TS2HES) differ:'// . 5X,'C1 = ',D21.14/5X,'C2 = ',D21.14/ . 5X,'C3 = ',D21.14) ENDIF EC = MAX(EC,ABS(FC(PX,PY)-C1)) ECX = MAX(ECX,ABS(FX(PX,PY)-CX)) ECY = MAX(ECY,ABS(FY(PX,PY)-CY)) ECXX = MAX(ECXX,ABS(FXX(PX,PY)-CXX)) ECXY = MAX(ECXY,ABS(FXY(PX,PY)-CXY)) ECYY = MAX(ECYY,ABS(FYY(PX,PY)-CYY)) 6 CONTINUE 7 CONTINUE C C Print errors and the ratio EC/EPS. C RC = EC/EPS WRITE (LOUT,100) WRITE (LOUT,110) EC, RC WRITE (LOUT,120) ECX WRITE (LOUT,130) ECY WRITE (LOUT,140) ECXX WRITE (LOUT,150) ECXY WRITE (LOUT,160) ECYY STOP 100 FORMAT (///1X,'Maximum absolute errors in the ', . 'interpolant C and partial'/ . 1X,'derivatives CX, ..., CYY relative ', . 'to machine precision EPS'// . 1X,10X,'Function',3X,'Max Error',3X, . 'Max Error/EPS'/) 110 FORMAT (1X,13X,'C',7X,E9.3,7X,F5.2) 120 FORMAT (1X,13X,'CX',6X,E9.3) 130 FORMAT (1X,13X,'CY',6X,E9.3) 140 FORMAT (1X,13X,'CXX',5X,E9.3) 150 FORMAT (1X,13X,'CXY',5X,E9.3) 160 FORMAT (1X,13X,'CYY',5X,E9.3) C C Error in TSHEP2. C 10 WRITE (LOUT,200) IER STOP 200 FORMAT (///1X,'*** Error in TSHEP2 -- IER =',I2, . ' ***') C C Error in TS2GRD. C 11 WRITE (LOUT,210) IER STOP 210 FORMAT (///1X,'*** Error in TS2GRD -- IER =',I2, . ' ***') END DOUBLE PRECISION FUNCTION STORE (X) DOUBLE PRECISION X C C*********************************************************** C C From TRIPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 03/18/90 C C This function forces its argument X to be stored in a C memory location, thus providing a means of determining C floating point number characteristics (such as the machine C precision) when it is necessary to avoid computation in C high precision registers. C C C On input: C C X = Value to be stored. C C X is not altered by this function. C C On output: C C STORE = Value of X after it has been stored and C possibly truncated or rounded to the single C precision word length. C C Modules required by STORE: None C C*********************************************************** C DOUBLE PRECISION Y COMMON/STCOM/Y C Y = X STORE = Y RETURN END SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << SHAR_EOF > 'res1' TS2TST: TSHEP2D Test Program Specify a node set number in the range 1 to 5: Specify the first test function (1 to 10): Specify the last test function ( 1 to 10): N = 100 Specify the number of nodes NC for the least squares fits. NC = 18 is recommended. NC GE 9. NC = Specify the number of nodes NW for the weights. NW = 32 is recommended. 1 LE NW LE MIN(40,N-1). NW = Specify the number of rows and columns NR in the uniform grid of cells used to locate nearest neighbors. NR = Sqrt(N/3) is recommended. 1 LE NR LE 10. NR = SHAR_EOF fi # end of overwriting check if test -f 'res1_2' then echo shar: will not over-write existing file "'res1_2'" else cat << SHAR_EOF > 'res1_2' TS2TST Output: TSHEP2D Interpolation Errors for N nodes and an NI by NI Uniform Grid of Interpolationion Points in the Unit Square Node set 1 N = 100 NI = 33 NC = 17 NW = 30 NR = 6 Function Max Error Mean Error SSE/SSM 1 0.0361 0.00411 0.000586 2 0.0477 0.00208 0.002289 3 0.0088 0.00085 0.000294 4 0.0179 0.00145 0.001170 5 0.0069 0.00048 0.000166 6 0.0428 0.00417 0.010234 7 0.7777 0.08287 0.017348 8 0.1724 0.01898 0.004223 9 10.9061 0.65699 0.002766 10 0.2334 0.02018 0.021557 Average SSE/SSM over the test functions = 0.006063 SHAR_EOF fi # end of overwriting check if test -f 'res2' then echo shar: will not over-write existing file "'res2'" else cat << SHAR_EOF > 'res2' Maximum absolute errors in the interpolant C and partial derivatives CX, ..., CYY relative to machine precision EPS Function Max Error Max Error/EPS C 0.777E-15 3.50 CX 0.195E-12 CY 0.279E-12 CXX 0.171E-09 CXY 0.463E-10 CYY 0.228E-09 SHAR_EOF fi # end of overwriting check if test -f 'src_792.f' then echo shar: will not over-write existing file "'src_792.f'" else cat << SHAR_EOF > 'src_792.f' SUBROUTINE TESTDT (K, N,X,Y) DOUBLE PRECISION X(100), Y(100) INTEGER K, N C C*********************************************************** C C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 03/28/97 C C This subroutine returns one of five sets of nodes used C for testing scattered data fitting methods. All five sets C approximately cover the unit square [0,1] X [0,1]: the C convex hulls of sets 1 and 3 extend slightly outside the C square but do not completely cover it, those of sets 2 and C 5 coincide with the unit square, and the convex hull of C set 4 is a large subset of the unit square. C C On input: C C K = Integer in the range 1 to 5 which determines the C choice of data set as follows: C C K = 1 - Franke's 100-node set C K = 2 - Franke's 33-node set C K = 3 - Lawson's 25-node set C K = 4 - Random 100-node set C K = 5 - Gridded 81-node set C C X,Y = Arrays of length at least N(K), where C N(1) = 100, N(2) = 33, N(3) = 25, C N(4) = 100, and N(5) = 81. C C Input parameters are not altered by this routine. C C On output: C C N = Number of nodes in set K, or 0 if K is outside C its valid range. C C X,Y = Nodal coordinates of node set K. C C Subprograms required by TESTDT: None C C*********************************************************** C DOUBLE PRECISION X1(100), Y1(100), X2(33), Y2(33), . X3(25), Y3(25), X4(100), Y4(100), . X5(81), Y5(81) INTEGER I C C Node set 1: Franke's 100-node set. C DATA (X1(I),Y1(I), I = 1,20)/ . 0.0227035, -0.0310206, 0.0539888, 0.1586742, . 0.0217008, 0.2576924, 0.0175129, 0.3414014, . 0.0019029, 0.4943596, -0.0509685, 0.5782854, . 0.0395408, 0.6993418, -0.0487061, 0.7470194, . 0.0315828, 0.9107649, -0.0418785, 0.9962890, . 0.1324189, 0.0501330, 0.1090271, 0.0918555, . 0.1254439, 0.2592973, 0.0934540, 0.3381592, . 0.0767578, 0.4171125, 0.1451874, 0.5615563, . 0.0626494, 0.6552235, 0.1452734, 0.7524066, . 0.0958668, 0.9146523, 0.0695559, 0.9632421/ DATA (X1(I),Y1(I), I = 21,40)/ . 0.2645602, 0.0292939, 0.2391645, 0.0602303, . 0.2088990, 0.2668783, 0.2767329, 0.3696044, . 0.1714726, 0.4801738, 0.2266781, 0.5940595, . 0.1909212, 0.6878797, 0.1867647, 0.8185576, . 0.2304634, 0.9046507, 0.2426219, 0.9805412, . 0.3663168, 0.0396955, 0.3857662, 0.0684484, . 0.3832392, 0.2389548, 0.3179087, 0.3124129, . 0.3466321, 0.4902989, 0.3776591, 0.5199303, . 0.3873159, 0.6445227, 0.3812917, 0.8203789, . 0.3795364, 0.8938079, 0.2803515, 0.9711719/ DATA (X1(I),Y1(I), I = 41,60)/ . 0.4149771, -0.0284618, 0.4277679, 0.1560965, . 0.4200010, 0.2262471, 0.4663631, 0.3175094, . 0.4855658, 0.3891417, 0.4092026, 0.5084949, . 0.4792578, 0.6324247, 0.4812279, 0.7511007, . 0.3977761, 0.8489712, 0.4027321, 0.9978728, . 0.5848691, -0.0271948, 0.5730076, 0.1272430, . 0.6063893, 0.2709269, 0.5013894, 0.3477728, . 0.5741311, 0.4259422, 0.6106955, 0.6084711, . 0.5990105, 0.6733781, 0.5380621, 0.7235242, . 0.6096967, 0.9242411, 0.5026188, 1.0308762/ DATA (X1(I),Y1(I), I = 61,80)/ . 0.6616928, 0.0255959, 0.6427836, 0.0707835, . 0.6396475, 0.2008336, 0.6703963, 0.3259843, . 0.7001181, 0.4890704, 0.6333590, 0.5096324, . 0.6908947, 0.6697880, 0.6895638, 0.7759569, . 0.6718889, 0.9366096, 0.6837675, 1.0064516, . 0.7736939, 0.0285374, 0.7635332, 0.1021403, . 0.7410424, 0.1936581, 0.8258981, 0.3235775, . 0.7306034, 0.4714228, 0.8086609, 0.6091595, . 0.8214531, 0.6685053, 0.7290640, 0.8022808, . 0.8076643, 0.8476790, 0.8170951, 1.0512371/ DATA (X1(I),Y1(I), I = 81,100)/ . 0.8424572, 0.0380499, 0.8684053, 0.0902048, . 0.8366923, 0.2083092, 0.9418461, 0.3318491, . 0.8478122, 0.4335632, 0.8599583, 0.5910139, . 0.9175700, 0.6307383, 0.8596328, 0.8144841, . 0.9279871, 0.9042310, 0.8512805, 0.9696030, . 1.0449820, -0.0120900, 0.9670631, 0.1334114, . 0.9857884, 0.2695844, 0.9676313, 0.3795281, . 1.0129299, 0.4396054, 0.9657040, 0.5044425, . 1.0019855, 0.6941519, 1.0359297, 0.7459923, . 1.0414677, 0.8682081, 0.9471506, 0.9801409/ C C Node set 2: Franke's 33-node set. C DATA (X2(I),Y2(I), I = 1,33)/ . 0.05, 0.45, 0.00, 0.50, . 0.00, 1.00, 0.00, 0.00, . 0.10, 0.15, 0.10, 0.75, . 0.15, 0.30, 0.20, 0.10, . 0.25, 0.20, 0.30, 0.35, . 0.35, 0.85, 0.50, 0.00, . 0.50, 1.00, 0.55, 0.95, . 0.60, 0.25, 0.60, 0.65, . 0.60, 0.85, 0.65, 0.70, . 0.70, 0.20, 0.70, 0.65, . 0.70, 0.90, 0.75, 0.10, . 0.75, 0.35, 0.75, 0.85, . 0.80, 0.40, 0.80, 0.65, . 0.85, 0.25, 0.90, 0.35, . 0.90, 0.80, 0.95, 0.90, . 1.00, 0.00, 1.00, 0.50, . 1.00, 1.00/ C C Node set 3: Lawson's 25-node set. C DATA (X3(I),Y3(I), I = 1,25)/ . 0.13750, 0.97500, 0.91250, 0.98750, . 0.71250, 0.76250, 0.22500, 0.83750, . -0.05000, 0.41250, 0.47500, 0.63750, . 0.05000, -0.05000, 0.45000, 1.03750, . 1.08750, 0.55000, 0.53750, 0.80000, . -0.03750, 0.75000, 0.18750, 0.57500, . 0.71250, 0.55000, 0.85000, 0.43750, . 0.70000, 0.31250, 0.27500, 0.42500, . 0.45000, 0.28750, 0.81250, 0.18750, . 0.45000, -0.03750, 1.00000, 0.26250, . 0.50000, 0.46250, 0.18750, 0.26250, . 0.58750, 0.12500, 1.05000, -0.06125, . 0.10000, 0.11250/ C C Node set 4: Random 100-node set. C DATA (X4(I),Y4(I), I = 1,20)/ . 0.0096326, 0.3083158, 0.0216348, 0.2450434, . 0.0298360, 0.8613847, 0.0417447, 0.0977864, . 0.0470462, 0.3648355, 0.0562965, 0.7156339, . 0.0646857, 0.5311312, 0.0740377, 0.9755672, . 0.0873907, 0.1781117, 0.0934832, 0.5452797, . 0.1032216, 0.1603881, 0.1110176, 0.7837139, . 0.1181193, 0.9982015, 0.1251704, 0.6910589, . 0.1327330, 0.1049580, 0.1439536, 0.8184662, . 0.1564861, 0.7086405, 0.1651043, 0.4456593, . 0.1786039, 0.1178342, 0.1886405, 0.3189021/ DATA (X4(I),Y4(I), I = 21,40)/ . 0.2016706, 0.9668446, 0.2099886, 0.7571834, . 0.2147003, 0.2016598, 0.2204141, 0.3232444, . 0.2343715, 0.4368583, 0.2409660, 0.8907869, . 0.2527740, 0.0647260, 0.2570839, 0.5692618, . 0.2733365, 0.2947027, 0.2853833, 0.4332426, . 0.2901755, 0.3347464, 0.2964854, 0.7436284, . 0.3019725, 0.1066265, 0.3125695, 0.8845357, . 0.3307163, 0.5158730, 0.3378504, 0.9425637, . 0.3439061, 0.4799701, 0.3529922, 0.1783069, . 0.3635507, 0.1146760, 0.3766172, 0.8225797/ DATA (X4(I),Y4(I), I = 41,60)/ . 0.3822429, 0.2270688, 0.3869838, 0.4073598, . 0.3973137, 0.8875080, 0.4170708, 0.7631616, . 0.4255588, 0.9972804, 0.4299218, 0.4959884, . 0.4372839, 0.3410421, 0.4705033, 0.2498120, . 0.4736655, 0.6409007, 0.4879299, 0.1058690, . 0.4940260, 0.5411969, 0.5055324, 0.0089792, . 0.5162593, 0.8784268, 0.5219219, 0.5515874, . 0.5348529, 0.4038952, 0.5483213, 0.1654023, . 0.5569571, 0.2965158, 0.5638611, 0.3660356, . 0.5784908, 0.0366554, 0.5863950, 0.9502420/ DATA (X4(I),Y4(I), I = 61,80)/ . 0.5929148, 0.2638101, 0.5987839, 0.9277386, . 0.6117561, 0.5377694, 0.6252296, 0.7374676, . 0.6331381, 0.4674627, 0.6399048, 0.9186109, . 0.6488972, 0.0416884, 0.6558537, 0.1291029, . 0.6677405, 0.6763676, 0.6814074, 0.8444238, . 0.6887812, 0.3273328, 0.6940896, 0.1893879, . 0.7061687, 0.0645923, 0.7160957, 0.0180147, . 0.7317445, 0.8904992, 0.7370798, 0.4160648, . 0.7462030, 0.4688995, 0.7566957, 0.2174508, . 0.7699998, 0.5734231, 0.7879347, 0.8853319/ DATA (X4(I),Y4(I), I = 81,100)/ . 0.7944014, 0.8018436, 0.8164468, 0.6388941, . 0.8192794, 0.8931002, 0.8368405, 0.1000558, . 0.8500993, 0.2789506, 0.8588255, 0.9082948, . 0.8646496, 0.3259159, 0.8792329, 0.8318747, . 0.8837536, 0.0508513, 0.8900077, 0.9708450, . 0.8969894, 0.5120548, 0.9044917, 0.2859716, . 0.9083947, 0.9581641, 0.9203972, 0.6183429, . 0.9347906, 0.3779934, 0.9434519, 0.4010423, . 0.9490328, 0.9478657, 0.9569571, 0.7425486, . 0.9772067, 0.8883287, 0.9983493, 0.5496750/ C C Node set 5: 9 by 9 uniform grid. C DATA (X5(I),Y5(I), I = 1,20)/ . 0.125, 0.000, 0.000, 0.125, . 0.000, 0.250, 0.000, 0.375, . 0.000, 0.500, 0.000, 0.625, . 0.000, 0.750, 0.000, 0.875, . 0.000, 1.000, 0.000, 0.000, . 0.125, 0.125, 0.125, 0.250, . 0.125, 0.375, 0.125, 0.500, . 0.125, 0.625, 0.125, 0.750, . 0.125, 0.875, 0.125, 1.000, . 0.250, 0.000, 0.250, 0.125/ DATA (X5(I),Y5(I), I = 21,40)/ . 0.250, 0.250, 0.250, 0.375, . 0.250, 0.500, 0.250, 0.625, . 0.250, 0.750, 0.250, 0.875, . 0.250, 1.000, 0.375, 0.000, . 0.375, 0.125, 0.375, 0.250, . 0.375, 0.375, 0.375, 0.500, . 0.375, 0.625, 0.375, 0.750, . 0.375, 0.875, 0.375, 1.000, . 0.500, 0.000, 0.500, 0.125, . 0.500, 0.250, 0.500, 0.375/ DATA (X5(I),Y5(I), I = 41,60)/ . 0.500, 0.500, 0.500, 0.625, . 0.500, 0.750, 0.500, 0.875, . 0.500, 1.000, 0.625, 0.000, . 0.625, 0.125, 0.625, 0.250, . 0.625, 0.375, 0.625, 0.500, . 0.625, 0.625, 0.625, 0.750, . 0.625, 0.875, 0.625, 1.000, . 0.750, 0.000, 0.750, 0.125, . 0.750, 0.250, 0.750, 0.375, . 0.750, 0.500, 0.750, 0.625/ DATA (X5(I),Y5(I), I = 61,81)/ . 0.750, 0.750, 0.750, 0.875, . 0.750, 1.000, 0.875, 0.000, . 0.875, 0.125, 0.875, 0.250, . 0.875, 0.375, 0.875, 0.500, . 0.875, 0.625, 0.875, 0.750, . 0.875, 0.875, 0.875, 1.000, . 1.000, 0.000, 1.000, 0.125, . 1.000, 0.250, 1.000, 0.375, . 1.000, 0.500, 1.000, 0.625, . 1.000, 0.750, 1.000, 0.875, . 1.000, 1.000/ C C Store node set K in (X,Y). C IF (K .EQ. 1) THEN DO 1 I = 1,100 X(I) = X1(I) Y(I) = Y1(I) 1 CONTINUE N = 100 ELSEIF (K .EQ. 2) THEN DO 2 I = 1,33 X(I) = X2(I) Y(I) = Y2(I) 2 CONTINUE N = 33 ELSEIF (K .EQ. 3) THEN DO 3 I = 1,25 X(I) = X3(I) Y(I) = Y3(I) 3 CONTINUE N = 25 ELSEIF (K .EQ. 4) THEN DO 4 I = 1,100 X(I) = X4(I) Y(I) = Y4(I) 4 CONTINUE N = 100 ELSEIF (K .EQ. 5) THEN DO 5 I = 1,81 X(I) = X5(I) Y(I) = Y5(I) 5 CONTINUE N = 81 ELSE N = 0 ENDIF RETURN END SUBROUTINE TSTFN1 (K,X,Y,IFLAG, F,FX,FY) INTEGER K, IFLAG DOUBLE PRECISION X, Y, F, FX, FY C C*********************************************************** C C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 10/14/98 C C This subroutine computes the value and, optionally, the C first partial derivatives of one of ten bivariate test C functions. The first six functions were chosen by Richard C Franke to test interpolation software (See the reference C below). The last four functions represent more chal- C lenging surface fitting problems. C C On input: C C K = Integer in the range 1 to 10 which determines C the choice of function as follows: C C K = 1 - Exponential C K = 2 - Cliff C K = 3 - Saddle C K = 4 - Gentle C K = 5 - Steep C K = 6 - Sphere C K = 7 - Trig C K = 8 - Synergistic Gaussian C K = 9 - Cloverleaf Asymmetric Peak/Valley C K = 10 - Cosine Peak C C Note that function 6 is only defined inside a circle of C radius 8/9 centered at (.5,.5). Thus, if (X-.5)**2 + C (Y-.5)**2 .GE. 64/81, the value (and partials if IFLAG=1) C are set to 0 for this function. Also, the first partial C derivatives of function 10 are not defined at (.5,.5) -- C again, zeros are returned. C C X,Y = Coordinates of the point at which the selected C function is to be evaluated. C C IFLAG = Derivative option indicator: C IFLAG = 0 if only a function value is C required. C IFLAG = 1 if both the function and its first C partial derivatives are to be C evaluated. C C Input parameters are not altered by this routine. C C On output: C C F = Value of function K at (X,Y). C C FX,FY = First partial derivatives of function K at C (X,Y) if IFLAG = 1, unaltered otherwise. C C Intrinsic functions called by TSTFN1: COS, EXP, SIN, C SQRT, TANH C C Reference: R. Franke, A Critical Comparison of Some C Methods for Interpolation of Scattered Data, C Naval Postgraduate School Technical Report, C NPS-53-79-003, 1979. C C*********************************************************** C DOUBLE PRECISION T1, T2, T3, T4 IF (K .LT. 1 .OR. K .GT. 10) RETURN GO TO (1,2,3,4,5,6,7,8,9,10), K C C Exponential: C 1 F = .75*EXP(-((9.*X-2.)**2 + (9.*Y-2.)**2)/4.) + . .75*EXP(-((9.*X+1.)**2)/49. - (9.*Y+1.)/10.) + . .5*EXP(-((9.*X-7.)**2 + (9.*Y-3.)**2)/4.) - . .2*EXP(-(9.*X-4.)**2 - (9.*Y-7.)**2) IF (IFLAG .NE. 1) RETURN T1 = EXP(-((9.*X-2.)**2 + (9.*Y-2.)**2)/4.) T2 = EXP(-((9.*X+1.)**2)/49. - (9.*Y+1.)/10.) T3 = EXP(-((9.*X-7.)**2 + (9.*Y-3.)**2)/4.) T4 = EXP(-(9.*X-4.)**2 - (9.*Y-7.)**2) FX = -3.375*(9.*X-2.)*T1 - (27./98.)*(9.*X+1.)*T2 . -2.25*(9.*X-7.)*T3 + 3.6*(9.*X-4.)*T4 FY = -3.375*(9.*Y-2.)*T1 - .675*T2 . -2.25*(9.*Y-3.)*T3 + 3.6*(9.*Y-7.)*T4 RETURN C C Cliff: C 2 F = (TANH(9.0*(Y-X)) + 1.0)/9.0 IF (IFLAG .NE. 1) RETURN T1 = 18.0*(Y-X) FX = -4.0/(EXP(T1) + 2.0 + EXP(-T1)) FY = -FX RETURN C C Saddle: C 3 F = (1.25 + COS(5.4*Y))/(6.0 + 6.0*(3.0*X-1.0)**2) IF (IFLAG .NE. 1) RETURN T1 = 5.4*Y T2 = 1.0 + (3.0*X-1.)**2 FX = -(3.0*X-1.0)*(1.25 + COS(T1))/(T2**2) FY = -.9*SIN(T1)/T2 RETURN C C Gentle: C 4 F = EXP(-5.0625*((X-.5)**2 + (Y-.5)**2))/3.0 IF (IFLAG .NE. 1) RETURN T1 = X - .5 T2 = Y - .5 T3 = -3.375*EXP(-5.0625*(T1**2 + T2**2)) FX = T1*T3 FY = T2*T3 RETURN C C Steep: C 5 F = EXP(-20.25*((X-.5)**2 + (Y-.5)**2))/3.0 IF (IFLAG .NE. 1) RETURN T1 = X - .5 T2 = Y - .5 T3 = -13.5*EXP(-20.25*(T1**2 + T2**2)) FX = T1*T3 FY = T2*T3 RETURN C C Sphere: C 6 T4 = 64.0 - 81.0*((X-.5)**2 + (Y-.5)**2) F = 0. IF (T4 .GE. 0.) F = SQRT(T4)/9.0 - .5 IF (IFLAG .NE. 1) RETURN T1 = X - .5 T2 = Y - .5 T3 = 0. IF (T4 .GT. 0.) T3 = -9.0/SQRT(T4) FX = T1*T3 FY = T2*T3 RETURN C C Trig: C 7 F = 2.0*COS(10.0*X)*SIN(10.0*Y) + SIN(10.0*X*Y) IF (IFLAG .NE. 1) RETURN T1 = 10.0*X T2 = 10.0*Y T3 = 10.0*COS(10.0*X*Y) FX = -20.0*SIN(T1)*SIN(T2) + T3*Y FY = 20.0*COS(T1)*COS(T2) + T3*X RETURN C C Gaussx(1,.5,.1) + Gaussy(.75,.5,.1) + Gaussx(1,.5,.1)* C Gaussy(.75,.5,.1), where Gaussx(a,b,c) is the Gaussian C function of x with amplitude a, center (mean) b, and C width (standard deviation) c. C 8 T1 = 5.0 - 10.0*X T2 = 5.0 - 10.0*Y T3 = EXP(-.5*T1*T1) T4 = EXP(-.5*T2*T2) F = T3 + .75*T4*(1.0+T3) IF (IFLAG .NE. 1) RETURN FX = T1*T3*(10.0 + 7.5*T4) FY = T2*T4*(7.5 + 7.5*T3) RETURN C C Cloverleaf Asymmetric Hill/Valley: C 9 T1 = EXP((10.0 - 20.0*X)/3.0) T2 = EXP((10.0 - 20.0*Y)/3.0) T3 = 1.0/(1.0 + T1) T4 = 1.0/(1.0 + T2) F = ((20.0/3.0)**3 * T1*T2)**2 * (T3*T4)**5 * . (T1-2.0*T3)*(T2-2.0*T4) IF (IFLAG .NE. 1) RETURN FX = ((20.0/3.0)*T1)**2 * ((20.0/3.0)*T3)**5 * . (2.0*T1-3.0*T3-5.0+12.0*T3*T3)*T2*T2*T4**5 * . (T2-2.0*T4) FY = ((20.0/3.0)*T1)**2 * ((20.0/3.0)*T3)**5 * . (2.0*T2-3.0*T4-5.0+12.0*T4*T4)*T2*T2*T4**5 * . (T1-2.0*T3) RETURN C C Cosine Peak: C 10 T1 = SQRT( (80.0*X - 40.0)**2 + (90.0*Y - 45.0)**2 ) T2 = EXP(-.04*T1) T3 = COS(.15*T1) F = T2*T3 IF (IFLAG .NE. 1) RETURN T4 = SIN(.15*T1) FX = 0. FY = 0. IF (T1 .EQ. 0.) RETURN T4 = SIN(.15*T1) FX = -T2*(12.0*T4 + 3.2*T3)*(80.0*X - 40.0)/T1 FY = -T2*(13.5*T4 + 3.6*T3)*(90.0*Y - 45.0)/T1 RETURN END SUBROUTINE TSTFN2 (K,X,Y,IFLAG, F,FX,FY,FXX,FXY,FYY) INTEGER K, IFLAG DOUBLE PRECISION X, Y, F, FX, FY, FXX, FXY, FYY C C*********************************************************** C C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 10/14/98 C C This subroutine computes the value and, optionally, the C first and/or second partial derivatives of one of ten C bivariate test functions. The first six functions were C chosen by Richard Franke to test interpolation software. C (See the reference below.) The last four functions repre- C sent more challenging surface fitting problems. C C On input: C C K = Integer in the range 1 to 10 which determines C the choice of function as follows: C C K = 1 - Exponential C K = 2 - Cliff C K = 3 - Saddle C K = 4 - Gentle C K = 5 - Steep C K = 6 - Sphere C K = 7 - Trig C K = 8 - Synergistic Gaussian C K = 9 - Cloverleaf Asymmetric Peak/Valley C K = 10 - Cosine Peak C C Note that function 6 is only defined inside a circle of C radius 8/9 centered at (.5,.5). Thus, if (X-.5)**2 + C (Y-.5)**2 .GE. 64/81, the value (and partials if IFLAG=1) C are set to 0 for this function. Also, the first partial C derivatives of function 10 are not defined at (.5,.5) -- C again, zeros are returned. C C X,Y = Coordinates of the point at which the selected C function is to be evaluated. C C IFLAG = Derivative option indicator: C IFLAG = 0 if only a function value is C required. C IFLAG = 1 if both the function and its first C partial derivatives are to be C evaluated. C IFLAG = 2 if the function, its first partial C derivatives, and its second partial C derivatives are to be evaluated. C C Input parameters are not altered by this routine. C C On output: C C F = Value of function K at (X,Y). C C FX,FY = First partial derivatives of function K at C (X,Y) if IFLAG >= 1, unaltered otherwise. C C FXX,FXY,FYY = Second partial derivatives of function C K at (X,Y) if IFLAG >= 2, unaltered C otherwise. C C Intrinsic functions called by TSTFN2: COS, EXP, SIN, C SQRT, TANH C C Reference: R. Franke, A Critical Comparison of Some C Methods for Interpolation of Scattered Data, C Naval Postgraduate School Technical Report, C NPS-53-79-003, 1979. C C*********************************************************** C DOUBLE PRECISION T1, T2, T3, T4, T5, T6 IF (K .LT. 1 .OR. K .GT. 10) RETURN GO TO (1,2,3,4,5,6,7,8,9,10), K C C Exponential: C 1 F = .75*EXP(-((9.*X-2.)**2 + (9.*Y-2.)**2)/4.) + . .75*EXP(-((9.*X+1.)**2)/49. - (9.*Y+1.)/10.) + . .5*EXP(-((9.*X-7.)**2 + (9.*Y-3.)**2)/4.) - . .2*EXP(-(9.*X-4.)**2 - (9.*Y-7.)**2) IF (IFLAG .LT. 1) RETURN T1 = EXP(-((9.*X-2.)**2 + (9.*Y-2.)**2)/4.) T2 = EXP(-((9.*X+1.)**2)/49. - (9.*Y+1.)/10.) T3 = EXP(-((9.*X-7.)**2 + (9.*Y-3.)**2)/4.) T4 = EXP(-(9.*X-4.)**2 - (9.*Y-7.)**2) FX = -3.375*(9.*X-2.)*T1 - (27./98.)*(9.*X+1.)*T2 . -2.25*(9.*X-7.)*T3 + 3.6*(9.*X-4.)*T4 FY = -3.375*(9.*Y-2.)*T1 - .675*T2 . -2.25*(9.*Y-3.)*T3 + 3.6*(9.*Y-7.)*T4 IF (IFLAG .LT. 2) RETURN FXX = 15.1875*((9.*X-2.)**2 - 2.)*T1 + . 60.75*((9.*X+1.)**2 - 24.5)*T2 + . 10.125*((9.*X-7.)**2 - 2.)*T3 - . 64.8*((9.*X-4.)**2 - .5)*T4 FXY = 15.1875*(9.*X-2.)*(9.*Y-2.)*T1 + . (243./980.)*(9.*X+1.)*T2 + . 10.125*(9.*X-7.)*(9.*Y-3.)*T3 - . 64.8*(9.*X-4.)*(9.*Y-7.)*T4 FYY = 15.1875*((9.*Y-2.)**2 - 2.)*T1 + . .6075*T2 + . 10.125*((9.*Y-3.)**2 - 2.)*T3 - . 64.8*((9.*Y-7.)**2 - .5)*T4 RETURN C C Cliff: C 2 F = (TANH(9.0*(Y-X)) + 1.0)/9.0 IF (IFLAG .LT. 1) RETURN T1 = 18.0*(Y-X) FX = -4.0/(EXP(T1) + 2.0 + EXP(-T1)) FY = -FX IF (IFLAG .LT. 2) RETURN FXX = 18.0*TANH(0.5*T1)*FX FXY = -FXX FYY = FXX RETURN C C Saddle: C 3 F = (1.25 + COS(5.4*Y))/(6.0 + 6.0*(3.0*X-1.0)**2) IF (IFLAG .LT. 1) RETURN T1 = 5.4*Y T2 = 1.0 + (3.0*X-1.)**2 FX = -(3.0*X-1.0)*(1.25 + COS(T1))/(T2**2) FY = -.9*SIN(T1)/T2 IF (IFLAG .LT. 2) RETURN FXX = 3.0*(1.25 + COS(T1))*(3.0*T2-4.0)/(T2**3) FXY = 5.4*(3.0*X-1.0)*SIN(T1)/(T2**2) FYY = -4.86*COS(T1)/T2 RETURN C C Gentle: C 4 F = EXP(-5.0625*((X-.5)**2 + (Y-.5)**2))/3.0 IF (IFLAG .LT. 1) RETURN T1 = X - .5 T2 = Y - .5 T3 = -3.375*EXP(-5.0625*(T1**2 + T2**2)) FX = T1*T3 FY = T2*T3 IF (IFLAG .LT. 2) RETURN FXX = (1.0 - 10.125*T1*T1)*T3 FXY = -10.125*T1*T2*T3 FYY = (1.0 - 10.125*T2*T2)*T3 RETURN C C Steep: C 5 F = EXP(-20.25*((X-.5)**2 + (Y-.5)**2))/3.0 IF (IFLAG .LT. 1) RETURN T1 = X - .5 T2 = Y - .5 T3 = -13.5*EXP(-20.25*(T1**2 + T2**2)) FX = T1*T3 FY = T2*T3 IF (IFLAG .LT. 2) RETURN FXX = (1.0 - 40.5*T1*T1)*T3 FXY = -40.5*T1*T2*T3 FYY = (1.0 - 40.5*T2*T2)*T3 RETURN C C Sphere: C 6 T4 = 64.0 - 81.0*((X-.5)**2 + (Y-.5)**2) F = 0. IF (T4 .GE. 0.) F = SQRT(T4)/9.0 - .5 IF (IFLAG .LT. 1) RETURN T1 = X - .5 T2 = Y - .5 T3 = 0. IF (T4 .GT. 0.) T3 = -9.0/SQRT(T4) FX = T1*T3 FY = T2*T3 IF (IFLAG .LT. 2) RETURN FXX = (1.0 + FX*FX)*T3 FXY = FX*FY FYY = (1.0 + FY*FY)*T3 RETURN C C Trig: C 7 F = 2.0*COS(10.0*X)*SIN(10.0*Y) + SIN(10.0*X*Y) IF (IFLAG .LT. 1) RETURN T1 = 10.0*X T2 = 10.0*Y T3 = 10.0*COS(10.0*X*Y) FX = -20.0*SIN(T1)*SIN(T2) + T3*Y FY = 20.0*COS(T1)*COS(T2) + T3*X IF (IFLAG .LT. 2) RETURN T4 = 100.0*SIN(10.0*X*Y) FXX = -200.0*COS(T1)*SIN(T2) - T4*Y*Y FXY = -200.0*SIN(T1)*COS(T2) + T3 - T4*X*Y FYY = -200.0*COS(T1)*SIN(T2) - T4*X*X RETURN C C Gaussx(1,.5,.1) + Gaussy(.75,.5,.1) + Gaussx(1,.5,.1)* C Gaussy(.75,.5,.1), where Gaussx(a,b,c) is the Gaussian C function of x with amplitude a, center (mean) b, and C width (standard deviation) c. C 8 T1 = 5.0 - 10.0*X T2 = 5.0 - 10.0*Y T3 = EXP(-.5*T1*T1) T4 = EXP(-.5*T2*T2) F = T3 + .75*T4*(1.0+T3) IF (IFLAG .LT. 1) RETURN FX = T1*T3*(10.0 + 7.5*T4) FY = T2*T4*(7.5 + 7.5*T3) IF (IFLAG .LT. 2) RETURN FXX = T3*(T1*T1-1.0)*(100.0 + 75.0*T4) FXY = 75.0*T1*T2*T3*T4 FYY = T4*(T2*T2-1.0)*(75.0 + 75.0*T3) RETURN C C Cloverleaf Asymmetric Hill/Valley: C 9 T1 = EXP((10.0 - 20.0*X)/3.0) T2 = EXP((10.0 - 20.0*Y)/3.0) T3 = 1.0/(1.0 + T1) T4 = 1.0/(1.0 + T2) T5 = 20.0/3.0 F = (T5**3 * T1*T2)**2 * (T3*T4)**5 * . (T1-2.0*T3)*(T2-2.0*T4) IF (IFLAG .LT. 1) RETURN T6 = (T5*T1*T2)**2 * (T5*T3*T4)**5 FX = T6 * (T2-2.0*T4) * . ((12.0*T3 - 3.0)*T3 + 2.0*T1 - 5.0) FY = T6 * (T1-2.0*T3) * . ((12.0*T4 - 3.0)*T4 + 2.0*T2 - 5.0) IF (IFLAG .LT. 2) RETURN FXX = T5*T6 * (T2-2.0*T4) * . (((-84.0*T3 + 78.0)*T3 + 23.0)*T3 + 4.0*T1-25.0) FXY = T5*T6 * . ((12.0*T4 - 3.0)*T4 + 2.0*T2 - 5.0) * . ((12.0*T3 - 3.0)*T3 + 2.0*T1 - 5.0) FYY = T5*T6 * (T1-2.0*T3) * . (((-84.0*T4 + 78.0)*T4 + 23.0)*T4 + 4.0*T2-25.0) RETURN C C Cosine Peak: C 10 T1 = SQRT( (80.0*X - 40.0)**2 + (90.0*Y - 45.0)**2 ) T2 = EXP(-.04*T1) T3 = COS(.15*T1) F = T2*T3 IF (IFLAG .LT. 1) RETURN T4 = SIN(.15*T1) FX = 0. FY = 0. IF (T1 .EQ. 0.) RETURN T4 = SIN(.15*T1) FX = -T2*(12.0*T4 + 3.2*T3)*(80.0*X - 40.0)/T1 FY = -T2*(13.5*T4 + 3.6*T3)*(90.0*Y - 45.0)/T1 IF (IFLAG .LT. 2) RETURN FXX = 0. FXY = 0. FYY = 0. IF (T1 .EQ. 0.) RETURN T5 = T2/(T1**3) FXX = T5*(T1*(76.8*T4 - 133.76*T3)*(80.0*X-40.0)**2 - . (960.0*T4 + 256.0*T3)*(90.0*Y-45.0)**2 ) FXY = T5*(T1*(86.4*T4 - 150.48*T3) + 1080.0*T4 + . 288.0*T3)*(80.0*X-40.0)*(90.0*Y-45.0) FYY = T5*(T1*(97.2*T4 - 169.29*T3)*(90.0*Y-45.0)**2 - . (1215.0*T4 + 324.0*T3)*(80.0*X-40.0)**2) RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << SHAR_EOF > 'src.f' SUBROUTINE TSHEP2 (N,X,Y,F,NC,NW,NR, LCELL,LNEXT,XMIN, . YMIN,DX,DY,SX,SY,RMAX,RW,A,IER) INTEGER N, NC, NW, NR, LCELL(NR,NR), LNEXT(N), IER DOUBLE PRECISION X(N), Y(N), F(N), XMIN, YMIN, DX, . DY, SX, SY, RMAX, RW(N), A(10,N) C C*********************************************************** C C From TSHEP2D C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 03/23/97 C C This subroutine computes a set of parameters defining a C C2 (twice continuously differentiable) bivariate function C C(X,Y) which interpolates data values F at a set of N C arbitrarily distributed points (X,Y) in the plane (nodes). C The interpolant C may be evaluated at an arbitrary point C by function TS2VAL, and its first partial derivatives are C computed by Subroutine TS2GRD. C C The interpolation scheme is a modified Shepard method C based on a cosine series: C C C = [W(1)*C(1)+W(2)*C(2)+..+W(N)*C(N)]/[W(1)+W(2)+..+W(N)] C C for bivariate functions W(k) and C(k). The nodal func- C tions are given by C C C(k)(x,y) = A(1,k) + A(2,k)*cos(px) + A(3,k)*cos(py) + C A(4,k)*cos(2*px) + A(5,k)*cos(px)*cos(py) + C A(6,k)*cos(2*py) + A(7,k)*cos(3*px) + C A(8,k)*cos(2*px)*cos(py) + C A(9,k)*cos(px)*cos(2*py) + A(10,k)*cos(3*py), C C for px = ((x-XMIN)/(XMAX-XMIN))*Pi and C py = ((y-YMIN)/(YMAX-YMIN))*Pi , C C where [XMIN,XMAX] X [YMIN,YMAX] is the smallest rectangle C that contains the nodes. The method exactly reproduces C true cubic functions of cos(px) and cos(py). C C The coefficients A(,k) of C(k) are obtained by a weighted C least squares fit to the closest NC data points (along C with point k) with weights related to inverse distance. C Note that the radius of influence for the least squares C fit is fixed for each k, but varies with k. C C The weights are taken to be C C W(k)(x,y) = ( (R(k)-D(k))+ / R(k)*D(k) )**3 , C C where (R(k)-D(k))+ = 0 if R(k) < D(k), and D(k)(x,y) is C the Euclidean distance between (x,y) and (X(k),Y(k)). The C radius of influence R(k) varies with k and is chosen so C that NW nodes are within the radius. Note that W(k) is C not defined at node (X(k),Y(k)), but C(x,y) has limit F(k) C as (x,y) approaches (X(k),Y(k)). C C On input: C C N = Number of nodes and data values. N .GE. 10. C C X,Y = Arrays of length N containing the Cartesian C coordinates of the nodes. C C F = Array of length N containing the data values C in one-to-one correspondence with the nodes. C C NC = Number of data points (in addition to point K) C to be used in the least squares fit for coeffi- C cients defining the nodal functions C(k). C Values found to be optimal for test data sets C ranged from 14 to 19. The recommended value is C NC = 18, and NC must be in the range 9 to C Min(40,N-1). C C NW = Number of nodes within (and defining) the radii C of influence R(k) which enter into the weights C W(k). For N sufficiently large, a recommended C value is NW = 32. In general, NW should be at C least 1.5*NC. 1 .LE. NW .LE. Min(40,N-1). C C NR = Number of rows and columns in the cell grid de- C fined in Subroutine STORE2. A rectangle con- C taining the nodes is partitioned into cells in C order to increase search efficiency. NR = C Sqrt(N/3) is recommended. NR .GE. 1. C C The above parameters are not altered by this routine. C C LCELL = Array of length .GE. NR**2. C C LNEXT = Array of length .GE. N. C C RW = Array of length .GE. N. C C A = Array of length .GE. 10N. C C On output: C C LCELL = NR by NR array of nodal indexes associated C with cells. Refer to Subroutine STORE2. C C LNEXT = Array of length N containing next-node C indexes. Refer to Subroutine STORE2. C C XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell C dimensions. Refer to Subroutine C STORE2. C C SX,SY = Scale factors for mapping [XMIN,XMAX] and C [YMIN,YMAX] to [0,Pi]: C SX = Pi/(XMAX-XMIN) and SY = Pi/(YMAX-YMIN). C C RMAX = Largest element in RW -- maximum radius R(k). C C RW = Array containing the the radii R(k) which enter C into the weights W(k). C C A = 10 by N array containing the coefficients for C nodal function C(k) in column k. C C Note that the output parameters described above are not C defined unless IER = 0. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N, NC, NW, or NR is outside its C valid range. C IER = 2 if duplicate nodes were encountered. C IER = 3 if all nodes are collinear. C C Modules required by TSHEP2: GETNP2, GIVENS, ROTATE, C SETUPC, STORE2 C C Intrinsic functions called by TSHEP2: ABS, ACOS, DBLE, C MAX, MIN, SQRT C C*********************************************************** C INTEGER LMX PARAMETER (LMX=40) INTEGER I, IERR, IP1, IRM1, IROW, J, JP1, K, LMAX, . LNP, NEQ, NN, NNC, NNW, NP, NPTS(LMX), . NCWMAX DOUBLE PRECISION B(11,11), C, D, DMIN, DTOL, FK, PI, . RC, RS(0:LMX), RSMX, RTOL, RWS, S, . STF, T, W, WMAX, WSF, XK, XMAX, XNP, . YK, YMAX, YNP C DATA DTOL/.01/, RTOL/1.D-5/, WSF/1.D10/ C C Local parameters: C C B = Transpose of the augmented regression matrix C C = First component of the plane rotation used to C zero the lower triangle of B**T -- computed C by Subroutine GIVENS C D = Distance between nodes K and NP C DMIN = Minimum of the magnitudes of the diagonal C elements of the regression matrix after C zeros are introduced below the diagonal C DTOL = Tolerance for detecting an ill-conditioned C system. The system is accepted when C DMIN .GE. DTOL. C FK = Data value at mode K -- F(K) C I = Index for A, B, and NPTS C IERR = Error flag for the call to Subroutine STORE2 C IP1 = I+1 C IRM1 = IROW-1 C IROW = Row index for B**T C J = Index for A and B C JP1 = J+1 C K = Nodal function index and column index for A C LMAX = Maximum number of NPTS elements C LMX = Maximum value of LMAX (value of LMAX for N C sufficiently large C LNP = Current length of NPTS C NEQ = Number of equations in the least squares fit C NN,NNC,NNW = Local copies of N, NC, and NW C NP = NPTS element C NPTS = Array containing the indexes of a sequence of C nodes to be used (along with node K) in the C least squares fit associated with C(k), or C to compute RW. The nodes are ordered by C distance from K, and the last element C (usually indexed by LNP) is used only to C determine RC, or RW(K) if NW > NC. C NCWMAX = Max(NC,NW) C PI = Pi C RC = Radius of influence which enters into the C weights W for C(K) C RS = Array containing squared distances between C nodes K and NPTS(LNP) -- used to compute C RC and RW(K) C RSMX = Maximum squared RW element encountered C RTOL = Tolerance for detecting a sufficiently large C relative change in consecutive RS elements. C If the change is not greater than RTOL, the C nodes are treated as being the same C distance from K C RWS = Current squared value of RW(K) C STF = Marquardt stabilization factor used to damp C out the last four solution components when C the system is ill-conditioned. C T = Temporary variable for accumulating a scalar C product in the back solve C W = Weight associated with node NP in the least C squares system for nodal function K: W = C (RC-D)/(RC*D) = 1/D - 1/RC C WMAX = Weight associated with node K in the least C squares system for nodal function K. A C large weight is used to force interpolation C of the data value F(K). WMAX is the C largest W value scaled by WSF. C WSF = Scale factor for WMAX C XK,YK = Coordinates of node K -- X(K), Y(K) C XMAX,YMAX = Maximum nodal coordinates C XNP,YNP = Coordinates of node NP transformed by the C affine transformation that maps [XMIN,XMAX] C X [YMIN,YMAX] to [0,Pi] X [0,Pi] C NN = N NNC = NC NNW = NW NCWMAX = MAX(NNC,NNW) LMAX = MIN(LMX,NN-1) IF (NNC .LT. 9 .OR. NNW .LT. 1 .OR. NCWMAX .GT. . LMAX .OR. NR .LT. 1) GO TO 21 C C Create the cell data structure, compute XMAX, YMAX, SX, C and SY, and initialize RSMX. C CALL STORE2 (NN,X,Y,NR, LCELL,LNEXT,XMIN,YMIN,DX,DY, . IERR) XMAX = XMIN + DBLE(NR)*DX YMAX = YMIN + DBLE(NR)*DY PI = ACOS(-1.D0) SX = PI/(XMAX-XMIN) SY = PI/(YMAX-YMIN) IF (IERR .NE. 0) GO TO 23 RSMX = 0. C C Outer loop on node K: C DO 15 K = 1,NN XK = X(K) YK = Y(K) FK = F(K) C C Mark node K to exclude it from the search for nearest C neighbors. C LNEXT(K) = -LNEXT(K) C C Initialize for loop on NPTS elements. C RWS = 0. RC = 0. LNP = 0 RS(0) = 0. C C Compute NPTS, LNP, RWS, NEQ, and RC. C 1 IF (LNP .EQ. LMAX) GO TO 2 LNP = LNP + 1 CALL GETNP2 (XK,YK,X,Y,NR,LCELL,LNEXT,XMIN,YMIN, . DX,DY, NP,RS(LNP)) IF (RS(LNP) .EQ. 0.) GO TO 22 NPTS(LNP) = NP IF ((RS(LNP)-RS(LNP-1))/RS(LNP) .LT. RTOL) GO TO 1 IF (RWS .EQ. 0. .AND. LNP .GT. NNW) . RWS = RS(LNP) IF (RC .EQ. 0. .AND. LNP .GT. NNC) THEN C C RC = 0 (not yet computed) and LNP > NC. RC = C Sqrt(RS(LNP)) is sufficiently large to (strictly) C include NC nodes. The least squares fit will C include NEQ = LNP equations for 9 .LE. NC .LT. NEQ C .LE. LMAX .LE. N-1. C NEQ = LNP RC = SQRT(RS(LNP)) ENDIF C C Bottom of loop -- test for termination. C IF (LNP .GT. NCWMAX) GO TO 3 GO TO 1 C C All LMAX nodes are included in NPTS. RWS and/or RC**2 is C (arbitrarily) taken to be 10 percent larger than the C distance to the last node included. C 2 IF (RWS .EQ. 0.) RWS = 1.1*RS(LMAX) IF (RC .EQ. 0.) THEN NEQ = LMAX + 1 RC = SQRT(1.1*RS(LMAX)) ENDIF C C Store RW(K) and update RSMX if necessary. C 3 RW(K) = SQRT(RWS) IF (RWS .GT. RSMX) RSMX = RWS C C A Q-R decomposition is used to solve the least squares C system. The transpose of the augmented regression C matrix is stored in B. C C Store the equation associated with node K in the first C row. C WMAX = WSF*(1./SQRT(RS(1)) - 1./RC) XNP = SX*(XK-XMIN) YNP = SY*(YK-YMIN) CALL SETUPC (XNP,YNP,FK,WMAX, B(1,1)) C C Set up the remaining equations and zero out the lower C triangle with Givens rotations. C I = 0 4 I = I + 1 NP = NPTS(I) IROW = MIN(I+1,11) D = SQRT(RS(I)) W = 1./D - 1./RC XNP = SX*(X(NP)-XMIN) YNP = SY*(Y(NP)-YMIN) CALL SETUPC (XNP,YNP,F(NP),W, B(1,IROW)) IRM1 = IROW-1 DO 5 J = 1,IRM1 JP1 = J + 1 CALL GIVENS (B(J,J),B(J,IROW),C,S) CALL ROTATE (11-J,C,S,B(JP1,J),B(JP1,IROW)) 5 CONTINUE IF (I+1 .LT. NEQ) GO TO 4 C C Test the system for ill-conditioning. C DMIN = MIN( ABS(B(1,1)),ABS(B(2,2)),ABS(B(3,3)), . ABS(B(4,4)),ABS(B(5,5)),ABS(B(6,6)), . ABS(B(7,7)),ABS(B(8,8)),ABS(B(9,9)), . ABS(B(10,10)) ) IF (DMIN .GE. DTOL) GO TO 11 IF (NEQ .GT. LMAX) GO TO 7 C C Increase RC and add another equation to the system to C improve the conditioning. The number of NPTS elements C is also increased if necessary. C 6 NEQ = NEQ + 1 IF (NEQ .EQ. LMAX+1) THEN RC = SQRT(1.1*RS(LMAX)) GO TO 4 ENDIF IF (NEQ .LE. LNP) THEN C C NEQ .LE. LNP. C IF ((RS(NEQ)-RS(NEQ-1))/RS(NEQ) .LT. RTOL) GO TO 6 RC = SQRT(RS(NEQ)) GO TO 4 ENDIF C C NEQ = LNP+1. Add an element to NPTS. C LNP = LNP + 1 CALL GETNP2 (XK,YK,X,Y,NR,LCELL,LNEXT,XMIN,YMIN, . DX,DY, NP,RS(LNP)) IF (NP .EQ. 0) GO TO 22 NPTS(LNP) = NP IF ( (RS(LNP)-RS(LNP-1))/RS(LNP) .LT. RTOL ) GO TO 6 RC = SQRT(RS(LNP)) GO TO 4 C C Stabilize the system by damping the last four basis C functions -- add multiples of the last four unit C vectors to the last four equations. C 7 STF = 1.0/RC DO 10 I = 7,10 B(I,11) = STF IP1 = I + 1 DO 8 J = IP1,11 B(J,11) = 0. 8 CONTINUE DO 9 J = I,10 JP1 = J + 1 CALL GIVENS (B(J,J),B(J,11),C,S) CALL ROTATE (11-J,C,S,B(JP1,J),B(JP1,11)) 9 CONTINUE 10 CONTINUE C C Test the damped system for ill-conditioning. C DMIN = MIN( ABS(B(5,5)),ABS(B(6,6)),ABS(B(7,7)), . ABS(B(8,8)),ABS(B(9,9)),ABS(B(10,10)) ) IF (DMIN .LT. DTOL) GO TO 23 C C Solve the order-10 triangular system for the coefficients. C 11 DO 13 I = 10,1,-1 T = 0. IF (I .NE. 10) THEN IP1 = I + 1 DO 12 J = IP1,10 T = T + B(J,I)*A(J,K) 12 CONTINUE ENDIF A(I,K) = (B(11,I)-T)/B(I,I) 13 CONTINUE C C Unmark K and the elements of NPTS. C LNEXT(K) = -LNEXT(K) DO 14 I = 1,LNP NP = NPTS(I) LNEXT(NP) = -LNEXT(NP) 14 CONTINUE 15 CONTINUE C C No errors encountered. C RMAX = SQRT(RSMX) IER = 0 RETURN C C N, NC, NW, or NR is outside its valid range. C 21 IER = 1 RETURN C C Duplicate nodes were encountered by GETNP2. C 22 IER = 2 RETURN C C No unique solution due to collinear nodes. C 23 IER = 3 RETURN END DOUBLE PRECISION FUNCTION TS2VAL (PX,PY,N,X,Y,F,NR, . LCELL,LNEXT,XMIN,YMIN,DX,DY,SX,SY, . RMAX,RW,A) INTEGER N, NR, LCELL(NR,NR), LNEXT(N) DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN, . DX, DY, SX, SY, RMAX, RW(N), A(10,N) C C*********************************************************** C C From TSHEP2D C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 02/19/97 C C This function returns the value C(PX,PY), where C is the C weighted sum of nodal functions defined in Subroutine C TSHEP2. TS2GRD may be called to compute a gradient C of C along with the value, and/or to test for errors. C TS2HES may be called to compute a value, first partial C derivatives, and second partial derivatives at a point. C C On input: C C PX,PY = Cartesian coordinates of the point P at C which C is to be evaluated. C C N = Number of nodes and data values defining C. C N .GE. 10. C C X,Y,F = Arrays of length N containing the nodes and C data values interpolated by C. C C NR = Number of rows and columns in the cell grid. C Refer to Subroutine STORE2. NR .GE. 1. C C LCELL = NR by NR array of nodal indexes associated C with cells. Refer to Subroutine STORE2. C C LNEXT = Array of length N containing next-node C indexes. Refer to Subroutine STORE2. C C XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell C dimensions. DX and DY must be C positive. Refer to Subroutine C STORE2. C C SX,SY = Scale factors for mapping [XMIN,XMAX] and C [YMIN,YMAX] to [0,Pi]: C SX = Pi/(XMAX-XMIN) and SY = Pi/(YMAX-YMIN). C C RMAX = Largest element in RW -- maximum radius R(k). C C RW = Array containing the the radii R(k) which enter C into the weights W(k) defining C. C C A = 10 by N array containing the coefficients for C nodal function C(k) in column k. C C Input parameters are not altered by this function. The C parameters other than PX and PY should be input unaltered C from their values on output from TSHEP2. This function C should not be called if a nonzero error flag was returned C by TSHEP2. C C On output: C C TS2VAL = Function value C(PX,PY) unless N, NR, DX, C DY, or RMAX is invalid, in which case no C value is returned. C C Modules required by TS2VAL: NONE C C Intrinsic functions called by TS2VAL: COS, INT, SQRT C C*********************************************************** C INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP DOUBLE PRECISION CX1, CX2, CX3, CY1, CY2, CY3, D, R, . SW, SWC, W, XP, YP C C Local parameters: C C CX1 = Cos(XP) C CX2 = Cos(2*XP) C CX3 = Cos(3*XP) C CY1 = Cos(YP) C CY2 = Cos(2*YP) C CY3 = Cos(3*YP) C D = Distance between P and node K C I = Cell row index in the range IMIN to IMAX C IMIN,IMAX = Range of cell row indexes of the cells C intersected by a disk of radius RMAX C centered at P C J = Cell column index in the range JMIN to JMAX C JMIN,JMAX = Range of cell column indexes of the cells C intersected by a disk of radius RMAX C centered at P C K = Index of a node in cell (I,J) C KP = Previous value of K in the sequence of nodes C in cell (I,J) C R = Radius of influence for node K C SW = Sum of weights W(K) C SWC = Sum of weighted nodal function values at P C W = Weight W(K) value at P: ((R-D)+/(R*D))**3, C where (R-D)+ = 0 if R < D C XP,YP = Coordinates of node NP transformed by the C affine transformation that maps [XMIN,XMAX] C X [YMIN,YMAX] to [0,Pi] X [0,Pi] C XP = SX*(PX-XMIN) YP = SY*(PY-YMIN) IF (N .LT. 10 .OR. NR .LT. 1 .OR. DX .LE. 0. .OR. . DY .LE. 0. .OR. RMAX .LT. 0.) RETURN C C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining C the range of the search for nodes whose radii include C P. The cells which must be searched are those inter- C sected by (or contained in) a circle of radius RMAX C centered at P. C IMIN = INT((PX-XMIN-RMAX)/DX) + 1 IMAX = INT((PX-XMIN+RMAX)/DX) + 1 IF (IMIN .LT. 1) IMIN = 1 IF (IMAX .GT. NR) IMAX = NR JMIN = INT((PY-YMIN-RMAX)/DY) + 1 JMAX = INT((PY-YMIN+RMAX)/DY) + 1 IF (JMIN .LT. 1) JMIN = 1 IF (JMAX .GT. NR) JMAX = NR C C The following is a test for no cells within the circle C of radius RMAX. C IF (IMIN .GT. IMAX .OR. JMIN .GT. JMAX) GO TO 6 C C Compute trig function values at (XP,YP). C CX1 = COS(XP) CY1 = COS(YP) CX2 = COS(2.0*XP) CY2 = COS(2.0*YP) CX3 = COS(3.0*XP) CY3 = COS(3.0*YP) C C Accumulate weight values in SW and weighted nodal function C values in SWC. The weights are W(K) = ((R-D)+/(R*D))**3 C for R = RW(K) and D = distance between P and node K. C SW = 0. SWC = 0. C C Outer loop on cells (I,J). C DO 4 J = JMIN,JMAX DO 3 I = IMIN,IMAX K = LCELL(I,J) IF (K .EQ. 0) GO TO 3 C C Inner loop on nodes K. C 1 D = SQRT((PX-X(K))**2 + (PY-Y(K))**2) R = RW(K) IF (D .GE. R) GO TO 2 IF (D .EQ. 0.) GO TO 5 W = (1.0/D - 1.0/R)**3 SW = SW + W SWC = SWC + W*(A(1,K) + A(2,K)*CX1 + A(3,K)*CY1 + . A(4,K)*CX2 + A(5,K)*CX1*CY1 + . A(6,K)*CY2 + A(7,K)*CX3 + . A(8,K)*CX2*CY1 + A(9,K)*CX1*CY2 + . A(10,K)*CY3) C C Bottom of loop on nodes in cell (I,J). C 2 KP = K K = LNEXT(KP) IF (K .NE. KP) GO TO 1 3 CONTINUE 4 CONTINUE C C SW = 0 iff P is not within the radius R(K) for any node K. C IF (SW .EQ. 0.) GO TO 6 TS2VAL = SWC/SW RETURN C C (PX,PY) = (X(K),Y(K)). C 5 TS2VAL = F(K) RETURN C C All weights are 0 at P. C 6 TS2VAL = 0. RETURN END SUBROUTINE TS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN, . YMIN,DX,DY,SX,SY,RMAX,RW,A, C,CX, . CY,IER) INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN, . DX, DY, SX, SY, RMAX, RW(N), A(10,N), . C, CX, CY C C*********************************************************** C C From TSHEP2D C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 02/19/97 C C This subroutine computes the value and gradient at P = C (PX,PY) of the interpolatory function C defined in Sub- C routine TSHEP2. C is a weighted sum of cosine series C nodal functions. C C On input: C C PX,PY = Cartesian coordinates of the point P at C which C and its partial derivatives are C to be evaluated. C C N = Number of nodes and data values defining C. C N .GE. 10. C C X,Y,F = Arrays of length N containing the nodes and C data values interpolated by C. C C NR = Number of rows and columns in the cell grid. C Refer to Subroutine STORE2. NR .GE. 1. C C LCELL = NR by NR array of nodal indexes associated C with cells. Refer to Subroutine STORE2. C C LNEXT = Array of length N containing next-node C indexes. Refer to Subroutine STORE2. C C XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell C dimensions. DX and DY must be C positive. Refer to Subroutine C STORE2. C C SX,SY = Scale factors for mapping [XMIN,XMAX] and C [YMIN,YMAX] to [0,Pi]: C SX = Pi/(XMAX-XMIN) and SY = Pi/(YMAX-YMIN). C C RMAX = Largest element in RW -- maximum radius R(k). C C RW = Array of length N containing the the radii R(k) C which enter into the weights W(k) defining C. C C A = 10 by N array containing the coefficients for C nodal function C(k) in column k. C C Input parameters are not altered by this subroutine. C The parameters other than PX and PY should be input C unaltered from their values on output from TSHEP2. This C subroutine should not be called if a nonzero error flag C was returned by TSHEP2. C C On output: C C C = Value of C at (PX,PY) unless IER .EQ. 1, in C which case no values are returned. C C CX,CY = First partial derivatives of C at (PX,PY) C unless IER .EQ. 1. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N, NR, DX, DY or RMAX is invalid. C IER = 2 if no errors were encountered but C (PX,PY) is not within the radius R(k) C for any node k (and thus C=CX=CY=0). C C Modules required by TS2GRD: None C C Intrinsic functions called by TS2GRD: COS, INT, SIN, SQRT C C*********************************************************** C INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP DOUBLE PRECISION CK, CKX, CKY, CX1, CX2, CX3, CY1, . CY2, CY3, D, DELX, DELY, R, SW, SWC, . SWCX, SWCY, SWS, SWX, SWY, SX1, SX2, . SX3, SY1, SY2, SY3, T, W, WX, WY, XP, . YP C C Local parameters: C C CK = Value of nodal function C(K) at P C CKX,CKY = Partial derivatives of C(K) with respect to C X and Y, respectively C CX1,CX2,CX3 = Cos(XP), Cos(2*XP), and Cos(3*XP), C respectively C CY1,CY2,CY3 = Cos(YP), Cos(2*YP), and Cos(3*YP), C respectively C D = Distance between P and node K C DELX = PX - X(K) C DELY = PY - Y(K) C I = Cell row index in the range IMIN to IMAX C IMIN,IMAX = Range of cell row indexes of the cells C intersected by a disk of radius RMAX C centered at P C J = Cell column index in the range JMIN to JMAX C JMIN,JMAX = Range of cell column indexes of the cells C intersected by a disk of radius RMAX C centered at P C K = Index of a node in cell (I,J) C KP = Previous value of K in the sequence of nodes C in cell (I,J) C R = Radius of influence for node K C SW = Sum of weights W(K) C SWC = Sum of weighted nodal function values at P C SWCX,SWCY = Partial derivatives of SWC with respect to X C and Y, respectively C SWS = SW**2 C SWX,SWY = Partial derivatives of SW with respect to X C and Y, respectively C SX1,SX2,SX3 = Sin(XP), Sin(2*XP), and Sin(3*XP), C respectively C SY1,SY2,SY3 = Sin(YP), Sin(2*YP), and Sin(3*YP), C respectively C T = Temporary variable C W = Weight W(K) value at P: ((R-D)+/(R*D))**3, C where (R-D)+ = 0 if R < D C WX,WY = Partial derivatives of W with respect to X C and Y, respectively C XP,YP = Coordinates of node NP transformed by the C affine transformation that maps [XMIN,XMAX] C X [YMIN,YMAX] to [0,Pi] X [0,Pi] C XP = SX*(PX-XMIN) YP = SY*(PY-YMIN) IF (N .LT. 10 .OR. NR .LT. 1 .OR. DX .LE. 0. .OR. . DY .LE. 0. .OR. RMAX .LT. 0.) GO TO 6 C C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining C the range of the search for nodes whose radii include C P. The cells which must be searched are those inter- C sected by (or contained in) a circle of radius RMAX C centered at P. C IMIN = INT((PX-XMIN-RMAX)/DX) + 1 IMAX = INT((PX-XMIN+RMAX)/DX) + 1 IF (IMIN .LT. 1) IMIN = 1 IF (IMAX .GT. NR) IMAX = NR JMIN = INT((PY-YMIN-RMAX)/DY) + 1 JMAX = INT((PY-YMIN+RMAX)/DY) + 1 IF (JMIN .LT. 1) JMIN = 1 IF (JMAX .GT. NR) JMAX = NR C C The following is a test for no cells within the circle C of radius RMAX. C IF (IMIN .GT. IMAX .OR. JMIN .GT. JMAX) GO TO 7 C C Compute trig function values at (XP,YP). C CX1 = COS(XP) CY1 = COS(YP) CX2 = COS(2.0*XP) CY2 = COS(2.0*YP) CX3 = COS(3.0*XP) CY3 = COS(3.0*YP) SX1 = SIN(XP) SY1 = SIN(YP) SX2 = SIN(2.0*XP) SY2 = SIN(2.0*YP) SX3 = SIN(3.0*XP) SY3 = SIN(3.0*YP) C C C = SWC/SW = Sum(W(K)*C(K))/Sum(W(K)), where the sum is C from K = 1 to N, C(K) is the nodal function value, C and W(K) = ((R-D)+/(R*D))**3 for radius R(K) and dist- C ance D(K). Thus C C CX = (SWCX*SW - SWC*SWX)/SW**2 and C CY = (SWCY*SW - SWC*SWY)/SW**2 C C where SWCX and SWX are partial derivatives with respect C to X of SWC and SW, respectively. SWCY and SWY are de- C fined similarly. C SW = 0. SWX = 0. SWY = 0. SWC = 0. SWCX = 0. SWCY = 0. C C Outer loop on cells (I,J). C DO 4 J = JMIN,JMAX DO 3 I = IMIN,IMAX K = LCELL(I,J) IF (K .EQ. 0) GO TO 3 C C Inner loop on nodes K. C 1 DELX = PX - X(K) DELY = PY - Y(K) D = SQRT(DELX*DELX + DELY*DELY) R = RW(K) IF (D .GE. R) GO TO 2 C CK = A(1,K) + A(2,K)*CX1 + A(3,K)*CY1 + . A(4,K)*CX2 + A(5,K)*CX1*CY1 + A(6,K)*CY2 + . A(7,K)*CX3 + A(8,K)*CX2*CY1 + . A(9,K)*CX1*CY2 + A(10,K)*CY3 CKX = -SX*(A(2,K)*SX1 + 2.0*A(4,K)*SX2 + . A(5,K)*SX1*CY1 + 3.0*A(7,K)*SX3 + . 2.0*A(8,K)*SX2*CY1 + A(9,K)*SX1*CY2) CKY = -SY*(A(3,K)*SY1 + A(5,K)*CX1*SY1 + . 2.0*A(6,K)*SY2 + A(8,K)*CX2*SY1 + . 2.0*A(9,K)*CX1*SY2 + 3.0*A(10,K)*SY3) C IF (D .EQ. 0.) GO TO 5 T = (1.0/D - 1.0/R) W = T**3 T = -3.0*T*T/(D**3) WX = DELX*T WY = DELY*T C SW = SW + W SWX = SWX + WX SWY = SWY + WY SWC = SWC + W*CK SWCX = SWCX + WX*CK + W*CKX SWCY = SWCY + WY*CK + W*CKY C C Bottom of loop on nodes in cell (I,J). C 2 KP = K K = LNEXT(KP) IF (K .NE. KP) GO TO 1 3 CONTINUE 4 CONTINUE C C SW = 0 iff P is not within the radius R(K) for any node K. C IF (SW .EQ. 0.) GO TO 7 C = SWC/SW SWS = SW*SW CX = (SWCX*SW - SWC*SWX)/SWS CY = (SWCY*SW - SWC*SWY)/SWS IER = 0 RETURN C C (PX,PY) = (X(K),Y(K)). C 5 C = F(K) CX = CKX CY = CKY IER = 0 RETURN C C Invalid input parameter. C 6 IER = 1 RETURN C C No cells contain a point within RMAX of P, or C SW = 0 and thus D .GE. RW(K) for all K. C 7 C = 0. CX = 0. CY = 0. IER = 2 RETURN END SUBROUTINE TS2HES (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN, . YMIN,DX,DY,SX,SY,RMAX,RW,A, C,CX, . CY,CXX,CXY,CYY,IER) INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN, . DX, DY, SX, SY, RMAX, RW(N), A(10,N), . C, CX, CY, CXX, CXY, CYY C C*********************************************************** C C From TSHEP2D C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 02/19/97 C C This subroutine computes the value, gradient, and C Hessian at P = (PX,PY) of the interpolatory function C C defined in Subroutine TSHEP2. C is a weighted sum of C cosine series nodal functions. C C On input: C C PX,PY = Cartesian coordinates of the point P at C which C and its partial derivatives are C to be evaluated. C C N = Number of nodes and data values defining C. C N .GE. 10. C C X,Y,F = Arrays of length N containing the nodes and C data values interpolated by C. C C NR = Number of rows and columns in the cell grid. C Refer to Subroutine STORE2. NR .GE. 1. C C LCELL = NR by NR array of nodal indexes associated C with cells. Refer to Subroutine STORE2. C C LNEXT = Array of length N containing next-node C indexes. Refer to Subroutine STORE2. C C XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell C dimensions. DX and DY must be C positive. Refer to Subroutine C STORE2. C C SX,SY = Scale factors for mapping [XMIN,XMAX] and C [YMIN,YMAX] to [0,Pi]: C SX = Pi/(XMAX-XMIN) and SY = Pi/(YMAX-YMIN). C C RMAX = Largest element in RW -- maximum radius R(k). C C RW = Array of length N containing the the radii R(k) C which enter into the weights W(k) defining C. C C A = 10 by N array containing the coefficients for C nodal function C(k) in column k. C C Input parameters are not altered by this subroutine. C The parameters other than PX and PY should be input C unaltered from their values on output from TSHEP2. This C subroutine should not be called if a nonzero error flag C was returned by TSHEP2. C C On output: C C C = Value of C at (PX,PY) unless IER .EQ. 1, in C which case no values are returned. C C CX,CY = First partial derivatives of C at (PX,PY) C unless IER .EQ. 1. C C CXX,CXY,CYY = Second partial derivatives of C at C (PX,PY) unless IER .EQ. 1. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N, NR, DX, DY or RMAX is invalid. C IER = 2 if no errors were encountered but C (PX,PY) is not within the radius R(k) C for any node k (and thus C = 0). C C Modules required by TS2HES: None C C Intrinsic functions called by TS2HES: COS, INT, SIN, SQRT C C*********************************************************** C INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP DOUBLE PRECISION CK, CKX, CKXX, CKXY, CKY, CKYY, CX1, . CX2, CX3, CY1, CY2, CY3, D, DELX, . DELY, DXSQ, DYSQ, R, SW, SWC, SWCX, . SWCXX, SWCXY, SWCY, SWCYY, SWS, . SWX, SWXX, SWXY, SWY, SWYY, SX1, SX2, . SX3, SY1, SY2, SY3, T1, T2, W, WX, . WXX, WXY, WY, WYY, XP, YP C C Local parameters: C C CK = Value of nodal function C(K) at P C CKX,CKY = Partial derivatives of C(K) with respect C to X and Y, respectively C CKXX,CKXY,CKYY = Second partial derivatives of CK C CX1,CX2,CX3 = Cos(XP), Cos(2*XP), and Cos(3*XP), C respectively C CY1,CY2,CY3 = Cos(YP), Cos(2*YP), and Cos(3*YP), C respectively C D = Distance between P and node K C DELX = PX - X(K) C DELY = PY - Y(K) C DXSQ,DYSQ = DELX**2, DELY**2 C I = Cell row index in the range IMIN to IMAX C IMIN,IMAX = Range of cell row indexes of the cells C intersected by a disk of radius RMAX C centered at P C J = Cell column index in the range JMIN to C JMAX C JMIN,JMAX = Range of cell column indexes of the cells C intersected by a disk of radius RMAX C centered at P C K = Index of a node in cell (I,J) C KP = Previous value of K in the sequence of C nodes in cell (I,J) C R = Radius of influence for node K C SW = Sum of weights W(K) C SWC = Sum of weighted nodal function values C at P C SWCX,SWCY = Partial derivatives of SWC with respect C to X and Y, respectively C SWCXX,SWCXY,SWCYY = Second partial derivatives of SWC C SWS = SW**2 C SWX,SWY = Partial derivatives of SW with respect to C X and Y, respectively C SWXX,SWXY,SWYY = Second partial derivatives of SW C SX1,SX2,SX3 = Sin(XP), Sin(2*XP), and Sin(3*XP), C respectively C SY1,SY2,SY3 = Sin(YP), Sin(2*YP), and Sin(3*YP), C respectively C T1,T2 = Temporary variables C W = Weight W(K) value at P: C ((R-D)+/(R*D))**3, where (R-D)+ = 0 C if R < D C WX,WY = Partial derivatives of W with respect to C X and Y, respectively C WXX,WXY,WYY = Second partial derivatives of W C XP,YP = Coordinates of node NP transformed by the C affine transformation that maps C [XMIN,XMAX] X [YMIN,YMAX] to [0,Pi] X C [0,Pi] C XP = SX*(PX-XMIN) YP = SY*(PY-YMIN) IF (N .LT. 10 .OR. NR .LT. 1 .OR. DX .LE. 0. .OR. . DY .LE. 0. .OR. RMAX .LT. 0.) GO TO 6 C C Set IMIN, IMAX, JMIN, and JMAX to cell indexes defining C the range of the search for nodes whose radii include C P. The cells which must be searched are those inter- C sected by (or contained in) a circle of radius RMAX C centered at P. C IMIN = INT((PX-XMIN-RMAX)/DX) + 1 IMAX = INT((PX-XMIN+RMAX)/DX) + 1 IF (IMIN .LT. 1) IMIN = 1 IF (IMAX .GT. NR) IMAX = NR JMIN = INT((PY-YMIN-RMAX)/DY) + 1 JMAX = INT((PY-YMIN+RMAX)/DY) + 1 IF (JMIN .LT. 1) JMIN = 1 IF (JMAX .GT. NR) JMAX = NR C C The following is a test for no cells within the circle C of radius RMAX. C IF (IMIN .GT. IMAX .OR. JMIN .GT. JMAX) GO TO 7 C C Compute trig function values at (XP,YP). C CX1 = COS(XP) CY1 = COS(YP) CX2 = COS(2.0*XP) CY2 = COS(2.0*YP) CX3 = COS(3.0*XP) CY3 = COS(3.0*YP) SX1 = SIN(XP) SY1 = SIN(YP) SX2 = SIN(2.0*XP) SY2 = SIN(2.0*YP) SX3 = SIN(3.0*XP) SY3 = SIN(3.0*YP) C C C = SWC/SW = Sum(W(K)*C(K))/Sum(W(K)), where the sum is C from K = 1 to N, C(K) is the nodal function value, C and W(K) = ((R-D)+/(R*D))**3 for radius R(K) and dist- C ance D(K). Thus C C CX = (SWCX*SW - SWC*SWX)/SW**2 and C CY = (SWCY*SW - SWC*SWY)/SW**2 C C where SWCX and SWX are partial derivatives with respect C to x of SWC and SW, respectively. SWCY and SWY are de- C fined similarly. The second partials are C C CXX = ( SW*(SWCXX - 2*SWX*CX) - SWC*SWXX )/SW**2 C CXY = ( SW*(SWCXY-SWX*CY-SWY*CX) - SWC*SWXY )/SW**2 C CYY = ( SW*(SWCYY - 2*SWY*CY) - SWC*SWYY )/SW**2 C C where SWCXX and SWXX are second partials with respect C to x, SWCXY and SWXY are mixed partials, and SWCYY and C SWYY are second partials with respect to y. C SW = 0. SWX = 0. SWY = 0. SWXX = 0. SWXY = 0. SWYY = 0. SWC = 0. SWCX = 0. SWCY = 0. SWCXX = 0. SWCXY = 0. SWCYY = 0. C C Outer loop on cells (I,J). C DO 4 J = JMIN,JMAX DO 3 I = IMIN,IMAX K = LCELL(I,J) IF (K .EQ. 0) GO TO 3 C C Inner loop on nodes K. C 1 DELX = PX - X(K) DELY = PY - Y(K) DXSQ = DELX*DELX DYSQ = DELY*DELY D = SQRT(DXSQ + DYSQ) R = RW(K) IF (D .GE. R) GO TO 2 C CK = A(1,K) + A(2,K)*CX1 + A(3,K)*CY1 + . A(4,K)*CX2 + A(5,K)*CX1*CY1 + A(6,K)*CY2 + . A(7,K)*CX3 + A(8,K)*CX2*CY1 + . A(9,K)*CX1*CY2 + A(10,K)*CY3 CKX = -SX*(A(2,K)*SX1 + 2.0*A(4,K)*SX2 + . A(5,K)*SX1*CY1 + 3.0*A(7,K)*SX3 + . 2.0*A(8,K)*SX2*CY1 + A(9,K)*SX1*CY2) CKY = -SY*(A(3,K)*SY1 + A(5,K)*CX1*SY1 + . 2.0*A(6,K)*SY2 + A(8,K)*CX2*SY1 + . 2.0*A(9,K)*CX1*SY2 + 3.0*A(10,K)*SY3) CKXX = -SX*SX*(A(2,K)*CX1 + 4.0*A(4,K)*CX2 + . A(5,K)*CX1*CY1 + 9.0*A(7,K)*CX3 + . 4.0*A(8,K)*CX2*CY1 + A(9,K)*CX1*CY2) CKXY = SX*SY*(A(5,K)*SX1*SY1 + 2.0*A(8,K)*SX2*SY1 + . 2.0*A(9,K)*SX1*SY2) CKYY = -SY*SY*(A(3,K)*CY1 + A(5,K)*CX1*CY1 + . 4.0*A(6,K)*CY2 + A(8,K)*CX2*CY1 + . 4.0*A(9,K)*CX1*CY2 + 9.0*A(10,K)*CY3) C IF (D .EQ. 0.) GO TO 5 T1 = (1.0/D - 1.0/R) W = T1**3 T2 = -3.0*T1*T1/(D**3) WX = DELX*T2 WY = DELY*T2 T1 = 3.0*T1*(2.0+3.0*D*T1)/(D**6) WXX = T1*DXSQ + T2 WXY = T1*DELX*DELY WYY = T1*DYSQ + T2 C SW = SW + W SWX = SWX + WX SWY = SWY + WY SWXX = SWXX + WXX SWXY = SWXY + WXY SWYY = SWYY + WYY SWC = SWC + W*CK SWCX = SWCX + WX*CK + W*CKX SWCY = SWCY + WY*CK + W*CKY SWCXX = SWCXX + W*CKXX + 2.0*WX*CKX + CK*WXX SWCXY = SWCXY + W*CKXY + WX*CKY + WY*CKX + CK*WXY SWCYY = SWCYY + W*CKYY + 2.0*WY*CKY + CK*WYY C C Bottom of loop on nodes in cell (I,J). C 2 KP = K K = LNEXT(KP) IF (K .NE. KP) GO TO 1 3 CONTINUE 4 CONTINUE C C SW = 0 iff P is not within the radius R(K) for any node K. C IF (SW .EQ. 0.) GO TO 7 C = SWC/SW SWS = SW*SW CX = (SWCX*SW - SWC*SWX)/SWS CY = (SWCY*SW - SWC*SWY)/SWS CXX = (SW*(SWCXX-2.0*SWX*CX) - SWC*SWXX)/SWS CXY = (SW*(SWCXY-SWY*CX-SWX*CY) - SWC*SWXY)/SWS CYY = (SW*(SWCYY-2.0*SWY*CY) - SWC*SWYY)/SWS IER = 0 RETURN C C (PX,PY) = (X(K),Y(K)). C 5 C = F(K) CX = CKX CY = CKY CXX = CKXX CXY = CKXY CYY = CKYY IER = 0 RETURN C C Invalid input parameter. C 6 IER = 1 RETURN C C No cells contain a point within RMAX of P, or C SW = 0 and thus D .GE. RW(K) for all K. C 7 C = 0. CX = 0. CY = 0. CXX = 0. CXY = 0. CYY = 0. IER = 2 RETURN END SUBROUTINE GETNP2 (PX,PY,X,Y,NR,LCELL,LNEXT,XMIN,YMIN, . DX,DY, NP,DSQ) INTEGER NR, LCELL(NR,NR), LNEXT(*), NP DOUBLE PRECISION PX, PY, X(*), Y(*), XMIN, YMIN, DX, . DY, DSQ C C*********************************************************** C C From TSHEP2D C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 02/03/97 C C Given a set of N nodes and the data structure defined in C Subroutine STORE2, this subroutine uses the cell method to C find the closest unmarked node NP to a specified point P. C NP is then marked by setting LNEXT(NP) to -LNEXT(NP). (A C node is marked if and only if the corresponding LNEXT ele- C ment is negative. The absolute values of LNEXT elements, C however, must be preserved.) Thus, the closest M nodes to C P may be determined by a sequence of M calls to this rou- C tine. Note that if the nearest neighbor to node K is to C be determined (PX = X(K) and PY = Y(K)), then K should be C marked before the call to this routine. C C The search is begun in the cell containing (or closest C to) P and proceeds outward in rectangular layers until all C cells which contain points within distance R of P have C been searched, where R is the distance from P to the first C unmarked node encountered (infinite if no unmarked nodes C are present). C C This code is essentially unaltered from the subroutine C of the same name in QSHEP2D. C C On input: C C PX,PY = Cartesian coordinates of the point P whose C nearest unmarked neighbor is to be found. C C X,Y = Arrays of length N, for N .GE. 2, containing C the Cartesian coordinates of the nodes. C C NR = Number of rows and columns in the cell grid. C Refer to Subroutine STORE2. NR .GE. 1. C C LCELL = NR by NR array of nodal indexes associated C with cells. Refer to Subroutine STORE2. C C LNEXT = Array of length N containing next-node C indexes (or their negatives). Refer to C Subroutine STORE2. C C XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell C dimensions. DX and DY must be C positive. Refer to Subroutine C STORE2. C C Input parameters other than LNEXT are not altered by C this routine. With the exception of (PX,PY) and the signs C of LNEXT elements, these parameters should be unaltered C from their values on output from Subroutine STORE2. C C On output: C C NP = Index (for X and Y) of the nearest unmarked C node to P, or 0 if all nodes are marked or NR C .LT. 1 or DX .LE. 0 or DY .LE. 0. LNEXT(NP) C .LT. 0 IF NP .NE. 0. C C DSQ = Squared Euclidean distance between P and node C NP, or 0 if NP = 0. C C Modules required by GETNP2: None C C Intrinsic functions called by GETNP2: ABS, INT, SQRT C C*********************************************************** C INTEGER I, I0, I1, I2, IMAX, IMIN, J, J0, J1, J2, . JMAX, JMIN, L, LMIN, LN LOGICAL FIRST DOUBLE PRECISION DELX, DELY, R, RSMIN, RSQ, XP, YP C C Local parameters: C C DELX,DELY = PX-XMIN, PY-YMIN C FIRST = Logical variable with value TRUE iff the C first unmarked node has yet to be C encountered C I,J = Cell indexes in the range [I1,I2] X [J1,J2] C I0,J0 = Indexes of the cell containing or closest C to P C I1,I2,J1,J2 = Range of cell indexes defining the layer C whose intersection with the range C [IMIN,IMAX] X [JMIN,JMAX] is currently C being searched C IMIN,IMAX = Cell row indexes defining the range of the C search C JMIN,JMAX = Cell column indexes defining the range of C the search C L,LN = Indexes of nodes in cell (I,J) C LMIN = Current candidate for NP C R = Distance from P to node LMIN C RSMIN = Squared distance from P to node LMIN C RSQ = Squared distance from P to node L C XP,YP = Local copy of PX,PY -- coordinates of P C XP = PX YP = PY C C Test for invalid input parameters. C IF (NR .LT. 1 .OR. DX .LE. 0. .OR. DY .LE. 0.) . GO TO 9 C C Initialize parameters. C FIRST = .TRUE. IMIN = 1 IMAX = NR JMIN = 1 JMAX = NR DELX = XP - XMIN DELY = YP - YMIN I0 = INT(DELX/DX) + 1 IF (I0 .LT. 1) I0 = 1 IF (I0 .GT. NR) I0 = NR J0 = INT(DELY/DY) + 1 IF (J0 .LT. 1) J0 = 1 IF (J0 .GT. NR) J0 = NR I1 = I0 I2 = I0 J1 = J0 J2 = J0 C C Outer loop on layers, inner loop on layer cells, excluding C those outside the range [IMIN,IMAX] X [JMIN,JMAX]. C 1 DO 6 J = J1,J2 IF (J .GT. JMAX) GO TO 7 IF (J .LT. JMIN) GO TO 6 DO 5 I = I1,I2 IF (I .GT. IMAX) GO TO 6 IF (I .LT. IMIN) GO TO 5 IF (J .NE. J1 .AND. J .NE. J2 .AND. I .NE. I1 . .AND. I .NE. I2) GO TO 5 C C Search cell (I,J) for unmarked nodes L. C L = LCELL(I,J) IF (L .EQ. 0) GO TO 5 C C Loop on nodes in cell (I,J). C 2 LN = LNEXT(L) IF (LN .LT. 0) GO TO 4 C C Node L is not marked. C RSQ = (X(L)-XP)**2 + (Y(L)-YP)**2 IF (.NOT. FIRST) GO TO 3 C C Node L is the first unmarked neighbor of P encountered. C Initialize LMIN to the current candidate for NP, and C RSMIN to the squared distance from P to LMIN. IMIN, C IMAX, JMIN, and JMAX are updated to define the smal- C lest rectangle containing a circle of radius R = C Sqrt(RSMIN) centered at P, and contained in [1,NR] X C [1,NR] (except that, if P is outside the rectangle C defined by the nodes, it is possible that IMIN > NR, C IMAX < 1, JMIN > NR, or JMAX < 1). FIRST is reset to C FALSE. C LMIN = L RSMIN = RSQ R = SQRT(RSMIN) IMIN = INT((DELX-R)/DX) + 1 IF (IMIN .LT. 1) IMIN = 1 IMAX = INT((DELX+R)/DX) + 1 IF (IMAX .GT. NR) IMAX = NR JMIN = INT((DELY-R)/DY) + 1 IF (JMIN .LT. 1) JMIN = 1 JMAX = INT((DELY+R)/DY) + 1 IF (JMAX .GT. NR) JMAX = NR FIRST = .FALSE. GO TO 4 C C Test for node L closer than LMIN to P. C 3 IF (RSQ .GE. RSMIN) GO TO 4 C C Update LMIN and RSMIN. C LMIN = L RSMIN = RSQ C C Test for termination of loop on nodes in cell (I,J). C 4 IF (ABS(LN) .EQ. L) GO TO 5 L = ABS(LN) GO TO 2 5 CONTINUE 6 CONTINUE C C Test for termination of loop on cell layers. C 7 IF (I1 .LE. IMIN .AND. I2 .GE. IMAX .AND. . J1 .LE. JMIN .AND. J2 .GE. JMAX) GO TO 8 I1 = I1 - 1 I2 = I2 + 1 J1 = J1 - 1 J2 = J2 + 1 GO TO 1 C C Unless no unmarked nodes were encountered, LMIN is the C closest unmarked node to P. C 8 IF (FIRST) GO TO 9 NP = LMIN DSQ = RSMIN LNEXT(LMIN) = -LNEXT(LMIN) RETURN C C Error: NR, DX, or DY is invalid or all nodes are marked. C 9 NP = 0 DSQ = 0. RETURN END SUBROUTINE GIVENS ( A,B, C,S) DOUBLE PRECISION A, B, C, S C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 09/01/88 C C This subroutine constructs the Givens plane rotation, C C ( C S) C G = ( ) , where C*C + S*S = 1, C (-S C) C C which zeros the second component of the vector (A,B)**T C (transposed). Subroutine ROTATE may be called to apply C the transformation to a 2 by N matrix. C C This routine is identical to subroutine SROTG from the C LINPACK BLAS (Basic Linear Algebra Subroutines). C C On input: C C A,B = Components of the vector defining the rota- C tion. These are overwritten by values R C and Z (described below) which define C and S. C C On output: C C A = Signed Euclidean norm R of the input vector: C R = +/-SQRT(A*A + B*B) C C B = Value Z such that: C C = SQRT(1-Z*Z) and S=Z if ABS(Z) .LE. 1, and C C = 1/Z and S = SQRT(1-C*C) if ABS(Z) > 1. C C C = +/-(A/R) or 1 if R = 0. C C S = +/-(B/R) or 0 if R = 0. C C Modules required by GIVENS: None C C Intrinsic functions called by GIVENS: ABS, SQRT C C*********************************************************** C DOUBLE PRECISION AA, BB, R, U, V C C Local parameters: C C AA,BB = Local copies of A and B C R = C*A + S*B = +/-SQRT(A*A+B*B) C U,V = Variables used to scale A and B for computing R C AA = A BB = B IF (ABS(AA) .LE. ABS(BB)) GO TO 1 C C ABS(A) > ABS(B). C U = AA + AA V = BB/U R = SQRT(.25 + V*V) * U C = AA/R S = V * (C + C) C C Note that R has the sign of A, C > 0, and S has C SIGN(A)*SIGN(B). C B = S A = R RETURN C C ABS(A) .LE. ABS(B). C 1 IF (BB .EQ. 0.) GO TO 2 U = BB + BB V = AA/U C C Store R in A. C A = SQRT(.25 + V*V) * U S = BB/A C = V * (S + S) C C Note that R has the sign of B, S > 0, and C has C SIGN(A)*SIGN(B). C B = 1. IF (C .NE. 0.) B = 1./C RETURN C C A = B = 0. C 2 C = 1. S = 0. RETURN END SUBROUTINE ROTATE (N,C,S, X,Y ) INTEGER N DOUBLE PRECISION C, S, X(N), Y(N) C C*********************************************************** C C From SRFPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 09/01/88 C C ( C S) C This subroutine applies the Givens rotation ( ) to C (-S C) C (X(1) ... X(N)) C the 2 by N matrix ( ) . C (Y(1) ... Y(N)) C C This routine is identical to subroutine SROT from the C LINPACK BLAS (Basic Linear Algebra Subroutines). C C On input: C C N = Number of columns to be rotated. C C C,S = Elements of the Givens rotation. Refer to C subroutine GIVENS. C C The above parameters are not altered by this routine. C C X,Y = Arrays of length .GE. N containing the compo- C nents of the vectors to be rotated. C C On output: C C X,Y = Arrays containing the rotated vectors (not C altered if N < 1). C C Modules required by ROTATE: None C C*********************************************************** C INTEGER I DOUBLE PRECISION XI, YI C DO 1 I = 1,N XI = X(I) YI = Y(I) X(I) = C*XI + S*YI Y(I) = -S*XI + C*YI 1 CONTINUE RETURN END SUBROUTINE SETUPC (XI,YI,ZI,W, ROW) DOUBLE PRECISION XI, YI, ZI, W, ROW(11) C C*********************************************************** C C From TSHEP2D C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 02/19/97 C C This subroutine sets up the I-th row of an augmented re- C gression matrix for a weighted least squares fit of a C cosine series function f(x,y) to a set of data values z. C C On input: C C XI,YI = Normalized nodal coordinates (in the range C 0 to Pi). C C ZI = Data value at node I. C C W = Weight associated with node I. C C The above parameters are not altered by this routine. C C ROW = Array of length 11. C C On output: C C ROW = Array containing a row of the augmented re- C gression matrix. C C Modules required by SETUPC: None C C Intrinsic function called by SETUPC: COS C C*********************************************************** C DOUBLE PRECISION CX1, CX2, CX3, CY1, CY2, CY3 C C Local parameters: C C CX1 = Cos(XI) C CX2 = Cos(2*XI) C CX3 = Cos(3*XI) C CY1 = Cos(YI) C CY2 = Cos(2*YI) C CY3 = Cos(3*YI) C CX1 = COS(XI) CY1 = COS(YI) CX2 = COS(2.0*XI) CY2 = COS(2.0*YI) CX3 = COS(3.0*XI) CY3 = COS(3.0*YI) C ROW(1) = W ROW(2) = W*CX1 ROW(3) = W*CY1 ROW(4) = W*CX2 ROW(5) = W*CX1*CY1 ROW(6) = W*CY2 ROW(7) = W*CX3 ROW(8) = W*CX2*CY1 ROW(9) = W*CX1*CY2 ROW(10) = W*CY3 ROW(11) = W*ZI RETURN END SUBROUTINE STORE2 (N,X,Y,NR, LCELL,LNEXT,XMIN,YMIN,DX, . DY,IER) INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER DOUBLE PRECISION X(N), Y(N), XMIN, YMIN, DX, DY C C*********************************************************** C C From TSHEP2D C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 03/28/97 C C Given a set of N arbitrarily distributed nodes in the C plane, this subroutine creates a data structure for a C cell-based method of solving closest-point problems. The C smallest rectangle containing the nodes is partitioned C into an NR by NR uniform grid of cells, and nodes are as- C sociated with cells. In particular, the data structure C stores the indexes of the nodes contained in each cell. C For a uniform random distribution of nodes, the nearest C node to an arbitrary point can be determined in constant C expected time. C C This code is essentially unaltered from the subroutine C of the same name in QSHEP2D. C C On input: C C N = Number of nodes. N .GE. 2. C C X,Y = Arrays of length N containing the Cartesian C coordinates of the nodes. C C NR = Number of rows and columns in the grid. The C cell density (average number of nodes per cell) C is D = N/(NR**2). A recommended value, based C on empirical evidence, is D = 3 -- NR = C Sqrt(N/3). NR .GE. 1. C C The above parameters are not altered by this routine. C C LCELL = Array of length .GE. NR**2. C C LNEXT = Array of length .GE. N. C C On output: C C LCELL = NR by NR cell array such that LCELL(I,J) C contains the index (for X and Y) of the C first node (node with smallest index) in C cell (I,J), or LCELL(I,J) = 0 if no nodes C are contained in the cell. The upper right C corner of cell (I,J) has coordinates (XMIN+ C I*DX,YMIN+J*DY). LCELL is not defined if C IER .NE. 0. C C LNEXT = Array of next-node indexes such that C LNEXT(K) contains the index of the next node C in the cell which contains node K, or C LNEXT(K) = K if K is the last node in the C cell for K = 1,...,N. (The nodes contained C in a cell are ordered by their indexes.) C If, for example, cell (I,J) contains nodes C 2, 3, and 5 (and no others), then LCELL(I,J) C = 2, LNEXT(2) = 3, LNEXT(3) = 5, and C LNEXT(5) = 5. LNEXT is not defined if C IER .NE. 0. C C XMIN,YMIN = Cartesian coordinates of the lower left C corner of the rectangle defined by the C nodes (smallest nodal coordinates) un- C less IER = 1. The upper right corner is C (XMAX,YMAX) for XMAX = XMIN + NR*DX and C YMAX = YMIN + NR*DY. C C DX,DY = Dimensions of the cells unless IER = 1. DX C = (XMAX-XMIN)/NR and DY = (YMAX-YMIN)/NR, C where XMIN, XMAX, YMIN, and YMAX are the C extrema of X and Y. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N < 2 or NR < 1. C IER = 2 if DX = 0 or DY = 0. C C Modules required by STORE2: None C C Intrinsic functions called by STORE2: DBLE, INT C C*********************************************************** C INTEGER I, J, K, L, NN, NNR DOUBLE PRECISION DELX, DELY, XMN, XMX, YMN, YMX C C Local parameters: C C DELX,DELY = Components of the cell dimensions -- local C copies of DX,DY C I,J = Cell indexes C K = Nodal index C L = Index of a node in cell (I,J) C NN = Local copy of N C NNR = Local copy of NR C XMN,XMX = Range of nodal X coordinates C YMN,YMX = Range of nodal Y coordinates C NN = N NNR = NR IF (NN .LT. 2 .OR. NNR .LT. 1) GO TO 5 C C Compute the dimensions of the rectangle containing the C nodes. C XMN = X(1) XMX = XMN YMN = Y(1) YMX = YMN DO 1 K = 2,NN IF (X(K) .LT. XMN) XMN = X(K) IF (X(K) .GT. XMX) XMX = X(K) IF (Y(K) .LT. YMN) YMN = Y(K) IF (Y(K) .GT. YMX) YMX = Y(K) 1 CONTINUE XMIN = XMN YMIN = YMN C C Compute cell dimensions and test for zero area. C DELX = (XMX-XMN)/DBLE(NNR) DELY = (YMX-YMN)/DBLE(NNR) DX = DELX DY = DELY IF (DELX .EQ. 0. .OR. DELY .EQ. 0.) GO TO 6 C C Initialize LCELL. C DO 3 J = 1,NNR DO 2 I = 1,NNR LCELL(I,J) = 0 2 CONTINUE 3 CONTINUE C C Loop on nodes, storing indexes in LCELL and LNEXT. C DO 4 K = NN,1,-1 I = INT((X(K)-XMN)/DELX) + 1 IF (I .GT. NNR) I = NNR J = INT((Y(K)-YMN)/DELY) + 1 IF (J .GT. NNR) J = NNR L = LCELL(I,J) LNEXT(K) = L IF (L .EQ. 0) LNEXT(K) = K LCELL(I,J) = K 4 CONTINUE C C No errors encountered. C IER = 0 RETURN C C Invalid input parameter. C 5 IER = 1 RETURN C C DX = 0 or DY = 0. C 6 IER = 2 RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0