subroutine lu6chk( mode, m, n, nslack, w, lena, luparm, parmlu, &
a, indc, indr, p, q, lenc, lenr, locc, locr, inform )
integer(ip), intent(in) :: mode, m, n, nslack, lena
integer(ip), intent(in) :: indc(lena), indr(lena), p(m), q(n), &
lenc(n), lenr(m), locc(n), locr(m)
real(rp), intent(in) :: a(lena)
integer(ip), intent(inout) :: inform
integer(ip), intent(inout) :: luparm(30)
real(rp), intent(inout) :: parmlu(30)
real(rp), intent(inout) :: w(n)
!-----------------------------------------------------------------
! lu6chk looks at the LU factorization A = L*U.
!
! If mode = 1, lu6chk is being called by lu1fac.
! (Other modes not yet implemented.)
! The important input parameters are
!
! lprint = luparm(2)
! luparm(6) = 1 if TRP
! keepLU = luparm(8)
! Utol1 = parmlu(4)
! Utol2 = parmlu(5)
!
! and the significant output parameters are
!
! inform = luparm(10)
! nsing = luparm(11)
! jsing = luparm(12)
! jumin = luparm(19)
! Lmax = parmlu(11)
! Umax = parmlu(12)
! DUmax = parmlu(13)
! DUmin = parmlu(14)
! and w(*).
!
! Lmax and Umax return the largest elements in L and U.
! DUmax and DUmin return the largest and smallest diagonals of U
! (excluding diagonals that are exactly zero).
!
! In general, w(j) is set to the maximum absolute element in
! the j-th column of U. However, if the corresponding diagonal
! of U is small in absolute terms or relative to w(j)
! (as judged by the parameters Utol1, Utol2 respectively),
! then w(j) is changed to - w(j).
!
! Thus, if w(j) is not positive, the j-th column of A
! appears to be dependent on the other columns of A.
! The number of such columns, and the position of the last one,
! are returned as nsing and jsing.
!
! Note that nrank is assumed to be set already, and is not altered.
! Typically, nsing will satisfy nrank + nsing = n, but if
! Utol1 and Utol2 are rather large, nsing > n - nrank may occur.
!
! If keepLU = 0,
! Lmax and Umax are already set by lu1fac.
! The diagonals of U are in the top of A.
! Only Utol1 is used in the singularity test to set w(*).
!
! inform = 0 if A appears to have full column rank (nsing = 0).
! inform = 1 otherwise (nsing > 0).
!
! 00 Jul 1987: Early version.
! 09 May 1988: f77 version.
! 11 Mar 2001: Allow for keepLU = 0.
! 17 Nov 2001: Briefer output for singular factors.
! 05 May 2002: Comma needed in format 1100 (via Kenneth Holmstrom).
! 06 May 2002: With keepLU = 0, diags of U are in natural order.
! They were not being extracted correctly.
! 23 Apr 2004: TRP can judge singularity better by comparing
! all diagonals to DUmax.
! 27 Jun 2004: (PEG) Allow write only if nout > 0 .
! 13 Dec 2011: First f90 version.
! 12 Dec 2015: nslack ensures slacks are kept with w(j) > 0.
!------------------------------------------------------------------
character(1) :: mnkey
logical :: keepLU, TRP
integer(ip) :: i, j, jsing, jumin, k, l, l1, l2, ldiagU, lenL, &
lprint, ndefic, nout, nrank, nsing
real(rp) :: aij, diag, DUmax, DUmin, Lmax, Umax, Utol1, Utol2
nout = luparm(1)
lprint = luparm(2)
TRP = luparm(6) == 1 ! Threshold Rook Pivoting
keepLU = luparm(8) /= 0
nrank = luparm(16)
lenL = luparm(23)
Utol1 = parmlu(4)
Utol2 = parmlu(5)
inform = 0
Lmax = zero
Umax = zero
nsing = 0
jsing = 0
jumin = 0
DUmax = zero
DUmin = 1.0d+30
! w(j) is already set by lu1slk.
! w(1:n) = zero
if (keepLU) then
!--------------------------------------------------------------
! Find Lmax.
!--------------------------------------------------------------
do l = lena + 1 - lenL, lena
Lmax = max( Lmax, abs(a(l)) )
end do
!--------------------------------------------------------------
! Find Umax and set w(j) = maximum element in j-th column of U.
!--------------------------------------------------------------
do k = nslack + 1, nrank ! 12 Dec 2015: Allow for nslack.
i = p(k)
l1 = locr(i)
l2 = l1 + lenr(i) - 1
do l = l1, l2
j = indr(l)
aij = abs( a(l) )
w(j) = max( w(j), aij )
Umax = max( Umax, aij )
end do
end do
parmlu(11) = Lmax
parmlu(12) = Umax
!--------------------------------------------------------------
! Find DUmax and DUmin, the extreme diagonals of U.
!--------------------------------------------------------------
do k = nslack + 1, nrank ! 12 Dec 2015: Allow for nslack.
j = q(k)
i = p(k)
l1 = locr(i)
diag = abs( a(l1) )
DUmax = max( DUmax, diag )
if (DUmin > diag) then
DUmin = diag
jumin = j
end if
end do
else
!--------------------------------------------------------------
! keepLU = 0.
! Only diag(U) is stored. Set w(*) accordingly.
! Find DUmax and DUmin, the extreme diagonals of U.
!--------------------------------------------------------------
ldiagU = lena - n
do k = nslack + 1, nrank ! 12 Dec 2015: Allow for nslack.
j = q(k)
! diag = abs( a(ldiagU + k) ) ! 06 May 2002: Diags
diag = abs( a(ldiagU + j) ) ! are in natural order
w(j) = diag
DUmax = max( DUmax, diag )
if (DUmin > diag) then
DUmin = diag
jumin = j
end if
end do
end if
!--------------------------------------------------------------
! Negate w(j) if the corresponding diagonal of U is
! too small in absolute terms or relative to the other elements
! in the same column of U.
!
! 23 Apr 2004: TRP ensures that diags are NOT small relative to
! other elements in their own column.
! Much better, we can compare all diags to DUmax.
! 13 Nov 2015: This causes slacks to replace slacks when DUmax
! is big. It seems better to leave Utol1 alone.
! 12 Dec 2015: Allow for nslack.
! DUmax now excludes slack rows, so we can
! reset Utol1 again for TRP.
!--------------------------------------------------------------
if (mode == 1 .and. TRP) then
Utol1 = max( Utol1, Utol2*DUmax )
end if
if (keepLU) then
do k = nslack + 1, n ! 12 Dec 2015: Allow for nslack.
j = q(k)
if (k > nrank) then
diag = zero
else
i = p(k)
l1 = locr(i)
diag = abs( a(l1) )
end if
if (diag <= Utol1 .or. diag <= Utol2*w(j)) then
nsing = nsing + 1
jsing = j
w(j) = - w(j)
end if
end do
else ! keepLU = 0
do k = nslack + 1, n ! 12 Dec 2015: Allow for nslack.
j = q(k)
diag = w(j)
if (diag <= Utol1) then
nsing = nsing + 1
jsing = j
w(j) = - w(j)
end if
end do
end if
!-----------------------------------------------------------------
! Set output parameters.
!-----------------------------------------------------------------
if (jumin == 0) DUmin = zero
luparm(11) = nsing
luparm(12) = jsing
luparm(19) = jumin
parmlu(13) = DUmax
parmlu(14) = DUmin
if (nsing > 0) then ! The matrix has been judged singular.
inform = 1
ndefic = n - nrank
if (nout > 0 .and. lprint >= 0) then
if (m > n) then
mnkey = '>'
else if (m == n) then
mnkey = '='
else
mnkey = '<'
end if
write(nout, 1100) mnkey, nrank, ndefic, nsing
end if
end if
! Exit.
luparm(10) = inform
return
1100 format(' Singular(m', a, 'n)', ' rank', i9, ' n-rank', i8, ' nsing', i9)
end subroutine lu6chk