newnot Subroutine

subroutine newnot ( break, coef, l, k, brknew, lnew, coefg )

  !*****************************************************************************80
  !
  !! NEWNOT returns LNEW+1 knots which are equidistributed on (A,B).
  !
  !  Discussion:
  !
  !    The knots are equidistributed on (A,B) = ( BREAK(1), BREAK(L+1) ) 
  !    with respect to a certain monotone function G related to the K-th root of 
  !    the K-th derivative of the piecewise polynomial function F whose 
  !    piecewise polynomial representation is contained in BREAK, COEF, L, K.
  !
  !  Method:
  !
  !    The K-th derivative of the given piecewise polynomial function F does 
  !    not exist, except perhaps as a linear combination of delta functions. 
  !
  !    Nevertheless, we construct a piecewise constant function H with 
  !    breakpoint sequence BREAK which is approximately proportional 
  !    to abs ( d**K(F) ).
  !
  !    Specifically, on (BREAK(I), BREAK(I+1)),
  !
  !          abs(jump at BREAK(I) of PC)     abs(jump at BREAK(I+1) of PC)
  !      H = ---------------------------  +  ----------------------------
  !          BREAK(I+1) - BREAK(I-1)         BREAK(I+2) - BREAK(I)
  !
  !    with PC the piecewise constant (K-1)st derivative of F.
  !
  !    Then, the piecewise linear function G is constructed as
  !
  !      G(X) = integral ( A <= Y <= X )  H(Y)**(1/K) dY,
  !
  !    and its piecewise polynomial coefficients are stored in COEFG.
  !
  !    Then BRKNEW is determined by
  !
  !      BRKNEW(I) = A + G**(-1)((I-1)*STEP), for I = 1:LNEW+1,
  !
  !    where STEP = G(B) / LNEW and (A,B) = ( BREAK(1), BREAK(L+1) ).
  !
  !    In the event that PC = d**(K-1)(F) is constant in ( A, B ) and
  !    therefore H = 0 identically, BRKNEW is chosen uniformly spaced.
  !
  !    If IPRINT is set positive, then the piecewise polynomial coefficients
  !    of G will be printed out.
  !
  !  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 ) BREAK(L+1), real ( kind = 8 ) COEF(K,L), 
  !    integer ( kind = 4 ) L, integer K, the piecewise polynomial representation
  !    of a certain function F of order K.  Specifically,
  !      d**(k-1) F(X) = COEF(K,I) for BREAK(I) <= X < BREAK(I+1).
  !
  !    Input, integer ( kind = 4 ) LNEW, the number of intervals into which the 
  !    interval (A,B) is to be divided by the new breakpoint sequence BRKNEW. 
  !
  !    Output, real ( kind = 8 ) BRKNEW(LNEW+1), the new breakpoint sequence.
  !
  !    Output, real ( kind = 8 ) COEFG(2,L), the coefficient part of the piecewise
  !    polynomial representation BREAK, COEFG, L, 2 for the monotone piecewise 
  !    linear function G with respect to which BRKNEW will be equidistributed.
  !
  implicit none

  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) lnew

  real ( kind = 8 ) break(l+1)
  real ( kind = 8 ) brknew(lnew+1)
  real ( kind = 8 ) coef(k,l)
  real ( kind = 8 ) coefg(2,l)
  real ( kind = 8 ) dif
  real ( kind = 8 ) difprv
  integer ( kind = 4 ) i
  integer ( kind = 4 ), save :: iprint = 0
  integer ( kind = 4 ) j
  real ( kind = 8 ) oneovk
  real ( kind = 8 ) step
  real ( kind = 8 ) stepi
  !
  !  If G is constant, BRKNEW is uniform.
  !
  if ( l <= 1 ) then
     call evnnot ( break, coef, l, k, brknew, lnew, coefg )
     return
  end if

  brknew(1) = break(1)
  brknew(lnew+1) = break(l+1)
  !
  !  Construct the continuous piecewise linear function G.
  !
  oneovk = 1.0D+00 / real ( k, kind = 8 )
  coefg(1,1) = 0.0D+00
  difprv = abs ( coef(k,2) - coef(k,1) ) / ( break(3) - break(1) )

  do i = 2, l
     dif = abs ( coef(k,i) - coef(k,i-1) ) / ( break(i+1) - break(i-1) )
     coefg(2,i-1) = ( dif + difprv )**oneovk
     coefg(1,i) = coefg(1,i-1) + coefg(2,i-1) * ( break(i) - break(i-1) )
     difprv = dif
  end do

  coefg(2,l) = ( 2.0D+00 * difprv )**oneovk
  !
  !  STEP = G(B) / LNEW.
  !
  step = ( coefg(1,l) + coefg(2,l) * ( break(l+1) - break(l) ) ) &
       / real ( lnew, kind = 8 )

  if ( 0 < iprint ) then
     write ( *, '(2x,e16.7)' ) step
     do i = 1, l
        write ( *, '(i5,2e16.5)' ) i, coefg(1:2,i)
     end do
  end if
  !
  !  If G is constant, BRKNEW is uniform.
  !
  if ( step <= 0.0D+00 ) then
     call evnnot ( break, coef, l, k, brknew, lnew, coefg )
     return
  end if
  !
  !  For I = 2,..., LNEW, construct  BRKNEW(I) = A + G**(-1)(STEPI),
  !  with STEPI = ( I - 1 ) * STEP.  
  !
  !  This requires inversion of the piecewise linear function G.  
  !
  !  For this, J is found so that
  !
  !    G(BREAK(J)) <= STEPI <= G(BREAK(J+1))
  !
  !  and then
  !
  !    BRKNEW(I) = BREAK(J) + ( STEPI - G(BREAK(J)) ) / DG(BREAK(J) ).
  !
  !  The midpoint is chosen if DG(BREAK(J)) = 0.
  !
  j = 1

  do i = 2, lnew

     stepi = real ( i - 1, kind = 8 ) * step

     do

        if ( j == l ) then
           exit
        end if

        if ( stepi <= coefg(1,j+1) ) then
           exit
        end if

        j = j + 1

     end do

     if ( coefg(2,j) /= 0.0D+00 ) then
        brknew(i) = break(j) + ( stepi - coefg(1,j) ) / coefg(2,j)
     else
        brknew(i) = ( break(j) + break(j+1) ) / 2.0D+00
     end if

  end do

  return
end subroutine newnot