lu1mar – AoC-2023

private subroutine lu1mar(m, n, lena, maxmn, TCP, aijtol, Ltol, maxcol, maxrow, ibest, jbest, mbest, a, indc, indr, p, q, lenc, lenr, locc, locr, iploc, iqloc)

! TPP ! Absolute test for Complete Pivoting

Arguments

Type	Intent		Name
integer(kind=ip),	intent(in)	::	m
integer(kind=ip),	intent(in)	::	n
integer(kind=ip),	intent(in)	::	lena
integer(kind=ip),	intent(in)	::	maxmn
logical,	intent(in)	::	TCP
real(kind=rp),	intent(in)	::	aijtol
real(kind=rp),	intent(in)	::	Ltol
integer(kind=ip),	intent(in)	::	maxcol
integer(kind=ip),	intent(in)	::	maxrow
integer(kind=ip),	intent(out)	::	ibest
integer(kind=ip),	intent(out)	::	jbest
integer(kind=ip),	intent(out)	::	mbest
real(kind=rp),	intent(in)	::	a(lena)
integer(kind=ip),	intent(in)	::	indc(lena)
integer(kind=ip),	intent(in)	::	indr(lena)
integer(kind=ip),	intent(in)	::	p(m)
integer(kind=ip),	intent(in)	::	q(n)
integer(kind=ip),	intent(in)	::	lenc(n)
integer(kind=ip),	intent(in)	::	lenr(m)
integer(kind=ip),	intent(in)	::	locc(n)
integer(kind=ip),	intent(in)	::	locr(m)
integer(kind=ip),	intent(in)	::	iploc(n)
integer(kind=ip),	intent(in)	::	iqloc(m)

Called by

Help

Source Code

  subroutine lu1mar( m    , n     , lena  , maxmn,          &
                     TCP  , aijtol, Ltol  , maxcol, maxrow, &
                     ibest, jbest , mbest ,                 &
                     a    , indc  , indr  , p     , q    ,  &
                     lenc , lenr  , locc  , locr  , iploc, iqloc )

    logical,       intent(in)    :: TCP
    integer(ip),   intent(in)    :: m, n, lena, maxmn, maxcol, maxrow
    real(rp),      intent(in)    :: aijtol, Ltol, a(lena)
    integer(ip),   intent(in)    :: indc(lena), indr(lena), p(m)    , q(n)    , &
                                    lenc(n)   , lenr(m)   , iploc(n), iqloc(m), &
                                    locc(n)   , locr(m)
    integer(ip),   intent(out)   :: ibest, jbest, mbest

    !------------------------------------------------------------------
    ! lu1mar  uses a Markowitz criterion to select a pivot element
    ! for the next stage of a sparse LU factorization,
    ! subject to a Threshold Partial Pivoting stability criterion (TPP)
    ! that bounds the elements of L.
    !
    ! 00 Jan 1986: Version documented in LUSOL paper:
    !              Gill, Murray, Saunders and Wright (1987),
    !              "Maintaining LU factors of a general sparse matrix",
    !              Linear algebra and its applications 88/89, 239-270.
    !
    ! 02 Feb 1989: Following Suhl and Aittoniemi (1987), the largest
    !              element in each column is now kept at the start of
    !              the column, i.e. in position locc(j) of a and indc.
    !              This should speed up the Markowitz searches.
    !
    ! 26 Apr 1989: Both columns and rows searched during spars1 phase.
    !              Only columns searched during spars2 phase.
    !              maxtie replaced by maxcol and maxrow.
    ! 05 Nov 1993: Initializing  "mbest = m * n"  wasn't big enough when
    !              m = 10, n = 3, and last column had 7 nonzeros.
    ! 09 Feb 1994: Realised that "mbest = maxmn * maxmn" might overflow.
    !              Changed to    "mbest = maxmn * 1000".
    ! 27 Apr 2000: On large example from Todd Munson,
    !              that allowed  "if (mbest .le. nz1**2) go to 900"
    !              to exit before any pivot had been found.
    !              Introduced kbest = mbest / nz1.
    !              Most pivots can be rejected with no integer(ip) multiply.
    !              True merit is evaluated only if it's as good as the
    !              best so far (or better).  There should be no danger
    !              of integer(ip) overflow unless A is incredibly
    !              large and dense.
    !
    ! 10 Sep 2000  TCP, aijtol added for Threshold Complete Pivoting.
    !
    ! 10 Jan 2010: First f90 version.
    ! 12 Dec 2011: Declare intent.
    !------------------------------------------------------------------

    integer(ip)            :: i, j, kbest, lc, lc1, lc2, len1,          &
                              lp, lp1, lp2, lq, lq1, lq2, lr, lr1, lr2, &
                              merit, ncol, nrow, nz, nz1
    real(rp)               :: abest, aij, amax, cmax, lbest
    real(rp),    parameter :: gamma  = 2.0

    ! gamma  is "gamma" in the tie-breaking rule TB4 in the LUSOL paper.

    !------------------------------------------------------------------
    ! Search cols of length nz = 1, then rows of length nz = 1,
    ! then   cols of length nz = 2, then rows of length nz = 2, etc.
    !------------------------------------------------------------------
    abest  = zero
    lbest  = zero
    ibest  = 0
    kbest  = maxmn + 1
    mbest  = -1
    ncol   = 0
    nrow   = 0
    nz1    = 0

