subroutine lu1or1( m, n, nelem, lena, small, &
a, indc, indr, lenc, lenr, &
Amax, numnz, lerr, inform )
integer(ip), intent(in) :: m, n, nelem, lena
real(rp), intent(in) :: small
real(rp), intent(inout) :: a(lena)
integer(ip), intent(inout) :: indc(lena), indr(lena)
integer(ip), intent(out) :: lerr, inform
integer(ip), intent(out) :: lenc(n), lenr(m)
!------------------------------------------------------------------
! lu1or1 organizes the elements of an m by n matrix A as
! follows. On entry, the parallel arrays a, indc, indr,
! contain nelem entries of the form aij, i, j,
! in any order. nelem must be positive.
!
! Entries not larger than the input parameter small are treated as
! zero and removed from a, indc, indr. The remaining entries are
! defined to be nonzero. numnz returns the number of such nonzeros
! and Amax returns the magnitude of the largest nonzero.
! The arrays lenc, lenr return the number of nonzeros in each
! column and row of A.
!
! inform = 0 on exit, except inform = 1 if any of the indices in
! indc, indr imply that the element aij lies outside the m by n
! dimensions of A.
!
! xx Feb 1985: Original version.
! 17 Oct 2000: a, indc, indr now have size lena to allow nelem = 0.
!
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent and local variables.
!------------------------------------------------------------------
integer(ip) :: i, j, l, numnz
real(rp) :: Amax
lenr(1:m) = 0
lenc(1:n) = 0
lerr = 0
Amax = zero
numnz = nelem
do l = nelem, 1, -1
if (abs(a(l)) > small) then
i = indc(l)
j = indr(l)
Amax = max( Amax, abs(a(l)) )
if (i < 1 .or. i > m) go to 910
if (j < 1 .or. j > n) go to 910
lenr(i) = lenr(i) + 1
lenc(j) = lenc(j) + 1
else
! Replace a negligible element by last element. Since
! we are going backwards, we know the last element is ok.
a(l) = a(numnz)
indc(l) = indc(numnz)
indr(l) = indr(numnz)
numnz = numnz - 1
end if
end do
inform = 0
return
910 lerr = l
inform = 1
return
end subroutine lu1or1