eqblok Subroutine

subroutine eqblok ( t, n, kpm, work1, work2, bloks, lenblk, integs, nbloks, b )

  !*****************************************************************************80
  !
  !! EQBLOK is to be called in COLLOC.
  !
  !  Method:
  !
  !    Each breakpoint interval gives rise to a block in the linear system.
  !    This block is determined by the K collocation equations in the interval
  !    with the side conditions, if any, in the interval interspersed
  !    appropriately, and involves the KPM B-splines having the interval in
  !    their support.  Correspondingly, such a block has NROW = K + ISIDEL
  !    rows, with ISIDEL = number of side conditions in this and the 
  !    previous intervals, and NCOL = KPM columns.
  !
  !    Further, because the interior knots have multiplicity K, we can
  !    carry out in SLVBLK K elimination steps in a block before pivoting
  !    might involve an equation from the next block.  In the last block,
  !    of course, all KPM elimination steps will be carried out in SLVBLK.
  !
  !    See the detailed comments in SLVBLK for further
  !    information about the almost block diagonal form used here.
  !
  !  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 ) T(N+KPM), the knot sequence.
  !
  !    Input, integer ( kind = 4 ) N, the dimension of the approximating spline 
  !    space, that is, the order of the linear system to be constructed.
  !
  !    Input, integer ( kind = 4 ) KPM, = K + M, the order of the approximating 
  !    spline.
  !
  !    Input, integer ( kind = 4 ) LENBLK, the maximum length of the array BLOKS,
  !    as allowed by the dimension statement in COLLOC.
  !
  !    Workspace, real ( kind = 8 ) WORK1(KPM,KPM), used in PUTIT.
  !
  !    Workspace, real ( kind = 8 ) WORK2(KPM,M+1), used in PUTIT.
  !
  !    Output, real ( kind = 8 ) BLOKS(*), the coefficient matrix of the
  !    linear system, stored in almost block diagonal form, of size
  !    KPM * sum ( INTEGS(1,1:NBLOKS) ).
  !
  !    Output, integer ( kind = 4 ) INTEGS(3,NBLOKS), describing the block 
  !    structure.
  !    INTEGS(1,I) = number of rows in block I;
  !    INTEGS(2,I) = number of columns in block I;
  !    INTEGS(3,I) = number of elimination steps which can be carried out in 
  !    block I before pivoting might bring in an equation from the next block.
  !
  !    Output, integer ( kind = 4 ) NBLOKS, the number of blocks, equals number
  !    of polynomial pieces.
  !
  !    Output, real ( kind = 8 ) B(*), the right hand side of the linear 
  !    system, stored corresponding to the almost block diagonal form, 
  !    of size sum ( INTEGS(1,1:NBLOKS) ).
  !
  implicit none

  integer ( kind = 4 ) kpm
  integer ( kind = 4 ) n

  real ( kind = 8 ) b(*)
  real ( kind = 8 ) bloks(*)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) index
  integer ( kind = 4 ) indexb
  integer ( kind = 4 ) integs(3,*)
  integer ( kind = 4 ) iside
  integer ( kind = 4 ) isidel
  integer ( kind = 4 ) itermx
  integer ( kind = 4 ) k
  integer ( kind = 4 ) left
  integer ( kind = 4 ) lenblk
  integer ( kind = 4 ) m
  integer ( kind = 4 ) nbloks
  integer ( kind = 4 ) nrow
  real ( kind = 8 ) rho
  real ( kind = 8 ) t(n+kpm)
  real ( kind = 8 ) work1(kpm,kpm)
  real ( kind = 8 ) work2(kpm,*)
  real ( kind = 8 ) xside

  save / other /
  save / side /

  common / other / itermx, k, rho(19)
  common / side / m, iside, xside(10)

  index = 1
  indexb = 1
  i = 0
  iside = 1

  do left = kpm, n, k

     i = i + 1
     !
     !  Determine INTEGS(:,I).
     !
     integs(2,i) = kpm

     if ( n <= left ) then

        integs(3,i) = kpm
        isidel = m
        !
        !  At this point, ISIDE - 1 gives the number of side conditions
        !  incorporated so far.  Adding to this the side conditions in the
        !  current interval gives the number ISIDEL.
        !
     else

        integs(3,i) = k

        isidel = iside - 1

        do

           if ( isidel == m ) then
              exit
           end if

           if ( t(left+1) <= xside(isidel+1) ) then
              exit
           end if

           isidel = isidel + 1

        end do

     end if

     nrow = k + isidel
     integs(1,i) = nrow
     !
     !  The detailed equations for this block are generated and put
     !  together in PUTIT.
     !
     if ( lenblk < index + nrow * kpm - 1 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'EQBLOK - Fatal error!'
        write ( *, '(a)' ) '  The dimension of BLOKS is too small.'
        write ( *, '(a,i8)' ) '  LENBLK = ', lenblk
        stop
     end if

     call putit ( t, kpm, left, work1, work2, bloks(index), nrow, b(indexb) )

     index = index + nrow * kpm
     indexb = indexb + nrow

  end do

  nbloks = i

  return
end subroutine eqblok