c
c*******************************************************************************
c
      subroutine LS ( A,      ldA,    m,      n,     b,      phi, 
c
     1                w,      alpha,  iprint, ipr,   mou,
c
     2                k,      jb,     jnb,    Q,     ldQ,   
c
     3                R,      ldR,    x,      resid, rnorm,  g,
c
     4                ntbch,  infor,  bmod,   rr,    xp,     z    )
c 
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      logical  ipr
c
      integer  ldA, m,   n,     iprint, mou, k, jb(n), jnb(n),  
     1         ldQ, ldR, ntbch, infor
c
      real  A(ldA,n), b(m),     phi,        w(n),  alpha, Q(ldQ,n),
     1      R(ldR,n), x(n),     resid(m+1), rnorm, g(n),  bmod(m+1),
     2      rr(m+1),  xp(n+1),  z(n+1) 
c
      external  ARRAY,  BASPR,  GSDEL,  GSDERR, GSERR,  GSINS,  R1MACH,
     1          SCOPY,  SDOT,   SNRM2,  SSCAL,  STRSL,  VECPR
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  A subroutine for solving the m x n constrained least squares problem
c
c             (  ||         ||2                              )
c         min (  || b - A x ||  :  w**T x = phi ,  x .ge. 0  ) ,
c          x  (  ||         ||                               )
c
c  by solving the (m+1) x n augmented, constrained least squares problem
c
c     (  alpha w**T  )    .  ( alpha phi )
c     (--------------) x  =  (-----------) ,  with  x .ge. 0 ,
c     (      A       )       (     b     )
c
c  with alpha chosen very large to make the first equation an exact equality.
c
c
c  It is also used to solve the m x n constrained problem
c
c             (  ||         ||2              )
c         min (  || b - A x ||  :  x .ge. 0  ) ,
c          x  (  ||         ||               )
c
c  by taking  alpha = 1.0 ,
c              w**T = last row of A ,
c               phi = last element of b ,
c  i.e.,          
c               (  A_1   )            (  b_1   )
c       if  A = (--------) , and  b = (--------) ,  then solve
c               ( a_m**T )            ( beta_m )
c
c       ( alpha a_m**T )    .  ( alpha beta_m )
c       (--------------) x  =  (--------------) ,  with  x .ge. 0 ,  alpha = 1.0
c       (     A_1      )       (     b_1      )
c
c-------------------------------------------------------------------------------
c 
c  INPUT
c  -----
c
c  A       Real matrix ( m x n)
c          Least squares model matrix
c
c  ldA     Integer
c          Row dimension of Fortran array A in the calling program
c
c  m       Integer
c          No. of rows in matrix A (and vector b)
c
c  n       Integer
c          No. of columns in matrix A (and no. elements in solution vector x)
c
c  b       Real vector (length = m)
c          Measured r.h.s. vector for least squares model
c
c  phi     Real scalar
c          Value of the linear equality constraint
c
c  w       Real vector (length = n)
c          Constraint vector for linear equality constraint
c
c  alpha   Real scalar   
c          Weight to make the first equation an exact equality
c
c  iprint  Integer
c          Print parameter determining amount of printed output, 
c          ranging from 0  (for minimal output)  to  5 (for copious output) 
c          0, 1, 2, or 3 are reasonable choices
c
c  ipr     Logical
c          Flag for output to a file as well as to standard output unit.
c          = .true.   means write duplicate output to a file for later printing,
c          = .false.  means write output only to standard output unit
c
c  mou     Integer
c          Unit number for duplicate output file (not used if ipr = .false.)
c
c-------------------------------------------------------------------------------
c
c  The method consists in recursively updating/downdating a Gram-Schmidt (QR)
c  factorization of an (m+1) x k basis submatrix A_b obtained from the 
c  (m+1) x n system matrix
c                          ( alpha w**T )
c                          ( ---------- )
c                          (     A      ) 
c  by deleting the columns corresponding to zero elements of the current 
c  estimate of the solution vector x.
c
c  The factorization can be written
c                                    A_b = Q R , 
c  where  
c                Q**T Q = I  , and  R  is k x k upper triangular.
c
c  It is not necessary to explicitly form the matrix A_b, but it is necessary
c  to form and maintain two indexing sets: 
c
c   ( jb(j), j=1,k )     Set of column indices of columns currently in basis,
c                        i.e., the indices of the non-zero elements of x,
c
c   ( jnb(j), j=1,n-k )  Set of column indices of columns not in the basis,
c                        i.e., the indices of the zero elements of x.
c
c  If the program is called with  k = 0  a completely new basis is constructed,
c  but if it is called with  k .gt. 0  it proceeds by updating/downdating the 
c  final basis from a preceding call.
c
c  The updating/downdating iteration is continued until the solution vector and
c  the gradient vector satisfy the saddle point theorem and thus are optimal.
c
c  The algorithm may cycle.  If cycling occurs when the initial basis set was
c  non-empty (k > 0), then the basis is cleared, k set to 0, and the algorithm
c  is restarted.  If cycling occurs when the initial basis set was empty (k = 0)
c  then the subroutine sets an information parameter and returns to the calling
c  program with all output parameters equal to their values when cycling was
c  detected.
c
c-------------------------------------------------------------------------------
c
c  INPUT and OUTPUT
c  ----------------
c
c  k       Integer ( k .le. n )
c          Number of columns in the basis
c
c  jb      Integer vector (max. length = n)
c          Column numbers of the basis columns
c
c  jnb     Integer vector (max. length = n)
c          Column numbers of columns not in the basis.
c
c  Q       Real matrix ( m+1 x n )
c          The m+1 x k matrix with orthogonal columns in the QR factorization
c          
c  ldQ     Integer ( ldQ .ge. m+1 )
c          Row dimension of Fortran array Q in the calling program
c
c  R       Real matrix ( n x n )
c          The k x k upper triangular matrix in the QR factorization
c
c  ldR     Integer ( ldR .ge. n )
c          Row dimension of Fortran array R in the calling program
c
c
c  OUTPUT
c  ------
c
c  x       Real vector (length = n)
c          The estimated solution vector
c
c  resid   Real vector (length = m+1)
c          Vector of residuals for the m+1 x n augmented least squares problem,
c          i.e.,
c                         ( alpha phi )   ( alpha w**T )
c                 resid = (-----------) - (------------) x
c                         (     b     )   (     A      )
c
c  rnorm   Real
c          Sum of squared residuals for original m x n (not augmented) problem,
c          i.e., 
c                 rnorm = || b - A x ||**2 .
c
c  g       Real vector (length = n)
c          Negative of the gradient vector evaluated at solution vector x,
c          i.e.,
c                 g = -2 ( alpha w , A**T ) resid .
c
c  ntbch   Integer
c          Number of basis changes used in constructing the solution vector x.
c
c  infor   Integer
c          Information parameter indicating the status of the calculation.
c          = 0  means no errors.  
c          = 1  means cycling occurred.  This usually means that the augmented
c               matrix has rank .lt. n  or is close to such a rank deficient 
c               matrix.  Perturbing the matrix slightly may enable the 
c               subroutine to terminate successfully.
c          = 2  means that termination occurred because all candidate vectors
c               for adding to the basis were immediately deleted again.  If the
c               maximum element of  g  is  **small**  and positive, then there
c               is no problem.  Otherwise, the matrix may have rank .lt. n.
c
c
c  WORKING ARRAYS
c  --------------
c
c  bmod    Real vector (length = m+1)
c          Used to form the augmented r.h.s. vector  ( alpha phi )
c                                                    ( --------- )
c                                                    (     b     )         
c  rr      Real vector (length m+1) 
c          Working vector used by subroutine GSINS. 
c
c  xp      Real vector (length = n+1)
c          Working vector used in forming candidate solution vectors.
c
c  z       Real vector (length = n+1)
c          Working vector used in forming candidate solution vectors.
c
c-------------------------------------------------------------------------------
c
c  SUBROUTINES CALLED
c  ------------------
c
c  R1MACH  Subroutine to return the machine constants,
c          e.g.,
c                 epsilon_mach = R1MACH ( 4 )
c
c  SCOPY   Basic linear algebra subroutines (BLAS)
c  SDOT    
c  SNRM2   see J.J. Dongarra, C.B. Moler, J.R. Bunch, and G.W. Stewart,
c  SSCAL       Linpack Users' Guide, SIAM (1979) Appendix A
c  STRSL
c
c  GSINS   Subroutine for updating the Gram-Schmidt factorization
c
c  GSERR   Subroutine to print error messages from GSINS
c
c  GSDEL   Subroutine for downdating the Gram-Schmidt factorization
c
c  GSDERR  Subroutine to print error messages from GSDEL
c
c  BASPR   Subroutine to print out the current basis indexing set (jb) 
c          and the current non-basis indexing set (jnb)
c
c  VECPR   Subroutine to print out a vector array
c
c  ARRAY   Subroutine to print out a matrix array 
c
c-------------------------------------------------------------------------------
c
c    This subroutine is a modified version of Hanson and Lawson's "NNLS".  For 
c    further information, see:
c
c    C.L. Lawson and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall,
c    (1974) Chapt. 23.
c
c    K.H. Haskell and R.J. Hanson, "An algorithm for linear least squares 
c    problems with equality and nonnegativity constraints," Math. Prog. 21 
c    (1981) pp 98-118.
c
c    D.P. O'Leary and B.W. Rust, "Confidence intervals for inequality
c    constrained least squares problems, with applications to ill-posed
c    problems," Siam J. on Sci. and Stat. Comput. 7 (1986) pp 473-489.
c
c-------------------------------------------------------------------------------
c
c  LOCAL VARIABLES
c  ---------------
c
      logical  xpos, zpos
