dwnlit Subroutine

   subroutine dwnlit (w, mdw, m, n, l, ipivot, itype, h, scale, &
                      rnorm, idope, dope, done)

   integer :: idope(*), ipivot(*), itype(*), l, m, mdw, n
   real(wp) :: dope(*), h(*), rnorm, scale(*), w(mdw,*)
   logical :: done

   real(wp) :: alsq, amax, eanorm, factor, hbar, rn, sparam(5), &
               t, tau
   integer :: i, i1, imax, ir, j, j1, jj, jp, krank, l1, lb, lend, me, &
              mend, niv, nsoln
   logical :: indep, recalc

   me     = idope(1)
   nsoln  = idope(2)
   l1     = idope(3)
   alsq   = dope(1)
   eanorm = dope(2)
   tau    = dope(3)
   lb     = min(m-1,l)
   recalc = .true.
   rnorm  = 0.0_wp
   krank  = 0

   ! We set FACTOR=1.0 so that the heavy weight ALAMDA will be
   ! included in the test for column independence.
   factor = 1.0_wp
   lend = l

   main : block

      do i=1,lb

         ! Set IR to point to the I-th row.
         ir = i
         mend = m
         call dwnlt1 (i, lend, m, ir, mdw, recalc, imax, hbar, h, scale, w)

         ! Update column SS and find pivot column.
         call dwnlt3 (i, imax, m, mdw, ipivot, h, w)

         do
            ! Perform column interchange.
            ! Test independence of incoming column.
            if (dwnlt2(me, mend, ir, factor, tau, scale, w(1,i))) then

               ! Eliminate I-th column below diagonal using modified Givens
               ! transformations applied to (A B).
               !
               ! When operating near the ME line, use the largest element
               ! above it as the pivot.
               do j=m,i+1,-1
                  jp = j-1
                  if (j==me+1) then
                     imax = me
                     amax = scale(me)*w(me,i)**2
                     do jp=j-1,i,-1
                        t = scale(jp)*w(jp,i)**2
                        if (t>amax) then
                           imax = jp
                           amax = t
                        endif
                     end do
                     jp = imax
                  endif
                  if (w(j,i)/=0.0_wp) then
                     call drotmg (scale(jp), scale(j), w(jp,i), w(j,i), &
                                 sparam)
                     w(j,i) = 0.0_wp
                     call drotm (n+1-i, w(jp,i+1), mdw, w(j,i+1), mdw, &
                                 sparam)
                  endif
               end do
               exit
            else if (lend>i) then
               ! Column I is dependent.  Swap with column LEND.
               ! Perform column interchange,
               ! and find column in remaining set with largest SS.
               call dwnlt3 (i, lend, m, mdw, ipivot, h, w)
               lend = lend - 1
               imax = idamax(lend-i+1, h(i), 1) + i - 1
               hbar = h(imax)
            else
               krank = i - 1
               exit main
            endif
         end do

      end do
      krank = l1

   end block main

   if (krank<me) then
      factor = alsq
      do i=krank+1,me
         call dcopy (l, [0.0_wp], 0, w(i,1), mdw)
      end do

      ! Determine the rank of the remaining equality constraint
      ! equations by eliminating within the block of constrained
      ! variables.  Remove any redundant constraints.

      recalc = .true.
      lb = min(l+me-krank, n)
      do i=l+1,lb
         ir = krank + i - l
         lend = n
         mend = me
         call dwnlt1 (i, lend, me, ir, mdw, recalc, imax, hbar, h, &
                      scale, w)

         ! Update col ss and find pivot col
         call dwnlt3 (i, imax, m, mdw, ipivot, h, w)

         ! Perform column interchange
         ! Eliminate elements in the I-th col.
         do j=me,ir+1,-1
            if (w(j,i)/=0.0_wp) then
              call drotmg (scale(j-1), scale(j), w(j-1,i), w(j,i), &
                           sparam)
               w(j,i) = 0.0_wp
               call drotm (n+1-i, w(j-1,i+1), mdw,w(j,i+1), mdw, &
                           sparam)
            endif
         end do

         ! I=column being eliminated.
         ! Test independence of incoming column.
         ! Remove any redundant or dependent equality constraints.
         if (.not.dwnlt2(me, mend, ir, factor,tau,scale,w(1,i))) then
            jj = ir
            do ir=jj,me
               call dcopy (n, [0.0_wp], 0, w(ir,1), mdw)
               rnorm = rnorm + (scale(ir)*w(ir,n+1)/alsq)*w(ir,n+1)
               w(ir,n+1) = 0.0_wp
               scale(ir) = 1.0_wp
               ! Reclassify the zeroed row as a least squares equation.
               itype(ir) = 1
            end do
            ! Reduce ME to reflect any discovered dependent equality
            ! constraints.
            me = jj - 1
            exit
         endif
      end do
   endif

   ! Try to determine the variables KRANK+1 through L1 from the
   ! least squares equations.  Continue the triangularization with
   ! pivot element W(ME+1,I).
   if (krank<l1) then
      recalc = .true.

      ! Set FACTOR=ALSQ to remove effect of heavy weight from
      ! test for column independence.
      factor = alsq
      do i=krank+1,l1

         ! Set IR to point to the ME+1-st row.
         ir = me+1
         lend = l
         mend = m
         call dwnlt1 (i, l, m, ir, mdw, recalc, imax, hbar, h, scale, w)

         ! Update column SS and find pivot column.
         call dwnlt3 (i, imax, m, mdw, ipivot, h, w)

         ! Perform column interchange.
         ! Eliminate I-th column below the IR-th element.
         do j=m,ir+1,-1
            if (w(j,i)/=0.0_wp) then
               call drotmg (scale(j-1), scale(j), w(j-1,i), w(j,i), sparam)
               w(j,i) = 0.0_wp
               call drotm (n+1-i, w(j-1,i+1),  mdw, w(j,i+1), mdw, sparam)
            endif
         end do

         ! Test if new pivot element is near zero.
         ! If so, the column is dependent.
         ! Then check row norm test to be classified as independent.
         t = scale(ir)*w(ir,i)**2
         indep = t > (tau*eanorm)**2
         if (indep) then
            rn = 0.0_wp
            do i1=ir,m
               do j1=i+1,n
                  rn = max(rn, scale(i1)*w(i1,j1)**2)
               end do
            end do
            indep = t > rn*tau**2
         endif

         ! If independent, swap the IR-th and KRANK+1-th rows to
         ! maintain the triangular form.  Update the rank indicator
         ! KRANK and the equality constraint pointer ME.
         if (.not.indep) exit
         call dswap(n+1, w(krank+1,1), mdw, w(ir,1), mdw)
         call dswap(1, scale(krank+1), 1, scale(ir), 1)

         ! Reclassify the least square equation as an equality
         ! constraint and rescale it.
         itype(ir) = 0
         t = sqrt(scale(krank+1))
         call dscal(n+1, t, w(krank+1,1), mdw)
         scale(krank+1) = alsq
         me = me+1
         krank = krank+1
      end do
   endif

   ! If pseudorank is less than L, apply Householder transformation.
   ! from right.
   if (krank<l) then
      do j=krank,1,-1
         call dh12 (1, j, krank+1, l, w(j,1), mdw, h(j), w, mdw, 1, &
                   j-1)
      end do
   endif

   niv = krank + nsoln - l
   if (l==n) done = .true.

   ! End of initial triangularization.
   idope(1) = me
   idope(2) = krank
   idope(3) = niv

   end subroutine dwnlit