lu7elm Subroutine

private subroutine lu7elm(m, n, jelm, v, lena, luparm, parmlu, lenL, lenU, lrow, nrank, a, indc, indr, p, q, lenr, locc, locr, inform, diag)

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: m
integer(kind=ip), intent(in) :: n
integer(kind=ip), intent(in) :: jelm
real(kind=rp), intent(in) :: v(m)
integer(kind=ip), intent(in) :: lena
integer(kind=ip), intent(inout) :: luparm(30)
real(kind=rp), intent(inout) :: parmlu(30)
integer(kind=ip), intent(inout) :: lenL
integer(kind=ip), intent(in) :: lenU
integer(kind=ip), intent(inout) :: lrow
integer(kind=ip), intent(in) :: nrank
real(kind=rp), intent(inout) :: a(lena)
integer(kind=ip), intent(inout) :: indc(lena)
integer(kind=ip), intent(inout) :: indr(lena)
integer(kind=ip), intent(inout) :: p(m)
integer(kind=ip), intent(in) :: q(n)
integer(kind=ip), intent(inout) :: lenr(m)
integer(kind=ip), intent(inout) :: locc(n)
integer(kind=ip), intent(inout) :: locr(m)
integer(kind=ip), intent(out) :: inform
real(kind=rp), intent(out) :: diag

Calls

proc~~lu7elm~~CallsGraph proc~lu7elm lusol::lu7elm proc~lu1rec lusol::lu1rec proc~lu7elm->proc~lu1rec

Called by

proc~~lu7elm~~CalledByGraph proc~lu7elm lusol::lu7elm proc~lu8rpc lusol::lu8rpc proc~lu8rpc->proc~lu7elm

Source Code

  subroutine lu7elm( m, n, jelm, v, lena, luparm, parmlu,     &
                     lenL, lenU, lrow, nrank,                 &
                     a, indc, indr, p, q, lenr, locc, locr,   &
                     inform, diag )

    integer(ip),   intent(in)    :: m, n, jelm, lena, nrank
    integer(ip),   intent(in)    :: lenU, q(n)   ! not used
    real(rp),      intent(in)    :: v(m)
    integer(ip),   intent(inout) :: luparm(30), lenL, lrow,       &
                                    indc(lena), indr(lena), p(m), &
                                    lenr(m), locc(n), locr(m)
    real(rp),      intent(inout) :: parmlu(30), a(lena)

    integer(ip),   intent(out)   :: inform
    real(rp),      intent(out)   :: diag

    !------------------------------------------------------------------
    ! lu7elm  eliminates the subdiagonal elements of a vector  v(*),
    ! where  L*v = y  for some vector y.
    ! If  jelm > 0,  y  has just become column  jelm  of the matrix  A.
    ! lu7elm  should not be called unless  m  is greater than  nrank.
    !
    ! inform = 0 if y contained no subdiagonal nonzeros to eliminate.
    ! inform = 1 if y contained at least one nontrivial subdiagonal.
    ! inform = 7 if there is insufficient storage.
    !
    ! 09 May 1988: First f77 version.
    !              No longer calls lu7for at end.  lu8rpc, lu8mod do so.
    ! 13 Dec 2011: First f90 version.
    ! 20 Dec 2015: ilast is now output by lu1rec.
    !------------------------------------------------------------------

    integer(ip)            :: i, ilast, imax, k, kmax, l, l1, l2, lmax, &
                              minfre, nfree, nrank1
    real(rp)               :: small, vi, vmax

    small  = parmlu(3)
    nrank1 = nrank + 1
    diag   = zero

    ! Compress row file if necessary.

    minfre = m - nrank
    nfree  = lena - lenL - lrow
    if (nfree < minfre) then
       call lu1rec( m, .true., luparm, lrow, ilast, &
                    lena, a, indr, lenr, locr )
       nfree  = lena - lenL - lrow
       if (nfree < minfre) go to 970
    end if

    ! Pack the subdiagonals of  v  into  L,  and find the largest.

    vmax   = zero
    kmax   = 0
    l      = lena - lenL + 1

    do k = nrank1, m
       i       = p(k)
       vi      = abs(v(i))
       if (vi <= small) cycle
       l       = l - 1
       a(l)    = v(i)
       indc(l) = i
       if (vmax >= vi ) cycle
       vmax    = vi
       kmax    = k
       lmax    = l
    end do

    if (kmax == 0) go to 900

    !------------------------------------------------------------------
    ! Remove  vmax  by overwriting it with the last packed  v(i).
    ! Then set the multipliers in  L  for the other elements.
    !------------------------------------------------------------------
    imax       = p(kmax)
    vmax       = a(lmax)
    a(lmax)    = a(l)
    indc(lmax) = indc(l)
    l1         = l + 1
    l2         = lena - lenL
    lenL       = lenL + (l2 - l)

    do l = l1, l2
       a(l)    = - a(l) / vmax
       indr(l) =   imax
    end do

    ! Move the row containing vmax to pivotal position nrank + 1.

    p(kmax  ) = p(nrank1)
    p(nrank1) = imax
    diag      = vmax

    !------------------------------------------------------------------
    ! If jelm is positive, insert  vmax  into a new row of  U.
    ! This is now the only subdiagonal element.
    !------------------------------------------------------------------

    if (jelm > 0) then
       lrow       = lrow + 1
       locr(imax) = lrow
       lenr(imax) = 1
       a(lrow)    = vmax
       indr(lrow) = jelm
    end if

    inform = 1
    go to 990

    ! No elements to eliminate.

900 inform = 0
    go to 990

    ! Not enough storage.

970 inform = 7

990 return

  end subroutine lu7elm