📄 shepard.cpp
字号:
/* SF = MARQUARDT STABILIZATION FACTOR USED TO DAMP */
/* OUT THE FIRST 3 SOLUTION COMPONENTS(SECOND */
/* PARTIALS OF THE QUADRATIC) WHEN THE SYSTEM */
/* IS ILL - CONDITIONED. AS SF INCREASES, THE */
/* FITTING FUNCTION APPROACHES A LINEAR */
/* SUM = SUM OF SQUARED EUCLIDEAN DISTANCES BETWEEN */
/* NODE K AND THE NODES USED IN THE LEAST */
/* SQUARES FIT(UNLESS ADDITIONAL NODES ARE */
/* ADDED FOR STABILITY) */
/* T = TEMPORARY VARIABLE FOR ACCUMULATING A SCALAR */
/* PRODUCT IN THE BACK SOLVE */
/* XK, YK = COORDINATES OF NODE K -- X(K), Y(K) */
/* XMN, YMN = LOCAL VARIABLES FOR XMIN AND YMIN */
nn = *n;
nnq = *nq;
nnw = *nw;
nnr = *nr;
nqwmax = max(nnq, nnw);
/* Computing MIN */
i__1 = 40, i__2 = nn - 1;
lmax = min(i__1, i__2);
if (5 > nnq || 1 > nnw || nqwmax > lmax || nnr < 1)
{
goto L20;
}
/* CREATE THE CELL DATA STRUCTURE, AND INITIALIZE RSMX. */
store2_(&nn, &x[1], &y[1], &nnr, &lcell[lcell_offset], &lnext[1], &xmn, &
ymn, &ddx, &ddy, &ierr);
if (ierr != 0)
{
goto L22;
}
rsmx = 0.f;
/* OUTER LOOP ON NODE K */
i__1 = nn;
for (k = 1; k <= i__1; ++k)
{
xk = x[k];
yk = y[k];
fk = f[k];
/* MARK NODE K TO EXCLUDE IT FROM THE SEARCH FOR NEAREST */
/* NEIGHBORS. */
lnext[k] = -lnext[k];
/* INITIALIZE FOR LOOP ON NPTS. */
rs = 0.f;
sum = 0.f;
rws = 0.f;
rq = 0.f;
lnp = 0;
/* COMPUTE NPTS, LNP, RWS, NEQ, RQ, AND AVSQ. */
L1:
sum += rs;
if (lnp == lmax)
{
goto L3;
}
++lnp;
rsold = rs;
getnp2_(&xk, &yk, &x[1], &y[1], &nnr, &lcell[lcell_offset], &lnext[1],
&xmn, &ymn, &ddx, &ddy, &np, &rs);
if (rs == 0.f)
{
goto L21;
}
npts[lnp - 1] = np;
if ((rs - rsold) / rs < rtol)
{
goto L1;
}
if (rws == 0.f && lnp > nnw)
{
rws = rs;
}
if (rq != 0.f || lnp <= nnq)
{
goto L2;
}
/* RQ = 0(NOT YET COMPUTED) AND LNP .GT. NQ. RQ = */
/* SQRT(RS) IS SUFFICIENTLY LARGE TO(STRICTLY) INCLUDE */
/* NQ NODES. THE LEAST SQUARES FIT WILL INCLUDE NEQ = */
/* LNP - 1 EQUATIONS FOR 5 .LE. NQ .LE. NEQ .LT. LMAX */
/* .LE. N - 1. */
neq = lnp - 1;
rq = sqrt(rs);
avsq = sum / (double) neq;
/* BOTTOM OF LOOP -- TEST FOR TERMINATION. */
L2:
if (lnp > nqwmax)
{
goto L4;
}
goto L1;
/* ALL LMAX NODES ARE INCLUDED IN NPTS. RWS AND/OR RQ**2 IS */
/*(ARBITRARILY) TAKEN TO BE 10 PERCENT LARGER THAN THE */
/* DISTANCE RS TO THE LAST NODE INCLUDED. */
L3:
if (rws == 0.f)
{
rws = rs * 1.1f;
}
if (rq != 0.f)
{
goto L4;
}
neq = lmax;
rq = sqrt(rs * 1.1f);
avsq = sum / (double) neq;
/* STORE RSQ(K), UPDATE RSMX IF NECESSARY, AND COMPUTE AV. */
L4:
rsq[k] = rws;
if (rws > rsmx)
{
rsmx = rws;
}
av = sqrt(avsq);
/* SET UP THE AUGMENTED REGRESSION MATRIX(TRANSPOSED) AS THE */
/* COLUMNS OF B, AND ZERO OUT THE LOWER TRIANGLE(UPPER */
/* TRIANGLE OF B) WITH GIVENS ROTATIONS -- QR DECOMPOSITION */
/* WITH ORTHOGONAL MATRIX Q NOT STORED. */
i__ = 0;
L5:
++i__;
np = npts[i__ - 1];
irow = min(i__, 6);
setup2_(&xk, &yk, &fk, &x[np], &y[np], &f[np], &av, &avsq, &rq, &b[
irow * 6 - 6]);
if (i__ == 1)
{
goto L5;
}
irm1 = irow - 1;
i__2 = irm1;
for (j = 1; j <= i__2; ++j)
{
jp1 = j + 1;
givens_(&b[j + j * 6 - 7], &b[j + irow * 6 - 7], &c__, &s);
/* L6: */
i__3 = 6 - j;
rotate_(&i__3, &c__, &s, &b[jp1 + j * 6 - 7], &b[jp1 + irow * 6 -
7]);
}
if (i__ < neq)
{
goto L5;
}
/* TEST THE SYSTEM FOR ILL - CONDITIONING. */
/* Computing MIN */
r__1 = dabs(b[0]), r__2 = dabs(b[7]), r__1 = min(r__1, r__2), r__2 =
dabs(b[14]), r__1 = min(r__1, r__2), r__2 = dabs(b[21]), r__1 =
min(r__1, r__2), r__2 = dabs(b[28]);
dmin__ = dmin(r__1, r__2);
if (dmin__ * rq >= dtol)
{
goto L13;
}
if (neq == lmax)
{
goto L10;
}
/* INCREASE RQ AND ADD ANOTHER EQUATION TO THE SYSTEM TO */
/* IMPROVE THE CONDITIONING. THE NUMBER OF NPTS ELEMENTS */
/* IS ALSO INCREASED IF NECESSARY. */
L7:
rsold = rs;
++neq;
if (neq == lmax)
{
goto L9;
}
if (neq == lnp)
{
goto L8;
}
/* NEQ .LT. LNP */
np = npts[neq];
/* Computing 2nd power */
r__1 = x[np] - xk;
/* Computing 2nd power */
r__2 = y[np] - yk;
rs = r__1 * r__1 + r__2 * r__2;
if ((rs - rsold) / rs < rtol)
{
goto L7;
}
rq = sqrt(rs);
goto L5;
/* ADD AN ELEMENT TO NPTS. */
L8:
++lnp;
getnp2_(&xk, &yk, &x[1], &y[1], &nnr, &lcell[lcell_offset], &lnext[1],
&xmn, &ymn, &ddx, &ddy, &np, &rs);
if (np == 0)
{
goto L21;
}
npts[lnp - 1] = np;
if ((rs - rsold) / rs < rtol)
{
goto L7;
}
rq = sqrt(rs);
goto L5;
L9:
rq = sqrt(rs * 1.1f);
goto L5;
/* STABILIZE THE SYSTEM BY DAMPING SECOND PARTIALS -- ADD */
/* MULTIPLES OF THE FIRST THREE UNIT VECTORS TO THE FIRST */
/* THREE EQUATIONS. */
L10:
for (i__ = 1; i__ <= 3; ++i__)
{
b[i__ + 29] = sf;
ip1 = i__ + 1;
for (j = ip1; j <= 6; ++j)
{
/* L11: */
b[j + 29] = 0.f;
}
for (j = i__; j <= 5; ++j)
{
jp1 = j + 1;
givens_(&b[j + j * 6 - 7], &b[j + 29], &c__, &s);
/* L12: */
i__3 = 6 - j;
rotate_(&i__3, &c__, &s, &b[jp1 + j * 6 - 7], &b[jp1 + 29]);
}
}
/* TEST THE STABILIZED SYSTEM FOR ILL - CONDITIONING. */
/* Computing MIN */
r__1 = dabs(b[0]), r__2 = dabs(b[7]), r__1 = min(r__1, r__2), r__2 =
dabs(b[14]), r__1 = min(r__1, r__2), r__2 = dabs(b[21]), r__1 =
min(r__1, r__2), r__2 = dabs(b[28]);
dmin__ = dmin(r__1, r__2);
if (dmin__ * rq < dtol)
{
goto L22;
}
/* SOLVE THE 5 BY 5 TRIANGULAR SYSTEM FOR THE COEFFICIENTS */
L13:
for (ib = 1; ib <= 5; ++ib)
{
i__ = 6 - ib;
t = 0.f;
if (i__ == 5)
{
goto L15;
}
ip1 = i__ + 1;
for (j = ip1; j <= 5; ++j)
{
/* L14: */
t += b[j + i__ * 6 - 7] * a[j + k * 5];
}
L15:
a[i__ + k * 5] = (b[i__ * 6 - 1] - t) / b[i__ + i__ * 6 - 7];
}
/* SCALE THE COEFFICIENTS TO ADJUST FOR THE COLUMN SCALING. */
for (i__ = 1; i__ <= 3; ++i__)
{
/* L16: */
a[i__ + k * 5] /= avsq;
}
a[k * 5 + 4] /= av;
a[k * 5 + 5] /= av;
/* UNMARK K AND THE ELEMENTS OF NPTS. */
lnext[k] = -lnext[k];
i__3 = lnp;
for (i__ = 1; i__ <= i__3; ++i__)
{
np = npts[i__ - 1];
/* L17: */
lnext[np] = -lnext[np];
}
/* L18: */
}
/* NO ERRORS ENCOUNTERED. */
*xmin = xmn;
*ymin = ymn;
*dx = ddx;
*dy = ddy;
*rmax = sqrt(rsmx);
*ier = 0;
return 0;
/* N, NQ, NW, OR NR IS OUT OF RANGE. */
L20:
*ier = 1;
return 0;
/* DUPLICATE NODES WERE ENCOUNTERED BY GETNP2. */
L21:
*ier = 2;
return 0;
/* NO UNIQUE SOLUTION DUE TO COLLINEAR NODES. */
L22:
*xmin = xmn;
*ymin = ymn;
*dx = ddx;
*dy = ddy;
*ier = 3;
return 0;
} /* qshep2_ */
double qs2val_(double *px, double *py, int *n, double *x, double *y, double *f,
int *nr, int *lcell, int *lnext, double *xmin, double *ymin,
double *dx, double *dy, double *rmax, double *rsq, double *a)
{
int lcell_dim1, lcell_offset, i__1, i__2;
double ret_val;
// double sqrt(double);
static int i__, j, k;
static double w, rd, ds;
static int kp;
static double rs, xp, yp, sw, rds, swq;
static int imin, jmin, imax, jmax;
static double delx, dely, dxsq, dysq;
/****************************************************************/
/* */
/* ROBERT RENKA */
/* UNIV. OF NORTH TEXAS */
/*(817) 565 - 2767 */
/* 10/28/87 */
/* */
/* THIS FUNCTION RETURNS THE VALUE Q(PX, PY) WHERE Q IS THE */
/* WEIGHTED SUM OF QUADRATIC NODAL FUNCTIONS DEFINED IN SUB- */
/* ROUTINE QSHEP2. QS2GRD MAY BE CALLED TO COMPUTE A GRADI- */
/* ENT OF Q ALONG WITH THE VALUE, AND/OR TO TEST FOR ERRORS. */
/* */
/* ON INPUT -- */
/* */
/* PX, PY = CARTESIAN COORDINATES OF THE POINT P AT */
/* WHICH Q IS TO BE EVALUATED. */
/* */
/* N = NUMBER OF NODES AND DATA VALUES DEFINING Q. */
/* N .GE. 6. */
/* */
/* X, Y, F = ARRAYS OF LENGTH N CONTAINING THE NODES AND */
/* DATA VALUES INTERPOLATED BY Q. */
/* */
/* NR = NUMBER OF ROWS AND COLUMNS IN THE CELL GRID. */
/* REFER TO STORE2. NR .GE. 1. */
/* */
/* LCELL = NR BY NR ARRAY OF NODAL INDICES ASSOCIATED */
/* WITH CELLS. REFER TO STORE2. */
/* */
/* LNEXT = ARRAY OF LENGTH N CONTAINING NEXT - NODE INDI- */
/* CES. REFER TO STORE2. */
/* */
/* XMIN, YMIN, DX, DY = MINIMUM NODAL COORDINATES AND CELL */
/* DIMENSIONS. DX AND DY MUST BE */
/* POSITIVE. REFER TO STORE2. */
/* */
/* RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ -- */
/* MAXIMUM RADIUS. */
/* */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -