https://github.com/cran/fields
Tip revision: 32c60b2ec8167f7d2b26f55147b4e380a8ad77b3 authored by Doug Nychka on 25 September 2011, 00:00:00 UTC
version 6.6.1
version 6.6.1
Tip revision: 32c60b2
dsetup.f
SUBROUTINE dSETUP(X,WGHT,Y,NPOINT,V,QTY,NMAX,itp,ierr)
C PUT DELX=X(.+1)-X(.) INTO V(.,4)
C PUT THE THREE BANDS OF THE MATRIX Q-TRANSP*D INTO
C V(.,1-3)
C PUT THE THREE BANDS OF (D*Q)-TRANSP*(D*Q) AT AND
C ABOVE THE DIAGONAL INTO V(.,5-7)
C HERE Q IS THE TRIDIAGONAL MATRIX OF ORDER (NPOINT
C -2,NPOINT) THAT SATISFIES Q-TRANSP*T=0 AND WGHT
C IS THE DIAGONAL MATRIX WHOSE DIAGONAL ENTRIES
C ARE THE SQUARE ROOTS OF THE WEIGHTS USED IN THE
C PENALIZED LEAST-SQUARES CRITERION
c
implicit double precision (a-h,o-z)
REAL*8 WGHT(NMAX),X(NMAX),y(NMAX)
REAL*8 QTY(NMAX),V(NMAX,7)
REAL*8 DIFF,PREV
INTEGER NPOINT,I,NPM1
c
NPM1=NPOINT -1
V(1,4)=X(2)-X(1)
if(v(1,4).eq.0.d0) then
ierr=5
return
endif
DO 11 I=2,NPM1
V(I,4)=X(I+1) - X(I)
if(v(I,4).eq.0.d0) then
ierr=5
return
endif
if(itp.eq.0) then
V(I,1)=WGHT(I-1)/V(I-1,4)
V(I,2)=-WGHT(I)/V(I,4) - WGHT(I)/V(I-1,4)
V(I,3)=WGHT(I+1)/V(I,4)
else
V(I,1)=1.d0/V(I-1,4)
V(I,2)=-1.d0/V(I,4) - 1.0/V(I-1,4)
V(I,3)=1.d0/V(I,4)
endif
11 continue
c
V(NPOINT,1)=0.d0
DO 12 I=2,NPM1
12 V(I,5)=V(I,1)**2 + V(I,2)**2 + V(I,3)**2
IF(NPM1 .LT. 3)GO TO 14
DO 13 I=3,NPM1
13 V(I-1,6)=V(I-1,2)*V(I,1) + V(I-1,3)*V(I,2)
14 V(NPM1,6)=0.d0
IF(NPM1 .LT. 4)GO TO 16
DO 15 I=4,NPM1
15 V(I-2,7)=V(I-2,3)*V(I,1)
16 V(NPM1-1,7)=0.d0
V(NPM1,7)=0.d0
c
C CONSTRUCT Q-TRANSP. * Y IN QTY
PREV=(Y(2) - Y(1))/V(1,4)
DO 21 I=2,NPM1
DIFF=(Y(I+1)-Y(I))/V(I,4)
QTY(I)=DIFF - PREV
21 PREV=DIFF
c
RETURN
END