lu1fac Subroutine

  subroutine lu1fac( m    , n    , nelem, lena , luparm, parmlu,       &
                     a    , indc , indr , p    , q     ,               &
                     lenc , lenr , locc , locr ,                       &
                     iploc, iqloc, ipinv, iqinv, w     , inform )

    integer(ip),   intent(in)    :: m, n, nelem, lena

    integer(ip),   intent(inout) :: luparm(30)
    integer(ip),   intent(inout) :: indc(lena), indr(lena),            &
                                    p(m)      , q(n)      ,            &
                                    lenc(n)   , lenr(m)   ,            &
                                    iploc(n)  , iqloc(m)  ,            &
                                    ipinv(m)  , iqinv(n)  ,            &
                                    locc(n)   , locr(m)
    real(rp),      intent(inout) :: parmlu(30), a(lena), w(n)

    integer(ip),   intent(out)   :: inform

    !------------------------------------------------------------------
    ! lu1fac computes a factorization A = L*U, where A is a sparse
    ! matrix with m rows and n columns, P*L*P' is lower triangular
    ! and P*U*Q is upper triangular for certain permutations P, Q
    ! (which are returned in the arrays p, q).
    ! Stability is ensured by limiting the size of the elements of L.
    !
    ! The nonzeros of A are input via the parallel arrays a, indc, indr,
    ! which should contain nelem entries of the form    aij,    i,    j
    ! in any order.  There should be no duplicate pairs         i,    j.
    !
    ! ******************************************************************
    ! *        Beware !!!   The row indices i must be in indc,         *
    ! *              and the column indices j must be in indr.         *
    ! *              (Not the other way round!)                        *
    ! ******************************************************************
    !
    ! It does not matter if some of the entries in a(*) are zero.
    ! Entries satisfying  abs( a(i) ) .le. parmlu(3)  are ignored.
    ! Other parameters in luparm and parmlu are described below.
    !
    ! The matrix A may be singular.  On exit, nsing = luparm(11) gives
    ! the number of apparent singularities.  This is the number of
    ! "small" diagonals of the permuted factor U, as judged by
    ! the input tolerances Utol1 = parmlu(4) and  Utol2 = parmlu(5).
    ! The diagonal element diagj associated with column j of A is
    ! "small" if
    !                 abs( diagj ) .le. Utol1
    ! or
    !                 abs( diagj ) .le. Utol2 * max( uj ),
    !
    ! where max( uj ) is the maximum element in the j-th column of U.
    ! The position of such elements is returned in w(*).  In general,
    ! w(j) = + max( uj ),  but if column j is a singularity,
    ! w(j) = - max( uj ).  Thus, w(j) .le. 0 if column j appears to be
    ! dependent on the other columns of A.
    !
    ! NOTE: lu1fac (like certain other sparse LU packages) does not
    ! treat dense columns efficiently.  This means it will be slow
    ! on "arrow matrices" of the form
    !                  A = (x       a)
    !                      (  x     b)
    !                      (    x   c)
    !                      (      x d)
    !                      (x x x x e)
    ! if the numerical values in the dense column allow it to be
    ! chosen LATE in the pivot order.
    !
    ! With TPP (Threshold Partial Pivoting), the dense column is
    ! likely to be chosen late.
    !
    ! With TCP (Threshold Complete Pivoting), if any of a,b,c,d
    ! is significantly larger than other elements of A, it will
    ! be chosen as the first pivot and the dense column will be
    ! eliminated, giving reasonably sparse factors.
    ! However, if element e is so big that TCP chooses it, the factors
    ! will become dense.  (It's hard to win on these examples!)
    !==================================================================
    !
    !
    ! Notes on the array names
    ! ------------------------
    !
    ! During the LU factorization, the sparsity pattern of the matrix
    ! being factored is stored twice: in a column list and a row list.
    !
    ! The column list is ( a, indc, locc, lenc )
    ! where
    !       a(*)    holds the nonzeros,
    !       indc(*) holds the indices for the column list,
    !       locc(j) points to the start of column j in a(*) and indc(*),
    !       lenc(j) is the number of nonzeros in column j.
    !
    ! The row list is    (    indr, locr, lenr )
    ! where
    !       indr(*) holds the indices for the row list,
    !       locr(i) points to the start of row i in indr(*),
    !       lenr(i) is the number of nonzeros in row i.
    !
    !
    ! At all stages of the LU factorization, p contains a complete
    ! row permutation.  At the start of stage k,  p(1), ..., p(k-1)
    ! are the first k-1 rows of the final row permutation P.
    ! The remaining rows are stored in an ordered list
    !                    ( p, iploc, ipinv )
    ! where
    !       iploc(nz) points to the start in p(*) of the set of rows
    !                 that currently contain nz nonzeros,
    !       ipinv(i)  points to the position of row i in p(*).
    !
    ! For example,
    !       iploc(1) = k   (and this is where rows of length 1 begin),
    !       iploc(2) = k+p  if there are p rows of length 1
    !                      (and this is where rows of length 2 begin).
    !
    ! Similarly for q, iqloc, iqinv.
    !==================================================================
    !
    !
    ! 00 Jun 1983  Original version.
    ! 00 Jul 1987  nrank  saved in luparm(16).
    ! 12 Apr 1989  ipinv, iqinv added as workspace.
    ! 26 Apr 1989  maxtie replaced by maxcol in Markowitz search.
    ! 16 Mar 1992  jumin  saved in luparm(19).
    ! 10 Jun 1992  lu1fad has to move empty rows and cols to the bottom
    !              (via lu1pq3) before doing the dense LU.
    ! 12 Jun 1992  Deleted dense LU (lu1ful, lu1vlu).
    ! 25 Oct 1993  keepLU implemented.
    ! 07 Feb 1994  Added new dense LU (lu1ful, lu1den).
    ! 21 Dec 1994  Bugs fixed in lu1fad (nrank) and lu1ful (ipvt).
    ! 08 Aug 1995  Use p instead of w as parameter to lu1or3 (for F90).
    ! 13 Sep 2000  TPP and TCP options implemented.
    ! 17 Oct 2000  Fixed troubles due to A = empty matrix (Todd Munson).
    ! 01 Dec 2000  Save Lmax, Umax, etc. after both lu1fad and lu6chk.
    !              lu1fad sets them when keepLU = false.
    !              lu6chk sets them otherwise, and includes items
    !              from the dense LU.
    ! 11 Mar 2001  lu6chk now looks at diag(U) when keepLU = false.
    ! 26 Apr 2002  New TCP implementation using heap routines to
    !              store largest element in each column.
    !              New workspace arrays Ha, Hj, Hk required.
    !              For compatibility, borrow space from a, indc, indr
    !              rather than adding new input parameters.
    ! 01 May 2002  lu1den changed to lu1DPP (dense partial  pivoting).
    !              lu1DCP implemented       (dense complete pivoting).
    !              Both TPP and TCP now switch to dense mode and end.
    !
    ! 10 Jan 2010: First f90 version.
    !---------------------------------------------------------------------
    !
    !
    !  INPUT PARAMETERS
    !
    !  m      (not altered) is the number of rows in A.
    !  n      (not altered) is the number of columns in A.
    !  nelem  (not altered) is the number of matrix entries given in
    !         the arrays a, indc, indr.
    !  lena   (not altered) is the dimension of  a, indc, indr.
    !         This should be significantly larger than nelem.
    !         Typically one should have
    !            lena > max( 2*nelem, 10*m, 10*n, 10000 )
    !         but some applications may need more.
    !         On machines with virtual memory it is safe to have
    !         lena "far bigger than necessary", since not all of the
    !         arrays will be used.
    !  a      (overwritten) contains entries   Aij  in   a(1:nelem).
    !  indc   (overwritten) contains the indices i in indc(1:nelem).
    !  indr   (overwritten) contains the indices j in indr(1:nelem).
    !
    !  luparm input parameters:                                Typical value
    !
    !  luparm( 1) = nout     File number for printed messages.         6
    !
    !  luparm( 2) = lprint   Print level.                              0
    !                   <  0 suppresses output.
    !                   =  0 gives error messages.
    !                  >= 10 gives statistics about the LU factors.
    !                  >= 50 gives debug output from lu1fac
    !                        (the pivot row and column and the
    !                        no. of rows and columns involved at
    !                        each elimination step).
    !
    !  luparm( 3) = maxcol   lu1fac: maximum number of columns         5
    !                        searched allowed in a Markowitz-type
    !                        search for the next pivot element.
    !                        For some of the factorization, the
    !                        number of rows searched is
    !                        maxrow = maxcol - 1.
    !
    !  luparm( 6) = 0    =>  TPP: Threshold Partial   Pivoting.        0
    !             = 1    =>  TRP: Threshold Rook      Pivoting.
    !             = 2    =>  TCP: Threshold Complete  Pivoting.
    !             = 3    =>  TSP: Threshold Symmetric Pivoting.
    !             = 4    =>  TDP: Threshold Diagonal  Pivoting.
    !                             (TDP not yet implemented).
    !                        TRP and TCP are more expensive than TPP but
    !                        more stable and better at revealing rank.
    !                        Take care with setting parmlu(1), especially
    !                        with TCP.
    !                        NOTE: TSP and TDP are for symmetric matrices
    !                        that are either definite or quasi-definite.
    !                        TSP is effectively TRP for symmetric matrices.
    !                        TDP is effectively TCP for symmetric matrices.
    !
    !  luparm( 8) = keepLU   lu1fac: keepLU = 1 means the numerical    1
    !                        factors will be computed if possible.
    !                        keepLU = 0 means L and U will be discarded
    !                        but other information such as the row and
    !                        column permutations will be returned.
    !                        The latter option requires less storage.
    !
    !  parmlu input parameters:                                Typical value
    !
    !  parmlu( 1) = Ltol1    Max Lij allowed during Factor.
    !                                                  TPP     10.0 or 100.0
    !                                                  TRP      4.0 or  10.0
    !                                                  TCP      5.0 or  10.0
    !                                                  TSP      4.0 or  10.0
    !                        With TRP and TCP (Rook and Complete Pivoting),
    !                        values less than 25.0 may be expensive
    !                        on badly scaled data.  However,
    !                        values less than 10.0 may be needed
    !                        to obtain a reliable rank-revealing
    !                        factorization.
    !  parmlu( 2) = Ltol2    Max Lij allowed during Updates.            10.0
    !                        during updates.
    !  parmlu( 3) = small    Absolute tolerance for       eps**0.8 = 3.0d-13
    !                        treating reals as zero.
    !  parmlu( 4) = Utol1    Absolute tol for flagging    eps**0.67= 3.7d-11
    !                        small diagonals of U.
    !  parmlu( 5) = Utol2    Relative tol for flagging    eps**0.67= 3.7d-11
    !                        small diagonals of U.
    !                        (eps = machine precision)
    !  parmlu( 6) = Uspace   Factor limiting waste space in  U.      3.0
    !                        In lu1fac, the row or column lists
    !                        are compressed if their length
    !                        exceeds Uspace times the length of
    !                        either file after the last compression.
    !  parmlu( 7) = dens1    The density at which the Markowitz      0.3
    !                        pivot strategy should search maxcol
    !                        columns and no rows.
    !                        (Use 0.3 unless you are experimenting
    !                        with the pivot strategy.)
    !  parmlu( 8) = dens2    the density at which the Markowitz      0.5
    !                        strategy should search only 1 column,
    !                        or (if storage is available)
    !                        the density at which all remaining
    !                        rows and columns will be processed
    !                        by a dense LU code.
    !                        For example, if dens2 = 0.1 and lena is
    !                        large enough, a dense LU will be used
    !                        once more than 10 per cent of the
    !                        remaining matrix is nonzero.
    !
    !
    !  OUTPUT PARAMETERS
    !
    !  a, indc, indr     contain the nonzero entries in the LU factors of A.
    !         If keepLU = 1, they are in a form suitable for use
    !         by other parts of the LUSOL package, such as lu6sol.
    !         U is stored by rows at the start of a, indr.
    !         L is stored by cols at the end   of a, indc.
    !         If keepLU = 0, only the diagonals of U are stored, at the
    !         end of a.
    !  p, q   are the row and column permutations defining the
    !         pivot order.  For example, row p(1) and column q(1)
    !         defines the first diagonal of U.
    !  lenc(1:numl0) contains the number of entries in nontrivial
    !         columns of L (in pivot order).
    !  lenr(1:m) contains the number of entries in each row of U
    !         (in original order).
    !  locc(1:n) = 0 (ready for the LU update routines).
    !  locr(1:m) points to the beginning of the rows of U in a, indr.
    !  iploc, iqloc, ipinv, iqinv  are undefined.
    !  w      indicates singularity as described above.
    !  inform = 0 if the LU factors were obtained successfully.
    !         = 1 if U appears to be singular, as judged by lu6chk.
    !         = 3 if some index pair indc(l), indr(l) lies outside
    !             the matrix dimensions 1:m , 1:n.
    !         = 4 if some index pair indc(l), indr(l) duplicates
    !             another such pair.
    !         = 7 if the arrays a, indc, indr were not large enough.
    !             Their length "lena" should be increase to at least
    !             the value "minlen" given in luparm(13).
    !         = 8 if there was some other fatal error.  (Shouldn't happen!)
    !         = 9 if no diagonal pivot could be found with TSP or TDP.
    !             The matrix must not be sufficiently definite
    !             or quasi-definite.
    !         =10 if there was some other fatal error.
    !
    !  luparm output parameters:
    !
    !  luparm(10) = inform   Return code from last call to any LU routine.
    !  luparm(11) = nsing    No. of singularities marked in the
    !                        output array w(*).
    !  luparm(12) = jsing    Column index of last singularity.
    !  luparm(13) = minlen   Minimum recommended value for  lena.
    !  luparm(14) = maxlen   ?
    !  luparm(15) = nupdat   No. of updates performed by the lu8 routines.
    !  luparm(16) = nrank    No. of nonempty rows of U.
    !  luparm(17) = ndens1   No. of columns remaining when the density of
    !                        the matrix being factorized reached dens1.
    !  luparm(18) = ndens2   No. of columns remaining when the density of
    !                        the matrix being factorized reached dens2.
    !  luparm(19) = jumin    The column index associated with DUmin.
    !  luparm(20) = numL0    No. of columns in initial  L.
    !  luparm(21) = lenL0    Size of initial  L  (no. of nonzeros).
    !  luparm(22) = lenU0    Size of initial  U.
    !  luparm(23) = lenL     Size of current  L.
    !  luparm(24) = lenU     Size of current  U.
    !  luparm(25) = lrow     Length of row file.
    !  luparm(26) = ncp      No. of compressions of LU data structures.
    !  luparm(27) = mersum   lu1fac: sum of Markowitz merit counts.
    !  luparm(28) = nUtri    lu1fac: triangular rows in U.
    !  luparm(29) = nLtri    lu1fac: triangular rows in L.
    !  luparm(30) = nslack   lu1fac: no. of unit vectors at start of U. (info only)
    !
    !
    !
    !  parmlu output parameters:
    !
    !  parmlu(10) = Amax     Maximum element in  A.
    !  parmlu(11) = Lmax     Maximum multiplier in current  L.
    !  parmlu(12) = Umax     Maximum element in current  U.
    !  parmlu(13) = DUmax    Maximum diagonal in  U.
    !  parmlu(14) = DUmin    Minimum diagonal in  U.
    !  parmlu(15) = Akmax    Maximum element generated at any stage
    !                        during TCP factorization.
    !  parmlu(16) = growth   TPP: Umax/Amax    TRP, TCP, TSP: Akmax/Amax
    !  parmlu(17) =
    !  parmlu(18) =
    !  parmlu(19) =
    !  parmlu(20) = resid    lu6sol: residual after solve with U or U'.
    !  ...
    !  parmlu(30) =
    !---------------------------------------------------------------------

    character(1)           :: mnkey
    character(2)           :: kPiv(0:3)
    integer(ip)            :: i, idummy, j, jsing, jumin,              &
                              k, l, l2, lena2, lenH, lenL,             &
                              lenLk, lenU, lenUk, lerr,                &
                              ll, llsave, lm, lmaxr, locH,             &
                              lprint, lPiv, lrow, ltopl,               &
                              lu, mersum, minlen, nbump,               &
                              ncp, ndens1, ndens2,                     &
                              nLtri, nmove, nout, nrank,               &
                              nsing, numl0, numnz, nslack, nUtri
    logical                :: keepLU, TCP, TPP, TRP, TSP
    real(rp)               :: Agrwth, Akmax, Amax, avgmer,             &
                              condU, delem, densty, dincr,             &
                              dm, dn, DUmax, DUmin, growth,            &
                              Lmax, Ltol, small, Ugrwth,               &
                              Umax

    ! Grab relevant input parameters.

    nout   = luparm(1)
    lprint = luparm(2)
    lPiv   = luparm(6)
    keepLU = luparm(8) /= 0

    Ltol   = parmlu(1)  ! Limit on size of Lij
    small  = parmlu(3)  ! Drop tolerance

    TPP    = lPiv == 0  ! Threshold Partial   Pivoting (normal).
    TRP    = lPiv == 1  ! Threshold Rook      Pivoting
    TCP    = lPiv == 2  ! Threshold Complete  Pivoting.
    TSP    = lPiv == 3  ! Threshold Symmetric Pivoting.
    kPiv(0)= 'PP'
    kPiv(1)= 'RP'
    kPiv(2)= 'CP'
    kPiv(3)= 'SP'

    ! Initialize output parameters.

    inform = 0
    minlen = nelem + 2*(m + n)
    numl0  = 0
    lenL   = 0
    lenU   = 0
    lrow   = 0
    mersum = 0
    nUtri  = m
    nLtri  = 0
    ndens1 = 0
    ndens2 = 0
    nrank  = 0
    nsing  = 0
    jsing  = 0
    jumin  = 0
    nslack = 0

    Amax   = zero
    Lmax   = zero
    Umax   = zero
    DUmax  = zero
    DUmin  = zero
    Akmax  = zero

    if (m > n) then
        mnkey  = '>'
    else if (m == n) then
        mnkey  = '='
    else
        mnkey  = '<'
    end if

    ! Float version of dimensions.

    dm     = m
    dn     = n
    delem  = nelem

    ! Initialize workspace parameters.

    luparm(26) = 0             ! ncp
    if (lena < minlen) go to 970

    !-------------------------------------------------------------------
    ! Organize the  aij's  in  a, indc, indr.
    ! lu1or1  deletes small entries, tests for illegal  i,j's,
    !         and counts the nonzeros in each row and column.
    ! lu1or2  reorders the elements of  A  by columns.
    ! lu1or3  uses the column list to test for duplicate entries
    !         (same indices  i,j).
    ! lu1or4  constructs a row list from the column list.
    !-------------------------------------------------------------------
    call lu1or1( m   , n    , nelem, lena , small,                     &
                 a   , indc , indr , lenc , lenr,                      &
                 Amax, numnz, lerr , inform )

    if (nout > 0  .and.  lprint >= 10) then
       densty = 100.0_rp * delem / (dm * dn)
       write(nout, 1000) m, mnkey, n, numnz, Amax, densty
    end if
    if (inform /= 0) go to 930

!!! nelem  = numnz     !!! Don't change nelem.
!!! nelem is now numnz below (it might be less than the input value).

    call lu1or2( n, numnz, lena, a, indc, indr, lenc, locc )
    call lu1or3( m, n, lena, indc, lenc, locc, p, lerr, inform )

    if (inform /= 0) go to 940

    call lu1or4( m, n, numnz, lena, indc, indr, lenc, lenr, locc, locr )

    !------------------------------------------------------------------
    ! Set up lists of rows and columns with equal numbers of nonzeros,
    ! using  indc(*)  as workspace.
    ! 12 Dec 2015: Always call lu1slk here now.
    ! This sets nslack and w(j) = 1.0 for slacks, else 0.0.
    !------------------------------------------------------------------
    call lu1pq1( m, n, lenr, p, iploc, ipinv, indc(numnz + 1) )
    call lu1pq1( n, m, lenc, q, iqloc, iqinv, indc(numnz + 1) )
    call lu1slk( m, n, lena, q, iqloc, a, indc, locc, nslack, w )
    luparm(30) = nslack

    !------------------------------------------------------------------
    ! For TCP, allocate Ha, Hj, Hk at the end of a, indc, indr.
    ! Then compute the factorization  A = L*U.
    !------------------------------------------------------------------
    lenH   = 0                ! Keep -Wmaybe-uninitialized happy.
    lena2  = 0                !
    locH   = 0                !
    lmaxr  = 0                !
    if (TPP .or. TSP) then
       lenH   = 1
       lena2  = lena
       locH   = lena
       lmaxr  = 1
    else if (TRP) then
       lenH   = 1             ! Dummy
       lena2  = lena  - m     ! Reduced length of      a
       locH   = lena          ! Dummy
       lmaxr  = lena2 + 1     ! Start of Amaxr      in a
    else if (TCP) then
       lenH   = n             ! Length of heap
       lena2  = lena  - lenH  ! Reduced length of      a, indc, indr
       locH   = lena2 + 1     ! Start of Ha, Hj, Hk in a, indc, indr
       lmaxr  = 1             ! Dummy
    end if

    call lu1fad( m     , n     , numnz , lena2 , luparm, parmlu,       &
                 a     , indc  , indr  , p     , q     ,               &
                 lenc  , lenr  , locc  , locr  ,                       &
                 iploc , iqloc , ipinv , iqinv , w     ,               &
                 lenH  ,a(locH), indc(locH), indr(locH), a(lmaxr),     &
                 inform, lenL  , lenU  , minlen, mersum,               &
                 nUtri , nLtri , ndens1, ndens2, nrank , nslack,       &
                 Lmax  , Umax  , DUmax , DUmin , Akmax )

    luparm(16) = nrank
    luparm(23) = lenL
    if (inform == 7) go to 970
    if (inform == 9) go to 985
    if (inform ==10) go to 981
    if (inform >  0) go to 980

    if ( keepLU ) then
       !---------------------------------------------------------------
       ! The LU factors are at the top of  a, indc, indr,
       ! with the columns of  L  and the rows of  U  in the order
       !
       ! ( free )   ... ( u3 ) ( l3 ) ( u2 ) ( l2 ) ( u1 ) ( l1 ).
       !
       ! Starting with ( l1 ) and ( u1 ), move the rows of  U  to the
       ! left and the columns of  L  to the right, giving
       !
       ! ( u1 ) ( u2 ) ( u3 ) ...   ( free )   ... ( l3 ) ( l2 ) ( l1 ).
       !
       ! Also, set  numl0 = the number of nonempty columns of L.
       !---------------------------------------------------------------
       lu     = 0
       ll     = lena  + 1
       lm     = lena2 + 1
       ltopl  = ll - lenL - lenU
       lrow   = lenU

       do k = 1, nrank
          i       =   p(k)
          lenUk   = - lenr(i)
          lenr(i) =   lenUk
          j       =   q(k)
          lenLk   = - lenc(j) - 1
          if (lenLk > 0) then
             numl0        = numl0 + 1
             iqloc(numl0) = lenLk
          end if

          if (lu + lenUk < ltopl) then
             !=========================================================
             ! There is room to move ( uk ).  Just right-shift ( lk ).
             !=========================================================
             do idummy = 1, lenLk
                ll       = ll - 1
                lm       = lm - 1
                a(ll)    = a(lm)
                indc(ll) = indc(lm)
                indr(ll) = indr(lm)
             end do
          else
             !=========================================================
             ! There is no room for ( uk ) yet.  We have to
             ! right-shift the whole of the remaining LU file.
             ! Note that ( lk ) ends up in the correct place.
             !=========================================================
             llsave = ll - lenLk
             nmove  = lm - ltopl

             do idummy = 1, nmove
                ll       = ll - 1
                lm       = lm - 1
                a(ll)    = a(lm)
                indc(ll) = indc(lm)
                indr(ll) = indr(lm)
             end do

             ltopl  = ll
             ll     = llsave
             lm     = ll
          end if

          !======================================================
          ! Left-shift ( uk ).
          !======================================================
          locr(i) = lu + 1
          l2      = lm - 1
          lm      = lm - lenUk

          do l = lm, l2
             lu       = lu + 1
             a(lu)    = a(l)
             indr(lu) = indr(l)
          end do
       end do

       !---------------------------------------------------------------
       ! Save the lengths of the nonempty columns of  L,
       ! and initialize  locc(j)  for the LU update routines.
       !---------------------------------------------------------------
       lenc(1:numl0) = iqloc(1:numl0)
       locc(1:n)     = 0

       !---------------------------------------------------------------
       ! Test for singularity.
       ! lu6chk  sets  nsing, jsing, jumin, Lmax, Umax, DUmax, DUmin
       ! (including entries from the dense LU).
       ! input      i1 = 1 means we're calling lu6chk from LUSOL.
       ! output inform = 1 if there are singularities (nsing > 0).
       ! 12 Dec 2015: nslack is now an input.
       !---------------------------------------------------------------
       call lu6chk( i1, m, n, nslack, w, lena, luparm, parmlu,         &
                    a, indc, indr, p, q,                               &
                    lenc, lenr, locc, locr, inform )
       nsing  = luparm(11)
       jsing  = luparm(12)
       jumin  = luparm(19)
       Lmax   = parmlu(11)
       Umax   = parmlu(12)
       DUmax  = parmlu(13)
       DUmin  = parmlu(14)

    else
       !---------------------------------------------------------------
       ! keepLU = 0.  L and U were not kept, just the diagonals of U.
       ! lu1fac will probably be called again soon with keepLU = .true.
       ! 11 Mar 2001: lu6chk revised.  We can call it with keepLU = 0,
       !              but we want to keep Lmax, Umax from lu1fad.
       ! 05 May 2002: Allow for TCP with new lu1DCP.  Diag(U) starts
       !              below lena2, not lena.  Need lena2 in next line.
       ! 12 Dec 2015: nslack is now an input.
       !---------------------------------------------------------------
       call lu6chk( i1, m, n, nslack, w, lena2, luparm, parmlu,        &
                    a, indc, indr, p, q,                               &
                    lenc, lenr, locc, locr, inform )
       nsing  = luparm(11)
       jsing  = luparm(12)
       jumin  = luparm(19)
       DUmax  = parmlu(13)
       DUmin  = parmlu(14)
    end if

    go to 990

    !------------
    ! Error exits.
    !------------
930 inform = 3
    if (lprint >= 0) write(nout, 1300) lerr, indc(lerr), indr(lerr)
    go to 990

940 inform = 4
    if (lprint >= 0) write(nout, 1400) lerr, indc(lerr), indr(lerr)
    go to 990

970 inform = 7
    if (lprint >= 0) write(nout, 1700) lena, minlen
    go to 990

980 inform = 8
    if (lprint >= 0) write(nout, 1800)
    go to 990

981 inform = 10
    go to 990

985 inform = 9
    if (lprint >= 0) write(nout, 1900)

    ! Store output parameters.

990 luparm(10) = inform
    luparm(11) = nsing
    luparm(12) = jsing
    luparm(13) = minlen
    luparm(15) = 0
    luparm(16) = nrank
    luparm(17) = ndens1
    luparm(18) = ndens2
    luparm(19) = jumin
    luparm(20) = numl0
    luparm(21) = lenL
    luparm(22) = lenU
    luparm(23) = lenL
    luparm(24) = lenU
    luparm(25) = lrow
    luparm(27) = mersum
    luparm(28) = nUtri
    luparm(29) = nLtri

    parmlu(10) = Amax
    parmlu(11) = Lmax
    parmlu(12) = Umax
    parmlu(13) = DUmax
    parmlu(14) = DUmin
    parmlu(15) = Akmax

    Agrwth = Akmax  / (Amax + 1.0e-20_rp)
    Ugrwth = Umax   / (Amax + 1.0e-20_rp)
    if ( TPP ) then
        growth = Ugrwth
    else ! TRP or TCP or TSP
        growth = Agrwth
    end if
    parmlu(16) = growth

    !------------------------------------------------------------------
    ! Print statistics for the LU factors.
    !------------------------------------------------------------------
    ncp    = luparm(26)
    condU  = DUmax / max( DUmin, 1.0e-20_rp )
    dincr  = lenL + lenU - nelem
    dincr  = dincr * 100.0_rp / max( delem, one )
    avgmer = mersum
    avgmer = avgmer / dm
    nbump  = m - nUtri - nLtri

    if (nout > 0  .and.  lprint >= 10) then
       if ( TPP ) then
          write(nout, 1100) avgmer, lenL, lenL+lenU, ncp, dincr,       &
                            nUtri, lenU, Ltol, Umax, Ugrwth,           &
                            nLtri, ndens1, Lmax

       else
          write(nout, 1120) kPiv(lPiv), avgmer,                        &
                            lenL, lenL+lenU, ncp, dincr,               &
                            nUtri, lenU, Ltol, Umax, Ugrwth,           &
                            nLtri, ndens1, Lmax, Akmax, Agrwth
       end if

       write(nout, 1200) nbump, ndens2, DUmax, DUmin, condU
    end if

    return

1000 format(' m', i12, ' ', a, 'n', i12, '  Elems', i9,                &
            '  Amax', es10.1, '  Density', f7.2)
1100 format(' Merit', f8.1, '  lenL', i9, '  L+U', i11,                &
            '  Cmpressns', i5, '  Incres', f8.2                        &
      /     ' Utri', i9, '  lenU', i9, '  Ltol', es10.2,               &
            '  Umax', es10.1, '  Ugrwth', es8.1                        &
      /     ' Ltri', i9, '  dense1', i7, '  Lmax', es10.2)
1120 format(' Mer', a2, f8.1, '  lenL', i9, '  L+U', i11,              &
            '  Cmpressns', i5, '  Incres', f8.2                        &
      /     ' Utri', i9, '  lenU', i9, '  Ltol', es10.2,               &
            '  Umax', es10.1, '  Ugrwth', es8.1                        &
      /     ' Ltri', i9, '  dense1', i7, '  Lmax', es10.2,             &
            '  Akmax', es9.1, '  Agrwth', es8.1)
1200 format(' bump', i9, '  dense2', i7, '  DUmax', es9.1,             &
            '  DUmin', es9.1, '  condU', es9.1)
1300 format(/ ' lu1fac  error...  entry  a(', i8, ')  has an illegal', &
              ' row or column index'                                   &
            //' indc, indr =', 2i8)
1400 format(/ ' lu1fac  error...  entry  a(', i8, ')  has the same',   &
              ' indices as an earlier entry'                           &
            //' indc, indr =', 2i8)
1700 format(/ ' lu1fac  error...  insufficient storage'                &
            //' Increase  lena  from', i10, '  to at least', i10)
1800 format(/ ' lu1fac  error...  fatal bug',                          &
              '   (sorry --- this should never happen)')
1900 format(/ ' lu1fac  error...  TSP used but',                       &
              ' diagonal pivot could not be found')

  end subroutine lu1fac