c
      integer  i,      ib,     info,   iw10,   j,      jd,  
     1         jmax,   jn,     js,     jss,    jx,
     2         kinitl, km1,    mbc,    mbcp10, mp1,    nbch,   nmk,
     3         nmkmrk, ntcycl
c
      real  absR,   alfa,   almax,  altest, bnorm,  eps,    gmax,
     1      Rcond,  rho,    rhofac, rhom,   rhomin, rhopro, rhoprv,
     2      Rmax,   Rmin,   rratio, twoeps, xpnrm   
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c                       I n i t i a l i z a t i o n
c
c
      if (    ipr    ) iw10 = 1
      if ( .not. ipr ) iw10 = 0
c
      eps = R1MACH( 4 )
      twoeps = 2.0*eps
c
      kinitl = k
      ntcycl = 0
c
      mbc = 3*(m + n)
c                       MAX. NO. basis changes before assuming cycling
      mbcp10 = mbc + 10
c                       WHEN CYCLING, print basis for 10 iterations, then stop,
      ntbch = 0
      nbch = 0
c
      if ( iprint .ge. 4 ) then 
                           print 4, phi, k
                if ( ipr ) write (mou,4) phi, k
      end if
    4 format (' '/' '/' Enter LS with   phi =', 1pe14.6, '   k =', i4 )
c
      mp1 = m + 1
c
      bmod(1) = alpha*phi
      call SCOPY ( m, b, 1, bmod(2), 1 )
c                                               ( alpha phi )
c                                        bmod = (-----------)
c                                               (     b     )
c
      if ( iprint .gt. 4 ) call VECPR ( bmod, mp1, 'bmod', iw10, mou )
c
      if ( alpha .ne. 1.0 ) bnorm = SNRM2( m,   bmod(2), 1 )
      if ( alpha .eq. 1.0 ) bnorm = SNRM2( mp1, bmod,    1 ) 
c                                                            bnorm = || b ||
c
c--------------
  50  continue
c-------------- TRANSFER BACK TO HERE with k = 0 after cycling is detected 
c               (from statements just below comments numbered 1160 and 1180).
c
      rhomin = R1MACH( 2 )
      rho = sqrt ( rhomin )
      rhoprv = 2.0*rho
      rratio = 0.0
c                   -- rhomin, rhoprv, and rratio will be used to test 
c                                                                 for cycling
      if ( k .eq. 0 )  then
c                          INITIALIZE BASIS if cycling or no prev. call to LS 
              nmk = n
              do 60 j = 1, nmk
                 jnb(j) = j
   60         continue
              call SCOPY ( n, 0.0, 0, xp, 1 )
              call SCOPY ( mp1, bmod, 1, resid, 1 )
              rho = SNRM2( mp1, bmod, 1 )
              rhoprv = 2 * rho
       else
c                      UPDATE BASIS from previous call to LS
c----------------------
  100         continue
c---------------------- TRANSFER BACK TO HERE (from statement just below
c                       comment numbered 150) until all basis columns
c                       correspond to pos. elements in solution vector
c
              call SCOPY ( mp1, bmod, 1, resid, 1 )
c            ----------------------------------------------------
              call GSINS ( Q, ldQ, mp1, k, resid, xp, rr, info )
c            ----------------------------------------------------
c                             CALCULATES PROJECTION of r.h.s. on cols. of Q.  
c                             Returns result, normalized to length 1, in resid,
c                             with norm in xp(kp1) and 
c                                ( xp(1) , ... , xp(k) )**T = Q**T (r.h.s.)
              rhoprv = rho
              rho = xp(k+1) 
              nmk = n - k
              if ( iprint .gt. 0 ) then
                      if ( info .ne. 0 )
