setupq.f

Go to the documentation of this file.

      subroutine setupq ( x, dx, y, npoint, v, qty )

c*********************************************************************72
c
cc SETUPQ is to be called in SMOOTH.
c
c  from  * a practical guide to splines *  by c. de boor    
c  to be called in  s m o o t h
c  put  delx = x(.+1) - x(.)  into  v(.,4),
c  put  the three bands of  q-transp*d  into  v(.,1-3), and
c  put the three bands of  (d*q)-transp*(d*q)  at and above the diagonal
c     into  v(.,5-7) .
c     here,  q is  the tridiagonal matrix of order (npoint-2,npoint)
c  with general row  1/delx(i) , -1/delx(i) - 1/delx(i+1) , 1/delx(i+1)
c  and   d  is the diagonal matrix  with general row  dx(i) .
c
      implicit none

      integer npoint,   i,npm1
      double precision dx(npoint),qty(npoint),v(npoint,7),
     &  x(npoint),y(npoint),
     &                                                       diff,prev
      npm1 = npoint - 1
      v(1,4) = x(2) - x(1)
      do 11 i=2,npm1
         v(i,4) = x(i+1) - x(i)
         v(i,1) = dx(i-1)/v(i-1,4)
         v(i,2) = - dx(i)/v(i,4) - dx(i)/v(i-1,4)
   11    v(i,3) = dx(i+1)/v(i,4)
       v(npoint,1) = 0.0d+00
      do 12 i=2,npm1
   12    v(i,5) = v(i,1)**2 + v(i,2)**2 + v(i,3)**2
      if (npm1 .lt. 3)                  go to 14
      do 13 i=3,npm1
   13    v(i-1,6) = v(i-1,2)*v(i,1) + v(i-1,3)*v(i,2)
   14 v(npm1,6) = 0.0d+00
      if (npm1 .lt. 4)                  go to 16
      do 15 i=4,npm1
   15    v(i-2,7) = v(i-2,3)*v(i,1)
   16 v(npm1-1,7) = 0.0d+00
      v(npm1,7) = 0.0d+00
construct  q-transp. * y  in  qty.
      prev = (y(2) - y(1))/v(1,4)
      do 21 i=2,npm1
         diff = (y(i+1)-y(i))/v(i,4)
         qty(i) = diff - prev
   21    prev = diff
                                        return
      end