subroutine lu1DCP( 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)
!------------------------------------------------------------------
! lu1DCP factors a dense m x n matrix A by Gaussian elimination,
! using Complete Pivoting (row and column interchanges) for
! stability.
! 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's dgefa, 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(*).
!
! 01 May 2002: First dense Complete Pivoting, derived from lu1DPP.
! 07 May 2002: Another break needed at end of first loop.
! 26 Mar 2006: Cosmetic mods while looking for "nsing" bug when m<n.
! nsing redefined (see below).
! Changed to implicit none.
!
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent and local variables.
! 03 Feb 2012: a(kp1:m,j) = t*a(kp1:m,k) + a(kp1:m,j) needs the last :m
! 21 Dec 2015: t = 0 caused divide by zero.
! Add test to exit if aijmax <= small.
!------------------------------------------------------------------
!
! 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
! that 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.
!------------------------------------------------------------------
real(rp) :: aijmax, ajmax, t
integer(ip) :: i, imax, j, jlast, jmax, jnew, &
k, kp1, l, last, lencol, rankU
rankU = 0
lencol = m + 1
last = n
!-----------------------------------------------------------------
! Start of elimination loop.
!-----------------------------------------------------------------
do k = 1, n
kp1 = k + 1
lencol = lencol - 1
! Find the biggest aij in row imax and column jmax.
aijmax = zero
imax = k
jmax = k
jlast = last
do j = k, jlast
10 l = jdamax( lencol, a(k:m,j), i1 ) + k - 1
ajmax = abs(a(l,j))
if (ajmax <= small) then
!========================================================
! Do column interchange, changing old column to zero.
! Reduce "last" and try again with same j.
!========================================================
jnew = q(last)
q(last) = q(j)
q(j) = jnew
do i = 1, k - 1
t = a(i,last)
a(i,last) = a(i,j)
a(i,j) = t
end do
do i = k, m
t = a(i,last)
a(i,last) = zero
a(i,j) = t
end do
last = last - 1
if (j <= last) go to 10 ! repeat
go to 200 ! break
end if
! Check if this column has biggest aij so far.
if (aijmax < ajmax) then
aijmax = ajmax
imax = l
jmax = j
end if
if (j >= last) go to 200 ! break
end do
200 ipvt(k) = imax
! 21 Dec 2015: Exit if aijmax is essentially zero.
if (aijmax <= small) go to 500
rankU = rankU + 1
if (jmax /= k) then ! Do column interchange (k and jmax).
jnew = q(jmax)
q(jmax) = q(k)
q(k) = jnew
do i = 1, m
t = a(i,jmax)
a(i,jmax) = a(i,k)
a(i,k) = t
end do
end if
if (k < m) then ! Do row interchange if necessary.
t = a(imax,k)
if (imax /= k) then
a(imax,k) = a(k,k)
a(k,k) = t
end if
!===========================================================
! Compute multipliers.
! Do row elimination with column indexing.
!===========================================================
t = - one / t
! call dscal ( m-k, t, a(kp1,k), i1 )
a(kp1:m,k) = t*a(kp1:m,k)
do j = kp1, last
t = a(imax,j)
if (imax /= k) then
a(imax,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
else
go to 500 ! break
end if
if (k >= last) go to 500 ! break
end do
! Set ipvt(*) for singular rows.
500 do k = last + 1, m
ipvt(k) = k
end do
nsing = n - rankU
end subroutine lu1DCP