NZS:do nz = 1, maxmn
       ! nz1    = nz - 1
       ! if (mbest .le. nz1**2) go to 900
       if (kbest <= nz1) exit NZS
       if (ibest >  0  ) then
          if (ncol >= maxcol) go to 200
       end if
       if (nz    >  m  ) go to 200

       !---------------------------------------------------------------
       ! Search the set of columns of length  nz.
       !---------------------------------------------------------------
       lq1    = iqloc(nz)
       lq2    = n
       if (nz < m) lq2 = iqloc(nz + 1) - 1

Cols:  do lq = lq1, lq2
          ncol   = ncol + 1
          j      = q(lq)
          lc1    = locc(j)
          lc2    = lc1 + nz1
          amax   = abs( a(lc1) )

          ! Test all aijs in this column.
          ! amax is the largest element (the first in the column).
          ! cmax is the largest multiplier if aij becomes pivot.

          if ( TCP ) then
             if (amax < aijtol) cycle Cols ! Nothing in whole column
          end if

Colj:     do lc = lc1, lc2
             i      = indc(lc)
             len1   = lenr(i) - 1
             ! merit  = nz1 * len1
           ! if (merit > mbest) cycle
             if (len1  > kbest) cycle Colj

             ! aij  has a promising merit.
             ! Apply the stability test.
             ! We require  aij  to be sufficiently large compared to
             ! all other nonzeros in column  j.  This is equivalent
             ! to requiring cmax to be bounded by Ltol.

             if (lc == lc1) then

                ! This is the maximum element, amax.
                ! Find the biggest element in the rest of the column
                ! and hence get cmax.  We know cmax .le. 1, but
                ! we still want it exactly in order to break ties.
                ! 27 Apr 2002: Settle for cmax = 1.

                aij    = amax
                cmax   = one

                ! cmax   = zero
                ! do 140 l = lc1 + 1, lc2
                ! cmax  = max( cmax, abs( a(l) ) )
                ! 140            continue
                ! cmax   = cmax / amax
             else

                ! aij is not the biggest element, so cmax >= 1.
                ! Bail out if cmax will be too big.

                aij    = abs( a(lc) )
                if ( TCP ) then ! Absolute test for Complete Pivoting
                   if (aij      < aijtol) cycle Colj
                else !!! TPP
                   if (aij*Ltol < amax  ) cycle Colj
                end if
                cmax   = amax / aij
             end if

             ! aij  is big enough.  Its maximum multiplier is cmax.

             merit  = nz1 * len1
             if (merit == mbest) then

                ! Break ties.
                ! (Initializing mbest < 0 prevents getting here if
                ! nothing has been found yet.)
                ! In this version we minimize cmax
                ! but if it is already small we maximize the pivot.

                if (lbest <= gamma  .and.  cmax <= gamma) then
                   if (abest >= aij ) cycle Colj
                else
                   if (lbest <= cmax) cycle Colj
                end if
             end if

             ! aij  is the best pivot so far.

             ibest  = i
             jbest  = j
             kbest  = len1
             mbest  = merit
             abest  = aij
             lbest  = cmax
             if (nz == 1) exit NZS
          end do Colj

          ! Finished with that column.

          if (ibest > 0) then
             if (ncol >= maxcol) exit Cols
          end if
       end do Cols

       !---------------------------------------------------------------
       ! Search the set of rows of length  nz.
       !---------------------------------------------------------------
