tautsp Subroutine

subroutine tautsp ( tau, gtau, ntau, gamma, s, break, coef, l, k, iflag )

  !*****************************************************************************80
  !
  !! TAUTSP constructs a cubic spline interpolant to given data.
  !
  !  Discussion:
  !
  !    If 0 < GAMMA, additional knots are introduced where needed to
  !    make the interpolant more flexible locally.  This avoids extraneous
  !    inflection points typical of cubic spline interpolation at knots to
  !    rapidly changing data.
  !
  !  Method:  
  !
  !    On the I-th interval, (TAU(I), TAU(I+1)), the interpolant is of the
  !    form:
  !
  !    (*)  F(U(X)) = A + B * U + C * H(U,Z) + D * H(1-U,1-Z),
  !
  !    with  
  !
  !      U = U(X) = ( X - TAU(I) ) / DTAU(I). 
  !
  !    Here,
  !
  !      Z(I) = ADDG(I+1) / ( ADDG(I) + ADDG(I+1) )
  !
  !    but if the denominator vanishes, we set Z(I) = 0.5
  !
  !    Also, we have
  !
  !      ADDG(J) = abs ( DDG(J) ), 
  !      DDG(J) = DG(J+1) - DG(J),
  !      DG(J) = DIVDIF(J) = ( GTAU(J+1) - GTAU(J) ) / DTAU(J)
  !
  !    and
  !
  !      H(U,Z) = ALPHA * U**3 
  !             + ( 1 - ALPHA ) * ( max ( ( ( U - ZETA ) / ( 1 - ZETA ) ), 0 )**3
  !
  !    with
  !
  !      ALPHA(Z) = ( 1 - GAMMA / 3 ) / ZETA
  !      ZETA(Z) = 1 - GAMMA * min ( ( 1 - Z ), 1/3 )
  !
  !    Thus, for 1/3 <= Z <= 2/3, F is just a cubic polynomial on
  !    the interval I.  Otherwise, it has one additional knot, at
  !
  !      TAU(I) + ZETA * DTAU(I).
  !
  !    As Z approaches 1, H(.,Z) has an increasingly sharp bend near 1,
  !    thus allowing F to turn rapidly near the additional knot.
  !
  !    In terms of F(J) = GTAU(J) and FSECND(J) = second derivative of F 
  !    at TAU(J), the coefficients for (*) are given as:
  !
  !      A = F(I) - D
  !      B = ( F(I+1) - F(I) ) - ( C - D )
  !      C = FSECND(I+1) * DTAU(I)**2 / HSECND(1,Z)
  !      D = FSECND(I) * DTAU(I)**2 / HSECND(1,1-Z)
  !
  !    Hence these can be computed once FSECND(1:NTAU) is fixed.
  !
  !    F is automatically continuous and has a continuous second derivative
  !    except when Z=0 or 1 for some I.  We determine FSECND from
  !    the requirement that the first derivative of F be continuous.
  !
  !    In addition, we require that the third derivative be continuous
  !    across TAU(2) and across TAU(NTAU-1).  This leads to a strictly
  !    diagonally dominant tridiagonal linear system for the FSECND(I)
  !    which we solve by Gauss elimination without pivoting.
  !
  !  Modified:
  !
  !    14 February 2007
  !
  !  Author:
  !
  !    Carl DeBoor
  !
  !  Reference:
  !
  !    Carl DeBoor,
  !    A Practical Guide to Splines,
  !    Springer, 2001,
  !    ISBN: 0387953663,
  !    LC: QA1.A647.v27.
  !
  !  Parameters:
  !
  !    Input, real ( kind = 8 ) TAU(NTAU), the sequence of data points.  
  !    TAU must be strictly increasing.
  !
  !    Input, real ( kind = 8 ) GTAU(NTAU), the corresponding sequence of
  !    function values.
  !
  !    Input, integer ( kind = 4 ) NTAU, the number of data points.  
  !    NTAU must be at least 4.
  !
  !    Input, real ( kind = 8 ) GAMMA, indicates whether additional flexibility
  !    is desired.
  !    GAMMA = 0.0, no additional knots;
  !    GAMMA in (0.0,3.0), under certain conditions on the given data at
  !    points I-1, I, I+1, and I+2, a knot is added in the I-th interval, 
  !    for I = 2,...,NTAU-2.  See description of method.  The interpolant 
  !    gets rounded with increasing gamma.  A value of 2.5 for GAMMA is typical.
  !    GAMMA in (3.0,6.0), same, except that knots might also be added in
  !    intervals in which an inflection point would be permitted.  A value 
  !    of 5.5 for GAMMA is typical.
  !
  !    Output, real ( kind = 8 ) BREAK(L), real ( kind = 8 ) COEF(K,L), 
  !    integer ( kind = 4 ) L, integer K, give the piecewise polynomial 
  !    representation of the interpolant.  Specifically, 
  !    for BREAK(i) <= X <= BREAK(I+1), the interpolant has the form:
  !      F(X) = COEF(1,I) +  DX    * ( 
  !             COEF(2,I) + (DX/2) * (
  !             COEF(3,I) + (DX/3) *
  !             COEF(4,I) ) )
  !    with  DX = X - BREAK(I) for I = 1,..., L.
  !
  !    Output, integer ( kind = 4 ) IFLAG, error flag.
  !    1, no error.
  !    2, input was incorrect.
  !
  !    Output, real ( kind = 8 ) S(NTAU,6).  The individual columns of this
  !    array contain the following quantities mentioned in the write up
  !    and below.
  !    S(.,1) = DTAU = TAU(.+1)-TAU;
  !    S(.,2) = DIAG = diagonal in linear system;
  !    S(.,3) = U = upper diagonal in linear system;
  !    S(.,4) = R = right hand side for linear system (initially)
  !           = FSECND = solution of linear system, namely the second
  !             derivatives of interpolant at TAU;
  !    S(.,5) = Z = indicator of additional knots;
  !    S(.,6) = 1/HSECND(1,X) with X = Z or 1-Z.
  !
  implicit none

  integer ( kind = 4 ) ntau

  real ( kind = 8 ) alph
  real ( kind = 8 ) alpha
  real ( kind = 8 ) break(*)
  real ( kind = 8 ) c
  real ( kind = 8 ) coef(4,*)
  real ( kind = 8 ) d
  real ( kind = 8 ) del
  real ( kind = 8 ) denom
  real ( kind = 8 ) divdif
  real ( kind = 8 ) entry
  real ( kind = 8 ) entry3
  real ( kind = 8 ) factor
  real ( kind = 8 ) factr2
  real ( kind = 8 ) gam
  real ( kind = 8 ) gamma
  real ( kind = 8 ) gtau(ntau)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iflag
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) method
  real ( kind = 8 ) onemg3
  real ( kind = 8 ) onemzt
  real ( kind = 8 ) ratio
  real ( kind = 8 ) s(ntau,6)
  real ( kind = 8 ) sixth
  real ( kind = 8 ) tau(ntau)
  real ( kind = 8 ) temp
  real ( kind = 8 ) x
  real ( kind = 8 ) z
  real ( kind = 8 ) zeta
  real ( kind = 8 ) zt2

  alph(x) = min ( 1.0D+00, onemg3 / x )
  !
  !  There must be at least 4 interpolation points.
  !
  if ( ntau < 4 ) then
     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) 'TAUTSP - Fatal error!'
     write ( *, '(a)' ) '  Input NTAU must be at least 4.'
     write ( *, '(a,i8)' ) '  NTAU = ', ntau
     iflag = 2
     stop
  end if
  !
  !  Construct delta TAU and first and second (divided) differences of data.
  !
  do i = 1, ntau - 1

     s(i,1) = tau(i+1) - tau(i)

     if ( s(i,1) <= 0.0D+00 ) then
        write ( *, '(a,i3,a,2e15.6,a)' ) &
             '  Point ', i, ' and the next ', tau(i), tau(i+1), ' are disordered.'
        iflag = 2
        return
     end if

     s(i+1,4) = ( gtau(i+1) - gtau(i) ) / s(i,1)

  end do

  do i = 2, ntau - 1
     s(i,4) = s(i+1,4) - s(i,4)
  end do
  !
  !  Construct system of equations for second derivatives at TAU.
  !
  !  At each interior data point, there is one continuity equation.
  !  At the first and last interior data point there is an additional 
  !  equation for a total of NTAU equations in NTAU unknowns.
  !
  i = 2
  s(2,2) = s(1,1) / 3.0D+00
  sixth = 1.0D+00 / 6.0D+00
  method = 2
  gam = gamma

  if ( gam <= 0.0D+00 ) then
     method = 1
  end if

  if ( 3.0D+00 < gam ) then
     method = 3
     gam = gam - 3.0D+00
  end if

  onemg3 = 1.0D+00 - gam / 3.0D+00
  !
  !  Loop over I.
  !
