subroutine lu1rec( n, reals, luparm, ltop, ilast, &
lena, a, ind, lenc, locc )
logical, intent(in) :: reals
integer(ip), intent(in) :: n, lena
integer(ip), intent(inout) :: ltop
integer(ip), intent(out) :: ilast
integer(ip), intent(inout) :: luparm(30), ind(lena), lenc(n), locc(n)
real(rp), intent(inout) :: a(lena)
!------------------------------------------------------------------
! lu1rec recovers space in the column or row lists.
! 00 Jun 1983: Original version of lu1rec followed John Reid's
! compression routine in LA05. It recovered space
! in ind(*) and optionally a(*) by eliminating entries
! with ind(l) = 0.
! The elements of ind(*) could not be negative.
! If len(i) was positive, entry i contained
! that many elements, starting at loc(i).
! Otherwise, entry i was eliminated.
!
! 23 Mar 2001: Realised we could have len(i) = 0 in rare cases!
! (Mostly during TCP when the pivot row contains
! a column of length 1 that couldn't be a pivot.)
! Revised storage scheme to
! keep entries with ind(l) > 0,
! squeeze out entries with -n <= ind(l) <= 0,
! and to allow len(i) = 0.
! Empty items are moved to the end of the compressed
! ind(*) and/or a(*) arrays are given one empty space.
! Items with len(i) < 0 are still eliminated.
!
! 27 Mar 2001: Decided to use only ind(l) > 0 and = 0 in lu1fad.
! Still have to keep entries with len(i) = 0.
!
! On exit:
! ltop is the length of useful entries in ind(*), a(*).
! ind(ltop+1) is "i=ilast" such that len(i), loc(i) belong to the
! last item in ind(*), a(*).
!
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent and local variables.
! 20 Dec 2015: ilast is output instead of ind(ltop+1).
!------------------------------------------------------------------
integer(ip) :: i, k, klast, l, leni, lprint, nempty, nout
nempty = 0
do i = 1, n
leni = lenc(i)
if (leni > 0) then
l = locc(i) + leni - 1
lenc(i) = ind(l)
ind(l) = - (n + i)
else if (leni == 0) then
nempty = nempty + 1
end if
end do
k = 0
klast = 0 ! Previous k
ilast = 0 ! Last entry moved.
do l = 1, ltop
i = ind(l)
if (i > 0) then
k = k + 1
ind(k) = i
if (reals) a(k) = a(l)
else if (i < -n) then ! This is the end of entry i.
i = - (i + n)
ilast = i
k = k + 1
ind(k) = lenc(i)
if (reals) a(k) = a(l)
locc(i) = klast + 1
lenc(i) = k - klast
klast = k
end if
end do
! Move any empty items to the end, adding 1 free entry for each.
if (nempty > 0) then
do i = 1, n
if (lenc(i) == 0) then
k = k + 1
locc(i) = k
ind(k) = 0
ilast = i
end if
end do
end if
nout = luparm(1)
lprint = luparm(2)
if (lprint >= 50) write(nout, 1000) ltop, k, reals, nempty
luparm(26) = luparm(26) + 1 ! ncp
! 20 Dec 2015: Return ilast itself instead of ind(ltop + 1).
ltop = k
! ind(ltop+1) = ilast
return
1000 format(' lu1rec. File compressed from', i10, ' to', i10, l3, ' nempty =', i8)
end subroutine lu1rec