subroutine lu1gau( m , melim , ncold , nspare, small , &
lpivc1, lpivc2, lfirst, lpivr2, lfree , minfre, &
ilast , jlast , lrow , lcol , lu , nfill , &
a , indc , indr , &
lenc , lenr , locc , locr , &
mark , al , markl , &
au , ifill , jfill )
integer(ip), intent(in) :: m, melim, ncold, nspare, &
lpivc1, lpivc2, lpivr2, lfree, minfre
integer(ip), intent(in) :: locr(*), mark(*)
real(rp), intent(in) :: small
real(rp), intent(in) :: al(melim), au(ncold)
integer(ip), intent(inout) :: ilast, jlast, lfirst, lrow, lcol, lu, nfill
real(rp), intent(inout) :: a(*)
integer(ip), intent(inout) :: locc(*), indc(*), indr(*), lenc(*), lenr(*), &
markl(melim), ifill(melim), jfill(ncold)
!------------------------------------------------------------------
! lu1gau does most of the work for each step of
! Gaussian elimination.
! A multiple of the pivot column is added to each other column j
! in the pivot row. The column list is fully updated.
! The row list is updated if there is room, but some fill-ins may
! remain, as indicated by ifill and jfill.
!
! Input:
! ilast is the row at the end of the row list.
! jlast is the column at the end of the column list.
! lfirst is the first column to be processed.
! lu + 1 is the corresponding element of U in au(*).
! nfill keeps track of pending fill-in.
! a(*) contains the nonzeros for each column j.
! indc(*) contains the row indices for each column j.
! al(*) contains the new column of L. A multiple of it is
! used to modify each column.
! mark(*) has been set to -1, -2, -3, ... in the rows
! corresponding to nonzero 1, 2, 3, ... of the col of L.
! au(*) contains the new row of U. Each nonzero gives the
! required multiple of the column of L.
!
! Workspace:
! markl(*) marks the nonzeros of L actually used.
! (A different mark, namely j, is used for each column.)
!
! Output:
! ilast New last row in the row list.
! jlast New last column in the column list.
! lfirst = 0 if all columns were completed,
! > 0 otherwise.
! lu returns the position of the last nonzero of U
! actually used, in case we come back in again.
! nfill keeps track of the total extra space needed in the
! row file.
! ifill(ll) counts pending fill-in for rows involved in the new
! column of L.
! jfill(lu) marks the first pending fill-in stored in columns
! involved in the new row of U.
!
! 16 Apr 1989: First version of lu1gau.
! 23 Apr 1989: lfirst, lu, nfill are now input and output
! to allow re-entry if elimination is interrupted.
! 23 Mar 2001: Introduced ilast, jlast.
! 27 Mar 2001: Allow fill-in "in situ" if there is already room
! up to but NOT INCLUDING the end of the
! row or column file.
! Seems safe way to avoid overwriting empty rows/cols
! at the end. (May not be needed though, now that we
! have ilast and jlast.)
!
! 10 Jan 2010: First f90 version.
! 28 Feb 2010: Declare intent and local variables.
!------------------------------------------------------------------
logical :: atend
integer(ip) :: i, j, k, l, l1, l2, last, lc, lc1, lc2, &
leni, lenj, ll, lr, lr1, lrep, &
ndone, ndrop, nfree
real(rp) :: aij, uj
do 600 lr = lfirst, lpivr2
j = indr(lr)
lenj = lenc(j)
nfree = lfree - lcol
if (nfree < minfre) go to 900
!---------------------------------------------------------------
! Inner loop to modify existing nonzeros in column j.
! The "do l = lc1, lc2" loop performs most of the arithmetic
! involved in the whole LU factorization.
! ndone counts how many multipliers were used.
! ndrop counts how many modified nonzeros are negligibly small.
!---------------------------------------------------------------
lu = lu + 1
uj = au(lu)
lc1 = locc(j)
lc2 = lc1 + lenj - 1
atend = j == jlast
ndone = 0
if (lenj == 0) go to 500
ndrop = 0
do l = lc1, lc2
i = indc(l)
ll = - mark(i)
if (ll > 0) then
ndone = ndone + 1
markl(ll) = j
a(l) = a(l) + al(ll) * uj
if (abs( a(l) ) <= small) then
ndrop = ndrop + 1
end if
end if
end do
!---------------------------------------------------------------
! Remove any negligible modified nonzeros from both
! the column file and the row file.
!---------------------------------------------------------------
if (ndrop == 0) go to 500
k = lc1
do l = lc1, lc2
i = indc(l)
if (abs( a(l) ) > small) then
a(k) = a(l)
indc(k) = i
k = k + 1
cycle
end if
! Delete the nonzero from the row file.
lenj = lenj - 1
lenr(i) = lenr(i) - 1
lr1 = locr(i)
last = lr1 + lenr(i)
do lrep = lr1, last
if (indr(lrep) == j) exit
end do
indr(lrep) = indr(last)
indr(last) = 0
if (i == ilast) lrow = lrow - 1
end do
! Free the deleted elements from the column file.
do l = k, lc2
indc(l) = 0
end do
if (atend) lcol = k - 1
!---------------------------------------------------------------
! Deal with the fill-in in column j.
!---------------------------------------------------------------
500 if (ndone == melim) go to 590
! See if column j already has room for the fill-in.
if (atend) go to 540
last = lc1 + lenj - 1
l1 = last + 1
l2 = last + (melim - ndone)
! 27 Mar 2001: Be sure it's not at or past end of the col file.
if (l2 >= lcol) go to 520
do l = l1, l2
if (indc(l) /= 0) go to 520
end do
go to 540
! We must move column j to the end of the column file.
! First, leave some spare room at the end of the
! current last column.
! 14 Jul 2015: (William Gandler) Fix deceptive loop
! do l = lcol + 1, lcol + nspare
! lcol = l
520 l1 = lcol + 1
l2 = lcol + nspare
do l = l1, l2
! lcol = l
indc(l) = 0 ! Spare space is free.
end do
lcol = l2
atend = .true.
jlast = j
l1 = lc1
lc1 = lcol + 1
locc(j) = lc1
do l = l1, last
lcol = lcol + 1
a(lcol) = a(l)
indc(lcol) = indc(l)
indc(l) = 0 ! Free space.
end do
!---------------------------------------------------------------
! Inner loop for the fill-in in column j.
! This is usually not very expensive.
!---------------------------------------------------------------
540 last = lc1 + lenj - 1
ll = 0
do lc = lpivc1, lpivc2
ll = ll + 1
if (markl(ll) == j ) cycle
aij = al(ll)*uj
if (abs(aij) <= small) cycle
lenj = lenj + 1
last = last + 1
a(last) = aij
i = indc(lc)
indc(last) = i
leni = lenr(i)
! Add 1 fill-in to row i if there is already room.
! 27 Mar 2001: Be sure it's not at or past the end
! of the row file.
l = locr(i) + leni
if (l < lrow .and. indr(l) <= 0) then
indr(l) = j
lenr(i) = leni + 1
else
! Row i does not have room for the fill-in.
! Increment ifill(ll) to count how often this has
! happened to row i. Also, add m to the row index
! indc(last) in column j to mark it as a fill-in that is
! still pending.
! If this is the first pending fill-in for row i,
! nfill includes the current length of row i
! (since the whole row has to be moved later).
! If this is the first pending fill-in for column j,
! jfill(lu) records the current length of column j
! (to shorten the search for pending fill-ins later).
if (ifill(ll) == 0) nfill = nfill + leni + nspare
if (jfill(lu) == 0) jfill(lu) = lenj
nfill = nfill + 1
ifill(ll) = ifill(ll) + 1
indc(last) = m + i
end if
end do
if ( atend ) lcol = last
! End loop for column j. Store its final length.
590 lenc(j) = lenj
600 end do
! Successful completion.
lfirst = 0
return
! Interruption. We have to come back in after the
! column file is compressed. Give lfirst a new value.
! lu and nfill will retain their current values.
900 lfirst = lr
return
end subroutine lu1gau