stxwx.f

   1 C Output from Public domain Ratfor, version 1.0
   2       subroutine stxwx(x,z,w,k,xknot,n,y,hs0,hs1,hs2,hs3)
   3
   4 c      implicit none
   5       integer k,n
   6       DOUBLE precision x(k),z(k),w(k), xknot(n+4),y(n),
   7      &     hs0(n),hs1(n),hs2(n),hs3(n)
   8 C local
   9       DOUBLE precision eps,vnikx(4,1),work(16)
  10       integer lenxk, i,j, ileft,mflag
  11 c
  12 c      integer interv
  13 c      external interv
  14
  15       lenxk=n+4
  16 C     Initialise the output vectors
  17       do 1 i=1,n
  18          y(i)=0d0
  19          hs0(i)=0d0
  20          hs1(i)=0d0
  21          hs2(i)=0d0
  22          hs3(i)=0d0
  23  1    continue
  24
  25 C Compute X' W^2 X -> hs0,hs1,hs2,hs3  and X' W^2 Z -> y
  26 C Note that here the weights w(i) == sqrt(wt[i])  where wt[] where original weights
  27       ileft=1
  28       eps= .1d-9
  29
  30       do 100 i=1,k
  31          call interv(xknot(1), n+1, x(i), ileft, mflag)
  32 C        if(mflag==-1) {write(6,'("Error in hess ",i2)')mflag;stop}
  33 C        if(mflag==-1) {return}
  34          if(mflag.eq. 1)then
  35             if(x(i).le.(xknot(ileft)+eps))then
  36                ileft=ileft-1
  37             else
  38                return
  39             endif
  40 C        else{write(6,'("Error in hess ",i2)')mflag;stop}}
  41          endif
  42          call bsplvd (xknot,lenxk,4,x(i),ileft,work,vnikx,1)
  43
  44          j= ileft-4+1
  45          y(j) = y(j)+w(i)**2*z(i)*vnikx(1,1)
  46          hs0(j)=hs0(j)+w(i)**2*vnikx(1,1)**2
  47          hs1(j)=hs1(j)+w(i)**2*vnikx(1,1)*vnikx(2,1)
  48          hs2(j)=hs2(j)+w(i)**2*vnikx(1,1)*vnikx(3,1)
  49          hs3(j)=hs3(j)+w(i)**2*vnikx(1,1)*vnikx(4,1)
  50          j= ileft-4+2
  51          y(j) = y(j)+w(i)**2*z(i)*vnikx(2,1)
  52          hs0(j)=hs0(j)+w(i)**2*vnikx(2,1)**2
  53          hs1(j)=hs1(j)+w(i)**2*vnikx(2,1)*vnikx(3,1)
  54          hs2(j)=hs2(j)+w(i)**2*vnikx(2,1)*vnikx(4,1)
  55          j= ileft-4+3
  56          y(j) = y(j)+w(i)**2*z(i)*vnikx(3,1)
  57          hs0(j)=hs0(j)+w(i)**2*vnikx(3,1)**2
  58          hs1(j)=hs1(j)+w(i)**2*vnikx(3,1)*vnikx(4,1)
  59          j= ileft-4+4
  60          y(j) = y(j)+w(i)**2*z(i)*vnikx(4,1)
  61          hs0(j)=hs0(j)+w(i)**2*vnikx(4,1)**2
  62  100  continue
  63
  64       return
  65       end