      subroutine GSPLVD(t,n,wcoef,left,index,dbiatx,dgiatx)
c
c     calculates all nontrivial cubic g-splines and their
c     derivatives for prescribed discontinuities in second derivatives
c     as linear combination of derivatives of local cubic B-splines
c
c*****INPUT*****
c           t      = extended knot sequence
c           n      = number of knots
c           wcoef  = prescribed jumps of secnd deriv.
c           left   = integer such that t(left).le.x.lt.t(left+1)
c           index  = integer;
c                           0 = generalized second deriv's
c                   otherwise = classical second deriv's
c           dbiatx = array (4,3) containing deriv's of local B-splines;
c                                (j-1)
c                    dbiatx(i,j)=D    B (x)
c                                      left+i-4
c*****OUTPUT*****
c           dgiatx = array (4,3) containing deriv's of G-splines;
c                                (j-1)
c                    dgiatx(i,j)=D    G (x)
c                                     left+i-4
c*****METHOD****
c     All gB-splines which do not vanish on (t(left),t(left+1))
c     are obtained as linear combinations of hermite cubic
c     B-splines. The coeff's in linear combination are discrete
c     cubic splines, which are calculated explicitly according
c     to the Theorem on discrete splines in thesis. The knot
c     sequence T and coef's WCOEF are expected in the form
c
c     t(1)=t(2)=t(3)=t(4)<t(5)<...<t(n-3)=t(n-2)=t(n-1)=t(n)
c        w(1) w(2) w(3) w(4)   ...     w(n-3) w(n-2) w(n-1).
c
c     The value of first (last) three coef's in WCOEF is
c     irrelevant and thus may be set to 1. The jumps in
c     second deriv's are determined so as to make the
c     the L_2 derivative of spline s continuous:
c
c        WCOEF(i+1)DDs(T_i+)-WCOEF(i)DDs(T_i-)=0 .
c
c     W A R N I N G...The restriction WCOEF(i)>0 imposed by
c     the theory of Chebyshev systems is nowhere checked for.
c     Unpredicted results (though not meanigless!) as well
c     as floating point zero divide may result otherwise.
c
      implicit double precision (a-h,o-z)
      dimension t(n),wcoef(1)
c     dimension wcoef(n-7)
      dimension dgiatx(4,3),dbiatx(4,3),ds(-2:2),gamma(3)
c
c     calcluating measures of subintervals with repect to WCOEF
c
      do, i=-2,2
          ds(i)=(t(left+i+1)-t(left+i))/wcoef(left+i)
      enddo
c
c     If one wishes to consider the gB-splines on a biinfinite
c     mesh, then all gB-splines have support of four intervals,
c     and the following two IF..ENDIF constructions should be
c     left out; the WCOEF sequence should then be modified
c     accordingly, that is, the supposition:
c           w(1)=w(2)=w(3) and w(n-2)=w(n-1)
c     is no longer required. Else, the IF...ENDIF enforces
c     deBoor's maxime: anything times zero is zero.
c
      if (left.eq.4) then
            ds(-2)=1.d0
            ds(-1)=0.d0
      endif
      if (left.eq.n-4) then
            ds(1)=0.d0
            ds(2)=1.d0
      endif
c
c     normalisation factors
c
      do, i=1,3
          gamma(i)=ds(i-3)/(ds(i-3)+ds(i-2))*(t(left+i-1)-t(left+i-3))+
     +             t(left+i-1)-t(left+i-2)+
     +             ds(i-1)/(ds(i-2)+ds(i-1))*(t(left+i)-t(left+i-2))
      enddo
c
c     calculating derivatives of G-splines as linear combination
c     of derivatives of local B-splines
c
      do, i=1,3
      dgiatx(1,i)=dbiatx(1,i)/gamma(1)*ds(0)/(ds(-1)+ds(0))*
     +                           (t(left+1)-t(left-1))
      temp       =ds(1)/(ds(0)+ds(1))*(t(left+2)-t(left))
      dgiatx(2,i)=dbiatx(1,i)/gamma(1)*( ds(-2)/(ds(-2)+ds(-1))*
     +                  (t(left)-t(left-2))+t(left)-t(left-1) )+
     +         dbiatx(2,i)/gamma(2)*( t(left+1)-t(left)+temp ) +
     +         dbiatx(3,i)/gamma(2)*temp
      temp       =ds(-1)/(ds(-1)+ds(0))*(t(left+1)-t(left-1))
      dgiatx(3,i)=dbiatx(2,i)/gamma(2)*temp+
     +         dbiatx(3,i)/gamma(2)*(temp+t(left+1)-t(left))+
     +         dbiatx(4,i)/gamma(3)*( t(left+2)-t(left+1)+ds(2)
     +                 /(ds(1)+ds(2))*(t(left+3)-t(left+1)) )
      dgiatx(4,i)=dbiatx(4,i)/gamma(3)*ds(0)/(ds(0)+ds(1))*
     +                 (t(left+2)-t(left))
      enddo
      if (index.eq.0) then
            do, j=1,4
                  dgiatx(j,3)=dgiatx(j,3)*wcoef(left)
            enddo
      endif
      return
      end
