SUBROUTINE CSHEP2 (N,X,Y,F,NC,NW,NR, LCELL,LNEXT,XMIN, . YMIN,DX,DY,RMAX,RW,A,IER) c*********************************************************************72 c cc CSHEP2 defines a C2 function F(X,Y) interpolating scattered data. C C From CSHEP2D C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 02/13/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 CS2VAL, and its first partial derivatives are C computed by Subroutine CS2GRD. C C The interpolation scheme is a modified Cubic Shepard C method: 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)*(x-X(k))**3 + C A(2,k)*(x-X(k))**2*(y-Y(k)) + C A(3,k)*(x-X(k))*(y-Y(k))**2 + C A(4,k)*(y-Y(k))**3 + A(5,k)*(x-X(k))**2 + C A(6,k)*(x-X(k))*(y-Y(k)) + A(7,k)*(y-Y(k))**2 C + A(8,k)*(x-X(k)) + A(9,k)*(y-Y(k)) + F(k) . C C Thus, C(k) is a cubic function which interpolates the data C value at node k. Its coefficients A(,k) are obtained by a C weighted least squares fit to the closest NC data points C with weights similar to W(k). Note that the radius of C influence for the least squares fit is fixed for each k, C 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 to be used in the least C squares fit for coefficients defining the nodal C functions C(k). Values found to be optimal for C test data sets ranged from 11 to 25. A recom- C mended value for general data sets is NC = 17. C For nodes lying on (or close to) a rectangular C grid, the recommended value is NC = 11. In any C case, NC must be in the range 9 to 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 = 30. In general, NW should be C about 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. 9N. 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 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 = 9 by N array containing the coefficients for C cubic 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 CSHEP2: GETNP2, GIVENS, ROTATE, C SETUP2, STORE2 C C Intrinsic functions called by CSHEP2: ABS, DBLE, MAX, C MIN, SQRT C INTEGER LMX PARAMETER (LMX=40) INTEGER I, IERR, IP1, IRM1, IROW, J, JP1, K, LMAX, . LNP, NEQ, NN, NNC, NNR, NNW, NP, NPTS(LMX), . NCWMAX DOUBLE PRECISION B(10,10), C, DDX, DDY, DMIN, DTOL, . FK, RC, RS, RSMX, RSOLD, RTOL, RWS, . S, SF, SFC, SFS, STF, SUM, T, XK, . XMN, YK, YMN C DATA RTOL/1.D-5/, DTOL/.01/ 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 DDX,DDY = Local variables for DX and DY 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*RC .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 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 C LNP = Current length of NPTS C NEQ = Number of equations in the least squares fit C NN,NNC,NNR = Local copies of N, NC, and NR C NNW = Local copy of NW C NP = NPTS element C NPTS = Array containing the indexes of a sequence of C nodes to be used in the least squares fit C or to compute RW. The nodes are ordered C by 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 RC = Radius of influence which enters into the C weights for C(K) (see Subroutine SETUP2) C RS = Squared distance between K and NPTS(LNP) -- C used to compute RC and RW(K) C RSMX = Maximum squared RW element encountered C RSOLD = Squared distance between K and NPTS(LNP-1) -- C used to compute a relative change in RS C between succeeding NPTS elements C RTOL = Tolerance for detecting a sufficiently large C relative change in RS. If the change is C not greater than RTOL, the nodes are C treated as being the same distance from K C RWS = Current squared value of RW(K) C S = Second component of the plane rotation deter- C mined by subroutine GIVENS C SF = Scale factor for the linear terms (columns 8 C and 9) in the least squares fit -- inverse C of the root-mean-square distance between K C and the nodes (other than K) in the least C squares fit C SFS = Scale factor for the quadratic terms (columns C 5, 6, and 7) in the least squares fit -- C SF*SF C SFC = Scale factor for the cubic terms (first 4 C columns) in the least squares fit -- SF**3 C STF = Marquardt stabilization factor used to damp C out the first 4 solution components (third C partials of the cubic) when the system is C ill-conditioned. As STF increases, the C fitting function approaches a quadratic C polynomial. C SUM = Sum of squared Euclidean distances between C node K and the nodes used in the least C squares fit (unless additional nodes are C added for stability) C T = Temporary variable for accumulating a scalar C product in the back solve C XK,YK = Coordinates of node K -- X(K), Y(K) C XMN,YMN = Local variables for XMIN and YMIN C INTEGER N, NC, NW, NR, LCELL(NR,NR), LNEXT(N), IER DOUBLE PRECISION X(N), Y(N), F(N), XMIN, YMIN, DX, . DY, RMAX, RW(N), A(9,N) NN = N NNC = NC NNW = NW NNR = NR NCWMAX = MAX(NNC,NNW) LMAX = MIN(LMX,NN-1) IF (NNC .LT. 9 .OR. NNW .LT. 1 .OR. NCWMAX .GT. . LMAX .OR. NNR .LT. 1) GO TO 21 C C Create the cell data structure, and initialize RSMX. C CALL STORE2 (NN,X,Y,NNR, LCELL,LNEXT,XMN,YMN,DDX,DDY, . IERR) IF (IERR .NE. 0) GO TO 23 RSMX = 0. C C Outer loop on node K: C DO 16 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. C RS = 0. SUM = 0. RWS = 0. RC = 0. LNP = 0 C C Compute NPTS, LNP, RWS, NEQ, RC, and SFS. C 1 SUM = SUM + RS IF (LNP .EQ. LMAX) GO TO 2 LNP = LNP + 1 RSOLD = RS CALL GETNP2 (XK,YK,X,Y,NNR,LCELL,LNEXT,XMN,YMN, . DDX,DDY, NP,RS) IF (RS .EQ. 0.) GO TO 22 NPTS(LNP) = NP IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 1 IF (RWS .EQ. 0. .AND. LNP .GT. NNW) RWS = RS IF (RC .EQ. 0. .AND. LNP .GT. NNC) THEN C C RC = 0 (not yet computed) and LNP > NC. RC = Sqrt(RS) C is sufficiently large to (strictly) include NC nodes. C The least squares fit will include NEQ = LNP - 1 C equations for 9 .LE. NC .LE. NEQ .LT. LMAX .LE. N-1. C NEQ = LNP - 1 RC = SQRT(RS) SFS = DBLE(NEQ)/SUM 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 RS to the last node included. C 2 IF (RWS .EQ. 0.) RWS = 1.1*RS IF (RC .EQ. 0.) THEN NEQ = LMAX RC = SQRT(1.1*RS) SFS = DBLE(NEQ)/SUM ENDIF C C Store RW(K), update RSMX if necessary, and compute SF C and SFC. C 3 RW(K) = SQRT(RWS) IF (RWS .GT. RSMX) RSMX = RWS SF = SQRT(SFS) SFC = SF*SFS 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 with columns (rows of B) defined C as follows: 1-4 are the cubic terms, 5-7 are the quad- C ratic terms, 8 and 9 are the linear terms, and the last C column is the right hand side. C C Set up the equations and zero out the lower triangle with C Givens rotations. C I = 0 4 I = I + 1 NP = NPTS(I) IROW = MIN(I,10) CALL SETUP2 (XK,YK,FK,X(NP),Y(NP),F(NP),SF,SFS, . SFC,RC, B(1,IROW)) IF (I .EQ. 1) GO TO 4 IRM1 = IROW-1 DO 5 J = 1,IRM1 JP1 = J + 1 CALL GIVENS (B(J,J),B(J,IROW),C,S) CALL ROTATE (10-J,C,S,B(JP1,J),B(JP1,IROW)) 5 CONTINUE IF (I .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)) ) IF (DMIN*RC .GE. DTOL) GO TO 11 IF (NEQ .EQ. 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 RSOLD = RS NEQ = NEQ + 1 IF (NEQ .EQ. LMAX) THEN RC = SQRT(1.1*RS) GO TO 4 ENDIF IF (NEQ .LT. LNP) THEN C C NEQ < LNP. C NP = NPTS(NEQ+1) RS = (X(NP)-XK)**2 + (Y(NP)-YK)**2 IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 6 RC = SQRT(RS) GO TO 4 ENDIF C C NEQ = LNP. Add an element to NPTS. C LNP = LNP + 1 CALL GETNP2 (XK,YK,X,Y,NNR,LCELL,LNEXT,XMN,YMN, . DDX,DDY, NP,RS) IF (NP .EQ. 0) GO TO 22 NPTS(LNP) = NP IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 6 RC = SQRT(RS) GO TO 4 C C Stabilize the system by damping third partials -- add C multiples of the first four unit vectors to the first C four equations. C 7 STF = 1.0/RC DO 10 I = 1,4 B(I,10) = STF IP1 = I + 1 DO 8 J = IP1,10 B(J,10) = 0. 8 CONTINUE DO 9 J = I,9 JP1 = J + 1 CALL GIVENS (B(J,J),B(J,10),C,S) CALL ROTATE (10-J,C,S,B(JP1,J),B(JP1,10)) 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)) ) IF (DMIN*RC .LT. DTOL) GO TO 23 C C Solve the 9 by 9 triangular system for the coefficients. C 11 DO 13 I = 9,1,-1 T = 0. IF (I .NE. 9) THEN IP1 = I + 1 DO 12 J = IP1,9 T = T + B(J,I)*A(J,K) 12 CONTINUE ENDIF A(I,K) = (B(10,I)-T)/B(I,I) 13 CONTINUE C C Scale the coefficients to adjust for the column scaling. C DO 14 I = 1,4 A(I,K) = A(I,K)*SFC 14 CONTINUE A(5,K) = A(5,K)*SFS A(6,K) = A(6,K)*SFS A(7,K) = A(7,K)*SFS A(8,K) = A(8,K)*SF A(9,K) = A(9,K)*SF C C Unmark K and the elements of NPTS. C LNEXT(K) = -LNEXT(K) DO 15 I = 1,LNP NP = NPTS(I) LNEXT(NP) = -LNEXT(NP) 15 CONTINUE 16 CONTINUE C C No errors encountered. C XMIN = XMN YMIN = YMN DX = DDX DY = DDY 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 XMIN = XMN YMIN = YMN DX = DDX DY = DDY IER = 3 RETURN END DOUBLE PRECISION FUNCTION CS2VAL (PX,PY,N,X,Y,F,NR, . LCELL,LNEXT,XMIN,YMIN,DX,DY,RMAX,RW,A) c*********************************************************************72 c cc CS2VAL evaluates an interpolant defined by CSHEP2. C C From CSHEP2D 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 This function returns the value C(PX,PY), where C is the C weighted sum of cubic nodal functions defined in Subrou- C tine CSHEP2. CS2GRD may be called to compute a gradient C of C along with the value, and/or to test for errors. C CS2HES 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 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 = 9 by N array containing the coefficients for C cubic 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 CSHEP2. This function C should not be called if a nonzero error flag was returned C by CSHEP2. C C On output: C C CS2VAL = 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 CS2VAL: NONE C C Intrinsic functions called by CS2VAL: INT, SQRT C INTEGER N, NR, LCELL(NR,NR), LNEXT(N) DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN, . DX, DY, RMAX, RW(N), A(9,N) INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP DOUBLE PRECISION D, DELX, DELY, R, SW, SWC, W, XP, YP C C Local parameters: C C D = Distance between P and node K C DELX = XP - X(K) C DELY = YP - 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 W = Weight W(K) value at P: ((R-D)+/(R*D))**3, C where (R-D)+ = 0 if R < D C XP,YP = Local copies of PX and PY -- coordinates of P C cs2val = 0.0 XP = PX YP = PY 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((XP-XMIN-RMAX)/DX) + 1 IMAX = INT((XP-XMIN+RMAX)/DX) + 1 IF (IMIN .LT. 1) IMIN = 1 IF (IMAX .GT. NR) IMAX = NR JMIN = INT((YP-YMIN-RMAX)/DY) + 1 JMAX = INT((YP-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 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 DELX = XP - X(K) DELY = YP - Y(K) D = SQRT(DELX*DELX + DELY*DELY) 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)*DELX+A(2,K)*DELY+ . A(5,K))*DELX + (A(3,K)*DELY+ . A(6,K))*DELY + A(8,K) )*DELX + . ( (A(4,K)*DELY+A(7,K))*DELY + . A(9,K) )*DELY + F(K) ) 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 CS2VAL = SWC/SW RETURN C C (PX,PY) = (X(K),Y(K)). C 5 CS2VAL = F(K) RETURN C C All weights are 0 at P. C 6 CS2VAL = 0. RETURN END SUBROUTINE CS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN, . YMIN,DX,DY,RMAX,RW,A, C,CX,CY,IER) c*********************************************************************72 c cc CS2GRD returns the value and gradients of a function defined by CSHEP2. C C From CSHEP2D 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 This subroutine computes the value and gradient at P = C (PX,PY) of the interpolatory function C defined in Sub- C routine CSHEP2. C is a weighted sum of cubic nodal C 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 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 = 9 by N array containing the coefficients for C cubic 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 CSHEP2. This C subroutine should not be called if a nonzero error flag C was returned by CSHEP2. 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 CS2GRD: None C C Intrinsic functions called by CS2GRD: INT, SQRT C INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN, . DX, DY, RMAX, RW(N), A(9,N), C, CX, . CY INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP DOUBLE PRECISION CK, CKX, CKY, D, DELX, DELY, R, SW, . SWC, SWCX, SWCY, SWS, SWX, SWY, T, W, . WX, WY, XP, YP C C Local parameters: C C CK = Value of cubic nodal function C(K) at P C CKX,CKY = Partial derivatives of C(K) with respect to X C and Y, respectively C D = Distance between P and node K C DELX = XP - X(K) C DELY = YP - 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 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 = Local copies of PX and PY -- coordinates of P C XP = PX YP = PY 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((XP-XMIN-RMAX)/DX) + 1 IMAX = INT((XP-XMIN+RMAX)/DX) + 1 IF (IMIN .LT. 1) IMIN = 1 IF (IMAX .GT. NR) IMAX = NR JMIN = INT((YP-YMIN-RMAX)/DY) + 1 JMAX = INT((YP-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 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 cubic 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 = XP - X(K) DELY = YP - Y(K) D = SQRT(DELX*DELX + DELY*DELY) R = RW(K) IF (D .GE. R) GO TO 2 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 T = A(2,K)*DELX + A(3,K)*DELY + A(6,K) CKY = ( 3.0*A(4,K)*DELY + A(3,K)*DELX + . 2.0*A(7,K) )*DELY + T*DELX + A(9,K) T = T*DELY + A(8,K) CKX = ( 3.0*A(1,K)*DELX + A(2,K)*DELY + . 2.0*A(5,K) )*DELX + T CK = ( (A(1,K)*DELX+A(5,K))*DELX + T )*DELX + . ( (A(4,K)*DELY+A(7,K))*DELY + A(9,K) )*DELY + . F(K) 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 = A(8,K) CY = A(9,K) 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 CS2HES (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN, . YMIN,DX,DY,RMAX,RW,A, C,CX,CY,CXX, . CXY,CYY,IER) c*********************************************************************72 c cc CS2HES returns value, gradients and Hessian of a function defined by CSHEP2. c C From CSHEP2D 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 This subroutine computes the value, gradient, and C Hessian at P = (PX,PY) of the interpolatory function C C defined in Subroutine CSHEP2. C is a weighted sum of C cubic 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 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 = 9 by N array containing the coefficients for C cubic 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 CSHEP2. This C subroutine should not be called if a nonzero error flag C was returned by CSHEP2. 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 CS2HES: None C C Intrinsic functions called by CS2HES: INT, SQRT C INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER DOUBLE PRECISION PX, PY, X(N), Y(N), F(N), XMIN, YMIN, . DX, DY, RMAX, RW(N), A(9,N), C, CX, . CY, CXX, CXY, CYY INTEGER I, IMAX, IMIN, J, JMAX, JMIN, K, KP DOUBLE PRECISION CK, CKX, CKXX, CKXY, CKY, CKYY, D, . DELX, DELY, DXSQ, DYSQ, R, SW, SWC, . SWCX, SWCXX, SWCXY, SWCY, SWCYY, SWS, . SWX, SWXX, SWXY, SWY, SWYY, T1, T2, . T3, T4, W, WX, WXX, WXY, WY, WYY, XP, . YP C C Local parameters: C C CK = Value of cubic nodal function C(K) at P C CKX,CKY = Partial derivatives of C(K) with respect to X C and Y, respectively C CKXX,CKXY,CKYY = Second partial derivatives of CK C D = Distance between P and node K C DELX = XP - X(K) C DELY = YP - 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 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 SWCXX,SWCXY,SWCYY = Second partial derivatives of SWC C SWS = SW**2 C SWX,SWY = Partial derivatives of SW with respect to X C and Y, respectively C SWXX,SWXY,SWYY = Second partial derivatives of SW C T1,T2,T3,T4 = Temporary variables 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 WXX,WXY,WYY = Second partial derivatives of W C XP,YP = Local copies of PX and PY -- coordinates of P C XP = PX YP = PY 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((XP-XMIN-RMAX)/DX) + 1 IMAX = INT((XP-XMIN+RMAX)/DX) + 1 IF (IMIN .LT. 1) IMIN = 1 IF (IMAX .GT. NR) IMAX = NR JMIN = INT((YP-YMIN-RMAX)/DY) + 1 JMAX = INT((YP-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 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 cubic 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 = XP - X(K) DELY = YP - Y(K) DXSQ = DELX*DELX DYSQ = DELY*DELY D = SQRT(DXSQ + DYSQ) R = RW(K) IF (D .GE. R) GO TO 2 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 T1 = A(1,K)*DELX + A(2,K)*DELY + A(5,K) T2 = T1 + T1 + A(1,K)*DELX T3 = A(4,K)*DELY + A(3,K)*DELX + A(7,K) T4 = T3 + T3 + A(4,K)*DELY CK = (T1*DELX + A(6,K)*DELY + A(8,K))*DELX + . (T3*DELY + A(9,K))*DELY + F(K) CKX = T2*DELX + (A(3,K)*DELY+A(6,K))*DELY + A(8,K) CKY = T4*DELY + (A(2,K)*DELX+A(6,K))*DELX + A(9,K) CKXX = T2 + 3.0*A(1,K)*DELX CKXY = 2.0*(A(2,K)*DELX + A(3,K)*DELY) + A(6,K) CKYY = T4 + 3.0*A(4,K)*DELY 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 = A(8,K) CY = A(9,K) CXX = 2.0*A(5,K) CXY = A(6,K) CYY = 2.0*A(7,K) 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) c*********************************************************************72 c cc GETNP2 gets the closest unmarked node to a specified point. C C From CSHEP2D 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 INTEGER NR, LCELL(NR,NR), LNEXT(*), NP DOUBLE PRECISION PX, PY, X(*), Y(*), XMIN, YMIN, DX, . DY, DSQ 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) c*********************************************************************72 c cc GIVENS constructs a Givens plane rotation. 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 DOUBLE PRECISION A, B, C, S 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 ) c*********************************************************************72 c cc ROTATE applies a Givens rotation. 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 INTEGER N DOUBLE PRECISION C, S, X(N), Y(N) 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 SETUP2 (XK,YK,ZK,XI,YI,ZI,S1,S2,S3,R, ROW) c*********************************************************************72 c cc SETUP2 sets a row of a regression matrix to fit a cubic function to data. C C From CSHEP2D 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 This subroutine sets up the I-th row of an augmented re- C gression matrix for a weighted least squares fit of a C cubic function f(x,y) to a set of data values z, where C f(XK,YK) = ZK. The first four columns (cubic terms) are C scaled by S3, the next three columns (quadratic terms) C are scaled by S2, and the eighth and ninth columns (lin- C ear terms) are scaled by S1. C C On input: C C XK,YK = Coordinates of node K. C C ZK = Data value at node K to be interpolated by f. C C XI,YI,ZI = Coordinates and data value at node I. C C S1,S2,S3 = Scale factors. C C R = Radius of influence about node K defining the C weight. C C The above parameters are not altered by this routine. C C ROW = Array of length 10. C C On output: C C ROW = Array containing a row of the augmented re- C gression matrix. C C Modules required by SETUP2: None C C Intrinsic function called by SETUP2: SQRT C DOUBLE PRECISION XK, YK, ZK, XI, YI, ZI, S1, S2, S3, . R, ROW(10) INTEGER I DOUBLE PRECISION D, DX, DXSQ, DY, DYSQ, W, W1, W2, W3 C C Local parameters: C C D = Distance between nodes K and I C DX = XI - XK C DXSQ = DX*DX C DY = YI - YK C DYSQ = DY*DY C I = DO-loop index C W = Weight associated with the row: (R-D)/(R*D) C (0 if D = 0 or D > R) C W1 = S1*W C W2 = S2*W C W3 = W3*W C DX = XI - XK DY = YI - YK DXSQ = DX*DX DYSQ = DY*DY D = SQRT(DXSQ + DYSQ) IF (D .LE. 0. .OR. D .GE. R) GO TO 1 W = (R-D)/R/D W1 = S1*W W2 = S2*W W3 = S3*W ROW(1) = DXSQ*DX*W3 ROW(2) = DXSQ*DY*W3 ROW(3) = DX*DYSQ*W3 ROW(4) = DYSQ*DY*W3 ROW(5) = DXSQ*W2 ROW(6) = DX*DY*W2 ROW(7) = DYSQ*W2 ROW(8) = DX*W1 ROW(9) = DY*W1 ROW(10) = (ZI - ZK)*W RETURN C C Nodes K and I coincide or node I is outside of the radius C of influence. Set ROW to the zero vector. C 1 DO 2 I = 1,10 ROW(I) = 0. 2 CONTINUE RETURN END SUBROUTINE STORE2 (N,X,Y,NR, LCELL,LNEXT,XMIN,YMIN,DX, . DY,IER) c*********************************************************************72 c cc STORE2 creates a closest-point data structure given points in the plane. c C From CSHEP2D 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 INTEGER N, NR, LCELL(NR,NR), LNEXT(N), IER DOUBLE PRECISION X(N), Y(N), XMIN, YMIN, DX, DY 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