c                              -----------------------------------------  
     1                          call GSERR ( info, n+1, k+1, ipr, mou )
c                              -----------------------------------------
                      if ( nbch .ge. mbc ) then
                                               print 160, nbch, rho
                                 if (ipr) write (mou,160) nbch, rho 
                      endif
              end if
c            ---------------------------------------
              call SSCAL ( mp1, xp(k+1), resid, 1 )
c            --------------------------------------- RESTORES LENGTH of resid
              if (iprint .ge. 4) 
     1                     call VECPR ( resid, mp1, 'rsid', iw10, mou )
c            ----------------------------------------
              call STRSL ( R, ldR, k, xp, 01, info )
c            ---------------------------------------- SOLVES TRIANGULAR SYSTEM
c                                                       R xp = Q**T (r.h.s.)
              if (iprint .ge. 4)
     1                         call VECPR ( xp, k, 'xp  ', iw10, mou ) 
              if ( (info .gt. 0) .and. (iprint .gt. 0) ) then 
                                                 print 170, info, k
                                   if (ipr) write (mou,170) info, k
              end if
c
c             -------- CHECK SOLUTION VECTOR for non-positive elements ---------
c
              xpos = .true.
              do 120 j=k,1,-1
                 if ( xp(j) .le. 0.0 ) then
c                                          REMOVE COLUMN j from basis
                         xpos = .false.
                         nbch = nbch + 1
c                       ------------------------------------------------
                         call GSDEL ( Q, ldQ, mp1, k, R, ldR, j, info )
c                       ------------------------------------------------
                         if ( (info .ne. 0) .and. (iprint .gt. 0) )
c                                       --------------------------------
     1                                   call GSDERR ( j, k, ipr, mou )
c                                       --------------------------------
                         k = k - 1
                         nmk = n - k
                         jnb(nmk) = jb(j)
                         do 110 ib=j,k
                            jb(ib) = jb(ib+1)
  110                    continue
                         if ( nbch .ge. mbc ) then
c                                                 ASSUME algorithm is cycling.
c                        ------------------------------------------------
                          call BASPR (jb, k, jnb, nmk, iprint, ipr, mou)
c                        ------------------------------------------------
c                              PRINT BASIS SET for 10 iterations,
c                              then stop trying.
c 115                                                  ------------
                               if ( nbch .ge. mbcp10 )  go to 1150
c                                                      ------------
                         end if
                         if (iprint .ge. 4) 
     1                    call BASPR (jb, k, jnb, nmk, iprint, ipr, mou)
                         if (iprint .ge. 5)
     1                    call ARRAY (Q, mp1, k, ldQ, 'Q   ', iw10, mou)
                 end if
  120         continue
c
c             HAVE NEW BASIS obtained by dropping some cols. of previous basis.
c             IF new basis contains cols. which give neg. solution components,
c             REPEAT ABOVE PROCEDURE until all cols. give pos. solution elements
c             or all cols. removed from basis.
c 150                                              -----------
              if ( (.not. xpos) .and. (k .gt. 0) )  go to 100
c                                                  -----------  
      end if
c
  160 format (' '/' nbch =', i5, 4x, 'rho =', 1pe14.6 )
  170 format (' '/' info =', i4, '  returned by STRSL to ',
     1                           'LS in calc. of (xp),  k =', i4 )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c             B a s i c   I t e r a t i o n   S t a r t s   H e r e
c
c--------------
  200 continue
c-------------- TRANSFER BACK TO HERE from bottom of iteration loop (from 
c               statements just below comments numbered 750 and 1125)
c
c----------------- INITIALIZE new feasible solution vector ---------------------
c
      call SCOPY ( n, 0.0, 0, x, 1 )
      do 210 j=1,k
         x( jb(j) ) = xp(j)
  210 continue
      if (iprint .ge. 4) call VECPR ( x, n, 'x   ', iw10, mou )
c
      if ( (rho .eq. rhomin) .and. (ntcycl .ge. 2) ) then
                                                     infor = 1
                              if ( iprint .gt. 0 ) print 230, k, rhomin
               if( ipr .and. (iprint .gt. 1) ) write(mou,230) k, rhomin
               call BASPR ( jb, k, jnb, nmk, iprint, ipr, mou )
c 220         ------------
               go to 1200
c             ------------ CYCLING HAS OCCURED TWICE.  Stop trying.
      end if
  230 format (/' k , rhomin = ', i3, 1pe15.6 )
