https://github.com/cran/fields
Tip revision: a82f141190b29d59547a12b9a44a7e2d1af8fb60 authored by Doug Nychka on 05 January 2007, 21:27:53 UTC
version 3.3.1
version 3.3.1
Tip revision: a82f141
dmaket.f
subroutine dmaket(m,n,dim,des,lddes,npoly,t,ldt,
* wptr,info,ptab,ldptab)
integer m,n,dim,lddes,npoly,ldt,wptr(dim),info,ptab(ldptab,dim)
double precision des(lddes,dim),t(ldt,*)
c
c Purpose: create t matrix and append s1 to it.
c
c On Entry:
c m order of the derivatives in the penalty
c n number of rows in des
c dim number of columns in des
c des(lddes,dim) variables to be splined
c lddes leading dimension of des as declared in the
c calling program
c ldt leading dimension of t as declared in the
c calling program
c
c npoly dimension of polynomial part of spline
c On Exit:
c t(ldt,npoly+ncov1) [t:s1]
c info error indication
c 0 : successful completion
c 1 : error in creation of t
c Work Arrays:
c wptr(dim) integer work vector
c
c Subprograms Called Directly:
c Blas - dcopy
c Other - mkpoly
c
c $Header: /fs/image/home/thoar/CVS.REPOS/fields/src/dmaket.f,v 1.1.1.1 2002/12/04 22:46:16 thoar Exp $
c
integer i,j,k,tt,nt,bptr,eptr
c
info = 0
c npoly = mkpoly(m,dim)
do 5 j=1,n
t(j,1)=1.0
5 continue
nt = 1
if (npoly .gt. 1) then
do 10 j=1,dim
nt = j + 1
wptr(j) = nt
ptab(nt,j)= ptab(nt,j) +1
call dcopy(n,des(1,j),1,t(1,nt),1)
10 continue
c
c get cross products of x's in null space for m>2
c
c WARNING: do NOT change next do loop unless you fully understand:
c This first gets x1*x1, x1*x2, x1*x3, then
c x2*x2, x2*x3, and finally x3*x3 for dim=3,n=3
c wptr(1) is always at the beginning of the current
c level of cross products, hence the end of the
c previous level which is used for the next.
c wptr(j) is at the start of xj * (previous level)
c
do 50 k=2,m-1
do 40 j=1,dim
bptr = wptr(j)
wptr(j) = nt + 1
eptr = wptr(1) - 1
do 30 tt=bptr,eptr
nt = nt + 1
do 21 jj= 1,dim
ptab(nt,jj)= ptab(tt,jj)
21 continue
ptab( nt,j)= 1+ ptab( nt,j)
do 20 i=1,n
t(i,nt) = des(i,j) * t(i,tt)
20 continue
30 continue
40 continue
50 continue
if (nt .ne. npoly) then
info = 1
return
endif
endif
c
end