subroutine lu1mRP( m , n , lena , maxmn, &
Ltol , maxcol, maxrow, &
ibest, jbest , mbest , &
a , indc , indr , p , q, &
lenc , lenr , locc , locr , &
iploc, iqloc , Amaxr )
integer(ip), intent(in) :: m, n, lena, maxmn, maxcol, maxrow
real(rp), intent(in) :: Ltol
integer(ip), intent(in) :: indc(lena), indr(lena), p(m) , q(n) , &
lenc(n) , lenr(m) , iploc(n), iqloc(m), &
locc(n) , locr(m)
real(rp), intent(in) :: a(lena) , Amaxr(m)
integer(ip), intent(out) :: ibest, jbest, mbest
!------------------------------------------------------------------
! lu1mRP uses a Markowitz criterion to select a pivot element
! for the next stage of a sparse LU factorization,
! subject to a Threshold Rook Pivoting stability criterion (TRP)
! that bounds the elements of L and U.
!
! 11 Jun 2002: First version of lu1mRP derived from lu1mar.
! 11 Jun 2002: Current version of lu1mRP.
!
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent.
!------------------------------------------------------------------
integer(ip) :: i, j, kbest, lc, lc1, lc2, len1, &
lp, lp1, lp2, lq, lq1, lq2, lr, lr1, lr2, &
merit, ncol, nrow, nz, nz1
real(rp) :: abest, aij, amax, atoli, atolj
!------------------------------------------------------------------
! Search cols of length nz = 1, then rows of length nz = 1,
! then cols of length nz = 2, then rows of length nz = 2, etc.
!------------------------------------------------------------------
abest = zero
ibest = 0
kbest = maxmn + 1
mbest = -1
ncol = 0
nrow = 0
nz1 = 0
NZS:do nz = 1, maxmn
! nz1 = nz - 1
! if (mbest .le. nz1**2) go to 900
if (kbest <= nz1) exit NZS
if (ibest > 0 ) then
if (ncol >= maxcol) go to 200
end if
if (nz > m ) go to 200
!---------------------------------------------------------------
! Search the set of columns of length nz.
!---------------------------------------------------------------
lq1 = iqloc(nz)
lq2 = n
if (nz < m) lq2 = iqloc(nz + 1) - 1
Cols: do lq = lq1, lq2
ncol = ncol + 1
j = q(lq)
lc1 = locc(j)
lc2 = lc1 + nz1
amax = abs( a(lc1) )
atolj = amax / Ltol ! Min size of pivots in col j
! Test all aijs in this column.
Colj: do lc = lc1, lc2
i = indc(lc)
len1 = lenr(i) - 1
! merit = nz1 * len1
! if (merit > mbest) cycle Colj
if (len1 > kbest) cycle Colj
! aij has a promising merit.
! Apply the Threshold Rook Pivoting stability test.
! First we require aij to be sufficiently large
! compared to other nonzeros in column j.
! Then we require aij to be sufficiently large
! compared to other nonzeros in row i.
aij = abs( a(lc) )
if (aij < atolj ) cycle Colj
if (aij*Ltol < Amaxr(i)) cycle Colj
! aij is big enough.
merit = nz1 * len1
if (merit == mbest) then
! Break ties.
! (Initializing mbest < 0 prevents getting here if
! nothing has been found yet.)
if (abest >= aij) cycle Colj
end if
! aij is the best pivot so far.
ibest = i
jbest = j
kbest = len1
mbest = merit
abest = aij
if (nz == 1) exit NZS
end do Colj
! Finished with that column.
if (ibest > 0) then
if (ncol >= maxcol) exit Cols
end if
end do Cols
!---------------------------------------------------------------
! Search the set of rows of length nz.
!---------------------------------------------------------------
! 200 if (mbest .le. nz*nz1) go to 900
200 if (kbest <= nz ) exit NZS
if (ibest > 0) then
if (nrow >= maxrow) go to 290
end if
if (nz > n) go to 290
lp1 = iploc(nz)
lp2 = m
if (nz < n) lp2 = iploc(nz + 1) - 1
Rows: do lp = lp1, lp2
nrow = nrow + 1
i = p(lp)
lr1 = locr(i)
lr2 = lr1 + nz1
atoli = Amaxr(i) / Ltol ! Min size of pivots in row i
Rowi: do lr = lr1, lr2
j = indr(lr)
len1 = lenc(j) - 1
! merit = nz1 * len1
! if (merit > mbest) cycle
if (len1 > kbest) cycle Rowi
! aij has a promising merit.
! Find where aij is in column j.
lc1 = locc(j)
lc2 = lc1 + len1
amax = abs( a(lc1) )
do lc = lc1, lc2
if (indc(lc) == i) exit
end do
! Apply the Threshold Rook Pivoting stability test.
! First we require aij to be sufficiently large
! compared to other nonzeros in row i.
! Then we require aij to be sufficiently large
! compared to other nonzeros in column j.
aij = abs( a(lc) )
if (aij < atoli) cycle Rowi
if (aij*Ltol < amax ) cycle Rowi
! aij is big enough.
merit = nz1 * len1
if (merit == mbest) then
! Break ties as before.
! (Initializing mbest < 0 prevents getting here if
! nothing has been found yet.)
if (abest >= aij ) cycle Rowi
end if
! aij is the best pivot so far.
ibest = i
jbest = j
kbest = len1
mbest = merit
abest = aij
if (nz == 1) exit NZS
end do Rowi ! This was loop 260
! Finished with that row.
if (ibest > 0) then
if (nrow >= maxrow) exit Rows
end if
end do Rows
! See if it's time to quit.
290 if (ibest > 0) then
if (nrow >= maxrow .and. ncol >= maxcol) exit NZS
end if
! Press on with next nz.
nz1 = nz
if (ibest > 0) kbest = mbest / nz1
end do NZS
end subroutine lu1mRP