c
c
c----------- DETERMINE GRADIENT VECTOR at new feasible solution ----------------
c
      if (iprint .gt. 4) call VECPR ( resid, mp1, 'rsid', iw10, mou )
c
      do 240 j=1,n
         g(j) = alpha*w(j)*resid(1) + SDOT( m, A(1,j), 1, resid(2), 1 )
  240 continue
c                  (g) = -(1/2)*(gradient) 
c
      if (iprint .ge. 4) call VECPR ( g, n, 'g   ', iw10, mou )
c
c       
c----------------- CHECK for CYCLING and for CONVERGENCE -----------------------
c
      if ( rho .gt. rhoprv ) then 
c                                 NEW SUM SQUARED RESIDUALS larger than previous
             rhofac = rho / rhoprv
c 250                                  ------------
             if ( rhofac .eq. rratio )  go to 1140
c                                      ------------ CYCLING, almost surely when
c                                                   same increase occurs twice
             rratio = rhofac
      end if
c
      if ( alpha .eq. 1.0 ) rhom = rho
      if ( alpha .ne. 1.0 ) rhom = SNRM2( m, resid(2), 1 )
c                                                           rhom = || b - A x ||
      if ( rhom .lt. bnorm*twoeps ) then
c                                                    || b - A x ||
c                                    CONVERGENCE if  ------------- < 2 eps_mach
c                                                      || b ||
              if ( iprint .gt. 0 ) then
                                   print 280, rhom, bnorm, twoeps
                      if (ipr) write(mou,280) rhom, bnorm, twoeps
              end if
              infor = 0
c 275        ------------
              go to 1200
c            ------------ STOP ITERATING because for all practical purposes
c                                           A x = b
      end if
  280 format ( /' *-- norm(A*x - b)/norm(b)  =', 1pe11.3,
     1                             '/', 1pe10.3, ' .lt.', 1pe10.3 / )
c
      nmkmrk = nmk
c                  NUMBER OF COLUMNS not yet tried as additions to current basis
c--------------
  300 continue
c-------------- TRANSFER BACK TO HERE (from statement just below comment
c               numbered 650) when candidate solution vector z has a non-pos. 
c               element corresponding to a column that would have been added 
c               to the basis.  
c
      if (nmkmrk .eq. 0) then
                         infor = 2 
c 350        ------------
              go to 1200
c            ------------
c                       CONVERGENCE if current basis contains all of the columns
      endif
c
c
c------- FIND MAXIMUM ELEMENT of (g), i.e., the min. elt. of gradient ----------
c
      jmax = 1
      jx = jnb(1)
      gmax = g(jx)
      do 400 j=2,nmkmrk
         jn = jnb(j) 
         if ( g(jn).gt.gmax ) then
                jmax = j
                jx = jn 
                gmax = g(jx)
         endif
  400 continue
      if (gmax .le. 0.0) then
                         infor = 0 
c 450       ------------
             go to 1200 
c           ------------
c                      CONVERGENCE if all non-basis gradient components are 
c                      non-positive, i.e., all non-basis elements of (g) > 0.0.
      endif
c
c
c------- DECREASE OBJECTIVE FUNCTION by adding column jx to the basis ----------
c         ( g(jx) = gmax > 0.0, so element jx of gradient is < 0.0 )
c
      nmk = nmk - 1 
      nmkmrk = nmkmrk - 1
      k = k + 1
      nbch = nbch + 1
      do 500 j=jmax,nmk
          jnb(j) = jnb(j+1)
  500 continue
      jb(k) = jx
      xp(k) = x(jx) 
c
      if ( nbch .ge. mbc ) then
                call BASPR ( jb, k, jnb, nmk, iprint, ipr, mou )
c 510                                   ------------
                if ( nbch .ge. mbcp10 )  go to 1150
c                                       ------------ CYCLING, almost surely
      end if
      if ( iprint .ge. 4 ) 
     1                call BASPR ( jb, k, jnb, nmk, iprint, ipr, mou )
c
      Q(1,k) = alpha*w(jx)
      call SCOPY( m, A(1,jx), 1, Q(2,k), 1 )
c    -------------------------------------------------------------
      call GSINS ( Q, ldQ, mp1, (k-1), Q(1,k), R(1,k), rr, info )
c    -------------------------------------------------------------
      if ( (info.ne.0) .and. (iprint.gt.0) ) 
c                               --------------------------------------
     1                           call GSERR ( info, jx, k, ipr, mou )
