1C Output from Public domain Ratfor, version 1.0 2 3C PURPOSE 4C Calculation of the cubic B-spline smoothness prior 5C for "usual" interior knot setup. 6C Uses BSPVD and INTRV in the CMLIB 7C sgm[0-3](nb) Symmetric matrix 'SIGMA' 8C whose (i,j)'th element contains the integral of 9C B''(i,.) B''(j,.) , i=1,2 ... nb and j=i,...nb. 10C Only the upper four diagonals are computed. 11 12 subroutine sgram(sg0,sg1,sg2,sg3,tb,nb) 13 14c implicit none 15C indices 16 integer nb 17 DOUBLE precision sg0(nb),sg1(nb),sg2(nb),sg3(nb), tb(nb+4) 18c ------------- 19 integer ileft,mflag, i,ii,jj, lentb 20 DOUBLE precision vnikx(4,3),work(16),yw1(4),yw2(4), wpt 21c 22 integer interv 23 external interv ! in ../../../appl/interv.c 24 25 lentb=nb+4 26C Initialise the sigma vectors 27 do i=1,nb 28 sg0(i)=0.d0 29 sg1(i)=0.d0 30 sg2(i)=0.d0 31 sg3(i)=0.d0 32 end do 33 34 ileft = 1 35 do i=1,nb 36C Calculate a linear approximation to the second derivative of the 37C non-zero B-splines over the interval [tb(i),tb(i+1)]. 38 ileft = interv(tb(1), nb+1,tb(i), 0,0, ileft, mflag) 39 40C Left end second derivatives 41 call bsplvd (tb,lentb,4,tb(i),ileft,work,vnikx,3) 42 43C Put values into yw1 44 do ii=1,4 45 yw1(ii) = vnikx(ii,3) 46 end do 47 48C Right end second derivatives 49 call bsplvd (tb,lentb,4,tb(i+1),ileft,work,vnikx,3) 50 51C Slope*(length of interval) in Linear Approximation to B'' 52 do ii=1,4 53 yw2(ii) = vnikx(ii,3) - yw1(ii) 54 end do 55 56C Calculate Contributions to the sigma vectors 57 wpt = tb(i+1) - tb(i) 58 if(ileft.ge.4) then 59 do ii=1,4 60 jj=ii 61 sg0(ileft-4+ii) = sg0(ileft-4+ii) + 62 & wpt*(yw1(ii)*yw1(jj)+ 63 & (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 64 & + yw2(ii)*yw2(jj)*0.3330d0) 65 jj=ii+1 66 if(jj.le.4) then 67 sg1(ileft+ii-4) = sg1(ileft+ii-4) + 68 & wpt* (yw1(ii)*yw1(jj) + 69 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 70 & +yw2(ii)*yw2(jj)*0.3330d0 ) 71 endif 72 jj=ii+2 73 if(jj.le.4) then 74 sg2(ileft+ii-4) = sg2(ileft+ii-4) + 75 & wpt* (yw1(ii)*yw1(jj) + 76 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 77 & +yw2(ii)*yw2(jj)*0.3330d0 ) 78 endif 79 jj=ii+3 80 if(jj.le.4) then 81 sg3(ileft+ii-4) = sg3(ileft+ii-4) + 82 & wpt* (yw1(ii)*yw1(jj) + 83 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 84 & +yw2(ii)*yw2(jj)*0.3330d0 ) 85 endif 86 end do 87 88 else if(ileft.eq.3) then 89 do ii=1,3 90 jj=ii 91 sg0(ileft-3+ii) = sg0(ileft-3+ii) + 92 & wpt* (yw1(ii)*yw1(jj) + 93 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 94 & +yw2(ii)*yw2(jj)*0.3330d0 ) 95 jj=ii+1 96 if(jj.le.3) then 97 sg1(ileft+ii-3) = sg1(ileft+ii-3) + 98 & wpt* (yw1(ii)*yw1(jj) + 99 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 100 & +yw2(ii)*yw2(jj)*0.3330d0 ) 101 endif 102 jj=ii+2 103 if(jj.le.3) then 104 sg2(ileft+ii-3) = sg2(ileft+ii-3) + 105 & wpt* (yw1(ii)*yw1(jj) + 106 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 107 & +yw2(ii)*yw2(jj)*0.3330d0 ) 108 endif 109 end do 110 111 else if(ileft.eq.2) then 112 do ii=1,2 113 jj=ii 114 sg0(ileft-2+ii) = sg0(ileft-2+ii) + 115 & wpt* (yw1(ii)*yw1(jj) + 116 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 117 & +yw2(ii)*yw2(jj)*0.3330d0 ) 118 jj=ii+1 119 if(jj.le.2) then 120 sg1(ileft+ii-2) = sg1(ileft+ii-2) + 121 & wpt* (yw1(ii)*yw1(jj) + 122 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 123 & +yw2(ii)*yw2(jj)*0.3330d0 ) 124 endif 125 end do 126 127 128 else if(ileft.eq.1) then 129 do ii=1,1 130 jj=ii 131 sg0(ileft-1+ii) = sg0(ileft-1+ii) + 132 & wpt* (yw1(ii)*yw1(jj) + 133 * (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0 134 & +yw2(ii)*yw2(jj)*0.3330d0 ) 135 end do 136 endif 137 138 end do 139 140 return 141 end 142