subroutine lu1ful( m , n , lena , lenD , lu1 , TPP, &
mleft , nleft, nrank, nrowu, &
lenL , lenU , nsing, &
keepLU, small, &
a , d , indc , indr , p , q, &
lenc , lenr , locc , ipinv, ipvt )
logical, intent(in) :: TPP, keepLU
integer(ip), intent(in) :: m, n, lena, lenD, lu1, &
mleft, nleft, nrank, nrowu
integer(ip), intent(in) :: locc(n)
real(rp), intent(in) :: small
integer(ip), intent(inout) :: lenL, lenU
integer(ip), intent(inout) :: indc(lena), indr(lena), p(m), q(n), &
lenc(n) , lenr(m)
real(rp), intent(inout) :: a(lena)
integer(ip), intent(out) :: ipvt(m), ipinv(m) ! workspace
integer(ip), intent(out) :: nsing ! not used outside
real(rp), intent(out) :: d(lenD)
!------------------------------------------------------------------
! lu1ful computes a dense (full) LU factorization of the
! mleft by nleft matrix that remains to be factored at the
! beginning of the nrowu-th pass through the main loop of lu1fad.
!
! 02 May 1989: First version.
! 05 Feb 1994: Column interchanges added to lu1DPP.
! 08 Feb 1994: ipinv reconstructed, since lu1pq3 may alter p.
!
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent and local variables.
!------------------------------------------------------------------
integer(ip) :: i, ibest, ipbase, j, jbest, k, l, l1, l2, &
la, lc, lc1, lc2, ld, ldbase, ldiagU, &
lkk, lkn, ll, lq, lu, nrowd, ncold
real(rp) :: ai, aj
!------------------------------------------------------------------
! If lu1pq3 moved any empty rows, reset ipinv = inverse of p.
!------------------------------------------------------------------
if (nrank < m) then
do l = 1, m
i = p(l)
ipinv(i) = l
end do
end if
!------------------------------------------------------------------
! Copy the remaining matrix into the dense matrix D.
!------------------------------------------------------------------
d(1:lenD) = zero
ipbase = nrowu - 1
ldbase = 1 - nrowu
do lq = nrowu, n
j = q(lq)
lc1 = locc(j)
lc2 = lc1 + lenc(j) - 1
do lc = lc1, lc2
i = indc(lc)
ld = ldbase + ipinv(i)
d(ld) = a(lc)
end do
ldbase = ldbase + mleft
end do
!------------------------------------------------------------------
! Call our favorite dense LU factorizer.
!------------------------------------------------------------------
if ( TPP ) then
call lu1DPP( d, mleft, mleft, nleft, small, nsing, ipvt, q(nrowu) )
else
call lu1DCP( d, mleft, mleft, nleft, small, nsing, ipvt, q(nrowu) )
end if
!------------------------------------------------------------------
! Move D to the beginning of A,
! and pack L and U at the top of a, indc, indr.
! In the process, apply the row permutation to p.
! lkk points to the diagonal of U.
!------------------------------------------------------------------
a(1:lenD) = d(1:lenD)
ldiagU = lena - n
lkk = 1
lkn = lenD - mleft + 1
lu = lu1
do k = 1, min( mleft, nleft )
l1 = ipbase + k
l2 = ipbase + ipvt(k)
if (l1 /= l2) then
i = p(l1)
p(l1) = p(l2)
p(l2) = i
end if
ibest = p(l1)
jbest = q(l1)
if ( keepLU ) then
!===========================================================
! Pack the next column of L.
!===========================================================
la = lkk
ll = lu
nrowd = 1
do i = k + 1, mleft
la = la + 1
ai = a(la)
if (abs( ai ) > small) then
nrowd = nrowd + 1
ll = ll - 1
a(ll) = ai
indc(ll) = p( ipbase + i )
indr(ll) = ibest
end if
end do
!===========================================================
! Pack the next row of U.
! We go backwards through the row of D
! so the diagonal ends up at the front of the row of U.
! Beware -- the diagonal may be zero.
!===========================================================
la = lkn + mleft
lu = ll
ncold = 0
do j = nleft, k, -1
la = la - mleft
aj = a(la)
if (abs( aj ) > small .or. j == k) then
ncold = ncold + 1
lu = lu - 1
a(lu) = aj
indr(lu) = q(ipbase + j)
end if
end do
lenr(ibest) = - ncold
lenc(jbest) = - nrowd
lenL = lenL + nrowd - 1
lenU = lenU + ncold
lkn = lkn + 1
else
!===========================================================
! Store just the diagonal of U, in natural order.
!===========================================================
a(ldiagU + jbest) = a(lkk)
end if
lkk = lkk + mleft + 1
end do
end subroutine lu1ful