c                               --------------------------------------
c
      if (iprint .ge. 5) then
                        call ARRAY ( Q, mp1, k, ldQ, 'Q   ', iw10, mou )
                        call ARRAY ( R, k,   k, ldR, 'R   ', iw10, mou )
      end if
      if (iprint .ge. 4)  then
             Rmax = 0.0
             Rmin = R1MACH( 2 )
             do 520 i=1,k
                absR = abs( R(i,i) )
                if ( absR .lt. Rmin ) Rmin = absR
                if ( absR .gt. Rmax ) Rmax = absR
  520        continue 
             Rcond = Rmin/Rmax
             print 525, Rcond
             if (ipr) write (mou,525) Rcond
      end if
  525 format (' '/' '/'  min |R(i,i)|', 
     1               /' --------------  .le. ', 1pe13.5,
     2               /'  max |R(i,i)|' )
c
c
c----- SOLVE FOR CANDIDATE SOLUTION VECTOR to be stored temporarily in (z) -----
c
c--------------
  600 continue
c-------------- TRANSFER BACK TO HERE after one or more columns eliminated
c               from the basis (transfer from statement just below comment
c               numbered 1130)
c
      call SCOPY ( mp1, bmod, 1, resid, 1 ) 
c    ---------------------------------------------------
      call GSINS ( Q, ldQ, mp1, k, resid, z, rr, info )
c    ---------------------------------------------------         ( alpha phi )
c                                            NOW z contains Q**T ( --------- )
c                                                                (     b     )
      if ( (info .ne. 0) .and. (iprint .gt. 0) )
c                             -----------------------------------------
     1                         call GSERR ( info, n+1, k+1, ipr, mou ) 
c                             -----------------------------------------
      rhopro = rhoprv
      rhoprv = rho
      rho = z(k+1)
c
      if ( (nbch .ge. mbc) .and. (iprint .gt. 0) ) then
                                                   print 160, nbch, rho
                                      if (ipr) write(mou,160) nbch, rho 
      end if
      if (iprint .ge. 4) then
                         call VECPR ( resid,  mp1, 'bprp', iw10, mou )
                         call VECPR ( z,      k,   'QT*b', iw10, mou )
                         call VECPR ( z(k+1), 1,   'rho ', iw10, mou )
      end if
c    --------------------------------------
      call SSCAL ( mp1, z(k+1), resid, 1 )
c    -------------------------------------- RESTORES LENGTH of resid
      if (iprint .ge. 4) 
     1                  call VECPR ( resid, mp1, 'rsid', iw10, mou )
c    ---------------------------------------
      call STRSL ( R, ldR, k, z, 01, info )
c    ---------------------------------------                    ( alpha phi )
c                          SOLVES TRIANGULAR SYSTEM  R z = Q**T ( --------- )
c                                                               (     b     )
      if ( (info .gt. 0) .and. (iprint .gt. 0) ) then 
                                                 print 620, info, k
                                   if (ipr) write (mou,620) info, k
      end if
  620 format (' '/' info =', i4, '  Returned by STRSL to ',
     1                             'LS in calc. of (z),  k =', i4 )
      if ( iprint .ge. 4 ) call VECPR ( z, k, 'z   ', iw10, mou )
c
      if ( z(k) .le. 0.0 ) then
c                               NEW SOLUTION ELEMENT z(k) is not positive
c                                     DROP COLUMN jb(k) from basis
                        nmk = nmk + 1
                        jnb(nmk) = jb(k)
                        k = k - 1
                        rho = rhoprv 
                        rhoprv = rhopro
c 650        -----------
              go to 300
c            ----------- GO BACK to search for new basis column candidate
      end if
c
c
c------------- CHECK WHETHER ALL ELEMENTS of (z) are positive ------------------
c
      zpos = .true.
      do 700 j=1,k
         if ( zpos ) jd = j
         if ( z(j) .le. 0.0 )  zpos = .false.
  700 continue
c
      if ( zpos ) then
c                                                      ( z )
c             ------------- ALL z(j) > 0.0, SET (xp) = (---) ---------------
c                                                      ( 0 )
              call SCOPY ( k, z, 1, xp, 1 ) 
              if ( k .ne. n )  call SCOPY ( nmk, 0.0, 0, xp(k+1), 1 )       
              if (iprint .ge. 4) call VECPR ( xp, k, 'xp  ', iw10, mou )
              if ( ((ntcycl .eq. 1) .or. (kinitl .eq. 0)) 
     1                           .and. (rho .lt. rhomin) )  rhomin = rho
c 750         -----------
               go to 200
c             -----------  GO BACK to calculate new gradient
      end if 