! 200  if (mbest .le. nz*nz1) go to 900
200    if (kbest <= nz    ) exit NZS
       if (ibest > 0) then
          if (nrow >= maxrow) go to 290
       end if
       if (nz > n) go to 290

       lp1    = iploc(nz)
       lp2    = m
       if (nz < n) lp2 = iploc(nz + 1) - 1

Rows:  do lp = lp1, lp2
          nrow   = nrow + 1
          i      = p(lp)
          lr1    = locr(i)
          lr2    = lr1 + nz1

Rowi:     do lr = lr1, lr2
             j      = indr(lr)
             len1   = lenc(j) - 1
             ! merit  = nz1 * len1
           ! if (merit > mbest) cycle
             if (len1  > kbest) cycle Rowi

             ! aij  has a promising merit.
             ! Find where  aij  is in column  j.

             lc1    = locc(j)
             lc2    = lc1 + len1
             amax   = abs( a(lc1) )
             do lc = lc1, lc2
                if (indc(lc) == i) exit
             end do

             ! Apply the same stability test as above.

             aij    = abs( a(lc) )
             if ( TCP ) then   !!! Absolute test for Complete Pivoting
                if (aij < aijtol) cycle Rowi
             end if

             if (lc == lc1) then

                ! This is the maximum element, amax.
                ! Find the biggest element in the rest of the column
                ! and hence get cmax.  We know cmax .le. 1, but
                ! we still want it exactly in order to break ties.
                ! 27 Apr 2002: Settle for cmax = 1.

                cmax   = one

                ! cmax   = zero
                !     do 240 l = lc1 + 1, lc2
                !        cmax  = max( cmax, abs( a(l) ) )
                ! 240 continue
                ! cmax   = cmax / amax
             else

                ! aij is not the biggest element, so cmax .ge. 1.
                ! Bail out if cmax will be too big.

                if ( TCP ) then
                   ! relax
                else
                   if (aij*Ltol < amax) cycle Rowi
                end if
                cmax   = amax / aij
             end if

             ! aij  is big enough.  Its maximum multiplier is cmax.

             merit  = nz1 * len1
             if (merit == mbest) then

                ! Break ties as before.
                ! (Initializing mbest < 0 prevents getting here if
                ! nothing has been found yet.)

                if (lbest <= gamma  .and.  cmax <= gamma) then
                   if (abest >= aij ) cycle Rowi
                else
                   if (lbest <= cmax) cycle Rowi
                end if
             end if

             ! aij  is the best pivot so far.

             ibest  = i
             jbest  = j
             kbest  = len1
             mbest  = merit
             abest  = aij
             lbest  = cmax
             if (nz == 1) exit NZS
          end do Rowi

          ! Finished with that row.

          if (ibest > 0) then
             if (nrow >= maxrow) exit Rows
          end if
       end do Rows

       ! See if it's time to quit.

290    if (ibest > 0) then
          if (nrow >= maxrow  .and.  ncol >= maxcol) exit NZS
       end if

       ! Press on with next nz.

       nz1    = nz
       if (ibest > 0) kbest = mbest / nz1
    end do NZS

  end subroutine lu1mar

lu1mar Subroutine

Contents

Source Code

private subroutine lu1mar(m, n, lena, maxmn, TCP, aijtol, Ltol, maxcol, maxrow, ibest, jbest, mbest, a, indc, indr, p, q, lenc, lenr, locc, locr, iploc, iqloc)

Arguments

Called by

Source Code