************************************************************************
*     File mi27lu.f and sn27lu.f.
*
*     This file contains most of the sparse LU package LUSOL
*     (the parts needed by MINOS, SQOPT and SNOPT).
*     The parts included are
*
*        27HEAD.f    These comments.
*        lusol1.f    Factor a given matrix A from scratch (lu1fac).
*        lusol2.f    Heap-management routines for lu1fac.
*        lusol6a.f   Solve with the current LU factors.
*        lusol7a.f   Utilities for all update routines.
*        lusol8a.f   Replace a column (Bartels-Golub update).
*
!***********************************************************************
!
!     File lusol1.f
!
!     lu1fac   lu1fad   lu1gau   lu1mar   lu1mRP   lu1mCP   lu1mSP
!     lu1pen   lu1mxc   lu1mxr   lu1or1   lu1or2   lu1or3   lu1or4
!     lu1pq1   lu1pq2   lu1pq3   lu1rec   lu1slk
!     lu1ful   lu1DPP   lu1DCP
!
! 26 Apr 2002: TCP implemented using heap data structure.
! 01 May 2002: lu1DCP implemented.
! 07 May 2002: lu1mxc must put 0.0 at top of empty columns.
! 09 May 2002: lu1mCP implements Markowitz with cols searched
!              in heap order.
!              Often faster (searching 20 or 40 cols) but more dense.
! 11 Jun 2002: TRP implemented.
!              lu1mRP implements Markowitz with Threshold Rook Pivoting.
!              lu1mxc maintains max col elements.  (Previously lu1max.)
!              lu1mxr maintains max row elements.
! 12 Jun 2002: lu1mCP seems too slow on big problems (e.g. memplus).
!              Disabled it for the moment.  (Use lu1mar + TCP.)
! 14 Dec 2002: TSP implemented.
!              lu1mSP implements Markowitz with
!              Threshold Symmetric Pivoting.
! 07 Mar 2003: character*1, character*2 changed to f90 form.
!              Comments changed from * in column to ! in column 1.
!              Comments kept within column 72 to avoid compiler warning.
! 19 Dec 2004: Hdelete(...) has new input argument Hlenin.
! 21 Dec 2004: Print Ltol and Lmax with e10.2 instead of e10.1.
! 26 Mar 2006: lu1fad: Ignore nsing from lu1ful.
!              lu1DPP: nsing redefined (but not used by lu1fad).
!              lu1DCP: nsing redefined (but not used by lu1fad).
! 30 Mar 2011: In lu1pen, loop 620 revised following advice from Ralf Östermark,
!              School of Business and Economics at Åbo Akademi University.
!              The previous version produced a core dump on the Cray XT.
! 06 Jun 2013: sn28lusol.f90 of 03 Apr 2013 contains updated f90 lu1mxr
!              to improve the efficiency of TRP.  Adapt this for lusol.f
!              = mi27lu.f and sn27lu.f.  The new lu1fad and lu1mxr need
!              3 extra work vectors that are implemented as automatic arrays
!              (not a feature of f77) but this should be ok as long as we
!              use an f90 compiler.
!              Ding Ma and Michael Saunders, Stanford University.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena)   , w(n)
      integer            indc(lena), indr(lena), ip(m)   , iq(n),
     &                   lenc(n)   , lenr(m)   ,
     &                   iploc(n)  , iqloc(m)  , ipinv(m), iqinv(n)
      integer            locc(n)   , locr(m)
!     ------------------------------------------------------------------
!     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 ip, iq).
!     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, ip contains a complete
!     row permutation.  At the start of stage k,  ip(1), ..., ip(k-1)
!     are the first k-1 rows of the final row permutation P.
!     The remaining rows are stored in an ordered list
!                          ( ip, iploc, ipinv )
!     where
!           iploc(nz) points to the start in ip(*) of the set of rows
!                     that currently contain nz nonzeros,
!           ipinv(i)  points to the position of row i in ip(*).
!
!     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 iq, 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 ip 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.
!
!     Systems Optimization Laboratory, Stanford University.
!  ---------------------------------------------------------------------
!
!
!  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.
!  ip, iq    are the row and column permutations defining the
!         pivot order.  For example, row ip(1) and column iq(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.
!
!  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) =
!
!
!
!  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          mnkey*1, kPiv(0:3)*2
      integer            lPiv
      logical            keepLU, TCP, TPP, TRP, TSP
      double precision   Lmax, Ltol

      double precision   zero         ,  one
      parameter        ( zero = 0.0d+0,  one = 1.0d+0 )

!     Grab relevant input parameters.

      nelem0 = nelem
      nout   = luparm(1)
      lprint = luparm(2)
      lPiv   = luparm(6)
      keepLU = luparm(8) .ne. 0

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

      TPP    = lPiv .eq. 0  ! Threshold Partial   Pivoting (normal).
      TRP    = lPiv .eq. 1  ! Threshold Rook      Pivoting
      TCP    = lPiv .eq. 2  ! Threshold Complete  Pivoting.
      TSP    = lPiv .eq. 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

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

      if (m .gt. n) then
         mnkey  = '>'
      else if (m .eq. 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 .lt. 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. gt. 0  .and.  lprint .ge. 10) then
         densty = 100.0d+0 * delem / (dm * dn)
         write(nout, 1000) m, mnkey, n, nelem, Amax, densty
      end if
      if (inform .ne. 0) go to 930

      nelem  = numnz

      call lu1or2( n, nelem, lena, a, indc, indr, lenc, locc )
      call lu1or3( m, n, lena, indc, lenc, locc, ip,
     &             lerr, inform )

      if (inform .ne. 0) go to 940

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

!     ------------------------------------------------------------------
!     Set up lists of rows and columns with equal numbers of nonzeros,
!     using  indc(*)  as workspace.
!     ------------------------------------------------------------------
      call lu1pq1( m, n, lenr, ip, iploc, ipinv, indc(nelem + 1) )
      call lu1pq1( n, m, lenc, iq, iqloc, iqinv, indc(nelem + 1) )