c                                                            xp(j)
c-------  z(jd) .le. 0.0, FIND THE z(j) WHICH MINIMIZES  --------------  -------
c                                                         xp(j) - z(j)
c 
c         ADJUST (xp) SO xp(j) = 0.0 AND REMOVE corresponding col. from
c         basis.  If more than one elt. of (z) had same associated min.
c         value, remove all the corresponding columns  from the basis.
c
c                               (    xp(j)     :               )
c              SET  almax = min ( ------------ : z(j) .le. 0.0 ) ,
c                               ( xp(j) - z(j) :               )
c
c         smallest value such that  (xp) = (xp) + almax*( (z) - (xp) )
c         zeroes an element of (xp).
c
      almax = xp(jd) / ( xp(jd) - z(jd) )
      do 800 j=jd+1,k
         if ( z(j) .le. 0.0 ) then
                              altest = xp(j) / ( xp(j) - z(j) )
                              if ( altest .lt. almax ) then
                                                       almax = altest
                                                       jd = j 
                              end if
         end if
  800 continue
      alfa = almax 
c
      if (iprint .ge. 4) then
                               print 900, alfa
                 if (ipr) write (mou,900) alfa
      end if
  900 format (' '/' '/' alfa =', 1pe13.5 ) 
c
c
c------------------ SET  (xp) = (xp) + alfa*( (z) - (xp) ) ---------------------
c
c  AND ELIMINATE COLUMN jd FROM BASIS (col. corresponding to alfa=almax) and 
c        all other cols. corresponding to elts. of (xp) driven to zero
c
      do 1000 j=1,k 
         xp(j) = xp(j) + almax*( z(j) - xp(j) )
 1000 continue
      if (iprint .ge. 4) call VECPR ( xp, k, 'xp  ', iw10, mou ) 
c
      jss = 1
      xpnrm = SNRM2( k, xp, 1 )
c
c--------------
 1100 continue
c-------------- TRANSFER BACK TO HERE (from statement just below comment
c               numbered 1115) to eliminate other columns whose corresponding 
c               elements of (xp) were driven to zero.
c
      nbch = nbch + 1
c    -------------------------------------------------
      call GSDEL ( Q, ldQ, mp1, k, R, ldR, jd, info )
c    ------------------------------------------------- DROPS COL. jd from basis
      if ( (info .ne. 0) .and. (iprint .gt. 0) )
c                                     --------------------------------- 
     1                                 call GSDERR ( jd, k, ipr, mou )
c                                     ---------------------------------
      k = k - 1
      nmk = n - k
      jnb(nmk) = jb(jd)
      do 1105 j=jd,k
         jb(j) = jb(j+1)
         xp(j) = xp(j+1)
 1105 continue
      xp(k+1) = 0.0
c
c     CHECK FOR CYCLING
c
      if ( nbch .ge. mbc ) then
            call BASPR ( jb, k, jnb, nmk, iprint, ipr, mou )
c1110                               ------------
            if ( nbch .ge. mbcp10 )  go to 1150
c                                   ------------ ALMOST SURELY CYCLING
      end if
      if ( (iprint .ge. 4) .and. (nbch .lt. mbc) )
     1                  call BASPR ( jb, k, jnb, nmk, iprint, ipr, mou )
      if ( iprint .ge. 5 ) 
     1                  call ARRAY ( Q, mp1, k, ldQ, 'Q   ', iw10, mou )
c
      js = jss
      km1 = k - 1
      do 1120 j=js,km1
         jd = j
         jss = jd
c1115                                ------------
         if ( xp(j) .le. eps*xpnrm )  go to 1100
c                                    ------------  xp(jd) WAS DRIVEN TO ZERO
 1120 continue
c
c     CHECK WHETHER BASIS SET IS EMPTY
c
      if ( k. eq. 0 ) then
                      rho = bnorm
c1125     -----------
           go to 200
c         -----------  BASIS SET EMPTY.  GO BACK to start building a new basis.
      end if
c
      rho = rhoprv 
      rhoprv = rhopro
c
c1130-----------
      go to 600
c    -----------  BASIS SET NOT EMPTY.  GO BACK to find new candidate solution.
c
c
c             E n d   o f   t h e   B a s i c   I t e r a t i o n
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c           C y c l i n g   H a n d l i n g   S t a r t s   H e r e
c
c--------------
 1140 continue
