subroutine lu7for( m, n, kfirst, klast, lena, luparm, parmlu, &
lenL, lenU, lrow, &
a, indc, indr, p, q, lenr, locc, locr, &
inform, diag )
integer(ip), intent(in) :: m, n, kfirst, klast, lena
integer(ip), intent(in) :: q(n)
integer(ip), intent(inout) :: luparm(30), lenL, lenU, lrow
integer(ip), intent(inout) :: indc(lena), indr(lena), &
p(m), lenr(m), locc(n), locr(m)
real(rp), intent(inout) :: parmlu(30), a(lena)
integer(ip), intent(out) :: inform
real(rp), intent(out) :: diag
!------------------------------------------------------------------
! lu7for (forward sweep) updates the LU factorization A = L*U
! when row iw = p(klast) of U is eliminated by a forward
! sweep of stabilized row operations, leaving p*U*q upper triangular.
!
! The row permutation p is updated to preserve stability and/or
! sparsity. The column permutation q is not altered.
!
! kfirst is such that row p(kfirst) is the first row involved
! in eliminating row iw. (Hence, kfirst marks the first nonzero
! in row iw in pivotal order.) If kfirst is unknown it may be
! input as 1.
!
! klast is such that row p(klast) is the row being eliminated.
! klast is not altered.
!
! lu7for should be called only if kfirst .le. klast.
! If kfirst = klast, there are no nonzeros to eliminate, but the
! diagonal element of row p(klast) may need to be moved to the
! front of the row.
!
! On entry, locc(*) must be zero.
!
! On exit:
! inform = 0 if row iw has a nonzero diagonal (could be small).
! inform = 1 if row iw has no diagonal.
! inform = 7 if there is not enough storage to finish the update.
!
! On a successful exit (inform le 1), locc(*) will again be zero.
!
! Jan 1985: Final f66 version.
! 09 May 1988: First f77 version.
! 13 Dec 2011: First f90 version.
! 20 Dec 2015: ilast is now output by lu1rec.
!------------------------------------------------------------------
logical :: swappd
integer(ip) :: ilast, iv, iw, j, jfirst, jlast, jv, &
k, kbegin, kstart, kstop, &
l, ldiag, lenv, lenw, lfirst, limit, &
lv, lv1, lv2, lv3, lw, lw1, lw2, &
minfre, nfree
real(rp) :: amult, Ltol, Uspace, small, vj, wj
Ltol = parmlu(2)
small = parmlu(3)
Uspace = parmlu(6)
kbegin = kfirst
swappd = .false.
! We come back here from below if a row interchange is performed.
100 iw = p(klast)
lenw = lenr(iw)
if (lenw == 0 ) go to 910
lw1 = locr(iw)
lw2 = lw1 + lenw - 1
jfirst = q(kbegin)
if (kbegin >= klast) go to 700
! Make sure there is room at the end of the row file
! in case row iw is moved there and fills in completely.
minfre = n + 1
nfree = lena - lenL - lrow
if (nfree < minfre) then
call lu1rec( m, .true., luparm, lrow, ilast, &
lena, a, indr, lenr, locr )
lw1 = locr(iw)
lw2 = lw1 + lenw - 1
nfree = lena - lenL - lrow
if (nfree < minfre) go to 970
end if
! Set markers on row iw.
do l = lw1, lw2
j = indr(l)
locc(j) = l
end do
!==================================================================
! Main elimination loop.
!==================================================================
kstart = kbegin
kstop = min( klast, n )
do k = kstart, kstop
jfirst = q(k)
lfirst = locc(jfirst)
if (lfirst == 0) go to 490
! Row iw has its first element in column jfirst.
wj = a(lfirst)
if (k == klast) go to 490
!---------------------------------------------------------------
! We are about to use the first element of row iv
! to eliminate the first element of row iw.
! However, we may wish to interchange the rows instead,
! to preserve stability and/or sparsity.
!---------------------------------------------------------------
iv = p(k)
lenv = lenr(iv)
lv1 = locr(iv)
vj = zero
if (lenv == 0 ) go to 150
if (indr(lv1) /= jfirst) go to 150
vj = a(lv1)
if ( swappd ) go to 200
if (Ltol * abs(wj) < abs(vj)) go to 200
if (Ltol * abs(vj) < abs(wj)) go to 150
if ( lenv <= lenw ) go to 200
!---------------------------------------------------------------
! Interchange rows iv and iw.
!---------------------------------------------------------------
150 p(klast) = iv
p(k) = iw
kbegin = k
swappd = .true.
go to 600
!---------------------------------------------------------------
! Delete the eliminated element from row iw
! by overwriting it with the last element.
!---------------------------------------------------------------
200 a(lfirst) = a(lw2)
jlast = indr(lw2)
indr(lfirst) = jlast
indr(lw2) = 0
locc(jlast) = lfirst
locc(jfirst) = 0
lenw = lenw - 1
lenU = lenU - 1
if (lrow == lw2) lrow = lrow - 1
lw2 = lw2 - 1
!---------------------------------------------------------------
! Form the multiplier and store it in the L file.
!---------------------------------------------------------------
if (abs(wj) <= small) go to 490
amult = - wj/vj
l = lena - lenL
a(l) = amult
indr(l) = iv
indc(l) = iw
lenL = lenL + 1
!---------------------------------------------------------------
! Add the appropriate multiple of row iv to row iw.
! We use two different inner loops. The first one is for the
! case where row iw is not at the end of storage.
!---------------------------------------------------------------
if (lenv == 1) go to 490
lv2 = lv1 + 1
lv3 = lv1 + lenv - 1
if (lw2 == lrow) go to 400
!...............................................................
! This inner loop will be interrupted only if
! fill-in occurs enough to bump into the next row.
!...............................................................
do lv = lv2, lv3
jv = indr(lv)
lw = locc(jv)
if (lw > 0) then ! No fill-in.
a(lw) = a(lw) + amult*a(lv)
if (abs(a(lw)) <= small) then ! Delete small computed element.
a(lw) = a(lw2)
j = indr(lw2)
indr(lw) = j
indr(lw2) = 0
locc(j) = lw
locc(jv) = 0
lenU = lenU - 1
lenw = lenw - 1
lw2 = lw2 - 1
end if
else ! Row iw doesn't have an element in column jv yet
! so there is a fill-in.
if (indr(lw2+1) /= 0) go to 360
lenU = lenU + 1
lenw = lenw + 1
lw2 = lw2 + 1
a(lw2) = amult * a(lv)
indr(lw2) = jv
locc(jv) = lw2
end if
end do
go to 490
! Fill-in interrupted the previous loop.
! Move row iw to the end of the row file.
360 lv2 = lv
locr(iw) = lrow + 1
do l = lw1, lw2
lrow = lrow + 1
a(lrow) = a(l)
j = indr(l)
indr(l) = 0
indr(lrow) = j
locc(j) = lrow
end do
lw1 = locr(iw)
lw2 = lrow
!...............................................................
! Inner loop with row iw at the end of storage.
!...............................................................
400 do lv = lv2, lv3
jv = indr(lv)
lw = locc(jv)
if (lw > 0) then ! No fill-in
a(lw) = a(lw) + amult*a(lv)
if (abs(a(lw)) <= small) then ! Delete small computed element
a(lw) = a(lw2)
j = indr(lw2)
indr(lw) = j
indr(lw2) = 0
locc(j) = lw
locc(jv) = 0
lenU = lenU - 1
lenw = lenw - 1
lw2 = lw2 - 1
end if
else ! Row iw doesn't have an element in column jv yet
! so there is a fill-in
lenU = lenU + 1
lenw = lenw + 1
lw2 = lw2 + 1
a(lw2) = amult * a(lv)
indr(lw2) = jv
locc(jv) = lw2
end if
end do
lrow = lw2
! The k-th element of row iw has been processed.
! Reset swappd before looking at the next element.
490 swappd = .false.
end do
!=================================================================
! End of main elimination loop.
!==================================================================
! Cancel markers on row iw.
600 lenr(iw) = lenw
if (lenw == 0) go to 910
do l = lw1, lw2
j = indr(l)
locc(j) = 0
end do
! Move the diagonal element to the front of row iw.
! At this stage, lenw > 0 and klast <= n.
700 do l = lw1, lw2
ldiag = l
if (indr(l) == jfirst) go to 730 ! not exit !!!
end do
go to 910
730 diag = a(ldiag)
a(ldiag) = a(lw1)
a(lw1) = diag
indr(ldiag) = indr(lw1)
indr(lw1) = jfirst
! If an interchange is needed, repeat from the beginning with the
! new row iw, knowing that the opposite interchange cannot occur.
if (swappd) go to 100
inform = 0
go to 950
! Singular
910 diag = zero
inform = 1
! Force a compression if the file for U is much longer than the
! no. of nonzeros in U (i.e. if lrow is much bigger than lenU).
! This should prevent memory fragmentation when there is far more
! memory than necessary (i.e. when lena is huge).
950 limit = int(Uspace*real(lenU)) + m + n + 1000
if (lrow > limit) then
call lu1rec( m, .true., luparm, lrow, ilast, &
lena, a, indr, lenr, locr )
end if
go to 990
! Not enough storage.
970 inform = 7
! Exit.
990 return
end subroutine lu7for