70 continue
  !
  !  Construct Z(I) and ZETA(I).
  !
  z = 0.5D+00

  if ( method == 1 ) then
     go to 100
  end if

  if ( method == 3 ) then
     go to 90
  end if

  if ( s(i,4) * s(i+1,4) < 0.0D+00 ) then
     go to 100
  end if

90 continue

  temp = abs ( s(i+1,4) )
  denom = abs ( s(i,4) ) + temp

  if ( denom /= 0.0D+00 ) then
     z = temp / denom
     if ( abs ( z - 0.5D+00 ) <= sixth ) then
        z = 0.5D+00
     end if
  end if

100 continue

  s(i,5) = z
  !
  !  Set up part of the I-th equation which depends on the I-th interval.
  !
  if ( z < 0.5D+00 ) then

     zeta = gam * z
     onemzt = 1.0D+00 - zeta
     zt2 = zeta**2
     alpha = alph(onemzt)
     factor = zeta / ( alpha * ( zt2 - 1.0D+00 ) + 1.0D+00 )
     s(i,6) = zeta * factor / 6.0D+00
     s(i,2) = s(i,2) + s(i,1) &
          * ( ( 1.0D+00 - alpha * onemzt ) * factor / 2.0D+00 - s(i,6) )
     !
     !  If Z = 0 and the previous Z = 1, then D(I) = 0.  
     !  Since then also U(I-1) = L(I+1) = 0, its value does not matter.  
     !  Reset D(I) = 1 to insure nonzero pivot in elimination.
     !
     if ( s(i,2) <= 0.0D+00 ) then
        s(i,2) = 1.0D+00
     end if

     s(i,3) = s(i,1) / 6.0D+00

  else if ( z == 0.5D+00 ) then

     s(i,2) = s(i,2) + s(i,1) / 3.0D+00
     s(i,3) = s(i,1) / 6.0D+00

  else if ( 0.5D+00 < z ) then

     onemzt = gam * ( 1.0D+00 - z )
     zeta = 1.0D+00 - onemzt
     alpha = alph(zeta)
     factor = onemzt / ( 1.0D+00 - alpha * zeta * ( 1.0D+00 + onemzt ) )
     s(i,6) = onemzt * factor / 6.0D+00
     s(i,2) = s(i,2) + s(i,1) / 3.0D+00
     s(i,3) = s(i,6) * s(i,1)

  end if

  if ( 2 < i ) then
     go to 190
  end if

  s(1,5) = 0.5D+00
  !
  !  The first two equations enforce continuity of the first and of
  !  the third derivative across TAU(2).
  !
  s(1,2) = s(1,1) / 6.0D+00
  s(1,3) = s(2,2)
  entry3 = s(2,3)

  if ( z < 0.5D+00 ) then

     factr2 = zeta * ( alpha * ( zt2 - 1.0D+00 ) + 1.0D+00 ) &
          / ( alpha * ( zeta * zt2 - 1.0D+00 ) + 1.0D+00 )

     ratio = factr2 * s(2,1) / s(1,2)
     s(2,2) = factr2 * s(2,1) + s(1,1)
     s(2,3) = - factr2 * s(1,1)

  else if ( z == 0.5D+00 ) then

     ratio = s(2,1) / s(1,2)
     s(2,2) = s(2,1) + s(1,1)
     s(2,3) = - s(1,1)

  else if ( 0.5D+00 < z ) then

     ratio = s(2,1) / s(1,2)
     s(2,2) = s(2,1) + s(1,1)
     s(2,3) = - s(1,1) * 6.0D+00 * alpha * s(2,6)

  end if
  !
  !  At this point, the first two equations read:
  !
  !              DIAG(1)*X1+U(1)*X2 + ENTRY3*X3 = R(2)
  !       -RATIO*DIAG(1)*X1+DIAG(2)*X2 + U(2)*X3 = 0.0
  !
  !  Eliminate first unknown from second equation.
  !
  s(2,2) = ratio * s(1,3) + s(2,2)
  s(2,3) = ratio * entry3 + s(2,3)
  s(1,4) = s(2,4)
  s(2,4) = ratio * s(1,4)

  go to 200

