factrb Subroutine

subroutine factrb ( w, ipivot, d, nrow, ncol, last, iflag )

  !*****************************************************************************80
  !
  !! FACTRB constructs a partial PLU factorization.
  !
  !  Discussion:
  !
  !    This factorization corresponds to steps 1 through LAST in Gauss 
  !    elimination for the matrix W of order ( NROW, NCOL ), using 
  !    pivoting of scaled rows.
  !
  !  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/output, real ( kind = 8 ) W(NROW,NCOL); on input, contains the
  !    matrix to be partially factored; on output, the partial factorization.
  !
  !    Output, integer ( kind = 4 ) IPIVOT(NROW), contains a record of the 
  !    pivoting strategy used; row IPIVOT(I) is used during the I-th elimination
  !    step, for I = 1, ..., LAST.
  !
  !    Workspace, real ( kind = 8 ) D(NROW), used to store the maximum entry
  !    in each row.
  !
  !    Input, integer ( kind = 4 ) NROW, the number of rows of W.
  !
  !    Input, integer ( kind = 4 ) NCOL, the number of columns of W.
  !
  !    Input, integer ( kind = 4 ) LAST, the number of elimination steps to 
  !    be carried out.
  !
  !    Input/output, integer ( kind = 4 ) IFLAG.  On output, equals the input 
  !    value times (-1)**(number of row interchanges during the factorization 
  !    process), in case no zero pivot was encountered.
  !    Otherwise, IFLAG = 0 on output.
  !
  implicit none

  integer ( kind = 4 ) ncol
  integer ( kind = 4 ) nrow

  real ( kind = 8 ) awikdi
  real ( kind = 8 ) colmax
  real ( kind = 8 ) d(nrow)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iflag
  integer ( kind = 4 ) ipivi
  integer ( kind = 4 ) ipivk
  integer ( kind = 4 ) ipivot(nrow)
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kp1
  integer ( kind = 4 ) last
  real ( kind = 8 ) ratio
  real ( kind = 8 ) rowmax
  real ( kind = 8 ) w(nrow,ncol)
  !
  !  Initialize IPIVOT and D.
  !
  do i = 1, nrow
     ipivot(i) = i
  end do

  do i = 1, nrow

     rowmax = maxval ( abs ( w(i,1:ncol) ) )

     if ( rowmax == 0.0D+00 ) then
        iflag = 0
        return
     end if

     d(i) = rowmax

  end do
  !
  !  Gauss elimination with pivoting of scaled rows, loop over K = 1,..., LAST.
  !
  k = 1
  !
  !  As pivot row for K-th step, pick among the rows not yet used,
  !  that is, from rows IPIVOT(K:NROW), the one whose K-th entry, compared 
  !  to the row size, is largest. 
  !
  !  If this row does not turn out to be row IPIVOT(K), redefine IPIVOT(K)
  !  appropriately and record this interchange by changing the sign
  !  of IFLAG.
  !
  do while ( k <= last )

     ipivk = ipivot(k)

     if ( k == nrow ) then
        if ( abs ( w(ipivk,nrow) ) + d(ipivk) <= d(ipivk) ) then
           iflag = 0
        end if
        return
     end if

     j = k
     kp1 = k + 1
     colmax = abs ( w(ipivk,k) ) / d(ipivk)
     !
     !  Find the largest pivot.
     !
     do i = kp1, nrow
        ipivi = ipivot(i)
        awikdi = abs ( w(ipivi,k) ) / d(ipivi)
        if ( colmax < awikdi ) then
           colmax = awikdi
           j = i
        end if
     end do

     if ( j /= k ) then
        ipivk = ipivot(j)
        ipivot(j) = ipivot(k)
        ipivot(k) = ipivk
        iflag = - iflag
     end if
     !
     !  If the pivot element is too small in absolute value, declare
     !  the matrix to be noninvertible and quit.
     !
     if ( abs ( w(ipivk,k) ) + d(ipivk) <= d(ipivk) ) then
        iflag = 0
        return
     end if
     !
     !  Otherwise, subtract the appropriate multiple of the pivot
     !  row from the remaining rows, that is, the rows IPIVOT(K+1:NROW),
     !  to make the K-th entry zero.
     !
     !  Save the multiplier in its place.
     !
     do i = kp1, nrow

        ipivi = ipivot(i)
        w(ipivi,k) = w(ipivi,k) / w(ipivk,k)

        ratio = - w(ipivi,k)
        w(ipivi,kp1:ncol) = ratio * w(ipivk,kp1:ncol) + w(ipivi,kp1:ncol)

     end do

     k = kp1

  end do

  return
end subroutine factrb