      subroutine GSPL(t,n,wcoef,left,biatx,giatx)
c
c     calculates all nontrivial cubic g-splines with
c     prescribed discontinuities in second derivatives
c     at the knots as linear combination 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           biatx= vector of 4 local non-vanishing B-splines
c
c*****OUTPUT*****
c
c           giatx= vector of 4 non-vanishing gB-splines
c
c*****OUTPUT*****
c
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 giatx(4),biatx(4),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 gB-splines as linear combination of B-splines
c
      giatx(1)=biatx(1)/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))
      giatx(2)=biatx(1)/gamma(1)*( ds(-2)/(ds(-2)+ds(-1))*
     +                  (t(left)-t(left-2))+t(left)-t(left-1) )+
     +         biatx(2)/gamma(2)*( t(left+1)-t(left)+temp ) +
     +         biatx(3)/gamma(2)*temp
      temp    =ds(-1)/(ds(-1)+ds(0))*(t(left+1)-t(left-1))
      giatx(3)=biatx(2)/gamma(2)*temp+
     +         biatx(3)/gamma(2)*(temp+t(left+1)-t(left))+
     +         biatx(4)/gamma(3)*( t(left+2)-t(left+1)+ds(2)
     +                 /(ds(1)+ds(2))*(t(left+3)-t(left+1)) )
      giatx(4)=biatx(4)/gamma(3)*ds(0)/(ds(0)+ds(1))*
     +                 (t(left+2)-t(left))
      return
      end