190 continue
  !
  !  The I-th equation enforces continuity of the first derivative
  !  across TAU(I).  It now reads:
  !
  !    - RATIO * DIAG(I-1) * X(I-1) + DIAG(I) * X(I) + U(I) * X(I+1) = R(I).
  !
  !  Eliminate (I-1)st unknown from this equation
  !
  s(i,2) = ratio * s(i-1,3) + s(i,2)
  s(i,4) = ratio * s(i-1,4) + s(i,4)
  !
  !  Set up the part of the next equation which depends on the I-th interval.
  !
200 continue

  if ( z < 0.5D+00 ) then

     ratio = - s(i,6) * s(i,1) / s(i,2)
     s(i+1,2) = s(i,1) / 3.0D+00

  else if ( z == 0.5D+00 ) then

     ratio = - ( s(i,1) / 6.0D+00 ) / s(i,2)
     s(i+1,2) = s(i,1) / 3.0D+00

  else if ( 0.5D+00 < z ) then

     ratio = - ( s(i,1) / 6.0D+00 ) / s(i,2)
     s(i+1,2) = s(i,1) &
          * ( ( 1.0D+00 - zeta * alpha ) * factor / 2.0D+00 - s(i,6) )

  end if
  !
  !  End of I loop.
  !
  i = i + 1

  if ( i < ntau - 1 ) then
     go to 70
  end if

  s(i,5) = 0.5D+00
  !
  !  The last two equations enforce continuity of third derivative and
  !  of first derivative across TAU(NTAU-1).
  !
  entry = ratio * s(i-1,3) + s(i,2) + s(i,1) / 3.0D+00
  s(i+1,2) = s(i,1) / 6.0D+00
  s(i+1,4) = ratio * s(i-1,4) + s(i,4)

  if ( z < 0.5D+00 ) then

     ratio = s(i,1) * 6.0D+00 * s(i-1,6) * alpha / s(i-1,2)
     s(i,2) = ratio * s(i-1,3) + s(i,1) + s(i-1,1)
     s(i,3) = - s(i-1,1)

  else if ( z == 0.5D+00 ) then

     ratio = s(i,1) / s(i-1,2)
     s(i,2) = ratio * s(i-1,3) + s(i,1) + s(i-1,1)
     s(i,3) = - s(i-1,1)

  else if ( 0.5D+00 < z ) then

     factr2 = onemzt * ( alpha * ( onemzt**2 - 1.0D+00 ) + 1.0D+00 ) &
          / ( alpha * ( onemzt**3 - 1.0D+00 ) + 1.0D+00 )

     ratio = factr2 * s(i,1) / s(i-1,2)
     s(i,2) = ratio * s(i-1,3) + factr2 * s(i-1,1) + s(i,1)
     s(i,3) = - factr2 * s(i-1,1)

  end if
  !
  !  At this point, the last two equations read:
  !
  !           DIAG(I)*XI+     U(I)*XI+1 = R(I)
  !    -RATIO*DIAG(I)*XI+DIAG(I+1)*XI+1 = R(I+1)
  !
  !  Eliminate XI from the last equation.
  !
  s(i,4) = ratio * s(i-1,4)
  ratio = - entry / s(i,2)
  s(i+1,2) = ratio * s(i,3) + s(i+1,2)
  s(i+1,4) = ratio * s(i,4) + s(i+1,4)
  !
  !  Back substitution.
  !
  s(ntau,4) = s(ntau,4) / s(ntau,2)

  do while ( 1 < i )

     s(i,4) = ( s(i,4) - s(i,3) * s(i+1,4) ) / s(i,2)
     i = i - 1

  end do

  s(1,4) = ( s(1,4) - s(1,3) * s(2,4) - entry3 * s(3,4) ) / s(1,2)
  !
  !  Construct polynomial pieces. 
  !
  break(1) = tau(1)
  l = 1

  do i = 1, ntau - 1

     coef(1,l) = gtau(i)
     coef(3,l) = s(i,4)
     divdif = ( gtau(i+1) - gtau(i) ) / s(i,1)
     z = s(i,5)

     if ( z == 0.0D+00 ) then

        coef(2,l) = divdif
        coef(3,l) = 0D+00
        coef(4,l) = 0.0D+00

     else if ( z < 0.5D+00 ) then

        zeta = gam * z
        onemzt = 1.0D+00 - zeta
        c = s(i+1,4) / 6.0D+00
        d = s(i,4) * s(i,6)
        l = l + 1
        del = zeta * s(i,1)
        break(l) = tau(i) + del
        zt2 = zeta**2
        alpha = alph(onemzt)
        factor = onemzt**2 * alpha
        coef(1,l) = gtau(i) + divdif * del &
             + s(i,1)**2 * ( d * onemzt * ( factor - 1.0D+00 ) &
             + c * zeta * ( zt2 - 1.0D+00 ) )
        coef(2,l) = divdif + s(i,1) * ( d * ( 1.0D+00 - 3.0D+00 * factor ) &
             + c * ( 3.0D+00 * zt2 - 1.0D+00 ) )
        coef(3,l) = 6.0D+00 * ( d * alpha * onemzt + c * zeta )
        coef(4,l) = 6.0D+00 * ( c - d * alpha ) / s(i,1)
        coef(4,l-1) = coef(4,l) &
             - 6.0D+00 * d * ( 1.0D+00 - alpha ) / ( del * zt2 )
        coef(2,l-1) = coef(2,l) - del * ( coef(3,l) &
             - ( del / 2.0D+00 ) * coef(4,l-1))

     else if ( z == 0.5D+00 ) then

        coef(2,l) = divdif &
             - s(i,1) * ( 2.0D+00 * s(i,4) + s(i+1,4) ) / 6.0D+00
        coef(4,l) = ( s(i+1,4) - s(i,4) ) / s(i,1)

     else if ( 0.5D+00 <= z ) then

        onemzt = gam * ( 1.0D+00 - z )

        if ( onemzt == 0.0D+00 ) then

           coef(2,l) = divdif
           coef(3,l) = 0D+00
           coef(4,l) = 0.0D+00

        else

           zeta = 1.0D+00 - onemzt
           alpha = alph(zeta)
           c = s(i+1,4) * s(i,6)
           d = s(i,4) / 6.0D+00
           del = zeta * s(i,1)
           break(l+1) = tau(i) + del
           coef(2,l) = divdif - s(i,1) * ( 2.0D+00 * d + c )
           coef(4,l) = 6.0D+00 * ( c * alpha - d ) / s(i,1)
           l = l + 1
           coef(4,l) = coef(4,l-1) + 6.0D+00 * ( 1.0D+00 - alpha ) * c &
                / ( s(i,1) * onemzt**3 )
           coef(3,l) = coef(3,l-1) + del * coef(4,l-1)
           coef(2,l) = coef(2,l-1) + del * ( coef(3,l-1) &
                + ( del / 2.0D+00 ) * coef(4,l-1) )
           coef(1,l) = coef(1,l-1) + del * ( coef(2,l-1) &
                + ( del / 2.0D+00 ) * ( coef(3,l-1) &
                + ( del / 3.0D+00 ) * coef(4,l-1) ) )

        end if

     end if

     l = l + 1
     break(l) = tau(i+1)

  end do

  l = l - 1
  k = 4
  iflag = 1

  return
end subroutine tautsp