!     ------------------------------------------------------------------
!     For TCP, allocate Ha, Hj, Hk at the end of a, indc, indr.
!     Then compute the factorization  A = L*U.
!     ------------------------------------------------------------------
      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     , nelem , lena2 , luparm, parmlu,
     &             a     , indc  , indr  , ip    , iq    ,
     &             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 ,
     &             Lmax  , Umax  , DUmax , DUmin , Akmax )

      luparm(16) = nrank
      luparm(23) = lenL
      if (inform .eq. 7) go to 970
      if (inform .eq. 9) go to 985
      if (inform .gt. 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       =   ip(k)
            lenUk   = - lenr(i)
            lenr(i) =   lenUk
            j       =   iq(k)
            lenLk   = - lenc(j) - 1
            if (lenLk .gt. 0) then
                numl0        = numl0 + 1
                iqloc(numl0) = lenLk
            end if

            if (lu + lenUk .lt. 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.
!        ---------------------------------------------------------------
         do k = 1, numl0
            lenc(k) = iqloc(k)
         end do

         do j = 1, n
            locc(j) = 0
         end do

!        ---------------------------------------------------------------
!        Test for singularity.
!        lu6chk  sets  nsing, jsing, jumin, Lmax, Umax, DUmax, DUmin
!        (including entries from the dense LU).
!        inform = 1  if there are singularities (nsing gt 0).
!        ---------------------------------------------------------------
         call lu6chk( 1, m, n, w, lena, luparm, parmlu,
     &                a, indc, indr, ip, iq,
     &                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.
!        ---------------------------------------------------------------
         call lu6chk( 1, m, n, w, lena2, luparm, parmlu,
     &                a, indc, indr, ip, iq,
     &                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 .ge. 0) write(nout, 1300) lerr, indc(lerr), indr(lerr)
      go to 990

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

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

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

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

!     Store output parameters.

  990 nelem      = nelem0
      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.0d-20)
      Ugrwth = Umax   / (Amax + 1.0d-20)
      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.0d-20 )
      dincr  = lenL + lenU - nelem
      dincr  = dincr * 100.0d+0 / max( delem, one )
      avgmer = mersum
      avgmer = avgmer / dm
      nbump  = m - nUtri - nLtri

      if (nout. gt. 0  .and.  lprint .ge. 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', 1p, e10.1, '  Density', 0p, f7.2)
 1100 format(' Merit', 0p, f8.1, '  lenL', i9, '  L+U', i11,
     &       '  Cmpressns', i5, '  Incres', 0p, f8.2
     & /     ' Utri', i9, '  lenU', i9, '  Ltol', 1p, e10.2,
     &       '  Umax', e10.1, '  Ugrwth', e8.1
     & /     ' Ltri', i9, '  dense1', i7, '  Lmax', e10.2)
 1120 format(' Mer', a2, 0p, f8.1, '  lenL', i9, '  L+U', i11,
     &       '  Cmpressns', i5, '  Incres', 0p, f8.2
     & /     ' Utri', i9, '  lenU', i9, '  Ltol', 1p, e10.2,
     &       '  Umax', e10.1, '  Ugrwth', e8.1
     & /     ' Ltri', i9, '  dense1', i7, '  Lmax', e10.2,
     &       '  Akmax', e9.1, '  Agrwth', e8.1)
 1200 format(' bump', i9, '  dense2', i7, '  DUmax', 1p, e9.1,
     &       '  DUmin', e9.1, '  condU', e9.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

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1fad( m     , n     , nelem , lena  , luparm, parmlu,
     &                   a     , indc  , indr  , ip    , iq    ,
     &                   lenc  , lenr  , locc  , locr  ,
     &                   iploc , iqloc , ipinv , iqinv , w     ,
     &                   lenH  , Ha    , Hj    , Hk    , Amaxr ,
     &                   inform, lenL  , lenU  , minlen, mersum,
     &                   nUtri , nLtri , ndens1, ndens2, nrank ,
     &                   Lmax  , Umax  , DUmax , DUmin , Akmax )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena), Amaxr(m), w(n)
      double precision   Ha(lenH)
      integer            indc(lena), indr(lena), ip(m), iq(n)
      integer            lenc(n)   , lenr(m)
      integer            locc(n)   , locr(m)
      integer            iploc(n)  , iqloc(m), ipinv(m), iqinv(n)
      integer            Hj(lenH)  , Hk(lenH)
      double precision   Lmax

!     ------------------------------------------------------------------
!     lu1fad  is a driver for the numerical phase of lu1fac.
!     At each stage it computes a column of  L  and a row of  U,
!     using a Markowitz criterion to select the pivot element,
!     subject to a stability criterion 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.
!                  To save time on highly triangular matrices, we wait
!                  until there are no further columns of length 1
!                  before setting and maintaining that property.
!
!     12 Apr 1989  ipinv and iqinv added (inverses of ip and iq)
!                  to save searching ip and iq for rows and columns
!                  altered in each elimination step.  (Used in lu1pq2)
!
!     19 Apr 1989  Code segmented to reduce its size.
!                  lu1gau does most of the Gaussian elimination work.
!                  lu1mar does just the Markowitz search.
!                  lu1mxc moves biggest elements to top of columns.
!                  lu1pen deals with pending fill-in in the row list.
!                  lu1pq2 updates the row and column permutations.
!
!     26 Apr 1989  maxtie replaced by maxcol, maxrow in the Markowitz
!                  search.  maxcol, maxrow change as density increases.
!
!     25 Oct 1993  keepLU implemented.
!
!     07 Feb 1994  Exit main loop early to finish off with a dense LU.
!                  densLU tells lu1fad whether to do it.
!     21 Dec 1994  Bug fixed.  nrank was wrong after the call to lu1ful.
!     12 Nov 1999  A parallel version of dcopy gave trouble in lu1ful
!                  during left-shift of dense matrix D within a(*).
!                  Fixed this unexpected problem here in lu1fad
!                  by making sure the first and second D don't overlap.
!
!     13 Sep 2000  TCP (Threshold Complete Pivoting) implemented.
!                  lu2max added
!                  (finds aijmax from biggest elems in each col).
!                  Utri, Ltri and Spars1 phases apply.
!                  No switch to Dense CP yet.  (Only TPP switches.)
!     14 Sep 2000  imax needed to remember row containing aijmax.
!     22 Sep 2000  For simplicity, lu1mxc always fixes
!                  all modified cols.
!                  (TPP spars2 used to fix just the first maxcol cols.)
!     08 Nov 2000: Speed up search for aijmax.
!                  Don't need to search all columns if the elimination
!                  didn't alter the col containing the current aijmax.
!     21 Nov 2000: lu1slk implemented for Utri phase with TCP
!                  to guard against deceptive triangular matrices.
!                  (Utri used to have aijtol >= 0.9999 to include
!                  slacks, but this allows other 1s to be accepted.)
!                  Utri now accepts slacks, but applies normal aijtol
!                  test to other pivots.
!     28 Nov 2000: TCP with empty cols must call lu1mxc and lu2max
!                  with ( lq1, n, ... ), not just ( 1, n, ... ).
!     23 Mar 2001: lu1fad bug with TCP.
!                  A col of length 1 might not be accepted as a pivot.
!                  Later it appears in a pivot row and temporarily
!                  has length 0 (when pivot row is removed
!                  but before the column is filled in).  If it is the
!                  last column in storage, the preceding col also thinks
!                  it is "last".  Trouble arises when the preceding col
!                  needs fill-in -- it overlaps the real "last" column.
!                  (Very rarely, same trouble might have happened if
!                  the drop tolerance caused columns to have length 0.)
!
!                  Introduced ilast to record the last row in row file,
!                             jlast to record the last col in col file.
!                  lu1rec returns ilast = indr(lrow + 1)
!                              or jlast = indc(lcol + 1).
!                  (Should be an output parameter, but didn't want to
!                  alter lu1rec's parameter list.)
!                  lu1rec also treats empty rows or cols safely.
!                  (Doesn't eliminate them!)
!
!     26 Apr 2002: Heap routines added for TCP.
!                  lu2max no longer needed.
!                  imax, jmax used only for printing.
!     01 May 2002: lu1DCP implemented (dense complete pivoting).
!                  Both TPP and TCP now switch to dense LU
!                  when density exceeds dens2.
!     06 May 2002: In dense mode, store diag(U) in natural order.
!     09 May 2002: lu1mCP implemented (Markowitz TCP via heap).
!     11 Jun 2002: lu1mRP implemented (Markowitz TRP).
!     28 Jun 2002: Fixed call to lu1mxr.
!     14 Dec 2002: lu1mSP implemented (Markowitz TSP).
!     15 Dec 2002: Both TPP and TSP can grab cols of length 1
!                  during Utri.
!     19 Dec 2004: Hdelete(...) has new input argument Hlenin.
!     26 Mar 2006: lu1fad returns nrank  = min( mrank, nrank )
!                  and ignores nsing from lu1ful
!     26 Mar 2006: Allow for empty columns before calling Hbuild.
!     03 Apr 2013: f90 lu1mxr recoded to improve efficiency of TRP.
!     06 Jun 2013: Adapted f90 lu1mxr for use in this f77 version.
!     Systems Optimization Laboratory, Stanford University.
!     ------------------------------------------------------------------

      logical            Utri, Ltri, spars1, spars2, dense
      logical            densLU, keepLU
      logical            TCP, TPP, TRP, TSP
      double precision   Lij, Ltol, small
      integer            Hlen, Hlenin, hops, h, lPiv

!     The following f90 automatic arrays were not allowed in f77
!     but they should be ok if we use an f90 compiler.
      integer            mark, markc(n), markr(m)

      integer            i1
      parameter        ( i1   = 1 )
      double precision   zero,           one
      parameter        ( zero = 0.0d+0,  one  = 1.0d+0 )

!     ------------------------------------------------------------------
!     Local variables
!     ---------------
!
!     lcol   is the length of the column file.  It points to the last
!            nonzero in the column list.
!     lrow   is the analogous quantity for the row file.
!     lfile  is the file length (lcol or lrow) after the most recent
!            compression of the column list or row list.
!     nrowd  and  ncold  are the number of rows and columns in the
!            matrix defined by the pivot column and row.  They are the
!            dimensions of the submatrix D being altered at this stage.
!     melim  and  nelim  are the number of rows and columns in the
!            same matrix D, excluding the pivot column and row.
!     mleft  and  nleft  are the number of rows and columns
!            still left to be factored.
!     nzchng is the increase in nonzeros in the matrix that remains
!            to be factored after the current elimination
!            (usually negative).
!     nzleft is the number of nonzeros still left to be factored.
!     nspare is the space we leave at the end of the last row or
!            column whenever a row or column is being moved to the end
!            of its file.  nspare = 1 or 2 might help reduce the
!            number of file compressions when storage is tight.
!
!     The row and column ordering permutes A into the form
!
!                        ------------------------
!                         \                     |
!                          \         U1         |
!                           \                   |
!                            --------------------
!                            |\
!                            | \
!                            |  \
!            P A Q   =       |   \
!                            |    \
!                            |     --------------
!                            |     |            |
!                            |     |            |
!                            | L1  |     A2     |
!                            |     |            |
!                            |     |            |
!                            --------------------
!
!     where the block A2 is factored as  A2 = L2 U2.
!     The phases of the factorization are as follows.
!
!     Utri   is true when U1 is being determined.
!            Any column of length 1 is accepted immediately (if TPP).
!
!     Ltri   is true when L1 is being determined.
!            lu1mar exits as soon as an acceptable pivot is found
!            in a row of length 1.
!
!     spars1 is true while the density of the (modified) A2 is less
!            than the parameter dens1 = parmlu(7) = 0.3 say.
!            lu1mar searches maxcol columns and maxrow rows,
!            where  maxcol = luparm(3),  maxrow = maxcol - 1.
!            lu1mxc is used to keep the biggest element at the top
!            of all remaining columns.
!
!     spars2 is true while the density of the modified A2 is less
!            than the parameter dens2 = parmlu(8) = 0.6 say.
!            lu1mar searches maxcol columns and no rows.
!            lu1mxc could fix up only the first maxcol cols (with TPP).
!            22 Sep 2000:  For simplicity, lu1mxc fixes all
!                          modified cols.
!
!     dense  is true once the density of A2 reaches dens2.
!            lu1mar searches only 1 column (the shortest).
!            lu1mxc could fix up only the first column (with TPP).
!            22 Sep 2000:  For simplicity, lu1mxc fixes all
!                          modified cols.
!
!     ------------------------------------------------------------------
!     To eliminate timings, comment out all lines containing "time".
!     ------------------------------------------------------------------

!     integer            eltime, mktime

!     call timer ( 'start ', 3 )
!     ntime  = (n / 4.0)


      nout   = luparm(1)
      lprint = luparm(2)
      maxcol = luparm(3)
      lPiv   = luparm(6)
      keepLU = luparm(8) .ne. 0

      TPP    = lPiv .eq. 0  ! Threshold Partial   Pivoting (normal).
      TRP    = lPiv .eq. 1  ! Threshold Rook      Pivoting
      TCP    = lPiv .eq. 2  ! Threshold Complete  Pivoting.
      TSP    = lPiv .eq. 3  ! Threshold Symmetric Pivoting.

      densLU = .false.
      maxrow = maxcol - 1
      ilast  = m                 ! Assume row m is last in the row file.
      jlast  = n                 ! Assume col n is last in the col file.
      lfile  = nelem
      lrow   = nelem
      lcol   = nelem
      minmn  = min( m, n )
      maxmn  = max( m, n )
      nzleft = nelem
      nspare = 1

      if ( keepLU ) then
         lu1    = lena   + 1
      else
!        Store only the diagonals of U in the top of memory.
         ldiagU = lena   - n
         lu1    = ldiagU + 1
      end if

      Ltol   = parmlu(1)
      small  = parmlu(3)
      uspace = parmlu(6)
      dens1  = parmlu(7)
      dens2  = parmlu(8)
      Utri   = .true.
      Ltri   = .false.
      spars1 = .false.
      spars2 = .false.
      dense  = .false.

!     Check parameters.

      Ltol   = max( Ltol, 1.0001d+0 )
      dens1  = min( dens1, dens2 )

!     Initialize output parameters.
!     lenL, lenU, minlen, mersum, nUtri, nLtri, ndens1, ndens2, nrank
!     are already initialized by lu1fac.

      Lmax   = zero
      Umax   = zero
      DUmax  = zero
      DUmin  = 1.0d+20
      if (nelem .eq. 0) Dumin = zero
      Akmax  = zero
      hops   = 0

      ! More initialization.

      if (TPP .or. TSP) then ! Don't worry yet about lu1mxc.
         aijmax = zero
         aijtol = zero
         Hlen   = 1

      else ! TRP or TCP
         ! Move biggest element to top of each column.
         ! Set w(*) to mark slack columns (unit vectors).

         call lu1mxc( i1, n, iq, a, indc, lenc, locc )
         call lu1slk( m, n, lena, iq, iqloc, a, locc, w )
      end if

      if (TRP) then
         ! Find biggest element in each row.
         ! 06 Jun 2013: Call new lu1mxr.

         mark = 0
         call lu1mxr( mark, i1, m, m, n, lena,
     &                a, indc, lenc, locc, indr, lenr, locr,
     &                ip, markc, markr, Amaxr )
      end if

      if (TCP) then
         ! Set Ha(1:Hlen) = biggest element in each column,
         !     Hj(1:Hlen) = corresponding column indices.

         Hlen  = 0
         do kk = 1, n
            Hlen     = Hlen + 1
            j        = iq(kk)
            if (lenc(j) .gt. 0) then
               lc    = locc(j)
               amax  = abs( a(lc) )
            else
               amax  = zero
            end if
            Ha(Hlen) = amax
            Hj(Hlen) = j
            Hk(j)    = Hlen
         end do

         ! Build the heap, creating new Ha, Hj and setting Hk(1:Hlen).

         call Hbuild( Ha, Hj, Hk, Hlen, Hlen, hops )
      end if

!     ------------------------------------------------------------------
!     Start of main loop.
!     ------------------------------------------------------------------
      mleft  = m + 1
      nleft  = n + 1

      do 800 nrowu = 1, minmn

!        mktime = (nrowu / ntime) + 4
!        eltime = (nrowu / ntime) + 9
         mleft  = mleft - 1
         nleft  = nleft - 1

!        Bail out if there are no nonzero rows left.

         if (iploc(1) .gt. m) go to 900

         ! For TCP, the largest Aij is at the top of the heap.

         if ( TCP ) then
            aijmax = Ha(1)      ! Marvelously easy !
            Akmax  = max( Akmax, aijmax )
            aijtol = aijmax / Ltol
         end if

!        ===============================================================
!        Find a suitable pivot element.
!        ===============================================================

         if ( Utri ) then
!           ------------------------------------------------------------
!           So far all columns have had length 1.
!           We are still looking for the (backward) triangular part of A
!           that forms the first rows and columns of U.
!           ------------------------------------------------------------

            lq1    = iqloc(1)
            lq2    = n
            if (m   .gt.   1) lq2 = iqloc(2) - 1

            if (lq1 .le. lq2) then  ! There are more cols of length 1.
               if (TPP .or. TSP) then
                  jbest  = iq(lq1)  ! Grab the first one.

               else ! TRP or TCP    ! Scan all columns of length 1.
                  jbest  = 0

                  do lq = lq1, lq2
                     j      = iq(lq)
                     if (w(j) .gt. zero) then ! Accept a slack
                        jbest  = j
                        go to 250
                     end if

                     lc     = locc(j)
                     amax   = abs( a(lc) )
                     if (TRP) then
                        i      = indc(lc)
                        aijtol = Amaxr(i) / Ltol
                     end if

                     if (amax .ge. aijtol) then
                        jbest  = j
                        go to 250
                     end if
                  end do
               end if

  250          if (jbest .gt. 0) then
                  lc     = locc(jbest)
                  ibest  = indc(lc)
                  mbest  = 0
                  go to 300
               end if
            end if

!           This is the end of the U triangle.
!           We will not return to this part of the code.
!           TPP and TSP call lu1mxc for the first time
!           (to move biggest element to top of each column).

            if (lprint .ge. 50) then
               write(nout, 1100) 'Utri ended.  spars1 = true'
            end if
            Utri   = .false.
            Ltri   = .true.
            spars1 = .true.
            nUtri  =  nrowu - 1
            if (TPP .or. TSP) then
               call lu1mxc( lq1, n, iq, a, indc, lenc, locc )
            end if
         end if

         if ( spars1 ) then
!           ------------------------------------------------------------
!           Perform a Markowitz search.
!           Search cols of length 1, then rows of length 1,
!           then   cols of length 2, then rows of length 2, etc.
!           ------------------------------------------------------------
!           call timer ( 'start ', mktime )

          ! if (TPP) then ! 12 Jun 2002: Next line disables lu1mCP below
            if (TPP .or. TCP) then
               call lu1mar( m    , n     , lena  , maxmn,
     &                      TCP  , aijtol, Ltol  , maxcol, maxrow,
     &                      ibest, jbest , mbest ,
     &                      a    , indc  , indr  , ip   , iq,
     &                      lenc , lenr  , locc  , locr ,
     &                      iploc, iqloc )

            else if (TRP) then
               call lu1mRP( m    , n     , lena  , maxmn,
     &                      Ltol , maxcol, maxrow,
     &                      ibest, jbest , mbest ,
     &                      a    , indc  , indr  , ip   , iq,
     &                      lenc , lenr  , locc  , locr ,
     &                      iploc, iqloc , Amaxr )

!           else if (TCP) then ! Disabled by test above
!              call lu1mCP( m    , n     , lena  , aijtol,
!    &              ibest, jbest , mbest ,
!    &              a    , indc  , indr  ,
!    &              lenc , lenr  , locc  ,
!    &              Hlen , Ha    , Hj    )

            else if (TSP) then
               call lu1mSP( m    , n     , lena  , maxmn,
     &                      Ltol , maxcol,
     &                      ibest, jbest , mbest ,
     &                      a    , indc  , iq    , locc , iqloc )
               if (ibest .eq. 0) go to 990
            end if

!           call timer ( 'finish', mktime )

            if ( Ltri ) then

!              So far all rows have had length 1.
!              We are still looking for the (forward) triangle of A
!              that forms the first rows and columns of L.

               if (mbest .gt. 0) then
                   Ltri   = .false.
                   nLtri  =  nrowu - 1 - nUtri
                   if (lprint .ge. 50) then
                      write(nout, 1100) 'Ltri ended.'
                   end if
               end if

            else

!              See if what's left is as dense as dens1.

               if (nzleft  .ge.  (dens1 * mleft) * nleft) then
                   spars1 = .false.
                   spars2 = .true.
                   ndens1 =  nleft
                   maxrow =  0
                   if (lprint .ge. 50) then
                      write(nout, 1100) 'spars1 ended.  spars2 = true'
                   end if
               end if
            end if

         else if ( spars2 .or. dense ) then
!           ------------------------------------------------------------
!           Perform a restricted Markowitz search,
!           looking at only the first maxcol columns.  (maxrow = 0.)
!           ------------------------------------------------------------
!           call timer ( 'start ', mktime )

          ! if (TPP) then ! 12 Jun 2002: Next line disables lu1mCP below
            if (TPP .or. TCP) then
               call lu1mar( m    , n     , lena  , maxmn,
     &              TCP  , aijtol, Ltol  , maxcol, maxrow,
     &              ibest, jbest , mbest ,
     &              a    , indc  , indr  , ip   , iq,
     &              lenc , lenr  , locc  , locr ,
     &              iploc, iqloc )

            else if (TRP) then
               call lu1mRP( m    , n     , lena  , maxmn,
     &              Ltol , maxcol, maxrow,
     &              ibest, jbest , mbest ,
     &              a    , indc  , indr  , ip   , iq,
     &              lenc , lenr  , locc  , locr ,
     &              iploc, iqloc , Amaxr )

!           else if (TCP) then ! Disabled by test above
!              call lu1mCP( m    , n     , lena  , aijtol,
!    &              ibest, jbest , mbest ,
!    &              a    , indc  , indr  ,
!    &              lenc , lenr  , locc  ,
!    &              Hlen , Ha    , Hj    )

            else if (TSP) then
               call lu1mSP( m    , n     , lena  , maxmn,
     &                      Ltol , maxcol,
     &                      ibest, jbest , mbest ,
     &                      a    , indc  , iq    , locc , iqloc )
               if (ibest .eq. 0) go to 985
            end if

!           call timer ( 'finish', mktime )

!           See if what's left is as dense as dens2.

            if ( spars2 ) then
               if (nzleft  .ge.  (dens2 * mleft) * nleft) then
                   spars2 = .false.
                   dense  = .true.
                   ndens2 =  nleft
                   maxcol =  1
                   if (lprint .ge. 50) then
                      write(nout, 1100) 'spars2 ended.  dense = true'
                   end if
               end if
            end if
         end if

!        ---------------------------------------------------------------
!        See if we can finish quickly.
!        ---------------------------------------------------------------
         if ( dense  ) then
            lenD   = mleft * nleft
            nfree  = lu1 - 1

            if (nfree .ge. 2 * lenD) then

!              There is room to treat the remaining matrix as
!              a dense matrix D.
!              We may have to compress the column file first.
!              12 Nov 1999: D used to be put at the
!                           beginning of free storage (lD = lcol + 1).
!                           Now put it at the end     (lD = lu1 - lenD)
!                           so the left-shift in lu1ful will not
!                           involve overlapping storage
!                           (fatal with parallel dcopy).
!
               densLU = .true.
               ndens2 = nleft
               lD     = lu1 - lenD
               if (lcol .ge. lD) then
                  call lu1rec( n, .true., luparm, lcol,
     &                         lena, a, indc, lenc, locc )
                  lfile  = lcol
                  jlast  = indc(lcol + 1)
               end if

               go to 900
            end if
         end if

!        ===============================================================
!        The best  aij  has been found.
!        The pivot row  ibest  and the pivot column  jbest
!        Define a dense matrix  D  of size  nrowd  by  ncold.
!        ===============================================================
  300    ncold  = lenr(ibest)
         nrowd  = lenc(jbest)
         melim  = nrowd  - 1
         nelim  = ncold  - 1
         mersum = mersum + mbest
         lenL   = lenL   + melim
         lenU   = lenU   + ncold
         if (lprint .ge. 50) then
            if (nrowu .eq. 1) then
               write(nout, 1100) 'lu1fad debug:'
            end if
            if ( TPP .or. TRP .or. TSP ) then
               write(nout, 1200) nrowu, ibest, jbest, nrowd, ncold
            else ! TCP
               jmax   = Hj(1)
               imax   = indc(locc(jmax))
               write(nout, 1200) nrowu, ibest, jbest, nrowd, ncold,
     &                           imax , jmax , aijmax
            end if
         end if

!        ===============================================================
!        Allocate storage for the next column of  L  and next row of  U.
!        Initially the top of a, indc, indr are used as follows:
!
!                   ncold       melim       ncold        melim
!
!        a      |...........|...........|ujbest..ujn|li1......lim|
!
!        indc   |...........|  lenr(i)  |  lenc(j)  |  markl(i)  |
!
!        indr   |...........| iqloc(i)  |  jfill(j) |  ifill(i)  |
!
!              ^           ^             ^           ^            ^
!              lfree   lsave             lu1         ll1          oldlu1
!
!        Later the correct indices are inserted:
!
!        indc   |           |           |           |i1........im|
!
!        indr   |           |           |jbest....jn|ibest..ibest|
!
!        ===============================================================
         if ( keepLU ) then
!           relax
         else
!           Always point to the top spot.
!           Only the current column of L and row of U will
!           take up space, overwriting the previous ones.
            lu1    = ldiagU + 1
         end if
         ll1    = lu1   - melim
         lu1    = ll1   - ncold
         lsave  = lu1   - nrowd
         lfree  = lsave - ncold

!        Make sure the column file has room.
!        Also force a compression if its length exceeds a certain limit.

         limit  = uspace * lfile  +  m  +  n  +  1000
         minfre = ncold  + melim
         nfree  = lfree  - lcol
         if (nfree .lt. minfre  .or.  lcol .gt. limit) then
            call lu1rec( n, .true., luparm, lcol,
     &                   lena, a, indc, lenc, locc )
            lfile  = lcol
            jlast  = indc(lcol + 1)
            nfree  = lfree - lcol
            if (nfree .lt. minfre) go to 970
         end if

!        Make sure the row file has room.

         minfre = melim + ncold
         nfree  = lfree - lrow
         if (nfree .lt. minfre  .or.  lrow .gt. limit) then
            call lu1rec( m, .false., luparm, lrow,
     &                   lena, a, indr, lenr, locr )
            lfile  = lrow
            ilast  = indr(lrow + 1)
            nfree  = lfree - lrow
            if (nfree .lt. minfre) go to 970
         end if

!        ===============================================================
!        Move the pivot element to the front of its row
!        and to the top of its column.
!        ===============================================================
         lpivr  = locr(ibest)
         lpivr1 = lpivr + 1
         lpivr2 = lpivr + nelim

         do 330 l = lpivr, lpivr2
            if (indr(l) .eq. jbest) go to 335
  330    continue

  335    indr(l)     = indr(lpivr)
         indr(lpivr) = jbest

         lpivc  = locc(jbest)
         lpivc1 = lpivc + 1
         lpivc2 = lpivc + melim

         do 340 l = lpivc, lpivc2
            if (indc(l) .eq. ibest) go to 345
  340    continue

  345    indc(l)     = indc(lpivc)
         indc(lpivc) = ibest
         abest       = a(l)
         a(l)        = a(lpivc)
         a(lpivc)    = abest

         if ( keepLU ) then
!           relax
         else
!           Store just the diagonal of U, in natural order.
!!!         a(ldiagU + nrowu) = abest ! This was in pivot order.
            a(ldiagU + jbest) = abest
         end if

         !==============================================================
         ! Delete pivot col from heap.
         ! Hk tells us where it is in the heap.
         !==============================================================
         if ( TCP ) then
            kbest  = Hk(jbest)
            Hlenin = Hlen
            call Hdelete( Ha, Hj, Hk, Hlenin, Hlen, n, kbest, h )
            hops   = hops + h
         end if

!        ===============================================================
!        Delete the pivot row from the column file
!        and store it as the next row of  U.
!        set  indr(lu) = 0     to initialize jfill ptrs on columns of D,
!             indc(lu) = lenj  to save the original column lengths.
!        ===============================================================
         a(lu1)    = abest
         indr(lu1) = jbest
         indc(lu1) = nrowd
         lu        = lu1

         diag      = abs( abest )
         Umax      = max(  Umax, diag )
         DUmax     = max( DUmax, diag )
         DUmin     = min( DUmin, diag )

         do 360 lr   = lpivr1, lpivr2
            lu       = lu + 1
            j        = indr(lr)
            lenj     = lenc(j)
            lenc(j)  = lenj - 1
            lc1      = locc(j)
            last     = lc1 + lenc(j)

            do 350 l = lc1, last
               if (indc(l) .eq. ibest) go to 355
  350       continue

  355       a(lu)      = a(l)
            indr(lu)   = 0
            indc(lu)   = lenj
            Umax       = max( Umax, abs( a(lu) ) )
            a(l)       = a(last)
            indc(l)    = indc(last)
            indc(last) = 0       ! Free entry
!???        if (j .eq. jlast) lcol = lcol - 1
  360    continue

!        ===============================================================
!        Delete the pivot column from the row file
!        and store the nonzeros of the next column of  L.
!        Set  indc(ll) = 0     to initialize markl(*) markers,
!             indr(ll) = 0     to initialize ifill(*) row fill-in cntrs,
!             indc(ls) = leni  to save the original row lengths,
!             indr(ls) = iqloc(i)    to save parts of  iqloc(*),
!             iqloc(i) = lsave - ls  to point to the nonzeros of  L
!                      = -1, -2, -3, ... in mark(*).
!        ===============================================================
         indc(lsave) = ncold
         if (melim .eq. 0) go to 700

         ll     = ll1 - 1
         ls     = lsave
         abest  = one / abest

         do 390 lc   = lpivc1, lpivc2
            ll       = ll + 1
            ls       = ls + 1
            i        = indc(lc)
            leni     = lenr(i)
            lenr(i)  = leni - 1
            lr1      = locr(i)
            last     = lr1 + lenr(i)

            do 380 l = lr1, last
               if (indr(l) .eq. jbest) go to 385
  380       continue

  385       indr(l)    = indr(last)
            indr(last) = 0       ! Free entry
!???        if (i .eq. ilast) lrow = lrow - 1

            a(ll)      = - a(lc) * abest
            Lij        = abs( a(ll) )
            Lmax       = max( Lmax, Lij )
!!!!! DEBUG
!           if (Lij .gt. Ltol) then
!              write( *  ,*) ' Big Lij!!!', nrowu
!              write(nout,*) ' Big Lij!!!', nrowu
!           end if

            indc(ll)   = 0
            indr(ll)   = 0
            indc(ls)   = leni
            indr(ls)   = iqloc(i)
            iqloc(i)   = lsave - ls
  390    continue

!        ===============================================================
!        Do the Gaussian elimination.
!        This involves adding a multiple of the pivot column
!        to all other columns in the pivot row.
!
!        Sometimes more than one call to lu1gau is needed to allow
!        compression of the column file.
!        lfirst  says which column the elimination should start with.
!        minfre  is a bound on the storage needed for any one column.
!        lu      points to off-diagonals of u.
!        nfill   keeps track of pending fill-in in the row file.
!        ===============================================================
         if (nelim .eq. 0) go to 700
         lfirst = lpivr1
         minfre = mleft + nspare
         lu     = 1
         nfill  = 0

! 400    call timer ( 'start ', eltime )
  400    call lu1gau( m     , melim , ncold , nspare, small ,
     &                lpivc1, lpivc2, lfirst, lpivr2, lfree , minfre,
     &                ilast , jlast , lrow  , lcol  , lu    , nfill ,
     &                a     , indc  , indr  ,
     &                lenc  , lenr  , locc  , locr  ,
     &                iqloc , a(ll1), indc(ll1),
     &                        a(lu1), indr(ll1), indr(lu1) )
!        call timer ( 'finish', eltime )

         if (lfirst .gt. 0) then

!           The elimination was interrupted.
!           Compress the column file and try again.
!           lfirst, lu and nfill have appropriate new values.

            call lu1rec( n, .true., luparm, lcol,
     &                   lena, a, indc, lenc, locc )
            lfile  = lcol
            jlast  = indc(lcol + 1)
            lpivc  = locc(jbest)
            lpivc1 = lpivc + 1
            lpivc2 = lpivc + melim
            nfree  = lfree - lcol
            if (nfree .lt. minfre) go to 970
            go to 400
         end if

!        ===============================================================
!        The column file has been fully updated.
!        Deal with any pending fill-in in the row file.
!        ===============================================================
         if (nfill .gt. 0) then

!           Compress the row file if necessary.
!           lu1gau has set nfill to be the number of pending fill-ins
!           plus the current length of any rows that need to be moved.

            minfre = nfill
            nfree  = lfree - lrow
            if (nfree .lt. minfre) then
               call lu1rec( m, .false., luparm, lrow,
     &                      lena, a, indr, lenr, locr )
               lfile  = lrow
               ilast  = indr(lrow + 1)
               lpivr  = locr(ibest)
               lpivr1 = lpivr + 1
               lpivr2 = lpivr + nelim
               nfree  = lfree - lrow
               if (nfree .lt. minfre) go to 970
            end if

!           Move rows that have pending fill-in to end of the row file.
!           Then insert the fill-in.

            call lu1pen( m     , melim , ncold , nspare, ilast,
     &                   lpivc1, lpivc2, lpivr1, lpivr2, lrow ,
     &                   lenc  , lenr  , locc  , locr  ,
     &                   indc  , indr  , indr(ll1), indr(lu1) )
         end if

!        ===============================================================
!        Restore the saved values of  iqloc.
!        Insert the correct indices for the col of L and the row of U.
!        ===============================================================
  700    lenr(ibest) = 0
         lenc(jbest) = 0

         ll          = ll1 - 1
         ls          = lsave

         do 710  lc  = lpivc1, lpivc2
            ll       = ll + 1
            ls       = ls + 1
            i        = indc(lc)
            iqloc(i) = indr(ls)
            indc(ll) = i
            indr(ll) = ibest
  710    continue

         lu          = lu1 - 1

         do 720  lr  = lpivr, lpivr2
            lu       = lu + 1
            indr(lu) = indr(lr)
  720    continue

!        ===============================================================
!        Free the space occupied by the pivot row
!        and update the column permutation.
!        Then free the space occupied by the pivot column
!        and update the row permutation.
!
!        nzchng is found in both calls to lu1pq2, but we use it only
!        after the second.
!        ===============================================================
         call lu1pq2( ncold, nzchng,
     &                indr(lpivr), indc( lu1 ), lenc, iqloc, iq, iqinv )

         call lu1pq2( nrowd, nzchng,
     &                indc(lpivc), indc(lsave), lenr, iploc, ip, ipinv )

         nzleft = nzleft + nzchng

!        ===============================================================
!        lu1mxr resets Amaxr(i) in each modified row i.
!        lu1mxc moves the largest aij to the top of each modified col j.
!        28 Jun 2002: Note that cols of L have an implicit diag of 1.0,
!                     so lu1mxr is called with ll1, not ll1+1, whereas
!                        lu1mxc is called with          lu1+1.
!        ===============================================================
         if (Utri .and. TPP) then
            ! Relax -- we're not keeping big elements at the top yet.

         else
            if (TRP  .and.  melim .gt. 0) then
!              call lu1mxr( ll1, ll, indc, Amaxr,
!    &                      a, indc, lenc, locc, indr, lenr, locr )
               ! Beware: The parts of ip that we need are in indc(ll1:ll)
               mark = mark + 1
               call lu1mxr( mark, ll1, ll, m, n, lena,            
     &                      a, indc, lenc, locc, indr, lenr, locr,
     &                      indc, markc, markr, Amaxr )
                          ! ^^^^  Here are the ip(k1:k2) needed by lu1mxr.
            end if

            if (nelim .gt. 0) then
               call lu1mxc( lu1+1, lu, indr, a, indc, lenc, locc )

               if (TCP) then ! Update modified columns in heap
                  do kk = lu1+1, lu
                     j    = indr(kk)
                     k    = Hk(j)
                     v    = abs( a(locc(j)) ) ! Biggest aij in column j
                     call Hchange( Ha, Hj, Hk, Hlen, n, k, v, j, h )
                     hops = hops + h
                  end do
               end if
            end if
         end if

!        ===============================================================
!        Negate lengths of pivot row and column so they will be
!        eliminated during compressions.
!        ===============================================================
         lenr(ibest) = - ncold
         lenc(jbest) = - nrowd

!        Test for fatal bug: row or column lists overwriting L and U.

         if (lrow .gt. lsave) go to 980
         if (lcol .gt. lsave) go to 980

!        Reset the file lengths if pivot row or col was at the end.

         if (ibest .eq. ilast) then
            lrow = locr(ibest)
         end if

         if (jbest .eq. jlast) then
            lcol = locc(jbest)
         end if
  800 continue

!     ------------------------------------------------------------------
!     End of main loop.
!     ------------------------------------------------------------------

!     ------------------------------------------------------------------
!     Normal exit.
!     Move empty rows and cols to the end of ip, iq.
!     Then finish with a dense LU if necessary.
!     ------------------------------------------------------------------
  900 inform = 0
      call lu1pq3( m, lenr, ip, ipinv, mrank )
      call lu1pq3( n, lenc, iq, iqinv, nrank )
      nrank  = min( mrank, nrank )

      if ( densLU ) then
!        call timer ( 'start ', 17 )
         call lu1ful( m     , n    , lena , lenD , lu1 , TPP,
     &                mleft , nleft, nrank, nrowu,
     &                lenL  , lenU , nsing,
     &                keepLU, small,
     &                a     , a(lD), indc , indr , ip  , iq,
     &                lenc  , lenr , locc , ipinv, locr )
!***     21 Dec 1994: Bug in next line.
!***     nrank  = nrank - nsing.  Changed to next line:
!***     nrank  = minmn - nsing

!***     26 Mar 2006: Previous line caused bug with m<n and nsing>0.
!        Don't mess with nrank any more.  Let end of lu1fac handle it.
!        call timer ( 'finish', 17 )
      end if

      minlen = lenL  +  lenU  +  2*(m + n)
      go to 990

!     Not enough space free after a compress.
!     Set  minlen  to an estimate of the necessary value of  lena.

  970 inform = 7
      minlen = lena  +  lfile  +  2*(m + n)
      go to 990

!     Fatal error.  This will never happen!
!    (Famous last words.)

  980 inform = 8
      go to 990

!     Fatal error with TSP.  Diagonal pivot not found.

  985 inform = 9

!     Exit.

  990 continue
!     call timer ( 'finish', 3 )
      return

 1100 format(/ 1x, a)
 1200 format(  ' nrowu', i7,
     &       '   i,jbest', 2i7,
     &       '   nrowd,ncold', 2i6,
     &       '   i,jmax', 2i7,
     &       '   aijmax', 1p, e10.2)

      end ! subroutine lu1fad

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1gau( m     , melim , ncold , nspare, small ,
     &                   lpivc1, lpivc2, lfirst, lpivr2, lfree , minfre,
     &                   ilast , jlast , lrow  , lcol  , lu    , nfill ,
     &                   a     , indc  , indr  ,
     &                   lenc  , lenr  , locc  , locr  ,
     &                   mark  , al    , markl ,
     &                           au    , ifill , jfill )

      implicit           double precision (a-h,o-z)
      double precision   a(*)        , al(melim)   , au(ncold)
      integer            indc(*)     , indr(*)     , lenc(*)  , lenr(*),
     &                   mark(*)     , markl(melim),
     &                   ifill(melim), jfill(ncold)
      integer            locc(*)     , locr(*)

!     ------------------------------------------------------------------
!     lu1gau does most of the work for each step of
!     Gaussian elimination.
!     A multiple of the pivot column is added to each other column j
!     in the pivot row.  The column list is fully updated.
!     The row list is updated if there is room, but some fill-ins may
!     remain, as indicated by ifill and jfill.
!
!
!  Input:
!     ilast    is the row    at the end of the row    list.
!     jlast    is the column at the end of the column list.
!     lfirst   is the first column to be processed.
!     lu + 1   is the corresponding element of U in au(*).
!     nfill    keeps track of pending fill-in.
!     a(*)     contains the nonzeros for each column j.
!     indc(*)  contains the row indices for each column j.
!     al(*)    contains the new column of L.  A multiple of it is
!              used to modify each column.
!     mark(*)  has been set to -1, -2, -3, ... in the rows
!              corresponding to nonzero 1, 2, 3, ... of the col of L.
!     au(*)    contains the new row of U.  Each nonzero gives the
!              required multiple of the column of L.
!
!  Workspace:
!     markl(*) marks the nonzeros of L actually used.
!              (A different mark, namely j, is used for each column.)
!
!  Output:
!     ilast     New last row    in the row    list.
!     jlast     New last column in the column list.
!     lfirst    = 0 if all columns were completed,
!               > 0 otherwise.
!     lu        returns the position of the last nonzero of U
!               actually used, in case we come back in again.
!     nfill     keeps track of the total extra space needed in the
!               row file.
!     ifill(ll) counts pending fill-in for rows involved in the new
!               column of L.
!     jfill(lu) marks the first pending fill-in stored in columns
!               involved in the new row of U.
!
!     16 Apr 1989: First version of lu1gau.
!     23 Apr 1989: lfirst, lu, nfill are now input and output
!                  to allow re-entry if elimination is interrupted.
!     23 Mar 2001: Introduced ilast, jlast.
!     27 Mar 2001: Allow fill-in "in situ" if there is already room
!                  up to but NOT INCLUDING the end of the
!                  row or column file.
!                  Seems safe way to avoid overwriting empty rows/cols
!                  at the end.  (May not be needed though, now that we
!                  have ilast and jlast.)
!     ------------------------------------------------------------------

      logical            atend

      do 600 lr = lfirst, lpivr2
         j      = indr(lr)
         lenj   = lenc(j)
         nfree  = lfree - lcol
         if (nfree .lt. minfre) go to 900

!        ---------------------------------------------------------------
!        Inner loop to modify existing nonzeros in column  j.
!        Loop 440 performs most of the arithmetic involved in the
!        whole LU factorization.
!        ndone counts how many multipliers were used.
!        ndrop counts how many modified nonzeros are negligibly small.
!        ---------------------------------------------------------------
         lu     = lu + 1
         uj     = au(lu)
         lc1    = locc(j)
         lc2    = lc1 + lenj - 1
         atend  = j .eq. jlast
         ndone  = 0
         if (lenj .eq. 0) go to 500

         ndrop  = 0

         do 440 l = lc1, lc2
            i        =   indc(l)
            ll       = - mark(i)
            if (ll .gt. 0) then
               ndone     = ndone + 1
               markl(ll) = j
               a(l)      = a(l)  +  al(ll) * uj
               if (abs( a(l) ) .le. small) then
                  ndrop  = ndrop + 1
               end if
            end if
  440    continue

!        ---------------------------------------------------------------
!        Remove any negligible modified nonzeros from both
!        the column file and the row file.
!        ---------------------------------------------------------------
         if (ndrop .eq. 0) go to 500
         k      = lc1

         do 480 l = lc1, lc2
            i        = indc(l)
            if (abs( a(l) ) .le. small) go to 460
            a(k)     = a(l)
            indc(k)  = i
            k        = k + 1
            go to 480

!           Delete the nonzero from the row file.

  460       lenj     = lenj    - 1
            lenr(i)  = lenr(i) - 1
            lr1      = locr(i)
            last     = lr1 + lenr(i)

            do 470 lrep = lr1, last
               if (indr(lrep) .eq. j) go to 475
  470       continue

  475       indr(lrep) = indr(last)
            indr(last) = 0
            if (i .eq. ilast) lrow = lrow - 1
  480    continue

!        Free the deleted elements from the column file.

         do 490  l  = k, lc2
            indc(l) = 0
  490    continue
         if (atend) lcol = k - 1

!        ---------------------------------------------------------------
!        Deal with the fill-in in column j.
!        ---------------------------------------------------------------
  500    if (ndone .eq. melim) go to 590

!        See if column j already has room for the fill-in.

         if (atend) go to 540
         last   = lc1  + lenj - 1
         l1     = last + 1
         l2     = last + (melim - ndone)
         ! 27 Mar 2001: Be sure it's not at or past end of the col file.
         if (l2 .ge. lcol) go to 520

         do 510 l = l1, l2
            if (indc(l) .ne. 0) go to 520
  510    continue
         go to 540

!        We must move column j to the end of the column file.
!        First, leave some spare room at the end of the
!        current last column.

  520    do 522  l  = lcol + 1, lcol + nspare
            lcol    = l
            indc(l) = 0     ! Spare space is free.
  522    continue

         atend   = .true.
         jlast   = j
         l1      = lc1
         lc1     = lcol + 1
         locc(j) = lc1

         do 525  l     = l1, last
            lcol       = lcol + 1
            a(lcol)    = a(l)
            indc(lcol) = indc(l)
            indc(l)    = 0      ! Free space.
  525    continue

!        ---------------------------------------------------------------
!        Inner loop for the fill-in in column j.
!        This is usually not very expensive.
!        ---------------------------------------------------------------
  540    last          = lc1 + lenj - 1
         ll            = 0

         do 560   lc   = lpivc1, lpivc2
            ll         = ll + 1
            if (markl(ll)  .eq. j    ) go to 560
            aij        = al(ll) * uj
            if (abs( aij ) .le. small) go to 560
            lenj       = lenj + 1
            last       = last + 1
            a(last)    = aij
            i          = indc(lc)
            indc(last) = i
            leni       = lenr(i)

!           Add 1 fill-in to row i if there is already room.
!           27 Mar 2001: Be sure it's not at or past the end
!                        of the row file.

            l          = locr(i) + leni
            if (     l  .ge. lrow) go to 550
            if (indr(l) .gt.    0) go to 550

            indr(l)    = j
            lenr(i)    = leni + 1
            go to 560

!           Row i does not have room for the fill-in.
!           Increment ifill(ll) to count how often this has
!           happened to row i.  Also, add m to the row index
!           indc(last) in column j to mark it as a fill-in that is
!           still pending.
!
!           If this is the first pending fill-in for row i,
!           nfill includes the current length of row i
!           (since the whole row has to be moved later).
!
!           If this is the first pending fill-in for column j,
!           jfill(lu) records the current length of column j
!           (to shorten the search for pending fill-ins later).

  550       if (ifill(ll) .eq. 0) nfill     = nfill + leni + nspare
            if (jfill(lu) .eq. 0) jfill(lu) = lenj
            nfill      = nfill     + 1
            ifill(ll)  = ifill(ll) + 1
            indc(last) = m + i
  560    continue

         if ( atend ) lcol = last

!        End loop for column  j.  Store its final length.

  590    lenc(j) = lenj
  600 continue

!     Successful completion.

      lfirst = 0
      return

!     Interruption.  We have to come back in after the
!     column file is compressed.  Give lfirst a new value.
!     lu and nfill will retain their current values.

  900 lfirst = lr
      return

      end ! subroutine lu1gau

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

      implicit           double precision (a-h,o-z)
      logical            TCP
      double precision   Ltol      , a(lena)
      integer            indc(lena), indr(lena), ip(m)   , iq(n)   ,
     &                   lenc(n)   , lenr(m)   , iploc(n), iqloc(m)
      integer            locc(n)   , locr(m)

!     ------------------------------------------------------------------
!     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 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 overflow unless A is incredibly
!                  large and dense.
!
!     10 Sep 2000  TCP, aijtol added for Threshold Complete Pivoting.
!
!     Systems Optimization Laboratory, Stanford University.
!     ------------------------------------------------------------------

      double precision   abest, lbest
      double precision   zero,           one,          gamma
      parameter        ( zero = 0.0d+0,  one = 1.0d+0, gamma  = 2.0d+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

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

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

         do 180 lq = lq1, lq2
            ncol   = ncol + 1
            j      = iq(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 .lt. aijtol) go to 180 ! Nothing in whole column
            end if

            do 160 lc = lc1, lc2
               i      = indc(lc)
               len1   = lenr(i) - 1
!              merit  = nz1 * len1
!              if (merit .gt. mbest) go to 160
               if (len1  .gt. kbest) go to 160

!              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 .eq. 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 .ge. 1.
!                 Bail out if cmax will be too big.

                  aij    = abs( a(lc) )
                  if ( TCP ) then ! Absolute test for Complete Pivoting
                     if (aij         .lt.  aijtol) go to 160
                  else !!! TPP
                     if (aij * Ltol  .lt.  amax  ) go to 160
                  end if
                  cmax   = amax / aij
               end if

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

               merit  = nz1 * len1
               if (merit .eq. 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 .le. gamma  .and.  cmax .le. gamma) then
                     if (abest .ge. aij ) go to 160
                  else
                     if (lbest .le. cmax) go to 160
                  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 .eq. 1) go to 900
  160       continue

!           Finished with that column.

            if (ibest .gt. 0) then
               if (ncol .ge. maxcol) go to 200
            end if
  180    continue

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

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

         do 280 lp = lp1, lp2
            nrow   = nrow + 1
            i      = ip(lp)
            lr1    = locr(i)
            lr2    = lr1 + nz1

            do 260 lr = lr1, lr2
               j      = indr(lr)
               len1   = lenc(j) - 1
!              merit  = nz1 * len1
!              if (merit .gt. mbest) go to 260
               if (len1  .gt. kbest) go to 260

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

               lc1    = locc(j)
               lc2    = lc1 + len1
               amax   = abs( a(lc1) )
               do 220 lc = lc1, lc2
                  if (indc(lc) .eq. i) go to 230
  220          continue

!              Apply the same stability test as above.

  230          aij    = abs( a(lc) )
               if ( TCP ) then   !!! Absolute test for Complete Pivoting
                  if (aij .lt. aijtol) go to 260
               end if

               if (lc .eq. 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  .lt.  amax) go to 260
                  end if
                  cmax   = amax / aij
               end if

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

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

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

                  if (lbest .le. gamma  .and.  cmax .le. gamma) then
                     if (abest .ge. aij ) go to 260
                  else
                     if (lbest .le. cmax) go to 260
                  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 .eq. 1) go to 900
  260       continue

!           Finished with that row.

            if (ibest .gt. 0) then
               if (nrow .ge. maxrow) go to 290
            end if
  280    continue

!        See if it's time to quit.

  290    if (ibest .gt. 0) then
            if (nrow .ge. maxrow  .and.  ncol .ge. maxcol) go to 900
         end if

!        Press on with next nz.

         nz1    = nz
         if (ibest .gt. 0) kbest  = mbest / nz1
      end do

  900 return

      end ! subroutine lu1mar

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1mRP( m    , n     , lena  , maxmn,
     &                   Ltol , maxcol, maxrow,
     &                   ibest, jbest , mbest ,
     &                   a    , indc  , indr  , ip   , iq,
     &                   lenc , lenr  , locc  , locr ,
     &                   iploc, iqloc , Amaxr )

      implicit
     &     none
      integer
     &     m, n, lena, maxmn, maxcol, maxrow, ibest, jbest, mbest
      double precision
     &     Ltol
      integer
     &     indc(lena), indr(lena), ip(m)   , iq(n)   ,
     &     lenc(n)   , lenr(m)   , iploc(n), iqloc(m),
     &     locc(n)   , locr(m)
      double precision
     &     a(lena)   , Amaxr(m)

!     ------------------------------------------------------------------
!     lu1mRP  uses a Markowitz criterion to select a pivot element
!     for the next stage of a sparse LU factorization,
!     subject to a Threshold Rook Pivoting stability criterion (TRP)
!     that bounds the elements of L and U.
!
!     11 Jun 2002: First version of lu1mRP derived from lu1mar.
!     11 Jun 2002: Current version of lu1mRP.
!     ------------------------------------------------------------------

      integer
     &     i, j, kbest, lc, lc1, lc2, len1,
     &     lp, lp1, lp2, lq, lq1, lq2, lr, lr1, lr2,
     &     merit, ncol, nrow, nz, nz1
      double precision
     &     abest, aij, amax, atoli, atolj

      double precision   zero
      parameter        ( zero = 0.0d+0 )

!     ------------------------------------------------------------------
!     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
      ibest  = 0
      kbest  = maxmn + 1
      mbest  = -1
      ncol   = 0
      nrow   = 0
      nz1    = 0

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

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

         do 180 lq = lq1, lq2
            ncol   = ncol + 1
            j      = iq(lq)
            lc1    = locc(j)
            lc2    = lc1 + nz1
            amax   = abs( a(lc1) )
            atolj  = amax / Ltol    ! Min size of pivots in col j

!           Test all aijs in this column.

            do 160 lc = lc1, lc2
               i      = indc(lc)
               len1   = lenr(i) - 1
!              merit  = nz1 * len1
!              if (merit .gt. mbest) go to 160
               if (len1  .gt. kbest) go to 160

!              aij  has a promising merit.
!              Apply the Threshold Rook Pivoting stability test.
!              First we require aij to be sufficiently large
!              compared to other nonzeros in column j.
!              Then  we require aij to be sufficiently large
!              compared to other nonzeros in row    i.

               aij    = abs( a(lc) )
               if (aij         .lt.  atolj   ) go to 160
               if (aij * Ltol  .lt.  Amaxr(i)) go to 160

!              aij  is big enough.

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

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

                  if (abest .ge. aij) go to 160
               end if

!              aij  is the best pivot so far.

               ibest  = i
               jbest  = j
               kbest  = len1
               mbest  = merit
               abest  = aij
               if (nz .eq. 1) go to 900
  160       continue

!           Finished with that column.

            if (ibest .gt. 0) then
               if (ncol .ge. maxcol) go to 200
            end if
  180    continue

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

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

         do 280 lp = lp1, lp2
            nrow   = nrow + 1
            i      = ip(lp)
            lr1    = locr(i)
            lr2    = lr1 + nz1
            atoli  = Amaxr(i) / Ltol   ! Min size of pivots in row i

            do 260 lr = lr1, lr2
               j      = indr(lr)
               len1   = lenc(j) - 1
!              merit  = nz1 * len1
!              if (merit .gt. mbest) go to 260
               if (len1  .gt. kbest) go to 260

!              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) .eq. i) go to 230
               end do

!              Apply the Threshold Rook Pivoting stability test.
!              First we require aij to be sufficiently large
!              compared to other nonzeros in row    i.
!              Then  we require aij to be sufficiently large
!              compared to other nonzeros in column j.

  230          aij    = abs( a(lc) )
               if (aij         .lt.  atoli) go to 260
               if (aij * Ltol  .lt.  amax ) go to 260

!              aij  is big enough.

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

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

                  if (abest .ge. aij ) go to 260
               end if

!              aij  is the best pivot so far.

               ibest  = i
               jbest  = j
               kbest  = len1
               mbest  = merit
               abest  = aij
               if (nz .eq. 1) go to 900
  260       continue

!           Finished with that row.

            if (ibest .gt. 0) then
               if (nrow .ge. maxrow) go to 290
            end if
  280    continue

!        See if it's time to quit.

  290    if (ibest .gt. 0) then
            if (nrow .ge. maxrow  .and.  ncol .ge. maxcol) go to 900
         end if

!        Press on with next nz.

         nz1    = nz
         if (ibest .gt. 0) kbest  = mbest / nz1
      end do

  900 return

      end ! subroutine lu1mRP

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1mCP( m     , n     , lena  , aijtol,
     &                   ibest , jbest , mbest ,
     &                   a     , indc  , indr  ,
     &                   lenc  , lenr  , locc  ,
     &                   Hlen  , Ha    , Hj    )

      integer            m, n, lena, ibest, jbest, mbest, Hlen
      integer            indc(lena), indr(lena),
     &                   lenc(n)   , lenr(m)   , locc(n), Hj(Hlen)
      double precision   aijtol
      double precision   a(lena)   , Ha(Hlen)

!     ------------------------------------------------------------------
!     lu1mCP  uses a Markowitz criterion to select a pivot element
!     for the next stage of a sparse LU factorization,
!     subject to a Threshold Complete Pivoting stability criterion (TCP)
!     that bounds the elements of L and U.
!
!     09 May 2002: First version of lu1mCP.
!                  It searches columns only, using the heap that
!                  holds the largest element in each column.
!     09 May 2002: Current version of lu1mCP.
!     ------------------------------------------------------------------

      integer
     &     j, kheap, lc, lc1, lc2, lenj, maxcol, ncol, nz1
      double precision
     &     abest, aij, amax, cmax, lbest

      double precision   zero,           one,          gamma
      parameter        ( zero = 0.0d+0,  one = 1.0d+0, gamma  = 2.0d+0 )

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

!     ------------------------------------------------------------------
!     Search up to maxcol columns stored at the top of the heap.
!     The very top column helps initialize mbest.
!     ------------------------------------------------------------------
      abest  = zero
      lbest  = zero
      ibest  = 0
      jbest  = Hj(1)               ! Column at the top of the heap
      lenj   = lenc(jbest)
      mbest  = lenj * Hlen         ! Bigger than any possible merit
      maxcol = 40                  ! ??? Big question
      ncol   = 0                   ! No. of columns searched

      do 500 kheap = 1, Hlen

         amax   = Ha(kheap)
         if (amax .lt. aijtol) go to 500

         ncol   = ncol + 1
         j      = Hj(kheap)

!        ---------------------------------------------------------------
!        This column has at least one entry big enough (the top one).
!        Search the column for other possibilities.
!        ---------------------------------------------------------------
         lenj   = lenc(j)
         nz1    = lenj - 1
         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.

         do 160 lc = lc1, lc2
            i      = indc(lc)
            len1   = lenr(i) - 1
            merit  = nz1 * len1
            if (merit .gt. mbest) go to 160

!           aij  has a promising merit.

            if (lc .eq. 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 .ge. 1.
!              Bail out if cmax will be too big.

               aij    = abs( a(lc) )
               if (aij   .lt.  aijtol    ) go to 160
               cmax   = amax / aij
            end if

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

            if (merit .eq. mbest) then

!              Break ties.
!              (Initializing mbest "too big" 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 .le. gamma  .and.  cmax .le. gamma) then
                  if (abest .ge. aij ) go to 160
               else
                  if (lbest .le. cmax) go to 160
               end if
            end if

!           aij  is the best pivot so far.

            ibest  = i
            jbest  = j
            mbest  = merit
            abest  = aij
            lbest  = cmax
            if (merit .eq. 0) go to 900 ! Col or row of length 1
  160    continue

         if (ncol .ge. maxcol) go to 900
  500 continue

  900 return

      end ! subroutine lu1mCP

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1mSP( m    , n     , lena  , maxmn,
     &                   Ltol , maxcol,
     &                   ibest, jbest , mbest ,
     &                   a    , indc  , iq    , locc , iqloc )

      implicit
     &     none
      integer
     &     m, n, lena, maxmn, maxcol, ibest, jbest, mbest
      double precision
     &     Ltol
      integer
     &     indc(lena), iq(n), iqloc(m), locc(n)
      double precision
     &     a(lena)

!     ------------------------------------------------------------------
!     lu1mSP  is intended for symmetric matrices that are either
!     definite or quasi-definite.
!     lu1mSP  uses a Markowitz criterion to select a pivot element for
!     the next stage of a sparse LU factorization of a symmetric matrix,
!     subject to a Threshold Symmetric Pivoting stability criterion
!     (TSP) restricted to diagonal elements to preserve symmetry.
!     This bounds the elements of L and U and should have rank-revealing
!     properties analogous to Threshold Rook Pivoting for unsymmetric
!     matrices.
!
!     14 Dec 2002: First version of lu1mSP derived from lu1mRP.
!                  There is no safeguard to ensure that A is symmetric.
!     14 Dec 2002: Current version of lu1mSP.
!     ------------------------------------------------------------------

      integer
     &     i, j, kbest, lc, lc1, lc2,
     &     lq, lq1, lq2, merit, ncol, nz, nz1
      double precision
     &     abest, aij, amax, atolj

      double precision   zero
      parameter        ( zero = 0.0d+0 )

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

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

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

         do 180 lq = lq1, lq2
            ncol   = ncol + 1
            j      = iq(lq)
            lc1    = locc(j)
            lc2    = lc1 + nz1
            amax   = abs( a(lc1) )
            atolj  = amax / Ltol    ! Min size of pivots in col j

!           Test all aijs in this column.
!           Ignore everything except the diagonal.

            do 160 lc = lc1, lc2
               i      = indc(lc)
               if (i .ne. j) go to 160     ! Skip off-diagonals.
!              merit  = nz1 * nz1
!              if (merit .gt. mbest) go to 160
               if (nz1   .gt. kbest) go to 160

!              aij  has a promising merit.
!              Apply the Threshold Partial Pivoting stability test
!              (which is equivalent to Threshold Rook Pivoting for
!              symmetric matrices).
!              We require aij to be sufficiently large
!              compared to other nonzeros in column j.

               aij    = abs( a(lc) )
               if (aij .lt. atolj  ) go to 160

!              aij  is big enough.

               merit  = nz1 * nz1
               if (merit .eq. mbest) then

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

                  if (abest .ge. aij) go to 160
               end if

!              aij  is the best pivot so far.

               ibest  = i
               jbest  = j
               kbest  = nz1
               mbest  = merit
               abest  = aij
               if (nz .eq. 1) go to 900
  160       continue

!           Finished with that column.

            if (ibest .gt. 0) then
               if (ncol .ge. maxcol) go to 200
            end if
  180    continue

!        See if it's time to quit.

  200    if (ibest .gt. 0) then
            if (ncol .ge. maxcol) go to 900
         end if

!        Press on with next nz.

         nz1    = nz
         if (ibest .gt. 0) kbest  = mbest / nz1
      end do

  900 return

      end ! subroutine lu1mSP

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1pen( m     , melim , ncold , nspare, ilast,
     &                   lpivc1, lpivc2, lpivr1, lpivr2, lrow ,
     &                   lenc  , lenr  , locc  , locr  ,
     &                   indc  , indr  , ifill , jfill )

      integer            indc(*)     , indr(*)     , lenc(*), lenr(*),
     &                   ifill(melim), jfill(ncold)
      integer            locc(*)     , locr(*)

!     ------------------------------------------------------------------
!     lu1pen deals with pending fill-in in the row file.
!     ifill(ll) says if a row involved in the new column of L
!               has to be updated.  If positive, it is the total
!               length of the final updated row.
!     jfill(lu) says if a column involved in the new row of U
!               contains any pending fill-ins.  If positive, it points
!               to the first fill-in in the column that has yet to be
!               added to the row file.
!
!     16 Apr 1989: First version of lu1pen.
!     23 Mar 2001: ilast used and updated.
!     30 Mar 2011: Loop 620 revised following advice from Ralf Östermark,
!                  School of Business and Economics at Åbo Akademi University
!                  after consulting Martti Louhivuori at the
!                  Centre of Scientific Computing in Helsinki.
!                  The previous version produced a core dump on the Cray XT.
!     ------------------------------------------------------------------

         ll     = 0

         do 650 lc  = lpivc1, lpivc2
            ll      = ll + 1
            if (ifill(ll) .eq. 0) go to 650

            ! Another row has pending fill.
            ! First, add some spare space at the end
            ! of the current last row.
            ! 30 Mar 2011: Better not to alter loop limits inside the loop.
            !              Some compilers might not freeze the limits
            !              before entering the loop (even though the
            !              Fortran standards say they should).

            do 620  l  = lrow + 1, lrow + nspare
            !  lrow    = l
               indr(l) = 0
  620       continue
            lrow   = lrow + nspare

            ! Now move row i to the end of the row file.

            i       = indc(lc)
            ilast   = i
            lr1     = locr(i)
            lr2     = lr1 + lenr(i) - 1
            locr(i) = lrow + 1

            do 630 lr = lr1, lr2
               lrow       = lrow + 1
               indr(lrow) = indr(lr)
               indr(lr)   = 0
  630       continue

            lrow    = lrow + ifill(ll)
  650    continue

!        Scan all columns of  D  and insert the pending fill-in
!        into the row file.

         lu     = 1

         do 680 lr = lpivr1, lpivr2
            lu     = lu + 1
            if (jfill(lu) .eq. 0) go to 680
            j      = indr(lr)
            lc1    = locc(j) + jfill(lu) - 1
            lc2    = locc(j) + lenc(j)   - 1

            do 670 lc = lc1, lc2
               i      = indc(lc) - m
               if (i .gt. 0) then
                  indc(lc)   = i
                  last       = locr(i) + lenr(i)
                  indr(last) = j
                  lenr(i)    = lenr(i) + 1
               end if
  670       continue
  680    continue

      end ! subroutine lu1pen

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1mxc( k1, k2, iq,
     &                   a, indc, lenc, locc )

      implicit           none
      integer            k1, k2
      integer            iq(k2), indc(*), lenc(*), locc(*)
      double precision   a(*)

!     ------------------------------------------------------------------
!     lu1mxc  moves the largest element in each of columns iq(k1:k2)
!     to the top of its column.
!     If k1 > k2, nothing happens.
!
!     06 May 2002: (and earlier)
!                  All columns k1:k2 must have one or more elements.
!     07 May 2002: Allow for empty columns.  The heap routines need to
!                  find 0.0 as the "largest element".
!     ------------------------------------------------------------------

      integer             i, j, k, l, lc, lc1, lc2, lenj
      double precision    amax
      double precision    zero
      parameter         ( zero = 0.0d+0 )

      do k = k1, k2
         j      = iq(k)
         lc1    = locc(j)
         lenj   = lenc(j)

         if (lenj .eq. 0) then
            a(lc1) = zero
         else

            ! The next 10 lines are equivalent to
            ! l      = idamax( lenc(j), a(lc1), 1 )  +  lc1 - 1
            ! >>>>>>>>
            lc2    = lc1 + lenc(j) - 1
            amax   = abs( a(lc1) )
            l      = lc1

            do lc = lc1+1, lc2
               if (amax .lt. abs( a(lc) )) then
                  amax   =  abs( a(lc) )
                  l      =  lc
               end if
            end do
            ! >>>>>>>>

            if (l .gt. lc1) then
               amax      = a(l)
               a(l)      = a(lc1)
               a(lc1)    = amax
               i         = indc(l)
               indc(l)   = indc(lc1)
               indc(lc1) = i
            end if
         end if
      end do

      end ! subroutine lu1mxc

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1mxr( mark, k1, k2, m, n, lena,
     &                   a, indc, lenc, locc, indr, lenr, locr,
     &                   p, markc, markr, Amaxr )

      implicit           none
      integer            mark, k1, k2, m, n, lena
      integer            indc(lena), lenc(n), locc(n),
     &                   indr(lena), lenr(m), locr(m), p(k2)
      integer            markc(n), markr(m)
      double precision   a(lena), Amaxr(m)

!------------------------------------------------------------------
! lu1mxr  finds the largest element in each of rows i = p(k1:k2)
! and stores it in each Amaxr(i).
! The nonzeros are stored column-wise in (a,indc,lenc,locc)
! and their structure is     row-wise in (  indr,lenr,locr).
!
! 11 Jun 2002: First version of lu1mxr.
!              Allow for empty columns.
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent.
! 03 Apr 2013: Recoded to improve efficiency.  Need new arrays
!              markc(n), markr(m) and local array cols(n).
!
!              First call:  mark = 0, k1 = 1, k2 = m.
!              Initialize all of markc(n), markr(m), Amaxr(m).
!              Columns are searched only once.
!              cols(n) is not used.
!
!              Later: mark := mark + 1 (greater than for previous call).
!              Cols involved in rows p(k1:k2) are searched only once.
!              cols(n) is local storage.
!              markc(:), markr(:) are marked (= mark) in some places.
!              For next call with new mark,
!              all of markc, markr will initially appear unmarked.
! 06 Jun 2013: Reverted to f77 for lusol.f, trusting that the
!              automatic array cols(n) is fine for f90 compilers.
!------------------------------------------------------------------

      integer             cols(n)
      integer             i, j, k, lc, lc1, lc2, lr, lr1, lr2, ncol
      double precision    zero
      parameter         ( zero = 0.0d+00 )

      if (mark == 0) then       ! First call: Find Amaxr(1:m) for original A.
         do i = 1, m
            markr(i) = 0
            Amaxr(i) = zero
         end do

         do j = 1, n
            markc(j) = 0
         end do

         do j = 1, n
            lc1   = locc(j)
            lc2   = lc1 + lenc(j) - 1
            do lc = lc1, lc2
               i  = indc(lc)
               Amaxr(i) = max( Amaxr(i), abs(a(lc)) )
            end do
         end do

      else                      ! Later calls: Find Amaxr(i) for rows i = p(k1:k2).

         ncol = 0
         do k = k1, k2          ! Search rows to find which cols are involved.
            i        = p(k)
            markr(i) = mark     ! Mark this row
            Amaxr(i) = zero    
            lr1   = locr(i)
            lr2   = lr1 + lenr(i) - 1
            do lr = lr1, lr2    ! Mark all unmarked cols in this row.
               j  = indr(lr)    ! Build up a list of which ones they are.
               if (markc(j) /= mark) then
                  markc(j)  = mark
                  ncol      = ncol + 1
                  cols(ncol)= j
               end if
            end do
         end do

         do k = 1, ncol         ! Search involved columns.
            j     = cols(k)
            lc1   = locc(j)
            lc2   = lc1 + lenc(j) - 1
            do lc = lc1, lc2
               i  = indc(lc)
               if (markr(i) == mark) then
                  Amaxr(i)  = max( Amaxr(i), abs(a(lc)) )
               end if
            end do
         end do
      end if

      end ! subroutine lu1mxr

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!      subroutine lu1mxr( k1, k2, ip, Amaxr,
!     &                   a, indc, lenc, locc, indr, lenr, locr )
!
!      implicit           none
!      integer            k1, k2
!      integer            ip(k2),
!     &                   indc(*), lenc(*), locc(*),
!     &                   indr(*), lenr(*), locr(*)
!      double precision   a(*), Amaxr(*)
!
!!     ------------------------------------------------------------------
!!     lu1mxr  finds the largest element in each of row ip(k1:k2)
!!     and stores it in Amaxr(*).  The nonzeros are stored column-wise
!!     in (a,indc,lenc,locc) and their structure is row-wise
!!     in (  indr,lenr,locr).
!!     If k1 > k2, nothing happens.
!!
!!     11 Jun 2002: First version of lu1mxr.
!!                  Allow for empty columns.
!!     19 Dec 2004: Declare Amaxr(*) instead of Amaxr(k2)
!!                  to stop grumbles from the NAG compiler.
!!                  (Mentioned by Mick Pont.)
!!     ------------------------------------------------------------------
!
!      integer             i, j, k, lc, lc1, lc2, lr, lr1, lr2
!      double precision    amax
!      double precision    zero
!      parameter         ( zero = 0.0d+0 )
!
!      do k = k1, k2
!         amax   = zero
!         i      = ip(k)
!
!         ! Find largest element in row i.
!
!         lr1    = locr(i)
!         lr2    = lr1 + lenr(i) - 1
!
!         do lr = lr1, lr2
!            j      = indr(lr)
!
!            ! Find where  aij  is in column  j.
!
!            lc1    = locc(j)
!            lc2    = lc1 + lenc(j) - 1
!            do lc = lc1, lc2
!               if (indc(lc) .eq. i) go to 230
!            end do
!
!  230       amax   = max( amax, abs( a(lc) ) )
!         end do
!
!         Amaxr(i) = amax
!      end do
!
!      end ! subroutine lu1mxr

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1or1( m, n, nelem, lena, small,
     &                   a, indc, indr, lenc, lenr,
     &                   Amax, numnz, lerr, inform )

      implicit           double precision (a-h,o-z)
      double precision   a(lena)
      integer            indc(lena), indr(lena)
      integer            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.
!     ------------------------------------------------------------------
      do 10 i = 1, m
         lenr(i) = 0
   10 continue

      do 20 j = 1, n
         lenc(j) = 0
   20 continue

      Amax   = 0.0d+0
      numnz  = nelem
      l      = nelem + 1

      do 100 ldummy = 1, nelem
         l      = l - 1
         if (abs( a(l) ) .gt. small) then
            i      = indc(l)
            j      = indr(l)
            Amax   = max( Amax, abs( a(l) ) )
            if (i .lt. 1  .or.  i .gt. m) go to 910
            if (j .lt. 1  .or.  j .gt. 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
  100 continue

      inform = 0
      return

  910 lerr   = l
      inform = 1
      return

      end ! subroutine lu1or1

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1or2( n, numa, lena, a, inum, jnum, len, loc )
      integer            inum(lena), jnum(lena), len(n)
      integer            loc(n)
      double precision   a(lena), ace, acep

!     ------------------------------------------------------------------
!     lu1or2  sorts a list of matrix elements  a(i,j)  into column
!     order, given  numa  entries  a(i,j),  i,  j  in the parallel
!     arrays  a, inum, jnum  respectively.  The matrix is assumed
!     to have  n  columns and an arbitrary number of rows.
!
!     On entry,  len(*)  must contain the length of each column.
!
!     On exit,  a(*) and inum(*)  are sorted,  jnum(*) = 0,  and
!     loc(j)  points to the start of column j.
!
!     lu1or2  is derived from  mc20ad,  a routine in the Harwell
!     Subroutine Library, author J. K. Reid.
!
!     xx Feb 1985: Original version.
!     17 Oct 2000: a, inum, jnum now have size lena to allow nelem = 0.
!     ------------------------------------------------------------------

!     Set  loc(j)  to point to the beginning of column  j.

      l = 1
      do 150 j  = 1, n
         loc(j) = l
         l      = l + len(j)
  150 continue

!     Sort the elements into column order.
!     The algorithm is an in-place sort and is of order  numa.

      do 230 i = 1, numa
!        Establish the current entry.
         jce     = jnum(i)
         if (jce .eq. 0) go to 230
         ace     = a(i)
         ice     = inum(i)
         jnum(i) = 0

!        Chain from current entry.

         do 200 j = 1, numa

!           The current entry is not in the correct position.
!           Determine where to store it.

            l        = loc(jce)
            loc(jce) = loc(jce) + 1

!           Save the contents of that location.

            acep = a(l)
            icep = inum(l)
            jcep = jnum(l)

!           Store current entry.

            a(l)    = ace
            inum(l) = ice
            jnum(l) = 0

!           If next current entry needs to be processed,
!           copy it into current entry.

            if (jcep .eq. 0) go to 230
            ace = acep
            ice = icep
            jce = jcep
  200    continue
  230 continue

!     Reset loc(j) to point to the start of column j.

      ja = 1
      do 250 j  = 1, n
         jb     = loc(j)
         loc(j) = ja
         ja     = jb
  250 continue

      end ! subroutine lu1or2

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1or3( m, n, lena,
     &                   indc, lenc, locc, iw,
     &                   lerr, inform )

      integer            indc(lena), lenc(n), iw(m)
      integer            locc(n)

!     ------------------------------------------------------------------
!     lu1or3  looks for duplicate elements in an  m by n  matrix  A
!     defined by the column list  indc, lenc, locc.
!     iw  is used as a work vector of length  m.
!
!     xx Feb 1985: Original version.
!     17 Oct 2000: indc, indr now have size lena to allow nelem = 0.
!     ------------------------------------------------------------------

      do 100 i = 1, m
         iw(i) = 0
  100 continue

      do 200 j = 1, n
         if (lenc(j) .gt. 0) then
            l1    = locc(j)
            l2    = l1 + lenc(j) - 1

            do 150 l = l1, l2
               i     = indc(l)
               if (iw(i) .eq. j) go to 910
               iw(i) = j
  150       continue
         end if
  200 continue

      inform = 0
      return

  910 lerr   = l
      inform = 1
      return

      end ! subroutine lu1or3

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1or4( m, n, nelem, lena,
     &                   indc, indr, lenc, lenr, locc, locr )

      integer            indc(lena), indr(lena), lenc(n), lenr(m)
      integer            locc(n), locr(m)

!     ------------------------------------------------------------------
!     lu1or4     constructs a row list  indr, locr
!     from a corresponding column list  indc, locc,
!     given the lengths of both columns and rows in  lenc, lenr.
!
!     xx Feb 1985: Original version.
!     17 Oct 2000: indc, indr now have size lena to allow nelem = 0.
!     ------------------------------------------------------------------

!     Initialize  locr(i)  to point just beyond where the
!     last component of row  i  will be stored.

      l      = 1
      do 10 i = 1, m
         l       = l + lenr(i)
         locr(i) = l
   10 continue

!     By processing the columns backwards and decreasing  locr(i)
!     each time it is accessed, it will end up pointing to the
!     beginning of row  i  as required.

      l2     = nelem
      j      = n + 1

      do 40 jdummy = 1, n
         j  = j - 1
         if (lenc(j) .gt. 0) then
            l1 = locc(j)

            do 30 l = l1, l2
               i        = indc(l)
               lr       = locr(i) - 1
               locr(i)  = lr
               indr(lr) = j
   30       continue

            l2     = l1 - 1
         end if
   40 continue

      end ! subroutine lu1or4

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1pq1( m, n, len, iperm, loc, inv, num )

      integer            len(m), iperm(m), loc(n), inv(m), num(n)

!     ------------------------------------------------------------------
!     lu1pq1  constructs a permutation  iperm  from the array  len.
!
!     On entry:
!     len(i)  holds the number of nonzeros in the i-th row (say)
!             of an m by n matrix.
!     num(*)  can be anything (workspace).
!
!     On exit:
!     iperm   contains a list of row numbers in the order
!             rows of length 0,  rows of length 1,..., rows of length n.
!     loc(nz) points to the first row containing  nz  nonzeros,
!             nz = 1, n.
!     inv(i)  points to the position of row i within iperm(*).
!     ------------------------------------------------------------------

!     Count the number of rows of each length.

      nzero  = 0
      do 10  nz  = 1, n
         num(nz) = 0
         loc(nz) = 0
   10 continue

      do 20  i   = 1, m
         nz      = len(i)
         if (nz .eq. 0) then
            nzero   = nzero   + 1
         else
            num(nz) = num(nz) + 1
         end if
   20 continue

!     Set starting locations for each length.

      l      = nzero + 1
      do 60  nz  = 1, n
         loc(nz) = l
         l       = l + num(nz)
         num(nz) = 0
   60 continue

!     Form the list.

      nzero  = 0
      do 100  i   = 1, m
         nz       = len(i)
         if (nz .eq. 0) then
            nzero    = nzero + 1
            iperm(nzero) = i
         else
            l        = loc(nz) + num(nz)
            iperm(l) = i
            num(nz)  = num(nz) + 1
         end if
  100 continue

!     Define the inverse of iperm.

      do 120 l  = 1, m
         i      = iperm(l)
         inv(i) = l
  120 continue

      end ! subroutine lu1pq1

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1pq2( nzpiv, nzchng,
     &                   indr , lenold, lennew, iqloc, iq, iqinv )

      integer            indr(nzpiv), lenold(nzpiv), lennew(*),
     &                   iqloc(*)   , iq(*)        , iqinv(*)

!     ===============================================================
!     lu1pq2 frees the space occupied by the pivot row,
!     and updates the column permutation iq.
!
!     Also used to free the pivot column and update the row perm ip.
!
!     nzpiv   (input)    is the length of the pivot row (or column).
!     nzchng  (output)   is the net change in total nonzeros.
!
!     14 Apr 1989  First version.
!     ===============================================================

      nzchng = 0

      do 200  lr  = 1, nzpiv
         j        = indr(lr)
         indr(lr) = 0
         nz       = lenold(lr)
         nznew    = lennew(j)

         if (nz .ne. nznew) then
            l        = iqinv(j)
            nzchng   = nzchng + (nznew - nz)

!           l above is the position of column j in iq  (so j = iq(l)).

            if (nz .lt. nznew) then

!              Column  j  has to move towards the end of  iq.

  110          next        = nz + 1
               lnew        = iqloc(next) - 1
               if (lnew .ne. l) then
                  jnew        = iq(lnew)
                  iq(l)       = jnew
                  iqinv(jnew) = l
               end if
               l           = lnew
               iqloc(next) = lnew
               nz          = next
               if (nz .lt. nznew) go to 110
            else

!              Column  j  has to move towards the front of  iq.

  120          lnew        = iqloc(nz)
               if (lnew .ne. l) then
                  jnew        = iq(lnew)
                  iq(l)       = jnew
                  iqinv(jnew) = l
               end if
               l           = lnew
               iqloc(nz)   = lnew + 1
               nz          = nz   - 1
               if (nz .gt. nznew) go to 120
            end if

            iq(lnew) = j
            iqinv(j) = lnew
         end if
  200 continue

      end ! subroutine lu1pq2

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1pq3( n, len, iperm, iw, nrank )

      integer            len(n), iperm(n)
      integer            iw(n)

!     ------------------------------------------------------------------
!     lu1pq3  looks at the permutation  iperm(*)  and moves any entries
!     to the end whose corresponding length  len(*)  is zero.
!
!     09 Feb 1994: Added work array iw(*) to improve efficiency.
!     ------------------------------------------------------------------

      nrank  = 0
      nzero  = 0

      do 10 k = 1, n
         i    = iperm(k)

         if (len(i) .eq. 0) then
            nzero        = nzero + 1
            iw(nzero)    = i
         else
            nrank        = nrank + 1
            iperm(nrank) = i
         end if
   10 continue

      do 20 k = 1, nzero
         iperm(nrank + k) = iw(k)
   20 continue

      end ! subroutine lu1pq3

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1rec( n, reals, luparm, ltop,
     &                   lena, a, ind, len, loc )

      logical            reals
      integer            luparm(30), ltop
      double precision   a(lena)
      integer            ind(lena), len(n)
      integer            loc(n)

!     ------------------------------------------------------------------
!     00 Jun 1983: Original version of lu1rec followed John Reid's
!                  compression routine in LA05.  It recovered
!                  space in ind(*) and optionally a(*)
!                  by eliminating entries with ind(l) = 0.
!                  The elements of ind(*) could not be negative.
!                  If len(i) was positive, entry i contained
!                  that many elements, starting at  loc(i).
!                  Otherwise, entry i was eliminated.
!
!     23 Mar 2001: Realised we could have len(i) = 0 in rare cases!
!                  (Mostly during TCP when the pivot row contains
!                  a column of length 1 that couldn't be a pivot.)
!                  Revised storage scheme to
!                     keep        entries with       ind(l) >  0,
!                     squeeze out entries with -n <= ind(l) <= 0,
!                  and to allow len(i) = 0.
!                  Empty items are moved to the end of the compressed
!                  ind(*) and/or a(*) arrays are given one empty space.
!                  Items with len(i) < 0 are still eliminated.
!
!     27 Mar 2001: Decided to use only ind(l) > 0 and = 0 in lu1fad.
!                  Still have to keep entries with len(i) = 0.
!
!     On exit:
!     ltop         is the length of useful entries in ind(*), a(*).
!     ind(ltop+1)  is "i" such that len(i), loc(i) belong to the last
!                  item in ind(*), a(*).
!     ------------------------------------------------------------------

      nempty = 0

      do 10 i = 1, n
         leni = len(i)
         if (leni .gt. 0) then
            l      = loc(i) + leni - 1
            len(i) = ind(l)
            ind(l) = - (n + i)
         else if (leni .eq. 0) then
            nempty = nempty + 1
         end if
   10 continue

      k      = 0
      klast  = 0    ! Previous k
      ilast  = 0    ! Last entry moved.

      do 20 l = 1, ltop
         i    = ind(l)
         if (i .gt. 0) then
            k      = k + 1
            ind(k) = i
            if (reals) a(k) = a(l)

         else if (i .lt. -n) then

!           This is the end of entry  i.

            i      = - (i + n)
            ilast  = i
            k      = k + 1
            ind(k) = len(i)
            if (reals) a(k) = a(l)
            loc(i) = klast + 1
            len(i) = k     - klast
            klast  = k
         end if
   20 continue

!     Move any empty items to the end, adding 1 free entry for each.

      if (nempty .gt. 0) then
         do i = 1, n
            if (len(i) .eq. 0) then
               k      = k + 1
               loc(i) = k
               ind(k) = 0
               ilast  = i
            end if
         end do
      end if

      nout   = luparm(1)
      lprint = luparm(2)
      if (lprint .ge. 50) write(nout, 1000) ltop, k, reals, nempty
      luparm(26) = luparm(26) + 1  ! ncp

!     Return ilast in ind(ltop + 1).

      ltop        = k
      ind(ltop+1) = ilast
      return

 1000 format(' lu1rec.  File compressed from', i10, '   to', i10, l3,
     &       '  nempty =', i8)

      end ! subroutine lu1rec

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1slk( m, n, lena, iq, iqloc, a, locc, w )

      implicit           none
      integer            m, n, lena
      integer            iq(n), iqloc(m), locc(n)
      double precision   a(lena), w(n)

!     ------------------------------------------------------------------
!     lu1slk  sets w(j) > 0 if column j is a unit vector.
!
!     21 Nov 2000: First version.  lu1fad needs it for TCP.
!                  Note that w(*) is nominally an integer array,
!                  but the only spare space is the double array w(*).
!     ------------------------------------------------------------------

      integer            j, lc1, lq, lq1, lq2

      do j = 1, n
         w(j) = 0.0d+0
      end do

      lq1    = iqloc(1)
      lq2    = n
      if (m .gt. 1) lq2 = iqloc(2) - 1

      do lq = lq1, lq2
         j      = iq(lq)
         lc1    = locc(j)
         if (abs( a(lc1) ) .eq. 1.0d+0) then
            w(j) = 1.0d+0
         end if
      end do

      end ! subroutine lu1slk

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1ful( m     , n    , lena , lenD , lu1 , TPP,
     &                   mleft , nleft, nrank, nrowu,
     &                   lenL  , lenU , nsing,
     &                   keepLU, small,
     &                   a     , d    , indc , indr , ip  , iq,
     &                   lenc  , lenr , locc , ipinv, ipvt )

      implicit           double precision (a-h,o-z)
      logical            TPP       , keepLU
      double precision   a(lena)   , d(lenD)
      integer            indc(lena), indr(lena), ip(m)   , iq(n),
     &                   lenc(n)   , lenr(m)   , ipinv(m)
      integer            locc(n)   , ipvt(m)

!     ------------------------------------------------------------------
!     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 ip.
!     ------------------------------------------------------------------

      parameter        ( zero = 0.0d+0 )


!     ------------------------------------------------------------------
!     If lu1pq3 moved any empty rows, reset ipinv = inverse of ip.
!     ------------------------------------------------------------------
      if (nrank .lt. m) then
         do 100 l    = 1, m
            i        = ip(l)
            ipinv(i) = l
  100    continue
      end if

!     ------------------------------------------------------------------
!     Copy the remaining matrix into the dense matrix D.
!     ------------------------------------------------------------------
!     call dload ( lenD, zero, d, 1 )
      do j = 1, lenD
         d(j) = zero
      end do
      ipbase = nrowu - 1
      ldbase = 1 - nrowu

      do 200 lq = nrowu, n
         j      = iq(lq)
         lc1    = locc(j)
         lc2    = lc1 + lenc(j) - 1

         do 150 lc = lc1, lc2
            i      = indc(lc)
            ld     = ldbase + ipinv(i)
            d(ld)  = a(lc)
  150    continue

         ldbase = ldbase + mleft
  200 continue

!     ------------------------------------------------------------------
!     Call our favorite dense LU factorizer.
!     ------------------------------------------------------------------
      if ( TPP ) then
         call lu1DPP( d, mleft, mleft, nleft, small, nsing,
     &                ipvt, iq(nrowu) )
      else
         call lu1DCP( d, mleft, mleft, nleft, small, nsing,
     &                ipvt, iq(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 ip.
!     lkk points to the diagonal of U.
!     ------------------------------------------------------------------
      call dcopy ( lenD, d, 1, a, 1 )

      ldiagU = lena   - n
      lkk    = 1
      lkn    = lenD  - mleft + 1
      lu     = lu1

      do 450  k = 1, min( mleft, nleft )
         l1     = ipbase + k
         l2     = ipbase + ipvt(k)
         if (l1 .ne. l2) then
            i      = ip(l1)
            ip(l1) = ip(l2)
            ip(l2) = i
         end if
         ibest  = ip( l1 )
         jbest  = iq( l1 )

         if ( keepLU ) then
!           ===========================================================
!           Pack the next column of L.
!           ===========================================================
            la     = lkk
            ll     = lu
            nrowd  = 1

            do 410  i = k + 1, mleft
               la     = la + 1
               ai     = a(la)
               if (abs( ai ) .gt. small) then
                  nrowd    = nrowd + 1
                  ll       = ll    - 1
                  a(ll)    = ai
                  indc(ll) = ip( ipbase + i )
                  indr(ll) = ibest
               end if
  410       continue

!           ===========================================================
!           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 420  j = nleft, k, -1
               la     = la - mleft
               aj     = a(la)
               if (abs( aj ) .gt. small  .or.  j .eq. k) then
                  ncold    = ncold + 1
                  lu       = lu    - 1
                  a(lu)    = aj
                  indr(lu) = iq( ipbase + j )
               end if
  420       continue

            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
  450 continue

      end ! subroutine lu1ful

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1DPP( a, lda, m, n, small, nsing,
     &                   ipvt, iq )

      implicit           none
      integer            lda, m, n, nsing
      integer            ipvt(m), iq(n)
      double precision   a(lda,n), small

!     ------------------------------------------------------------------
!     lu1DPP factors a dense m x n matrix A by Gaussian elimination,
!     using row interchanges for stability, as in dgefa from LINPACK.
!     This version also uses column interchanges if all elements in a
!     pivot column are smaller than (or equal to) "small".  Such columns
!     are changed to zero and permuted to the right-hand end.
!
!     As in LINPACK, ipvt(*) keeps track of pivot rows.
!     Rows of U are interchanged, but we don't have to physically
!     permute rows of L.  In contrast, column interchanges are applied
!     directly to the columns of both L and U, and to the column
!     permutation vector iq(*).
!
!     02 May 1989: First version derived from dgefa
!                  in LINPACK (version dated 08/14/78).
!     05 Feb 1994: Generalized to treat rectangular matrices
!                  and use column interchanges when necessary.
!                  ipvt is retained, but column permutations are applied
!                  directly to iq(*).
!     21 Dec 1994: Bug found via example from Steve Dirkse.
!                  Loop 100 added to set ipvt(*) for singular rows.
!     26 Mar 2006: nsing redefined (see below).
!                  Changed to implicit none.
!     ------------------------------------------------------------------
!
!     On entry:
!
!        a       Array holding the matrix A to be factored.
!
!        lda     The leading dimension of the array  a.
!
!        m       The number of rows    in  A.
!
!        n       The number of columns in  A.
!
!        small   A drop tolerance.  Must be zero or positive.
!
!     On exit:
!
!        a       An upper triangular matrix and the multipliers
!                which were used to obtain it.
!                The factorization can be written  A = L*U  where
!                L  is a product of permutation and unit lower
!                triangular matrices and  U  is upper triangular.
!
!        nsing   Number of singularities detected.
!                26 Mar 2006: nsing redefined to be more meaningful.
!                Users may define rankU = n - nsing and regard
!                U as upper-trapezoidal, with the first rankU columns
!                being triangular and the rest trapezoidal.
!                It would be better to return rankU, but we still
!                return nsing for compatibility (even though lu1fad
!                no longer uses it).
!
!        ipvt    Records the pivot rows.
!
!        iq      A vector to which column interchanges are applied.
!     ------------------------------------------------------------------

      double precision    t
      integer             idamax, i, j, k, kp1, l, last, lencol, rankU
      double precision    zero         ,  one
      parameter         ( zero = 0.0d+0,  one = 1.0d+0 )


      rankU  = 0
      k      = 1
      last   = n

!     ------------------------------------------------------------------
!     Start of elimination loop.
!     ------------------------------------------------------------------
   10 kp1    = k + 1
      lencol = m - k + 1

      ! Find l, the pivot row.

      l       = idamax( lencol, a(k,k), 1 ) + k - 1
      ipvt(k) = l

      if (abs( a(l,k) ) .le. small) then
         !==============================================================
         ! Do column interchange, changing old pivot column to zero.
         ! Reduce "last" and try again with same k.
         !==============================================================
         j        = iq(last)
         iq(last) = iq(k)
         iq(k)    = j

         do i = 1, k - 1
            t         = a(i,last)
            a(i,last) = a(i,k)
            a(i,k)    = t
         end do

         do i = k, m
            t         = a(i,last)
            a(i,last) = zero
            a(i,k)    = t
         end do

         last     = last - 1
         if (k .le. last) go to 10

      else
         rankU  = rankU + 1
         if (k .lt. m) then
            !===========================================================
            !Do row interchange if necessary.
            !===========================================================
            if (l .ne. k) then
               t      = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
            end if

            !===========================================================
            ! Compute multipliers.
            ! Do row elimination with column indexing.
            !===========================================================
            t = - one / a(k,k)
            call dscal ( m-k, t, a(kp1,k), 1 )

            do j = kp1, last
               t    = a(l,j)
               if (l .ne. k) then
                  a(l,j) = a(k,j)
                  a(k,j) = t
               end if
               call daxpy ( m-k, t, a(kp1,k), 1, a(kp1,j), 1 )
            end do

            k = k + 1
            if (k .le. last) go to 10
         end if
      end if

      ! Set ipvt(*) for singular rows.

      do 100 k = last + 1, m
         ipvt(k) = k
  100 continue

      nsing  = n - rankU

      end ! subroutine lu1DPP

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu1DCP( a, lda, m, n, small, nsing,
     &                   ipvt, iq )

      implicit           none
      integer            lda, m, n, nsing
      integer            ipvt(m), iq(n)
      double precision   a(lda,n), small

!     ------------------------------------------------------------------
!     lu1DCP factors a dense m x n matrix A by Gaussian elimination,
!     using Complete Pivoting (row and column interchanges) for
!     stability.
!     This version also uses column interchanges if all elements in a
!     pivot column are smaller than (or equal to) "small".  Such columns
!     are changed to zero and permuted to the right-hand end.
!
!     As in LINPACK's dgefa, ipvt(*) keeps track of pivot rows.
!     Rows of U are interchanged, but we don't have to physically
!     permute rows of L.  In contrast, column interchanges are applied
!     directly to the columns of both L and U, and to the column
!     permutation vector iq(*).
!
!     01 May 2002: First dense Complete Pivoting, derived from lu1DPP.
!     07 May 2002: Another break needed at end of first loop.
!     26 Mar 2006: Cosmetic mods while looking for "nsing" bug when m<n.
!                  nsing redefined (see below).
!                  Changed to implicit none.
!     ------------------------------------------------------------------
!
!     On entry:
!
!        a       Array holding the matrix A to be factored.
!
!        lda     The leading dimension of the array  a.
!
!        m       The number of rows    in  A.
!
!        n       The number of columns in  A.
!
!        small   A drop tolerance.  Must be zero or positive.
!
!     On exit:
!
!        a       An upper triangular matrix and the multipliers
!                which were used to obtain it.
!                The factorization can be written  A = L*U  where
!                L  is a product of permutation and unit lower
!                triangular matrices and  U  is upper triangular.
!
!        nsing   Number of singularities detected.
!                26 Mar 2006: nsing redefined to be more meaningful.
!                Users may define rankU = n - nsing and regard
!                U as upper-trapezoidal, with the first rankU columns
!                being triangular and the rest trapezoidal.
!                It would be better to return rankU, but we still
!                return nsing for compatibility (even though lu1fad
!                no longer uses it).
!
!        ipvt    Records the pivot rows.
!
!        iq      A vector to which column interchanges are applied.
!     ------------------------------------------------------------------

      double precision    aijmax, ajmax, t
      integer             idamax, i, imax, j, jlast, jmax, jnew,
     &                    k, kp1, l, last, lencol, rankU
      double precision    zero         ,  one
      parameter         ( zero = 0.0d+0,  one = 1.0d+0 )


      rankU  = 0
      lencol = m + 1
      last   = n

      !-----------------------------------------------------------------
      ! Start of elimination loop.
      !-----------------------------------------------------------------
      do k = 1, n
         kp1    = k + 1
         lencol = lencol - 1

         ! Find the biggest aij in row imax and column jmax.

         aijmax = zero
         imax   = k
         jmax   = k
         jlast  = last

         do j = k, jlast
   10       l      = idamax( lencol, a(k,j), 1 ) + k - 1
            ajmax  = abs( a(l,j) )

            if (ajmax .le. small) then
               !========================================================
               ! Do column interchange, changing old column to zero.
               ! Reduce  "last"  and try again with same j.
               !========================================================
               jnew     = iq(last)
               iq(last) = iq(j)
               iq(j)    = jnew

               do i = 1, k - 1
                  t         = a(i,last)
                  a(i,last) = a(i,j)
                  a(i,j)    = t
               end do

               do i = k, m
                  t         = a(i,last)
                  a(i,last) = zero
                  a(i,j)    = t
               end do

               last   = last - 1
               if (j .le. last) go to 10 ! repeat
               go to 200                 ! break
            end if

            ! Check if this column has biggest aij so far.

            if (aijmax .lt. ajmax) then
                aijmax  =   ajmax
                imax    =   l
                jmax    =   j
            end if

            if (j .ge. last) go to 200   ! break
         end do

  200    ipvt(k) = imax

         if (jmax .ne. k) then
            !==========================================================
            ! Do column interchange (k and jmax).
            !==========================================================
            jnew     = iq(jmax)
            iq(jmax) = iq(k)
            iq(k)    = jnew

            do i = 1, m
               t         = a(i,jmax)
               a(i,jmax) = a(i,k)
               a(i,k)    = t
            end do
         end if

         if (k .lt. m) then
            !===========================================================
            ! Do row interchange if necessary.
            !===========================================================
            t         = a(imax,k)
            if (imax .ne. k) then
               a(imax,k) = a(k,k)
               a(k,k)    = t
            end if

            !===========================================================
            ! Compute multipliers.
            ! Do row elimination with column indexing.
            !===========================================================
            t      = - one / t
            call dscal ( m-k, t, a(kp1,k), 1 )

            do j = kp1, last
               t         = a(imax,j)
               if (imax .ne. k) then
                  a(imax,j) = a(k,j)
                  a(k,j)    = t
               end if
               call daxpy ( m-k, t, a(kp1,k), 1, a(kp1,j), 1 )
            end do

         else
            go to 500               ! break
         end if

         if (k .ge. last) go to 500 ! break
      end do

      ! Set ipvt(*) for singular rows.

  500 do k = last + 1, m
         ipvt(k) = k
      end do

      nsing  = n - rankU

      end ! subroutine lu1DCP
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
*     File  lusol2.f
*
*     Hbuild   Hchange  Hdelete  Hdown    Hinsert  Hup
*
*     Heap-management routines for LUSOL's lu1fac.
*     May be useful for other applications.
*
* 11 Feb 2002: MATLAB version derived from "Algorithms" by R. Sedgewick.
* 03 Mar 2002: F77    version derived from MATLAB version.
* 07 May 2002: Safeguard input parameters k, N, Nk.
*              We don't want them to be output!
* 19 Dec 2004: Hdelete: Nin is new input parameter for length of Hj, Ha.
* 19 Dec 2004: Current version of lusol2.f.
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!     For LUSOL, the heap structure involves three arrays of length N.
!     N        is the current number of entries in the heap.
!     Ha(1:N)  contains the values that the heap is partially sorting.
!              For LUSOL they are double precision values -- the largest
!              element in each remaining column of the updated matrix.
!              The biggest entry is in Ha(1), the top of the heap.
!     Hj(1:N)  contains column numbers j.
!              Ha(k) is the biggest entry in column j = Hj(k).
!     Hk(1:N)  contains indices within the heap.  It is the
!              inverse of Hj(1:N), so  k = Hk(j)  <=>  j = Hj(k).
!              Column j is entry k in the heap.
!     hops     is the number of heap operations,
!              i.e., the number of times an entry is moved
!              (the number of "hops" up or down the heap).
!     Together, Hj and Hk let us find values inside the heap
!     whenever we want to change one of the values in Ha.
!     For other applications, Ha may need to be some other data type,
!     like the keys that sort routines operate on.
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Hbuild( Ha, Hj, Hk, N, Nk, hops )

      implicit
     &     none
      integer
     &     N, Nk, hops, Hj(N), Hk(Nk)
      double precision
     &     Ha(N)

*     ==================================================================
*     Hbuild initializes the heap by inserting each element of Ha.
*     Input:  Ha, Hj.
*     Output: Ha, Hj, Hk, hops.
*
*     01 May 2002: Use k for new length of heap, not k-1 for old length.
*     05 May 2002: Use kk in call to stop loop variable k being altered.
*                  (Actually Hinsert no longer alters that parameter.)
*     07 May 2002: ftnchek wants us to protect Nk, Ha(k), Hj(k) too.
*     07 May 2002: Current version of Hbuild.
*     ==================================================================

      integer
     &     h, jv, k, kk, Nkk
      double precision
     &     v

      Nkk  = Nk
      hops = 0
      do k = 1, N
         kk    = k
         v     = Ha(k)
         jv    = Hj(k)
         call Hinsert( Ha, Hj, Hk, kk, Nkk, v, jv, h )
         hops  = hops + h
      end do

      end ! subroutine Hbuild

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Hchange( Ha, Hj, Hk, N, Nk, k, v, jv, hops )

      implicit
     &     none
      integer
     &     N, Nk, k, jv, hops, Hj(N), Hk(Nk)
      double precision
     &     v, Ha(N)

*     ==================================================================
*     Hchange changes Ha(k) to v in heap of length N.
*
*     01 May 2002: Need Nk for length of Hk.
*     07 May 2002: Protect input parameters N, Nk, k.
*     07 May 2002: Current version of Hchange.
*     ==================================================================

      integer
     &     kx, Nx, Nkx
      double precision
     &     v1

      Nx     = N
      Nkx    = Nk
      kx     = k
      v1     = Ha(k)
      Ha(k)  = v
      Hj(k)  = jv
      Hk(jv) = k
      if (v1 .lt. v) then
         call Hup   ( Ha, Hj, Hk, Nx, Nkx, kx, hops )
      else
         call Hdown ( Ha, Hj, Hk, Nx, Nkx, kx, hops )
      end if

      end ! subroutine Hchange

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Hdelete( Ha, Hj, Hk, Nin, N, Nk, k, hops )

      implicit
     &     none
      integer
     &     N, Nin, Nk, k, hops, Hj(Nin), Hk(Nk)
      double precision
     &     Ha(Nin)

*     ==================================================================
*     Hdelete deletes Ha(k) from heap of length N.
*
*     03 Apr 2002: Current version of Hdelete.
*     01 May 2002: Need Nk for length of Hk.
*     07 May 2002: Protect input parameters N, Nk, k.
*     19 Dec 2004: Nin is new input parameter for length of Hj, Ha.
*     19 Dec 2004: Current version of Hdelete.
*     ==================================================================

      integer
     &     jv, kx, Nkx, Nx
      double precision
     &     v

      kx    = k
      Nkx   = Nk
      Nx    = N
      v     = Ha(N)
      jv    = Hj(N)
      N     = N - 1
      hops  = 0
      if (k .le. N) then
         call Hchange( Ha, Hj, Hk, Nx, Nkx, kx, v, jv, hops )
      end if

      end ! subroutine Hdelete

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Hdown ( Ha, Hj, Hk, N, Nk, kk, hops )

      implicit
     &     none
      integer
     &     N, Nk, kk, hops, Hj(N), Hk(Nk)
      double precision
     &     Ha(N)

*     ==================================================================
*     Hdown  updates heap by moving down tree from node k.
*
*     01 May 2002: Need Nk for length of Hk.
*     05 May 2002: Change input paramter k to kk to stop k being output.
*     05 May 2002: Current version of Hdown.
*     ==================================================================

      integer
     &     j, jj, jv, k, N2
      double precision
     &     v

      k     = kk
      hops  = 0
      v     = Ha(k)
      jv    = Hj(k)
      N2    = N/2

!     while 1
  100    if (k .gt. N2   ) go to 200   ! break
         hops   = hops + 1
         j      = k+k
         if (j .lt. N) then
            if (Ha(j) .lt. Ha(j+1)) j = j+1
         end if
         if (v .ge. Ha(j)) go to 200   ! break
         Ha(k)  = Ha(j)
         jj     = Hj(j)
         Hj(k)  = jj
         Hk(jj) =  k
         k      =  j
         go to 100
!     end while

  200 Ha(k)  =  v
      Hj(k)  = jv
      Hk(jv) =  k

      end ! subroutine Hdown

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Hinsert( Ha, Hj, Hk, N, Nk, v, jv, hops )

      implicit
     &     none
      integer
     &     N, Nk, jv, hops, Hj(N), Hk(Nk)
      double precision
     &     v, Ha(N)

*     ==================================================================
*     Hinsert inserts (v,jv) into heap of length N-1
*     to make heap of length N.
*
*     03 Apr 2002: First version of Hinsert.
*     01 May 2002: Require N to be final length, not old length.
*                  Need Nk for length of Hk.
*     07 May 2002: Protect input parameters N, Nk.
*     07 May 2002: Current version of Hinsert.
*     ==================================================================

      integer
     &     kk, Nkk, Nnew

      Nnew     = N
      Nkk      = Nk
      kk       = Nnew
      Ha(Nnew) =  v
      Hj(Nnew) = jv
      Hk(jv)   = Nnew
      call Hup   ( Ha, Hj, Hk, Nnew, Nkk, kk, hops )

      end ! subroutine Hinsert

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine Hup   ( Ha, Hj, Hk, N, Nk, kk, hops )

      implicit
     &     none
      integer
     &     N, Nk, kk, hops, Hj(N), Hk(Nk)
      double precision
     &     Ha(N)

*     ==================================================================
*     Hup updates heap by moving up tree from node k.
*
*     01 May 2002: Need Nk for length of Hk.
*     05 May 2002: Change input paramter k to kk to stop k being output.
*     05 May 2002: Current version of Hup.
*     ==================================================================

      integer
     &     j, jv, k, k2
      double precision
     &     v

      k     = kk
      hops  = 0
      v     = Ha(k)
      jv    = Hj(k)
!     while 1
  100    if (k .lt.  2    ) go to 200   ! break
         k2    = k/2
         if (v .lt. Ha(k2)) go to 200   ! break
         hops  = hops + 1
         Ha(k) = Ha(k2)
         j     = Hj(k2)
         Hj(k) =  j
         Hk(j) =  k
         k     = k2
         go to 100
!     end while

  200 Ha(k)  =  v
      Hj(k)  = jv
      Hk(jv) =  k

      end ! subroutine Hup
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!     File  lusol6a.f
!
!     lu6sol   lu6L     lu6Lt     lu6U     Lu6Ut   lu6LD
!     lu6chk
!
! 26 Apr 2002: lu6 routines put into a separate file.
! 15 Dec 2002: lu6sol modularized via lu6L, lu6Lt, lu6U, lu6Ut.
!              lu6LD implemented to allow solves with LDL' or L|D|L'.
! 23 Apr 2004: lu6chk modified.  TRP can judge singularity better
!              by comparing all diagonals to DUmax.
! 27 Jun 2004: lu6chk.  Allow write only if nout .gt. 0 .
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu6sol( mode, m, n, v, w,
     &                   lena, luparm, parmlu,
     &                   a, indc, indr, ip, iq,
     &                   lenc, lenr, locc, locr,
     &                   inform )

      implicit
     &     none
      integer
     &     luparm(30), mode, m, n, lena, inform,
     &     indc(lena), indr(lena), ip(m), iq(n),
     &     lenc(n), lenr(m), locc(n), locr(m)
      double precision
     &     parmlu(30), a(lena), v(m), w(n)

!-----------------------------------------------------------------------
!     lu6sol  uses the factorization  A = L U  as follows:
!
!     mode
!      1    v  solves   L v = v(input).   w  is not touched.
!      2    v  solves   L'v = v(input).   w  is not touched.
!      3    w  solves   U w = v.          v  is not altered.
!      4    v  solves   U'v = w.          w  is destroyed.
!      5    w  solves   A w = v.          v  is altered as in 1.
!      6    v  solves   A'v = w.          w  is destroyed.
!
!     If mode = 3,4,5,6, v and w must not be the same arrays.
!
!     If lu1fac has just been used to factorize a symmetric matrix A
!     (which must be definite or quasi-definite), the factors A = L U
!     may be regarded as A = LDL', where D = diag(U).  In such cases,
!
!     mode
!      7    v  solves   A v = L D L'v = v(input).   w  is not touched.
!      8    v  solves       L |D| L'v = v(input).   w  is not touched.
!
!     ip(*), iq(*)      hold row and column numbers in pivotal order.
!     lenc(k)           is the length of the k-th column of initial L.
!     lenr(i)           is the length of the i-th row of U.
!     locc(*)           is not used.
!     locr(i)           is the start  of the i-th row of U.
!
!     U is assumed to be in upper-trapezoidal form (nrank by n).
!     The first entry for each row is the diagonal element
!     (according to the permutations  ip, iq).  It is stored at
!     location locr(i) in a(*), indr(*).
!
!     On exit, inform = 0 except as follows.
!     If mode = 3,4,5,6 and if U (and hence A) is singular, then
!     inform = 1 if there is a nonzero residual in solving the system
!     involving U.  parmlu(20) returns the norm of the residual.
!
!       July 1987: Early version.
!     09 May 1988: f77 version.
!     27 Apr 2000: Abolished the dreaded "computed go to".
!                  But hard to change other "go to"s to "if then else".
!     15 Dec 2002: lu6L, lu6Lt, lu6U, lu6Ut added to modularize lu6sol.
!-----------------------------------------------------------------------

      if      (mode .eq. 1) then             ! Solve  L v(new) = v.
         call lu6L  (
     &        inform, m, n, v,
     &        lena, luparm, parmlu,
     &        a, indc, indr, lenc )

      else if (mode .eq. 2) then             ! Solve  L'v(new) = v.
         call lu6Lt (
     &        inform, m, n, v,
     &        lena, luparm, parmlu,
     &        a, indc, indr, lenc )

      else if (mode .eq. 3) then             ! Solve  U w = v.
         call lu6U  (
     &        inform, m, n, v, w,
     &        lena, luparm, parmlu,
     &        a, indr, ip, iq, lenr, locr )

      else if (mode .eq. 4) then             ! Solve  U'v = w.
         call lu6Ut (
     &        inform, m, n, v, w,
     &        lena, luparm, parmlu,
     &        a, indr, ip, iq, lenr, locr )

      else if (mode .eq. 5) then             ! Solve  A w      = v
         call lu6L  (                        ! via    L v(new) = v
     &        inform, m, n, v,
     &        lena, luparm, parmlu,
     &        a, indc, indr, lenc )
         call lu6U  (                        ! and    U w = v(new).
     &        inform, m, n, v, w,
     &        lena, luparm, parmlu,
     &        a, indr, ip, iq, lenr, locr )

      else if (mode .eq. 6) then             ! Solve  A'v = w
         call lu6Ut (                        ! via    U'v = w
     &        inform, m, n, v, w,
     &        lena, luparm, parmlu,
     &        a, indr, ip, iq, lenr, locr )
         call lu6Lt (                        ! and    L'v(new) = v.
     &        inform, m, n, v,
     &        lena, luparm, parmlu,
     &        a, indc, indr, lenc )

      else if (mode .eq. 7) then
         call lu6LD (                        ! Solve  LDv(bar) = v
     &        inform, 1, m, n, v,
     &        lena, luparm, parmlu,
     &        a, indc, indr, lenc, locr )
         call lu6Lt (                        ! and    L'v(new) = v(bar).
     &        inform, m, n, v,
     &        lena, luparm, parmlu,
     &        a, indc, indr, lenc )

      else if (mode .eq. 8) then
         call lu6LD (                        ! Solve  L|D|v(bar) = v
     &        inform, 2, m, n, v,
     &        lena, luparm, parmlu,
     &        a, indc, indr, lenc, locr )
         call lu6Lt (                        ! and    L'v(new) = v(bar).
     &        inform, m, n, v,
     &        lena, luparm, parmlu,
     &        a, indc, indr, lenc )
      end if

      end ! subroutine lu6sol

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu6L  (
     &     inform, m, n, v,
     &     lena, luparm, parmlu,
     &     a, indc, indr, lenc )

      implicit
     &     none
      integer
     &     inform, m, n, lena, luparm(30),
     &     indc(lena), indr(lena), lenc(n)
      double precision
     &     parmlu(30), a(lena), v(m)

!     ------------------------------------------------------------------
!     lu6L   solves   L v = v(input).
!
!     15 Dec 2002: First version derived from lu6sol.
!     15 Dec 2002: Current version.
!     ------------------------------------------------------------------

      integer
     &     i, ipiv, j, k, l, l1, ldummy, len, lenL, lenL0, numL, numL0
      double precision
     &     small, vpiv

      numL0  = luparm(20)
      lenL0  = luparm(21)
      lenL   = luparm(23)
      small  = parmlu(3)
      inform = 0
      l1     = lena + 1

      do k = 1, numL0
         len   = lenc(k)
         l     = l1
         l1    = l1 - len
         ipiv  = indr(l1)
         vpiv  = v(ipiv)

         if (abs( vpiv ) .gt. small) then
            !***** This loop could be coded specially.
            do ldummy = 1, len
               l    = l - 1
               j    = indc(l)
               v(j) = v(j)  +  a(l) * vpiv
            end do
         end if
      end do

      l      = lena - lenL0 + 1
      numL   = lenL - lenL0

      !***** This loop could be coded specially.

      do ldummy = 1, numL
         l      = l - 1
         i      = indr(l)
         if (abs( v(i) ) .gt. small) then
            j    = indc(l)
            v(j) = v(j)  +  a(l) * v(i)
         end if
      end do

!     Exit.

      luparm(10) = inform

      end ! subroutine lu6L

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu6Lt (
     &     inform, m, n, v,
     &     lena, luparm, parmlu,
     &     a, indc, indr, lenc )

      implicit
     &     none
      integer
     &     inform, m, n, lena, luparm(30),
     &     indc(lena), indr(lena), lenc(n)
      double precision
     &     parmlu(30), a(lena), v(m)

!     ------------------------------------------------------------------
!     lu6Lt  solves   L'v = v(input).
!
!     15 Dec 2002: First version derived from lu6sol.
!     15 Dec 2002: Current version.
!     ------------------------------------------------------------------

      integer
     &     i, ipiv, j, k, l, l1, l2, len, lenL, lenL0, numL0
      double precision
     &     small, sum

!     ------------------------------------------------------------------
      double precision   zero
      parameter        ( zero = 0.0d+0 )
!     ------------------------------------------------------------------

      numL0  = luparm(20)
      lenL0  = luparm(21)
      lenL   = luparm(23)
      small  = parmlu(3)
      inform = 0
      l1     = lena - lenL + 1
      l2     = lena - lenL0

      !***** This loop could be coded specially.
      do l = l1, l2
         j     = indc(l)
         if (abs( v(j) ) .gt. small) then
            i     = indr(l)
            v(i)  = v(i)  +  a(l) * v(j)
         end if
      end do

      do k = numL0, 1, -1
         len   = lenc(k)
         sum   = zero
         l1    = l2 + 1
         l2    = l2 + len

         !***** This loop could be coded specially.
         do l = l1, l2
            j     = indc(l)
            sum   = sum  +  a(l) * v(j)
         end do

         ipiv    = indr(l1)
         v(ipiv) = v(ipiv) + sum
      end do

!     Exit.

      luparm(10) = inform

      end ! subroutine lu6Lt

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu6U  (
     &     inform, m, n, v, w,
     &     lena, luparm, parmlu,
     &     a, indr, ip, iq, lenr, locr )

      implicit
     &     none
      integer
     &     inform, m, n, lena, luparm(30),
     &     indr(lena), ip(m), iq(n), lenr(m), locr(m)
      double precision
     &     parmlu(30), a(lena), v(m), w(n)

!     ------------------------------------------------------------------
!     lu6U   solves   U w = v.          v  is not altered.
!
!     15 Dec 2002: First version derived from lu6sol.
!     15 Dec 2002: Current version.
!     ------------------------------------------------------------------

      integer
     &     i, j, k, klast, l, l1, l2, l3, nrank, nrank1
      double precision
     &     resid, small, t

!     ------------------------------------------------------------------
      double precision   zero
      parameter        ( zero = 0.0d+0 )
!     ------------------------------------------------------------------

      nrank  = luparm(16)
      small  = parmlu(3)
      inform = 0
      nrank1 = nrank + 1
      resid  = zero

!     Find the first nonzero in v(1:nrank), counting backwards.

      do klast = nrank, 1, -1
         i      = ip(klast)
         if (abs( v(i) ) .gt. small) go to 320
      end do

  320 do k = klast + 1, n
         j     = iq(k)
         w(j)  = zero
      end do

!     Do the back-substitution, using rows 1:klast of U.

      do k  = klast, 1, -1
         i      = ip(k)
         t      = v(i)
         l1     = locr(i)
         l2     = l1 + 1
         l3     = l1 + lenr(i) - 1

         !***** This loop could be coded specially.
         do l = l2, l3
            j     = indr(l)
            t     = t  -  a(l) * w(j)
         end do

         j      = iq(k)
         if (abs( t ) .le. small) then
            w(j)  = zero
         else
            w(j)  = t / a(l1)
         end if
      end do

!     Compute residual for overdetermined systems.

      do k = nrank1, m
         i     = ip(k)
         resid = resid  +  abs( v(i) )
      end do

!     Exit.

      if (resid .gt. zero) inform = 1
      luparm(10) = inform
      parmlu(20) = resid

      end ! subroutine lu6U

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu6Ut (
     &     inform, m, n, v, w,
     &     lena, luparm, parmlu,
     &     a, indr, ip, iq, lenr, locr )

      implicit
     &     none
      integer
     &     inform, m, n, lena, luparm(30),
     &     indr(lena), ip(m), iq(n), lenr(m), locr(m)
      double precision
     &     parmlu(30), a(lena), v(m), w(n)

!     ------------------------------------------------------------------
!     lu6Ut  solves   U'v = w.          w  is destroyed.
!
!     15 Dec 2002: First version derived from lu6sol.
!     15 Dec 2002: Current version.
!     ------------------------------------------------------------------

      integer
     &     i, j, k, l, l1, l2, nrank, nrank1
      double precision
     &     resid, small, t

!     ------------------------------------------------------------------
      double precision   zero
      parameter        ( zero = 0.0d+0 )
!     ------------------------------------------------------------------

      nrank  = luparm(16)
      small  = parmlu(3)
      inform = 0
      nrank1 = nrank + 1
      resid  = zero

      do k = nrank1, m
         i     = ip(k)
         v(i)  = zero
      end do

!     Do the forward-substitution, skipping columns of U(transpose)
!     when the associated element of w(*) is negligible.

      do 480 k = 1, nrank
         i      = ip(k)
         j      = iq(k)
         t      = w(j)
         if (abs( t ) .le. small) then
            v(i) = zero
            go to 480
         end if

         l1     = locr(i)
         t      = t / a(l1)
         v(i)   = t
         l2     = l1 + lenr(i) - 1
         l1     = l1 + 1

         !***** This loop could be coded specially.
         do l = l1, l2
            j     = indr(l)
            w(j)  = w(j)  -  t * a(l)
         end do
  480 continue

!     Compute residual for overdetermined systems.

      do k = nrank1, n
         j     = iq(k)
         resid = resid  +  abs( w(j) )
      end do

!     Exit.

      if (resid .gt. zero) inform = 1
      luparm(10) = inform
      parmlu(20) = resid

      end ! subroutine lu6Ut

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu6LD (
     &     inform, mode, m, n, v,
     &     lena, luparm, parmlu,
     &     a, indc, indr, lenc, locr )

      implicit
     &     none
      integer
     &     luparm(30), inform, mode, m, n, lena,
     &     indc(lena), indr(lena), lenc(n), locr(m)
      double precision
     &     parmlu(30), a(lena), v(m)

!-----------------------------------------------------------------------
!     lu6LD  assumes lu1fac has computed factors A = LU of a
!     symmetric definite or quasi-definite matrix A,
!     using Threshold Symmetric Pivoting (TSP),   luparm(6) = 3,
!     or    Threshold Diagonal  Pivoting (TDP),   luparm(6) = 4.
!     It also assumes that no updates have been performed.
!     In such cases,  U = D L', where D = diag(U).
!     lu6LDL returns v as follows:
!
!     mode
!      1    v  solves   L D v = v(input).
!      2    v  solves   L|D|v = v(input).
!
!     15 Dec 2002: First version of lu6LD.
!     15 Dec 2002: Current version.
!-----------------------------------------------------------------------

      ! Solve L D v(new) = v  or  L|D|v(new) = v, depending on mode.
      ! The code for L is the same as in lu6L,
      ! but when a nonzero entry of v arises, we divide by
      ! the corresponding entry of D or |D|.

      integer
     &     ipiv, j, k, l, l1, ldummy, len, numL0
      double precision
     &     diag, small, vpiv

      numL0  = luparm(20)
      small  = parmlu(3)
      inform = 0
      l1     = lena + 1

      do k = 1, numL0
         len   = lenc(k)
         l     = l1
         l1    = l1 - len
         ipiv  = indr(l1)
         vpiv  = v(ipiv)

         if (abs( vpiv ) .gt. small) then
            !***** This loop could be coded specially.
            do ldummy = 1, len
               l    = l - 1
               j    = indc(l)
               v(j) = v(j)  +  a(l) * vpiv
            end do

            ! Find diag = U(ipiv,ipiv) and divide by diag or |diag|.

            l    = locr(ipiv)
            diag = A(l)
            if (mode .eq. 2) diag = abs( diag )
            v(ipiv) = vpiv / diag
         end if
      end do

      end ! subroutine lu6LD

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu6chk( mode, m, n, w,
     &                   lena, luparm, parmlu,
     &                   a, indc, indr, ip, iq,
     &                   lenc, lenr, locc, locr,
     &                   inform )

      implicit
     &     none
      integer
     &     mode, m, n, lena, inform,
     &     luparm(30), indc(lena), indr(lena), ip(m), iq(n),
     &     lenc(n), lenr(m), locc(n), locr(m)
      double precision
     &     parmlu(30), a(lena), w(n)

!     ------------------------------------------------------------------
!     lu6chk  looks at the LU factorization  A = L*U.
!
!     If mode = 1, lu6chk is being called by lu1fac.
!     (Other modes not yet implemented.)
!     The important input parameters are
!
!                    lprint = luparm(2)
!                             luparm(6) = 1 if TRP
!                    keepLU = luparm(8)
!                    Utol1  = parmlu(4)
!                    Utol2  = parmlu(5)
!
!     and the significant output parameters are
!
!                    inform = luparm(10)
!                    nsing  = luparm(11)
!                    jsing  = luparm(12)
!                    jumin  = luparm(19)
!                    Lmax   = parmlu(11)
!                    Umax   = parmlu(12)
!                    DUmax  = parmlu(13)
!                    DUmin  = parmlu(14)
!                    and      w(*).
!
!     Lmax  and Umax  return the largest elements in L and U.
!     DUmax and DUmin return the largest and smallest diagonals of U
!                     (excluding diagonals that are exactly zero).
!
!     In general, w(j) is set to the maximum absolute element in
!     the j-th column of U.  However, if the corresponding diagonal
!     of U is small in absolute terms or relative to w(j)
!     (as judged by the parameters Utol1, Utol2 respectively),
!     then w(j) is changed to - w(j).
!
!     Thus, if w(j) is not positive, the j-th column of A
!     appears to be dependent on the other columns of A.
!     The number of such columns, and the position of the last one,
!     are returned as nsing and jsing.
!
!     Note that nrank is assumed to be set already, and is not altered.
!     Typically, nsing will satisfy      nrank + nsing = n,  but if
!     Utol1 and Utol2 are rather large,  nsing > n - nrank   may occur.
!
!     If keepLU = 0,
!     Lmax  and Umax  are already set by lu1fac.
!     The diagonals of U are in the top of A.
!     Only Utol1 is used in the singularity test to set w(*).
!
!     inform = 0  if  A  appears to have full column rank  (nsing = 0).
!     inform = 1  otherwise  (nsing .gt. 0).
!
!     00 Jul 1987: Early version.
!     09 May 1988: f77 version.
!     11 Mar 2001: Allow for keepLU = 0.
!     17 Nov 2001: Briefer output for singular factors.
!     05 May 2002: Comma needed in format 1100 (via Kenneth Holmstrom).
!     06 May 2002: With keepLU = 0, diags of U are in natural order.
!                  They were not being extracted correctly.
!     23 Apr 2004: TRP can judge singularity better by comparing
!                  all diagonals to DUmax.
!     27 Jun 2004: (PEG) Allow write only if nout .gt. 0 .
!     ------------------------------------------------------------------

      character
     &     mnkey
      logical
     &     keepLU, TRP
      integer
     &     i, j, jsing, jumin, k, l, l1, l2, ldiagU, lenL, lprint,
     &     ndefic, nout, nrank, nsing
      double precision
     &     aij, diag, DUmax, DUmin, Lmax, Umax, Utol1, Utol2

      double precision   zero
      parameter        ( zero = 0.0d+0 )

      nout   = luparm(1)
      lprint = luparm(2)
      TRP    = luparm(6) .eq. 1  ! Threshold Rook Pivoting
      keepLU = luparm(8) .ne. 0
      nrank  = luparm(16)
      lenL   = luparm(23)
      Utol1  = parmlu(4)
      Utol2  = parmlu(5)

      inform = 0
      Lmax   = zero
      Umax   = zero
      nsing  = 0
      jsing  = 0
      jumin  = 0
      DUmax  = zero
      DUmin  = 1.0d+30

      do j = 1, n
         w(j) = zero
      end do


      if (keepLU) then
         !--------------------------------------------------------------
         ! Find  Lmax.
         !--------------------------------------------------------------
         do l = lena + 1 - lenL, lena
            Lmax  = max( Lmax, abs(a(l)) )
         end do

         !--------------------------------------------------------------
         ! Find Umax and set w(j) = maximum element in j-th column of U.
         !--------------------------------------------------------------
         do k = 1, nrank
            i     = ip(k)
            l1    = locr(i)
            l2    = l1 + lenr(i) - 1

            do l = l1, l2
               j     = indr(l)
               aij   = abs( a(l) )
               w(j)  = max( w(j), aij )
               Umax  = max( Umax, aij )
            end do
         end do

         parmlu(11) = Lmax
         parmlu(12) = Umax

         !--------------------------------------------------------------
         ! Find DUmax and DUmin, the extreme diagonals of U.
         !--------------------------------------------------------------
         do k = 1, nrank
            j      = iq(k)
            i      = ip(k)
            l1     = locr(i)
            diag   = abs( a(l1) )
            DUmax  = max( DUmax, diag )
            if (DUmin .gt. diag) then
               DUmin  =   diag
               jumin  =   j
            end if
         end do

      else
         !--------------------------------------------------------------
         ! keepLU = 0.
         ! Only diag(U) is stored.  Set w(*) accordingly.
         ! Find DUmax and DUmin, the extreme diagonals of U.
         !--------------------------------------------------------------
         ldiagU = lena - n

         do k = 1, nrank
            j      = iq(k)
          !!diag   = abs( a(ldiagU + k) ) ! 06 May 2002: Diags
            diag   = abs( a(ldiagU + j) ) ! are in natural order
            w(j)   = diag
            DUmax  = max( DUmax, diag )
            if (DUmin .gt. diag) then
               DUmin  =   diag
               jumin  =   j
            end if
         end do
      end if


      !--------------------------------------------------------------
      ! Negate w(j) if the corresponding diagonal of U is
      ! too small in absolute terms or relative to the other elements
      ! in the same column of  U.
      !
      ! 23 Apr 2004: TRP ensures that diags are NOT small relative to
      !              other elements in their own column.
      !              Much better, we can compare all diags to DUmax.
      !--------------------------------------------------------------
      if (mode .eq. 1  .and.  TRP) then
         Utol1 = max( Utol1, Utol2*DUmax )
      end if

      if (keepLU) then
         do k = 1, n
            j     = iq(k)
            if (k .gt. nrank) then
               diag   = zero
            else
               i      = ip(k)
               l1     = locr(i)
               diag   = abs( a(l1) )
            end if

            if (diag .le. Utol1  .or.  diag .le. Utol2*w(j)) then
               nsing  =   nsing + 1
               jsing  =   j
               w(j)   = - w(j)
            end if
         end do

      else ! keepLU = 0

         do k = 1, n
            j      = iq(k)
            diag   = w(j)

            if (diag .le. Utol1) then
               nsing  =   nsing + 1
               jsing  =   j
               w(j)   = - w(j)
            end if
         end do
      end if


      !-----------------------------------------------------------------
      ! Set output parameters.
      !-----------------------------------------------------------------
      if (jumin .eq. 0) DUmin = zero
      luparm(11) = nsing
      luparm(12) = jsing
      luparm(19) = jumin
      parmlu(13) = DUmax
      parmlu(14) = DUmin

      if (nsing .gt. 0) then  ! The matrix has been judged singular.
         inform = 1
         ndefic = n - nrank
         if (nout .gt. 0  .and.  lprint .ge. 0) then
            if (m .gt. n) then
               mnkey  = '>'
            else if (m .eq. n) then
               mnkey  = '='
            else
               mnkey  = '<'
            end if
            write(nout, 1100) mnkey, nrank, ndefic, nsing
         end if
      end if

!     Exit.

      luparm(10) = inform
      return

 1100 format(' Singular(m', a, 'n)',
     &       '  rank', i9, '  n-rank', i8, '  nsing', i9)

      end ! subroutine lu6chk

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
************************************************************************
*
*     File  lusol7a.f
*
*     lu7add   lu7cyc   lu7elm   lu7for   lu7rnk   lu7zap
*
*     Utilities for LUSOL's update routines.
*     lu7for is the most important -- the forward sweep.
*
* 01 May 2002: Derived from LUSOL's original lu7a.f file.
* 01 May 2002: Current version of lusol7a.f.
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu7add( m, n, jadd, v,
     $                   lena, luparm, parmlu,
     $                   lenL, lenU, lrow, nrank,
     $                   a, indr, ip, lenr, locr,
     $                   inform, klast, vnorm )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena), v(m)
      integer            indr(lena), ip(m), lenr(m)
      integer            locr(m)

*     ------------------------------------------------------------------
*     lu7add  inserts the first nrank elements of the vector v(*)
*     as column  jadd  of  U.  We assume that  U  does not yet have any
*     entries in this column.
*     Elements no larger than  parmlu(3)  are treated as zero.
*     klast  will be set so that the last row to be affected
*     (in pivotal order) is row  ip(klast).
*
*     09 May 1988: First f77 version.
*     ------------------------------------------------------------------

      parameter        ( zero = 0.0d+0 )

      small  = parmlu(3)
      vnorm  = zero
      klast  = 0

      do 200 k  = 1, nrank
         i      = ip(k)
         if (abs( v(i) ) .le. small) go to 200
         klast  = k
         vnorm  = vnorm  +  abs( v(i) )
         leni   = lenr(i)

*        Compress row file if necessary.

         minfre = leni + 1
         nfree  = lena - lenL - lrow
         if (nfree .lt. minfre) then
            call lu1rec( m, .true., luparm, lrow, lena,
     $                   a, indr, lenr, locr )
            nfree  = lena - lenL - lrow
            if (nfree .lt. minfre) go to 970
         end if

*        Move row  i  to the end of the row file,
*        unless it is already there.
*        No need to move if there is a gap already.

         if (leni .eq. 0) locr(i) = lrow + 1
         lr1    = locr(i)
         lr2    = lr1 + leni - 1
         if (lr2    .eq.   lrow) go to 150
         if (indr(lr2+1) .eq. 0) go to 180
         locr(i) = lrow + 1

         do 140 l = lr1, lr2
            lrow       = lrow + 1
            a(lrow)    = a(l)
            j          = indr(l)
            indr(l)    = 0
            indr(lrow) = j
  140    continue

  150    lr2     = lrow
         lrow    = lrow + 1

*        Add the element of  v.

  180    lr2       = lr2 + 1
         a(lr2)    = v(i)
         indr(lr2) = jadd
         lenr(i)   = leni + 1
         lenU      = lenU + 1
  200 continue

*     Normal exit.

      inform = 0
      go to 990

*     Not enough storage.

  970 inform = 7

  990 return

      end ! subroutine lu7add

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu7cyc( kfirst, klast, ip )

      integer            ip(klast)

*     ------------------------------------------------------------------
*     lu7cyc performs a cyclic permutation on the row or column ordering
*     stored in ip, moving entry kfirst down to klast.
*     If kfirst .ge. klast, lu7cyc should not be called.
*     Sometimes klast = 0 and nothing should happen.
*
*     09 May 1988: First f77 version.
*     ------------------------------------------------------------------

      if (kfirst .lt. klast) then
         ifirst = ip(kfirst)

         do 100 k = kfirst, klast - 1
            ip(k) = ip(k + 1)
  100    continue

         ip(klast) = ifirst
      end if

      end ! subroutine lu7cyc

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena), v(m)
      integer            indc(lena), indr(lena), ip(m), iq(n), lenr(m)
      integer            locc(n), locr(m)

*     ------------------------------------------------------------------
*     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.
*     ------------------------------------------------------------------

      parameter        ( zero = 0.0d+0 )

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

*     Compress row file if necessary.

      minfre = m - nrank
      nfree  = lena - lenL - lrow
      if (nfree .ge. minfre) go to 100
      call lu1rec( m, .true., luparm, lrow, lena, a, indr, lenr, locr )
      nfree  = lena - lenL - lrow
      if (nfree .lt. minfre) go to 970

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

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

      do 200 k = nrank1, m
         i       = ip(k)
         vi      = abs( v(i) )
         if (vi .le. small) go to 200
         l       = l - 1
         a(l)    = v(i)
         indc(l) = i
         if (vmax .ge. vi ) go to 200
         vmax    = vi
         kmax    = k
         lmax    = l
  200 continue

      if (kmax .eq. 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       = ip(kmax)
      vmax       = a(lmax)
      a(lmax)    = a(l)
      indc(lmax) = indc(l)
      l1         = l + 1
      l2         = lena - lenL
      lenL       = lenL + (l2 - l)

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

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

      ip(kmax  ) = ip(nrank1)
      ip(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 .gt. 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

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu7for( m, n, kfirst, klast,
     $                   lena, luparm, parmlu,
     $                   lenL, lenU, lrow,
     $                   a, indc, indr, ip, iq, lenr, locc, locr,
     $                   inform, diag )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena)
      integer            indc(lena), indr(lena), ip(m), iq(n), lenr(m)
      integer            locc(n), locr(m)

*     ------------------------------------------------------------------
*     lu7for  (forward sweep) updates the LU factorization  A = L*U
*     when row  iw = ip(klast)  of  U  is eliminated by a forward
*     sweep of stabilized row operations, leaving  ip * U * iq  upper
*     triangular.
*
*     The row permutation  ip  is updated to preserve stability and/or
*     sparsity.  The column permutation  iq  is not altered.
*
*     kfirst  is such that row  ip(kfirst)  is the first row involved
*     in eliminating row  iw.  (Hence,  kfirst  marks the first nonzero
*     in row  iw  in pivotal order.)  If  kfirst  is unknown it may be
*     input as  1.
*
*     klast   is such that row  ip(klast)  is the row being eliminated.
*     klast   is not altered.
*
*     lu7for  should be called only if  kfirst .le. klast.
*     If  kfirst = klast,  there are no nonzeros to eliminate, but the
*     diagonal element of row  ip(klast)  may need to be moved to the
*     front of the row.
*
*     On entry,  locc(*)  must be zero.
*
*     On exit:
*     inform = 0  if row iw has a nonzero diagonal (could be small).
*     inform = 1  if row iw has no diagonal.
*     inform = 7  if there is not enough storage to finish the update.
*
*     On a successful exit (inform le 1),  locc(*)  will again be zero.
*
*        Jan 1985: Final f66 version.
*     09 May 1988: First f77 version.
*     ------------------------------------------------------------------

      parameter        ( zero = 0.0d+0 )

      double precision   Ltol
      logical            swappd

      Ltol   = parmlu(2)
      small  = parmlu(3)
      uspace = parmlu(6)
      kbegin = kfirst
      swappd = .false.

*     We come back here from below if a row interchange is performed.

  100 iw     = ip(klast)
      lenw   = lenr(iw)
      if (lenw   .eq.   0  ) go to 910
      lw1    = locr(iw)
      lw2    = lw1 + lenw - 1
      jfirst = iq(kbegin)
      if (kbegin .ge. klast) go to 700

*     Make sure there is room at the end of the row file
*     in case row  iw  is moved there and fills in completely.

      minfre = n + 1
      nfree  = lena - lenL - lrow
      if (nfree .lt. minfre) then
         call lu1rec( m, .true., luparm, lrow, lena,
     $                a, indr, lenr, locr )
         lw1    = locr(iw)
         lw2    = lw1 + lenw - 1
         nfree  = lena - lenL - lrow
         if (nfree .lt. minfre) go to 970
      end if

*     Set markers on row  iw.

      do 120 l = lw1, lw2
         j       = indr(l)
         locc(j) = l
  120 continue


*     ==================================================================
*     Main elimination loop.
*     ==================================================================
      kstart = kbegin
      kstop  = min( klast, n )

      do 500 k  = kstart, kstop
         jfirst = iq(k)
         lfirst = locc(jfirst)
         if (lfirst .eq. 0) go to 490

*        Row  iw  has its first element in column  jfirst.

         wj     = a(lfirst)
         if (k .eq. klast) go to 490

*        ---------------------------------------------------------------
*        We are about to use the first element of row  iv
*               to eliminate the first element of row  iw.
*        However, we may wish to interchange the rows instead,
*        to preserve stability and/or sparsity.
*        ---------------------------------------------------------------
         iv     = ip(k)
         lenv   = lenr(iv)
         lv1    = locr(iv)
         vj     = zero
         if (lenv      .eq.   0   ) go to 150
         if (indr(lv1) .ne. jfirst) go to 150
         vj     = a(lv1)
         if (            swappd               ) go to 200
         if (Ltol * abs( wj )  .lt.  abs( vj )) go to 200
         if (Ltol * abs( vj )  .lt.  abs( wj )) go to 150
         if (            lenv  .le.  lenw     ) go to 200

*        ---------------------------------------------------------------
*        Interchange rows  iv  and  iw.
*        ---------------------------------------------------------------
  150    ip(klast) = iv
         ip(k)     = iw
         kbegin    = k
         swappd    = .true.
         go to 600

*        ---------------------------------------------------------------
*        Delete the eliminated element from row  iw
*        by overwriting it with the last element.
*        ---------------------------------------------------------------
  200    a(lfirst)    = a(lw2)
         jlast        = indr(lw2)
         indr(lfirst) = jlast
         indr(lw2)    = 0
         locc(jlast)  = lfirst
         locc(jfirst) = 0
         lenw         = lenw - 1
         lenU         = lenU - 1
         if (lrow .eq. lw2) lrow = lrow - 1
         lw2          = lw2  - 1

*        ---------------------------------------------------------------
*        Form the multiplier and store it in the  L  file.
*        ---------------------------------------------------------------
         if (abs( wj ) .le. small) go to 490
         amult   = - wj / vj
         l       = lena - lenL
         a(l)    = amult
         indr(l) = iv
         indc(l) = iw
         lenL    = lenL + 1

*        ---------------------------------------------------------------
*        Add the appropriate multiple of row  iv  to row  iw.
*        We use two different inner loops.  The first one is for the
*        case where row  iw  is not at the end of storage.
*        ---------------------------------------------------------------
         if (lenv .eq. 1) go to 490
         lv2    = lv1 + 1
         lv3    = lv1 + lenv - 1
         if (lw2 .eq. lrow) go to 400

*        ...............................................................
*        This inner loop will be interrupted only if
*        fill-in occurs enough to bump into the next row.
*        ...............................................................
         do 350 lv = lv2, lv3
            jv     = indr(lv)
            lw     = locc(jv)
            if (lw .gt. 0) then

*              No fill-in.

               a(lw)  = a(lw)  +  amult * a(lv)
               if (abs( a(lw) ) .le. small) then

*                 Delete small computed element.

                  a(lw)     = a(lw2)
                  j         = indr(lw2)
                  indr(lw)  = j
                  indr(lw2) = 0
                  locc(j)   = lw
                  locc(jv)  = 0
                  lenU      = lenU - 1
                  lenw      = lenw - 1
                  lw2       = lw2  - 1
               end if
            else

*              Row  iw  doesn't have an element in column  jv  yet
*              so there is a fill-in.

               if (indr(lw2+1) .ne. 0) go to 360
               lenU      = lenU + 1
               lenw      = lenw + 1
               lw2       = lw2  + 1
               a(lw2)    = amult * a(lv)
               indr(lw2) = jv
               locc(jv)  = lw2
            end if
  350    continue

         go to 490

*        Fill-in interrupted the previous loop.
*        Move row  iw  to the end of the row file.

  360    lv2      = lv
         locr(iw) = lrow + 1

         do 370 l = lw1, lw2
            lrow       = lrow + 1
            a(lrow)    = a(l)
            j          = indr(l)
            indr(l)    = 0
            indr(lrow) = j
            locc(j)    = lrow
  370    continue

         lw1    = locr(iw)
         lw2    = lrow

*        ...............................................................
*        Inner loop with row  iw  at the end of storage.
*        ...............................................................
  400    do 450 lv = lv2, lv3
            jv     = indr(lv)
            lw     = locc(jv)
            if (lw .gt. 0) then

*              No fill-in.

               a(lw)  = a(lw)  +  amult * a(lv)
               if (abs( a(lw) ) .le. small) then

*                 Delete small computed element.

                  a(lw)     = a(lw2)
                  j         = indr(lw2)
                  indr(lw)  = j
                  indr(lw2) = 0
                  locc(j)   = lw
                  locc(jv)  = 0
                  lenU      = lenU - 1
                  lenw      = lenw - 1
                  lw2       = lw2  - 1
               end if
            else

*              Row  iw  doesn't have an element in column  jv  yet
*              so there is a fill-in.

               lenU      = lenU + 1
               lenw      = lenw + 1
               lw2       = lw2  + 1
               a(lw2)    = amult * a(lv)
               indr(lw2) = jv
               locc(jv)  = lw2
            end if
  450    continue

         lrow   = lw2

*        The  k-th  element of row  iw  has been processed.
*        Reset  swappd  before looking at the next element.

  490    swappd = .false.
  500 continue

*     ==================================================================
*     End of main elimination loop.
*     ==================================================================

*     Cancel markers on row  iw.

  600 lenr(iw) = lenw
      if (lenw .eq. 0) go to 910
      do 620 l = lw1, lw2
         j       = indr(l)
         locc(j) = 0
  620 continue

*     Move the diagonal element to the front of row  iw.
*     At this stage,  lenw gt 0  and  klast le n.

  700 do 720 l = lw1, lw2
         ldiag = l
         if (indr(l) .eq. jfirst) go to 730
  720 continue
      go to 910

  730 diag        = a(ldiag)
      a(ldiag)    = a(lw1)
      a(lw1)      = diag
      indr(ldiag) = indr(lw1)
      indr(lw1)   = jfirst

*     If an interchange is needed, repeat from the beginning with the
*     new row  iw,  knowing that the opposite interchange cannot occur.

      if (swappd) go to 100
      inform = 0
      go to 950

*     Singular.

  910 diag   = zero
      inform = 1

*     Force a compression if the file for  U  is much longer than the
*     no. of nonzeros in  U  (i.e. if  lrow  is much bigger than  lenU).
*     This should prevent memory fragmentation when there is far more
*     memory than necessary  (i.e. when  lena  is huge).

  950 limit  = uspace * lenU + m + n + 1000
      if (lrow .gt. limit) then
         call lu1rec( m, .true., luparm, lrow, lena,
     $                a, indr, lenr, locr )
      end if
      go to 990

*     Not enough storage.

  970 inform = 7

*     Exit.

  990 return

      end ! subroutine lu7for

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu7rnk( m, n, jsing,
     $                   lena, luparm, parmlu,
     $                   lenL, lenU, lrow, nrank,
     $                   a, indc, indr, ip, iq, lenr, locc, locr,
     $                   inform, diag )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena)
      integer            indc(lena), indr(lena), ip(m), iq(n), lenr(m)
      integer            locc(n), locr(m)

*     ------------------------------------------------------------------
*     lu7rnk (check rank) assumes U is currently nrank by n
*     and determines if row nrank contains an acceptable pivot.
*     If not, the row is deleted and nrank is decreased by 1.
*
*     jsing is an input parameter (not altered).  If jsing is positive,
*     column jsing has already been judged dependent.  A substitute
*     (if any) must be some other column.
*
*     -- Jul 1987: First version.
*     09 May 1988: First f77 version.
*     ------------------------------------------------------------------

      parameter        ( zero = 0.0d+0 )

      Utol1    = parmlu(4)
      diag     = zero

*     Find Umax, the largest element in row nrank.

      iw       = ip(nrank)
      lenw     = lenr(iw)
      if (lenw .eq. 0) go to 400
      l1       = locr(iw)
      l2       = l1 + lenw - 1
      Umax     = zero
      lmax     = l1

      do 100 l = l1, l2
         if (Umax .lt. abs( a(l) )) then
             Umax   =  abs( a(l) )
             lmax   =  l
         end if
  100 continue

*     Find which column that guy is in (in pivotal order).
*     Interchange him with column nrank, then move him to be
*     the new diagonal at the front of row nrank.

      diag   = a(lmax)
      jmax   = indr(lmax)

      do 300 kmax = nrank, n
         if (iq(kmax) .eq. jmax) go to 320
  300 continue

  320 iq(kmax)  = iq(nrank)
      iq(nrank) = jmax
      a(lmax)   = a(l1)
      a(l1)     = diag
      indr(lmax)= indr(l1)
      indr(l1)  = jmax

*     See if the new diagonal is big enough.

      if (Umax .le. Utol1) go to 400
      if (jmax .eq. jsing) go to 400

*     ------------------------------------------------------------------
*     The rank stays the same.
*     ------------------------------------------------------------------
      inform = 0
      return

*     ------------------------------------------------------------------
*     The rank decreases by one.
*     ------------------------------------------------------------------
  400 inform = -1
      nrank  = nrank - 1
      if (lenw .gt. 0) then

*        Delete row nrank from U.

         lenU     = lenU - lenw
         lenr(iw) = 0
         do 420 l = l1, l2
            indr(l) = 0
  420    continue

         if (l2 .eq. lrow) then

*           This row was at the end of the data structure.
*           We have to reset lrow.
*           Preceding rows might already have been deleted, so we
*           have to be prepared to go all the way back to 1.

            do 450 l = 1, l2
               if (indr(lrow) .gt. 0) go to 900
               lrow  = lrow - 1
  450       continue
         end if
      end if

  900 return

      end ! subroutine lu7rnk

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu7zap( m, n, jzap, kzap,
     $                   lena, lenU, lrow, nrank,
     $                   a, indr, ip, iq, lenr, locr )

      implicit           double precision (a-h,o-z)
      double precision   a(lena)
      integer            indr(lena), ip(m), iq(n), lenr(m)
      integer            locr(m)

*     ------------------------------------------------------------------
*     lu7zap  eliminates all nonzeros in column  jzap  of  U.
*     It also sets  kzap  to the position of  jzap  in pivotal order.
*     Thus, on exit we have  iq(kzap) = jzap.
*
*     -- Jul 1987: nrank added.
*     10 May 1988: First f77 version.
*     ------------------------------------------------------------------

      do 100 k  = 1, nrank
         i      = ip(k)
         leni   = lenr(i)
         if (leni .eq. 0) go to 90
         lr1    = locr(i)
         lr2    = lr1 + leni - 1
         do 50 l = lr1, lr2
            if (indr(l) .eq. jzap) go to 60
   50    continue
         go to 90

*        Delete the old element.

   60    a(l)      = a(lr2)
         indr(l)   = indr(lr2)
         indr(lr2) = 0
         lenr(i)   = leni - 1
         lenU      = lenU - 1

*        Stop if we know there are no more rows containing  jzap.

   90    kzap   = k
         if (iq(k) .eq. jzap) go to 800
  100 continue

*     nrank must be smaller than n because we haven't found kzap yet.

      do 200 k = nrank+1, n
         kzap  = k
         if (iq(k) .eq. jzap) go to 800
  200 continue

*     See if we zapped the last element in the file.

  800 if (lrow .gt. 0) then
         if (indr(lrow) .eq. 0) lrow = lrow - 1
      end if

      end ! subroutine lu7zap
************************************************************************
*
*     File  lusol8a.f
*
*     lu8rpc
*
*     Sparse LU update: Replace Column
*     LUSOL's sparse implementation of the Bartels-Golub update.
*
* 01 May 2002: Derived from LUSOL's original lu8a.f file.
* 01 May 2002: Current version of lusol8a.f.
* 15 Sep 2004: Test nout. gt. 0 to protect write statements.
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu8rpc( mode1, mode2, m, n, jrep, v, w,
     $                   lena, luparm, parmlu,
     $                   a, indc, indr, ip, iq,
     $                   lenc, lenr, locc, locr,
     $                   inform, diag, vnorm )

      implicit           double precision(a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena), v(m), w(n)
      integer            indc(lena), indr(lena), ip(m), iq(n)
      integer            lenc(n), lenr(m)
      integer            locc(n), locr(m)

*     ------------------------------------------------------------------
*     lu8rpc  updates the LU factorization  A = L*U  when column  jrep
*     is replaced by some vector  a(new).
*
*     lu8rpc  is an implementation of the Bartels-Golub update,
*     designed for the case where A is rectangular and/or singular.
*     L is a product of stabilized eliminations (m x m, nonsingular).
*     P U Q is upper trapezoidal (m x n, rank nrank).
*
*     If  mode1 = 0,  the old column is taken to be zero
*                     (so it does not have to be removed from  U).
*
*     If  mode1 = 1,  the old column need not have been zero.
*
*     If  mode2 = 0,  the new column is taken to be zero.
*                     v(*)  is not used or altered.
*
*     If  mode2 = 1,  v(*)  must contain the new column  a(new).
*                     On exit,  v(*)  will satisfy  L*v = a(new).
*
*     If  mode2 = 2,  v(*)  must satisfy  L*v = a(new).
*
*     The array  w(*)  is not used or altered.
*
*     On entry, all elements of  locc  are assumed to be zero.
*     On a successful exit (inform ne 7), this will again be true.
*
*     On exit:
*     inform = -1  if the rank of U decreased by 1.
*     inform =  0  if the rank of U stayed the same.
*     inform =  1  if the rank of U increased by 1.
*     inform =  2  if the update seemed to be unstable
*                  (diag much bigger than vnorm).
*     inform =  7  if the update was not completed (lack of storage).
*     inform =  8  if jrep is not between 1 and n.
*
*     -- Jan 1985: Original F66 version.
*     -- Jul 1987: Modified to maintain U in trapezoidal form.
*     10 May 1988: First f77 version.
*     16 Oct 2000: Added test for instability (inform = 2).
*     ------------------------------------------------------------------

      logical            singlr
      parameter        ( zero = 0.0d+0 )

      nout   = luparm(1)
      lprint = luparm(2)
      nrank  = luparm(16)
      lenL   = luparm(23)
      lenU   = luparm(24)
      lrow   = luparm(25)
      Utol1  = parmlu(4)
      Utol2  = parmlu(5)
      nrank0 = nrank
      diag   = zero
      vnorm  = zero
      if (jrep .lt. 1) go to 980
      if (jrep .gt. n) go to 980

*     ------------------------------------------------------------------
*     If mode1 = 0, there are no elements to be removed from  U
*     but we still have to set  krep  (using a backward loop).
*     Otherwise, use lu7zap to remove column  jrep  from  U
*     and set  krep  at the same time.
*     ------------------------------------------------------------------
      if (mode1 .eq. 0) then
         krep   = n + 1

   10    krep   = krep - 1
         if (iq(krep) .ne. jrep) go to 10
      else
         call lu7zap( m, n, jrep, krep,
     $                lena, lenU, lrow, nrank,
     $                a, indr, ip, iq, lenr, locr )
      end if

*     ------------------------------------------------------------------
*     Insert a new column of u and find klast.
*     ------------------------------------------------------------------

      if (mode2 .eq. 0) then
         klast  = 0
      else
         if (mode2 .eq. 1) then

*           Transform v = a(new) to satisfy  L*v = a(new).

            call lu6sol( 1, m, n, v, w, lena, luparm, parmlu,
     $                   a, indc, indr, ip, iq,
     $                   lenc, lenr, locc, locr, inform )
         end if

*        Insert into  U  any nonzeros in the top of  v.
*        row  ip(klast)  will contain the last nonzero in pivotal order.
*        Note that  klast  will be in the range  ( 0, nrank ).

         call lu7add( m, n, jrep, v,
     $                lena, luparm, parmlu,
     $                lenL, lenU, lrow, nrank,
     $                a, indr, ip, lenr, locr,
     $                inform, klast, vnorm )
         if (inform .eq. 7) go to 970
      end if

*     ------------------------------------------------------------------
*     In general, the new column causes U to look like this:
*
*                 krep        n                 krep  n
*
*                ....a.........          ..........a...
*                 .  a        .           .        a  .
*                  . a        .            .       a  .
*                   .a        .             .      a  .
*        P U Q =     a        .    or        .     a  .
*                    b.       .               .    a  .
*                    b .      .                .   a  .
*                    b  .     .                 .  a  .
*                    b   ......                  ..a...  nrank
*                    c                             c
*                    c                             c
*                    c                             c     m
*
*     klast points to the last nonzero "a" or "b".
*     klast = 0 means all "a" and "b" entries are zero.
*     ------------------------------------------------------------------

      if (mode2 .eq. 0) then
         if (krep .gt. nrank) go to 900
      else if (nrank .lt. m) then

*        Eliminate any "c"s (in either case).
*        Row nrank + 1 may end up containing one nonzero.

         call lu7elm( m, n, jrep, v,
     $                lena, luparm, parmlu,
     $                lenL, lenU, lrow, nrank,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )
         if (inform .eq. 7) go to 970

         if (inform .eq. 1) then

*           The nonzero is apparently significant.
*           Increase nrank by 1 and make klast point to the bottom.

            nrank = nrank + 1
            klast = nrank
         end if
      end if

      if (nrank .lt. n) then

*        The column rank is low.
*
*        In the first case, we want the new column to end up in
*        position nrank, so the trapezoidal columns will have a chance
*        later on (in lu7rnk) to pivot in that position.
*
*        Otherwise the new column is not part of the triangle.  We
*        swap it into position nrank so we can judge it for singularity.
*        lu7rnk might choose some other trapezoidal column later.

         if (krep .lt. nrank) then
            klast     = nrank
         else
            iq(krep ) = iq(nrank)
            iq(nrank) = jrep
            krep      = nrank
         end if
      end if

*     ------------------------------------------------------------------
*     If krep .lt. klast, there are some "b"s to eliminate:
*
*                  krep
*
*                ....a.........
*                 .  a        .
*                  . a        .
*                   .a        .
*        P U Q =     a        .  krep
*                    b.       .
*                    b .      .
*                    b  .     .
*                    b   ......  nrank
*
*     If krep .eq. klast, there are no "b"s, but the last "a" still
*     has to be moved to the front of row krep (by lu7for).
*     ------------------------------------------------------------------

      if (krep .le. klast) then

*        Perform a cyclic permutation on the current pivotal order,
*        and eliminate the resulting row spike.  krep becomes klast.
*        The final diagonal (if any) will be correctly positioned at
*        the front of the new krep-th row.  nrank stays the same.

         call lu7cyc( krep, klast, ip )
         call lu7cyc( krep, klast, iq )

         call lu7for( m, n, krep, klast,
     $                lena, luparm, parmlu,
     $                lenL, lenU, lrow,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )
         if (inform .eq. 7) go to 970
         krep   = klast

*        Test for instability (diag much bigger than vnorm).

         singlr = vnorm .lt. Utol2 * abs( diag )
         if ( singlr ) go to 920
      end if

*     ------------------------------------------------------------------
*     Test for singularity in column krep (where krep .le. nrank).
*     ------------------------------------------------------------------

      diag   = zero
      iw     = ip(krep)
      singlr = lenr(iw) .eq. 0

      if (.not. singlr) then
         l1     = locr(iw)
         j1     = indr(l1)
         singlr = j1 .ne. jrep

         if (.not. singlr) then
            diag   = a(l1)
            singlr = abs( diag ) .le. Utol1          .or.
     $               abs( diag ) .le. Utol2 * vnorm
         end if
      end if

      if ( singlr  .and.  krep .lt. nrank ) then

*        Perform cyclic permutations to move column jrep to the end.
*        Move the corresponding row to position nrank
*        then eliminate the resulting row spike.

         call lu7cyc( krep, nrank, ip )
         call lu7cyc( krep, n    , iq )

         call lu7for( m, n, krep, nrank,
     $                lena, luparm, parmlu,
     $                lenL, lenU, lrow,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )
         if (inform .eq. 7) go to 970
      end if

*     Find the best column to be in position nrank.
*     If singlr, it can't be the new column, jrep.
*     If nothing satisfactory exists, nrank will be decreased.

      if ( singlr  .or.  nrank .lt. n ) then
         jsing  = 0
         if ( singlr ) jsing = jrep

         call lu7rnk( m, n, jsing,
     $                lena, luparm, parmlu,
     $                lenL, lenU, lrow, nrank,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )
      end if

*     ------------------------------------------------------------------
*     Set inform for exit.
*     ------------------------------------------------------------------

  900 if (nrank .eq. nrank0) then
         inform =  0
      else if (nrank .lt. nrank0) then
         inform = -1
         if (nrank0 .eq. n) then
            if (nout. gt. 0  .and.  lprint .ge. 0)
     &           write(nout, 1100) jrep, diag
         end if
      else
         inform =  1
      end if
      go to 990

*     Instability.

  920 inform = 2
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1200) jrep, diag
      go to 990

*     Not enough storage.

  970 inform = 7
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1700) lena
      go to 990

*     jrep  is out of range.

  980 inform = 8
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1800) m, n, jrep

*     Exit.

  990 luparm(10) = inform
      luparm(15) = luparm(15) + 1
      luparm(16) = nrank
      luparm(23) = lenL
      luparm(24) = lenU
      luparm(25) = lrow
      return

 1100 format(/ ' lu8rpc  warning.  Singularity after replacing column.',
     $       '    jrep =', i8, '    diag =', 1p, e12.2 )
 1200 format(/ ' lu8rpc  warning.  Instability after replacing column.',
     $       '    jrep =', i8, '    diag =', 1p, e12.2 )
 1700 format(/ ' lu8rpc  error...  Insufficient storage.',
     $         '    lena =', i8)
 1800 format(/ ' lu8rpc  error...  jrep  is out of range.',
     $         '    m =', i8, '    n =', i8, '    jrep =', i8)

      end ! subroutine lu8rpc
