1 C Output from Public domain Ratfor, version 1.0
4 C Calculation of the cubic B-spline smoothness prior
5 C for "usual" interior knot setup.
6 C Uses BSPVD and INTRV in the CMLIB
7 C sgm[0-3](nb) Symmetric matrix
8 C whose (i,j)'th element contains the integral of
9 C B''(i,.) B''(j,.) , i=1,2 ... nb and j=i,...nb.
10 C Only the upper four diagonals are computed.
12 subroutine sgram(sg0,sg1,sg2,sg3,tb,nb)
17 DOUBLE precision sg0(nb),sg1(nb),sg2(nb),sg3(nb), tb(nb+4)
19 integer ileft,mflag, i,ii,jj, lentb
20 DOUBLE precision vnikx(4,3),work(16),yw1(4),yw2(4), wpt
26 C Initialise the sigma vectors
36 C Calculate a linear approximation to the
37 C second derivative of the non-zero B-splines
38 C over the interval [tb(i),tb(i+1)].
39 C call intrv(tb(1),(nb+1),tb(i),ilo,ileft,mflag)
40 call interv(tb(1), nb+1,tb(i), ileft, mflag)
41 C Left end second derivatives
42 C call bspvd (tb,4,3,tb(i),ileft,4,vnikx,work)
43 call bsplvd (tb,lentb,4,tb(i),ileft,work,vnikx,3)
49 C Right end second derivatives
50 C call bspvd (tb,4,3,tb(i+1),ileft,4,vnikx,work)
51 call bsplvd (tb,lentb,4,tb(i+1),ileft,work,vnikx,3)
53 C Slope*(length of interval) in Linear Approximation to B''
55 yw2(ii) = vnikx(ii,3) - yw1(ii)
59 C Calculate Contributions to the sigma vectors
63 sg0(ileft-4+ii) = sg0(ileft-4+ii) +
64 & wpt*(yw1(ii)*yw1(jj)+
65 & (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
66 & + yw2(ii)*yw2(jj)*0.3330d0)
69 sg1(ileft+ii-4) = sg1(ileft+ii-4) +
70 & wpt* (yw1(ii)*yw1(jj) +
71 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
72 & +yw2(ii)*yw2(jj)*0.3330d0 )
76 sg2(ileft+ii-4) = sg2(ileft+ii-4) +
77 & wpt* (yw1(ii)*yw1(jj) +
78 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
79 & +yw2(ii)*yw2(jj)*0.3330d0 )
83 sg3(ileft+ii-4) = sg3(ileft+ii-4) +
84 & wpt* (yw1(ii)*yw1(jj) +
85 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
86 & +yw2(ii)*yw2(jj)*0.3330d0 )
90 else if(ileft.eq.3)then
93 sg0(ileft-3+ii) = sg0(ileft-3+ii) +
94 & wpt* (yw1(ii)*yw1(jj) +
95 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
96 & +yw2(ii)*yw2(jj)*0.3330d0 )
99 sg1(ileft+ii-3) = sg1(ileft+ii-3) +
100 & wpt* (yw1(ii)*yw1(jj) +
101 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
102 & +yw2(ii)*yw2(jj)*0.3330d0 )
106 sg2(ileft+ii-3) = sg2(ileft+ii-3) +
107 & wpt* (yw1(ii)*yw1(jj) +
108 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
109 & +yw2(ii)*yw2(jj)*0.3330d0 )
113 else if(ileft.eq.2)then
116 sg0(ileft-2+ii) = sg0(ileft-2+ii) +
117 & wpt* (yw1(ii)*yw1(jj) +
118 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
119 & +yw2(ii)*yw2(jj)*0.3330d0 )
122 sg1(ileft+ii-2) = sg1(ileft+ii-2) +
123 & wpt* (yw1(ii)*yw1(jj) +
124 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
125 & +yw2(ii)*yw2(jj)*0.3330d0 )
129 else if(ileft.eq.1)then
132 sg0(ileft-1+ii) = sg0(ileft-1+ii) +
133 & wpt* (yw1(ii)*yw1(jj) +
134 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
135 & +yw2(ii)*yw2(jj)*0.3330d0 )