subroutine lu1mSP( m , n , lena , maxmn, &
Ltol , maxcol, &
ibest, jbest , mbest , &
a , indc , q , locc , iqloc )
integer(ip), intent(in) :: m, n, lena, maxmn, maxcol
real(rp), intent(in) :: Ltol, a(lena)
integer(ip), intent(in) :: indc(lena), q(n), iqloc(m), locc(n)
integer(ip), intent(out) :: ibest, jbest, mbest
!------------------------------------------------------------------
! lu1mSP is intended for symmetric matrices that are either
! definite or quasi-definite.
! lu1mSP uses a Markowitz criterion to select a pivot element for
! the next stage of a sparse LU factorization of a symmetric matrix,
! subject to a Threshold Symmetric Pivoting stability criterion
! (TSP) restricted to diagonal elements to preserve symmetry.
! This bounds the elements of L and U and should have rank-revealing
! properties analogous to Threshold Rook Pivoting for unsymmetric
! matrices.
!
! 14 Dec 2002: First version of lu1mSP derived from lu1mRP.
! There is no safeguard to ensure that A is symmetric.
! 14 Dec 2002: Current version of lu1mSP.
!
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent.
!------------------------------------------------------------------
integer(ip) :: i, j, kbest, lc, lc1, lc2, &
lq, lq1, lq2, merit, ncol, nz, nz1
real(rp) :: abest, aij, amax, atolj
!------------------------------------------------------------------
! Search cols of length nz = 1, then cols of length nz = 2, etc.
!------------------------------------------------------------------
abest = zero
ibest = 0
kbest = maxmn + 1
mbest = -1
ncol = 0
nz1 = 0
NZS:do nz = 1, maxmn
! nz1 = nz - 1
! if (mbest <= nz1**2) exit
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.
! Ignore everything except the diagonal.
Colj: do lc = lc1, lc2
i = indc(lc)
if (i /= j) cycle Colj ! Skip off-diagonals.
! merit = nz1 * nz1
! if (merit > mbest) cycle
if (nz1 > kbest) cycle Colj
! aij has a promising merit.
! Apply the Threshold Partial Pivoting stability test
! (which is equivalent to Threshold Rook Pivoting for
! symmetric matrices).
! We require aij to be sufficiently large
! compared to other nonzeros in column j.
aij = abs( a(lc) )
if (aij < atolj ) cycle Colj
! aij is big enough.
merit = nz1 * nz1
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 = nz1
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
! See if it's time to quit.
200 if (ibest > 0) then
if (ncol >= maxcol) exit NZS
end if
! Press on with next nz.
nz1 = nz
if (ibest > 0) kbest = mbest / nz1
end do NZS
end subroutine lu1mSP