sgram.f

   1 C Output from Public domain Ratfor, version 1.0
   2
   3 C PURPOSE
   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.
  11
  12       subroutine sgram(sg0,sg1,sg2,sg3,tb,nb)
  13
  14 c      implicit none
  15 C indices
  16       integer nb
  17       DOUBLE precision sg0(nb),sg1(nb),sg2(nb),sg3(nb), tb(nb+4)
  18 c     -------------
  19       integer ileft,mflag, i,ii,jj, lentb
  20       DOUBLE precision vnikx(4,3),work(16),yw1(4),yw2(4), wpt
  21 c
  22 c      integer interv
  23 c      external interv
  24
  25       lentb=nb+4
  26 C Initialise the sigma vectors
  27       do 1 i=1,nb
  28          sg0(i)=0.d0
  29          sg1(i)=0.d0
  30          sg2(i)=0.d0
  31          sg3(i)=0.d0
  32  1    continue
  33
  34       ileft = 1
  35       do 2 i=1,nb
  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)
  44 C     Put values into yw1
  45          do 4 ii=1,4
  46             yw1(ii) = vnikx(ii,3)
  47  4       continue
  48
  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)
  52
  53 C     Slope*(length of interval) in Linear Approximation to B''
  54          do    6 ii=1,4
  55             yw2(ii) = vnikx(ii,3) - yw1(ii)
  56  6       continue
  57
  58          wpt = tb(i+1) - tb(i)
  59 C     Calculate Contributions to the sigma vectors
  60          if(ileft.ge.4) then
  61             do 10 ii=1,4
  62                jj=ii
  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)
  67                jj=ii+1
  68                if(jj.le.4)then
  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 )
  73                endif
  74                jj=ii+2
  75                if(jj.le.4)then
  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 )
  80                endif
  81                jj=ii+3
  82                if(jj.le.4)then
  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 )
  87                endif
  88  10         continue
  89
  90          else if(ileft.eq.3)then
  91             do    20 ii=1,3
  92                jj=ii
  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 )
  97                jj=ii+1
  98                if(jj.le.3)then
  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 )
 103                endif
 104                jj=ii+2
 105                if(jj.le.3)then
 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 )
 110                endif
 111  20         continue
 112
 113          else if(ileft.eq.2)then
 114             do    28 ii=1,2
 115                jj=ii
 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 )
 120                jj=ii+1
 121                if(jj.le.2)then
 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 )
 126                endif
 127  28         continue
 128
 129          else if(ileft.eq.1)then
 130             do 34 ii=1,1
 131                jj=ii
 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 )
 136  34         continue
 137
 138          endif
 139  2    continue
 140
 141       return
 142       end