C Output from Public domain Ratfor, version 1.0
      subroutine stxwx(x,z,w,k,xknot,n,y,hs0,hs1,hs2,hs3)

c      implicit none
      integer k,n
      DOUBLE precision x(k),z(k),w(k), xknot(n+4),y(n),
     &     hs0(n),hs1(n),hs2(n),hs3(n)
C local
      DOUBLE precision eps,vnikx(4,1),work(16)
      integer lenxk, i,j, ileft,mflag
c
c      integer interv
c      external interv

      lenxk=n+4
C     Initialise the output vectors
      do 1 i=1,n
         y(i)=0d0
         hs0(i)=0d0
         hs1(i)=0d0
         hs2(i)=0d0
         hs3(i)=0d0
 1    continue

C Compute X' W^2 X -> hs0,hs1,hs2,hs3  and X' W^2 Z -> y
C Note that here the weights w(i) == sqrt(wt[i])  where wt[] where original weights
      ileft=1
      eps= .1d-9

      do 100 i=1,k
         call interv(xknot(1), n+1, x(i), ileft, mflag)
C        if(mflag==-1) {write(6,'("Error in hess ",i2)')mflag;stop}
C        if(mflag==-1) {return}
         if(mflag.eq. 1)then
            if(x(i).le.(xknot(ileft)+eps))then
               ileft=ileft-1
            else
               return
            endif
C        else{write(6,'("Error in hess ",i2)')mflag;stop}}
         endif
         call bsplvd (xknot,lenxk,4,x(i),ileft,work,vnikx,1)

         j= ileft-4+1
         y(j) = y(j)+w(i)**2*z(i)*vnikx(1,1)
         hs0(j)=hs0(j)+w(i)**2*vnikx(1,1)**2
         hs1(j)=hs1(j)+w(i)**2*vnikx(1,1)*vnikx(2,1)
         hs2(j)=hs2(j)+w(i)**2*vnikx(1,1)*vnikx(3,1)
         hs3(j)=hs3(j)+w(i)**2*vnikx(1,1)*vnikx(4,1)
         j= ileft-4+2
         y(j) = y(j)+w(i)**2*z(i)*vnikx(2,1)
         hs0(j)=hs0(j)+w(i)**2*vnikx(2,1)**2
         hs1(j)=hs1(j)+w(i)**2*vnikx(2,1)*vnikx(3,1)
         hs2(j)=hs2(j)+w(i)**2*vnikx(2,1)*vnikx(4,1)
         j= ileft-4+3
         y(j) = y(j)+w(i)**2*z(i)*vnikx(3,1)
         hs0(j)=hs0(j)+w(i)**2*vnikx(3,1)**2
         hs1(j)=hs1(j)+w(i)**2*vnikx(3,1)*vnikx(4,1)
         j= ileft-4+4
         y(j) = y(j)+w(i)**2*z(i)*vnikx(4,1)
         hs0(j)=hs0(j)+w(i)**2*vnikx(4,1)**2
 100  continue

      return
      end