c
c*******************************************************************************
c
      subroutine MUSQ ( A,      maxrow, nr,    nc,     bmiAp,  wt,
c
     1                  r0,     gamma,  jb,    jnb,    WRK1,   WRK2,
c
     2                  maxrp1, work1,  work9, work10, work11, work12,
c
     3                  maxwrk, iprint, ipr,   mou,    iw10,   
c
     4                  simint, normal, zmin,  rsdwrk,  
c
     5                  elmin,  amusq, jmin,   kmin )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      external  SCOPY
c
      logical  ipr, simint, normal
c
      dimension  A(maxrow,nc),      wt(nc), 
c
     1           WRK1(maxrp1,nc),   WRK2(maxrp1,nc),  work1(maxwrk),
c
     2           work9(maxwrk),     work10(maxwrk),   work11(maxwrk),
c
     3           work12(maxwrk),    rsdwrk(maxrp1)
c
      dimension  zmin(nc), bmiAp(nr)
c
      integer  jb(nc), jnb(nc), jmin(nc)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  -------------
c  Given inputs:
c  -------------
c
c       A = nr x nc kernel matrix  (assumed to be already scaled)
c
c   bmiAp = The nr-vector  b - A*p  
c
c           where  b = scaled obs. vector in  b = Ax* + eta ,  eta ~ N(0, S**2)
c                  p = lower bound vector for x*
c
c      r0 = minimum sum of squared residuals for the problem
c
c                         {                          }
c                r0 = min { (b - A*x)**T * (b - A*x) }
c                      x  {                          }
c
c   gamma = confidence level parameter for desired conficence intervals
c
c                     { kappa (normal distr.), for one-at-a-time intervals
c             gamma = {
c                     { mu   (Chi-sq. distr.), for simultaneous intervals
c
c  --------------------                
c  This subr. computes:
c  --------------------  
c
c    zmin = solution of constrained least squares problem
c
c                    {                                  |          }
c       L_min =  min { (bmiAp - A*z)**T * (bmiAp - A*z) | z .ge. 0 }
c                 z  {                                  |          }
c  and
c  
c   amusq = mu**2 = L_min + gamma**2 , for one-at-a-time intervals
c
c  or
c
c   amusq = mu**2 = r0   +  gamma**2 , for simultaneous intervals
c
c  It uses subr. LS to do it by setting the equality constraint to 
c          the last row of  A  and  bmiAp  with weight = 1.0 .
c
c  It saves the basis column numbers corresponding to  zmin  
c                                       in the array  jmin(j), j = 1, kmin. 
c        & the non-basis column numbers in the array  jmin(j), j = kmin+1, n
c
c  -----------
c  It returns:
c  -----------
c
c    zmin = the solution vector (size = nc)
c
c  rsdwrk = residual vector (size = nr)
c
c   elmin = L_min
c
c   amusq = described above
c
c    jmin = described above
c
c    kmin = number of columns in the basis set (see above description of jmin)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  Solve the non-negatively constrained least squares problem
c
      kk = 0
      do 100 j=1,nc
       wt(j) = A(nr,j)
       jnb(j) = j
  100 continue
      row1wt = 1.0
c
c   --------------------------------------------------------------
      call LS ( A,      maxrow, nr-1,   nc,     bmiAp,  bmiAp(nr),
c
     1          wt,     row1wt, iprint, ipr,    mou,   
c
     2          kk,     jb,     jnb,    WRK1,   nr+1,   
c
     3          WRK2,   nr+1,   zmin,   rsdwrk, rnorm,  work11,
c
     4          nbch,   icycle, work12, work1,  work10, work9   )
c   --------------------------------------------------------------    
c
      if ( icycle .eq. 1 ) then
                           print 101
                  if (ipr) write (mou,101)
c
                  kk = 0
c
      end if
c
  101 format (/' '/' Cycling occurred in initial constrained point ',
     1             'estimation problem ( calc. of zmin = xmin - p )'/)
c
c    ---------------------------------------------
      call VECPR ( zmin, nc, 'z_min', iw10, mou )
c    ---------------------------------------------
c
      if ( iprint .gt. 0 ) then
c           -----------------------------------------------------------
             call VECPR ( rsdwrk, nr, '(b-A*p) - A*z_min', iw10, mou )
c           -----------------------------------------------------------
      end if
c
      if (iprint .ge. 4)   then
c           ---------------------------------------------------
             call VECPR ( work11, nc, 'grad', iw10, mou )
c           ---------------------------------------------------
      endif
c
c++++++++++++++++++++ Set up the basis and non-basis sets  +++++++++++++++++++++
c
      kmin = kk
c
      call SCOPY ( kk, jb, 1, jmin, 1 )
c
      if ( kk .lt. nc ) call SCOPY ( nc-kk, jnb, 1, jmin(kk+1), 1 )
c
      elmin = rnorm*rnorm + rsdwrk(1)*rsdwrk(1)
c
c          -- Note--   elmin = L_min
c
      if ( iprint .gt. 0 ) 
c                 -----------------------------------------------------
     1             call BASPR ( jb, kk, jnb, nc-kk, iprint, ipr, mou )
c                 -----------------------------------------------------
c
c+++++++++++ Set value of  mu**2  for desired confidence intervals +++++++++++++
c
      if ( simint ) then
                        amusq = r0    + gamma**2
              else
                        amusq = elmin + gamma**2
      end if
c
                   print 70, r0, elmin, gamma, amusq
      if (ipr) write(mou,70) r0, elmin, gamma, amusq
   70 format (' '/' '/' --------------------------------------'
     1        /' '/' ', 8x, '   r0 =',1pe21.14
     2        /' '/' ', 8x, 'L_min =',1pe21.14
     3        /' '/' ', 8x, 'gamma =',1pe21.14
     4        /' '/' ', 8x, 'mu**2 =',1pe21.14 /' '/
     5        ' ----------------------------------------' )
c
      if ( simint ) then
                  if ( normal ) then
                                         print 75
                                if (ipr) write (mou,75)
                        else
                                         print 76
                                if (ipr) write (mou,76)
                  end if
      end if
c
   75 format (' '/' '/' Will calculate simultaneous intervals ', 
     1                     'with   mu**2 = (gamma_m)**2' )
c
   76 format (' '/' '/' Will calculate simultaneous intervals ', 
     1                     'with   mu**2 = r0 + gamma**2' /' '/' ',
     2               46x, 'gamma**2 = n / ( 1.0 - (alpha)**(1/n) )' ) 
c
c
c++++ For simultaneous intervals, check for intersection with pos. orthant +++++
c
      if ( simint .and. (amusq .lt. elmin) ) then
c
                                   print 108, elmin
                     if (ipr) write (mou,108) elmin
c
                                   print 106
                     if (ipr) write (mou,106)
c
c            ++++++ 
              stop
c            ++++++ 
c
      end if
c
  108 format (/ ' '/' '/' '/9x,' L_min =',1pe21.14 )
c
  106 format (/' '/' '/' *************************** Warning ***', 
     1           '*******************'/' '/' '/' The problem is ',
     2           'inconsistent -- mu**2 .lt.  L_min'/' '/' '/ 
     3           ' ************************************************', 
     4           '**********' )
c
c
c+++++++ For one-at-a-time intervals, reset   mu**2 =  L_min + kappa**2 ++++++++
c
      if ( .not. simint ) then
                                        print 107
                          if (ipr) write (mou,107)
      end if
c
  107 format (/' '/' '/' ****************************************'
     1            /' '/'  Have set  mu**2 =  L_min + kappa**2 '
     2            /' '/'            for one-at-a-time intervals'
     3            /' '/' ****************************************' ) 
c
c
      return
c
      end
