subroutine lu6LD ( inform, mode, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc, locr )
integer(ip), intent(in) :: mode, m, n, lena
integer(ip), intent(in) :: indc(lena), indr(lena), lenc(n), locr(m)
real(rp), intent(in) :: a(lena)
integer(ip), intent(inout) :: luparm(30)
real(rp), intent(inout) :: parmlu(30), v(m)
integer(ip), intent(out) :: inform
!-------------------------------------------------------------------
! lu6LD assumes lu1fac has computed factors A = LU of a
! symmetric definite or quasi-definite matrix A,
! using Threshold Symmetric Pivoting (TSP), luparm(6) = 3,
! or Threshold Diagonal Pivoting (TDP), luparm(6) = 4.
! It also assumes that no updates have been performed.
! In such cases, U = D L', where D = diag(U).
! lu6LDL returns v as follows:
!
! mode
! 1 v solves L D v = v(input).
! 2 v solves L|D|v = v(input).
!
! 15 Dec 2002: First version of lu6LD.
! 15 Dec 2002: Current version.
! 13 Dec 2011: First f90 version.
!-----------------------------------------------------------------------
!
! Solve L D v(new) = v or L|D|v(new) = v, depending on mode.
! The code for L is the same as in lu6L,
! but when a nonzero entry of v arises, we divide by
! the corresponding entry of D or |D|.
integer(ip) :: ipiv, j, k, l, l1, ldummy, len, numL0
real(rp) :: diag, small, vpiv
numL0 = luparm(20)
small = parmlu(3)
inform = 0
l1 = lena + 1
do k = 1, numL0
len = lenc(k)
l = l1
l1 = l1 - len
ipiv = indr(l1)
vpiv = v(ipiv)
if (abs(vpiv) > small) then
!***** This loop could be coded specially.
do ldummy = 1, len
l = l - 1
j = indc(l)
v(j) = v(j) + a(l)*vpiv
end do
! Find diag = U(ipiv,ipiv) and divide by diag or |diag|.
l = locr(ipiv)
diag = a(l)
if (mode == 2) diag = abs(diag)
v(ipiv) = vpiv/diag
end if
end do
end subroutine lu6LD