subroutine lu1DPP( a, lda, m, n, small, nsing, ipvt, q )
integer(ip), intent(in) :: lda, m, n
real(rp), intent(in) :: small
integer(ip), intent(inout) :: q(n)
real(rp), intent(inout) :: a(lda,n)
integer(ip), intent(out) :: nsing ! not used outside
integer(ip), intent(out) :: ipvt(m)
!------------------------------------------------------------------
! lu1DPP factors a dense m x n matrix A by Gaussian elimination,
! using row interchanges for stability, as in dgefa from LINPACK.
! This version also uses column interchanges if all elements in a
! pivot column are smaller than (or equal to) "small". Such columns
! are changed to zero and permuted to the right-hand end.
!
! As in LINPACK, ipvt(*) keeps track of pivot rows.
! Rows of U are interchanged, but we don't have to physically
! permute rows of L. In contrast, column interchanges are applied
! directly to the columns of both L and U, and to the column
! permutation vector q(*).
!
! 02 May 1989: First version derived from dgefa
! in LINPACK (version dated 08/14/78).
! 05 Feb 1994: Generalized to treat rectangular matrices
! and use column interchanges when necessary.
! ipvt is retained, but column permutations are applied
! directly to q(*).
! 21 Dec 1994: Bug found via example from Steve Dirkse.
! Loop 100 added to set ipvt(*) for singular rows.
! 26 Mar 2006: nsing redefined (see below).
! Changed to implicit none.
!
! 10 Jan 2010: First f90 version. Need to do more f90-ing.
! 12 Dec 2011: Declare intent and local variables.
! 03 Feb 2012: a is intent(inout), not (out).
! a(kp1:m,j) = t*a(kp1:m,k) + a(kp1:m,j) needs the last :m
!------------------------------------------------------------------
!
! On entry:
!
! a Array holding the matrix A to be factored.
! lda The leading dimension of the array a.
! m The number of rows in A.
! n The number of columns in A.
! small A drop tolerance. Must be zero or positive.
!
! On exit:
!
! a An upper triangular matrix and the multipliers
! which were used to obtain it.
! The factorization can be written A = L*U where
! L is a product of permutation and unit lower
! triangular matrices and U is upper triangular.
! nsing Number of singularities detected.
! 26 Mar 2006: nsing redefined to be more meaningful.
! Users may define rankU = n - nsing and regard
! U as upper-trapezoidal, with the first rankU columns
! being triangular and the rest trapezoidal.
! It would be better to return rankU, but we still
! return nsing for compatibility (even though lu1fad
! no longer uses it).
! ipvt Records the pivot rows.
! q A vector to which column interchanges are applied.
!------------------------------------------------------------------
integer(ip) :: i, j, k, kp1, l, last, lencol, rankU
real(rp) :: t
rankU = 0
k = 1
last = n
!------------------------------------------------------------------
! Start of elimination loop.
!------------------------------------------------------------------
10 kp1 = k + 1
lencol = m - k + 1
! Find l, the pivot row.
l = jdamax( lencol, a(k:m,k), i1 ) + k - 1
ipvt(k) = l
if (abs( a(l,k) ) <= small) then
!==============================================================
! Do column interchange, changing old pivot column to zero.
! Reduce "last" and try again with same k.
!==============================================================
j = q(last)
q(last) = q(k)
q(k) = j
do i = 1, k - 1
t = a(i,last)
a(i,last) = a(i,k)
a(i,k) = t
end do
do i = k, m
t = a(i,last)
a(i,last) = zero
a(i,k) = t
end do
last = last - 1
if (k <= last) go to 10
else
rankU = rankU + 1
if (k < m) then
!===========================================================
! Do row interchange if necessary.
!===========================================================
if (l /= k) then
t = a(l,k)
a(l,k) = a(k,k)
a(k,k) = t
end if
!===========================================================
! Compute multipliers.
! Do row elimination with column indexing.
!===========================================================
t = - one / a(k,k)
! call dscal ( m-k, t, a(kp1,k), i1 )
a(kp1:m,k) = t*a(kp1:m,k)
do j = kp1, last
t = a(l,j)
if (l /= k) then
a(l,j) = a(k,j)
a(k,j) = t
end if
! call daxpy ( m-k, t, a(kp1,k), i1, a(kp1,j), i1 )
a(kp1:m,j) = t*a(kp1:m,k) + a(kp1:m,j)
end do
k = k + 1
if (k <= last) go to 10
end if
end if
! Set ipvt(*) for singular rows.
do k = last + 1, m
ipvt(k) = k
end do
nsing = n - rankU
end subroutine lu1DPP