subroutine lu1mCP( m , n , lena , aijtol, &
ibest , jbest , mbest , &
a , indc , indr , &
lenc , lenr , locc , &
Hlen , Ha , Hj )
integer(ip), intent(in) :: m, n, lena, Hlen
integer(ip), intent(in) :: indc(lena), indr(lena), &
lenc(n) , lenr(m) , locc(n), Hj(Hlen)
real(rp), intent(in) :: aijtol
real(rp), intent(in) :: a(lena) , Ha(Hlen)
integer(ip), intent(out) :: ibest, jbest, mbest
!------------------------------------------------------------------
! lu1mCP uses a Markowitz criterion to select a pivot element
! for the next stage of a sparse LU factorization,
! subject to a Threshold Complete Pivoting stability criterion (TCP)
! that bounds the elements of L and U.
!
! 09 May 2002: First version of lu1mCP.
! It searches columns only, using the heap that
! holds the largest element in each column.
! 09 May 2002: Current version of lu1mCP.
!
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent.
!------------------------------------------------------------------
integer(ip) :: i, j, kheap, lc, lc1, lc2, len1, lenj, &
maxcol, merit, ncol, nz1
real(rp) :: abest, aij, amax, cmax, lbest
real(rp), parameter :: gamma = 2.0
! gamma is "gamma" in the tie-breaking rule TB4 in the LUSOL paper.
!------------------------------------------------------------------
! Search up to maxcol columns stored at the top of the heap.
! The very top column helps initialize mbest.
!------------------------------------------------------------------
abest = zero
lbest = zero
ibest = 0
jbest = Hj(1) ! Column at the top of the heap
lenj = lenc(jbest)
mbest = lenj * Hlen ! Bigger than any possible merit
maxcol = 40 ! ??? Big question
ncol = 0 ! No. of columns searched
Cols: do kheap = 1, Hlen
amax = Ha(kheap)
if (amax < aijtol) cycle Cols
ncol = ncol + 1
j = Hj(kheap)
!---------------------------------------------------------------
! This column has at least one entry big enough (the top one).
! Search the column for other possibilities.
!---------------------------------------------------------------
lenj = lenc(j)
nz1 = lenj - 1
lc1 = locc(j)
lc2 = lc1 + nz1
!--amax = abs( a(lc1) )
! Test all aijs in this column.
! amax is the largest element (the first in the column).
! cmax is the largest multiplier if aij becomes pivot.
Colj: do lc = lc1, lc2
i = indc(lc)
len1 = lenr(i) - 1
merit = nz1 * len1
if (merit > mbest) cycle Colj
! aij has a promising merit.
if (lc == lc1) then
! This is the maximum element, amax.
! Find the biggest element in the rest of the column
! and hence get cmax. We know cmax .le. 1, but
! we still want it exactly in order to break ties.
! 27 Apr 2002: Settle for cmax = 1.
aij = amax
cmax = one
! cmax = zero
! do 140 l = lc1 + 1, lc2
! cmax = max( cmax, abs( a(l) ) )
! 140 continue
! cmax = cmax / amax
else
! aij is not the biggest element, so cmax .ge. 1.
! Bail out if cmax will be too big.
aij = abs( a(lc) )
if (aij < aijtol) cycle Colj
cmax = amax / aij
end if
! aij is big enough. Its maximum multiplier is cmax.
if (merit == mbest) then
! Break ties.
! (Initializing mbest "too big" prevents getting here if
! nothing has been found yet.)
! In this version we minimize cmax
! but if it is already small we maximize the pivot.
if (lbest <= gamma .and. cmax <= gamma) then
if (abest >= aij ) cycle Colj
else
if (lbest <= cmax) cycle Colj
end if
end if
! aij is the best pivot so far.
ibest = i
jbest = j
mbest = merit
abest = aij
lbest = cmax
if (merit == 0) exit Cols ! Col or row of length 1
end do Colj
if (ncol >= maxcol) exit Cols
end do Cols
end subroutine lu1mCP