c--------------  TRANSFER HERE (from statement just below comment numbered 250) 
c                after second time sum of sqared residuals exceeds the previous 
c                value by the same factor
c
      if ( iprint .gt. 0 ) then
                           print 1142, rho, rhoprv
              if (ipr) write(mou,1142) rho, rhoprv
      end if
 1142 format (/' LS cycling with  rho / rhoprv  =',
     1            1pe13.6, ' / ', 1pe12.6, '  repeating exactly'/ )
c
c--------------
 1150 continue
c--------------  TRANSFER HERE (from statements just below comments numbered 
c                115, 510 and 1110) when number of basis changes equals or 
c                exceeds maximum number allowed.
c
      ntcycl = ntcycl + 1
      if ( kinitl .eq. 0 )  ntcycl = ntcycl + 1
      if ( ntcycl .eq. 1 ) then
                           if ( iprint .gt. 0 ) then
                                                print 1165
                                   if (ipr) write(mou,1165)
                           end if
                           k = 0
                           ntbch = ntbch + nbch
                           nbch = 0
c1160         ----------
               go to 50
c             ----------  GO BACK to build a new basis
      end if
 1165 format (/' Current problem probably cycling.',
     1                      '  Clear basis set and try again'/)
c
c----- PROGRAM REACHES THIS POINT only after it has cycled more than once ------
c
      if ( iprint .gt. 0 ) then
                               print 1170, rhomin
                  if (ipr) write(mou,1170) rhomin
      end if
 1170 format (/' Current problem probably cycling again.',
     1                '  Stop trying this one'//' rhomin =', 1pe14.6/ )
c
      if ( (ntcycl .le. 2) .and. (rho .ne. rhomin) ) then
c                                                    PROGRAM HAS CYCLED TWICE,
c                                          BUT sum of squared residuals DIFFERS
c                                   FROM value obtained the last time it cycled
                        k = 0
                        ntbch = ntbch + nbch 
                        nbch = 0
c1180        ----------
              go to 50
c            ----------  GO BACK to build a new basis.
      end if
c
c---- PROGRAM REACHES THIS POINT only if it has cycled more than twice or ------
c           cycled twice to same value of sum of squared residuals.
c                  IN EITHER CASE, stop trying this problem.
      infor = 1
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c       C a l c u l a t e  sqrt( sum of residuals**2 )  a n d  R e t u r n
c
c--------------
 1200 continue
c--------------  TRANSFER HERE (from statements just below comments numbered 
c                275, 350, and 450) when convergence detected,
c                or (from statement just below comment numbered 220 and from 
c                directly above) for repeated cycling. 
c
      ntbch = ntbch + nbch
c
      rnorm = SNRM2( m, resid(2), 1 )
c
c    --------
      return
c    --------
c
      end 
c
c
c*******************************************************************************
c
c
      subroutine BASPR ( jb, k, jnb, nmk, iprint, ipr, mou )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      logical ipr
c
      integer  k, jb(k), nmk, jnb(nmk), iprint, mou
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c        THIS SUBROUTINE PRINTS THE BASIS AND NON-BASIS INDEXING SETS
c
c  INPUT   (all unchanged by this subroutine)
c  -----
c
c  jb      Integer vector
c          Column numbers of the basis columns
c
c  k       Integer 
c          Number of columns in the basis
c
c  jnb     Integer vector 
c          Column numbers of columns not in the basis
c
c  nmk     Integer
c          Number of columns not in the basis
c
c  iprint  Integer
c          Print parameter determining amount of printed output, 
c          ranging from 0  (for minimal output)  to  5 (for copious output) 
c          0, 1, 2, or 3 are reasonable choices
c
c  ipr     Logical
c          Flag for output to a file as well as to standard output unit.
c          = .true.   means write duplicate output to a file for later printing,
c          = .false.  means write output only to standard output unit
c
c  mou     Integer
c          Unit number for duplicate output file (not used if ipr = .false.)
c
c-------------------------------------------------------------------------------
c
c  LOCAL VARIABLE
c  --------------
c
      integer  ij
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++c
       if ( iprint .ge. 3 )  then
             print 106
             print 107, ( jb(ij), ij=1,k ) 
             print 108
             print 107, ( jnb(ij), ij=1,nmk )
       endif
c
       if ( (iprint .ge. 2) .and. ipr )  then
             write (mou,106)
             write (mou,107) ( jb(ij), ij=1,k )
             write (mou,108)
             write (mou,107) ( jnb(ij), ij=1,nmk )
       endif
c
  106   format (/' Basis indexing set')
  107   format (' ',10i5)
  108   format (/' Non-basis indexing set'/)
c
c    --------
      return
c    --------
c
      end 
