c
c*******************************************************************************
c
      subroutine BRACLS ( A,      ldA,    n,      mA,
c
     1                    WIND,   ldwind, nw,     xmin,
c
     2                    amusq,  b,      phlo,   phup,
c
     3                    epconv, iprint, ipr,    mou,
c
     4                    philo,  phiup,  cratlo, cratup,
c
     5                    phiest, nbctlo, nbctup, ntotbc,
c
     6                    work,   iwork,  elmin,  jmin,   kmin  )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  Programs written by
c                          Bert W. Rust   and   Dianne P. O'Leary
c                          Applied and Computational Mathematics Division
c                          Center for Computing and Applied Mathematics
c                          National Institue for Standards and Technology
c
c  This program computes confidence intervals for a set of linear functions
c
c                  phi = w**T * x
c
c  where
c                A * x = b + e
c
c  with
c
c       1.  A:   mA x n  known
c
c       2.  w:   n x 1   known
c
c       3.  x:   n x 1   unknown
c
c       4.  b:   mA x 1  known
c
c       5.  e:   mA x 1  normally distributed, mean 0, standard
c                        deviation  the identity matrix.
c
c       6.  x,  for physical reasons, is known to be nonnegative.
c
c  It is the last assumption that makes the problem nonstandard.
c
c  This program computes 100*alpha per cent confidence intervals for
c  phi, possibly for several vectors  w .
c
c  mu**2 = L_min + gamma**2 , 
c
c          where               
c                       {                          |          }
c          L_min =  min { (b - A x)**T * (b - A x) | x .ge. 0 } ,
c                    x  {                          |          }
c
c    *      gamma  is a number which determines the confidence level.
c
c  If    integral of chi-square distribution function with rank(A)
c        degrees of freedom from  0  to  gamma  equals  alpha,
c  Then  there is probability  alpha  that all of the true values
c        of the  phi's  lie within the computed bounds. 
c
c  If    integral of  normal(0,1) distribution function from
c        -gamma  to  gamma  equals  alpha,
c  Then  for each phi, there is probability  alpha  that the true
c        value lies within the computed bounds.
c
c  We will call the vectors  w  "window vectors".  they are stored in the 
c  matrix  WIND .
c
c  The algorithm used is a root-finder  BRAKET  coupled with a quadratic 
c  programming algorithm  LS .  The upper and lower bounds are defined by
c
c  { phup }   { max } {        |                                            }
c  {      } = {     } { w**T x | (b -  A x)**T (b -  A x) <= mu**2 , x >= 0 }
c  { phlo }   { min } {        |                                            }
c                x    
c
c  This problem can be shown to be equivalent to finding the two
c  solutions of
c
c        L(phi) - mu**2 = 0
c
c  where
c                     {                        |                        }
c        L(phi) = min { (b - A x)**T (b - A x) | w**T x = phi ,  x >= 0 }
c                  x  {                        |                        }
c
c   To calculate a confidence interval on a component x(j) of the
c   solution vector, set
c       w = the  j-th  unit vector
c         = all zeroes except for  one  in the  j-th  position.
c
c   For more information on the algorithm, see
c
c      Dianne P. O'Leary and Bert W. Rust, "Confidence intervals for inequality 
c         constrained least squares problems, with applications to ill-posed 
c         problems,"  SIAM J. on Scientific and Statistical Computing, Vol. 7, 
c         April 1986, pp. 473-489.
c
c      Bert W. Rust and Dianne P. O'Leary, "Confidence intervals for discrete 
c         approximations to ill-posed problems," Journal of Computational and
c         Graphical Statistics, Vol. 3, No. 1, pp. 67-96.
c
c
c  Input parameters: 
c  -----------------
c
c     A       real(ldA,n)
c             The constraint matrix.
c             If  A  has rank .lt. n ,  the program may fail to work.
c
c     ldA     integer
c             Leading dimension of  A  in the calling program.
c
c     n       integer
c             Number of variables  x.
c
c     mA      integer
c             Number of constraints (i.e., number of initialized rows of  A).
c
c     WIND    real(ldwind,nw)
c             This array stores the window vectors  w   as columns.
c
c     ldwind  integer
c             Leading dimension of  WIND  in the calling program.
c
c     nw      integer
c             Number of columns in the window matrix  WIND.
c
c     xmin    real(n)
c             Solution to the non-negatively constrained least squares problem.
c             If it is a zero vector, it will be computed herein.
c
c     amusq   real
c             mu**2 = Parameter defining the confidence level.
c
c     b       real(mA)
c             The data vector for the constraints.
c
c     phlo         real(nw)
c                  initial lower bounds for each window
c
c     phup         real(nw)
c                  initial upper bounds for each window.
c
c     epconv       real
c                  convergence criterion: 
c                  the iteration will stop when a value of
c                  the function is within  epconv  of  amusq
c                  ************check this****************** 
c
c     iprint       integer
c                  specifies the amount of output to the default
c                  output file (see ipr).
c                  if  iprint = 0, 1, 2, or 4, and  ipr  is true,
c                  then the amount of output to file 10 will be
c                  one level higher.
c
c                  =0  means error messages only. 
c                  =1  means minimal printing.
c                  =2  means a summary for each window.
c                  =3  means a summary for each iteration
c                      on each window.
c                  =4,5  means copious printing, used for
c                      debugging.
c
c     ipr         logical
c                 = .true.  if output should be sent to the 
c                           default output file (printer or terminal) 
c                           and more complete output to  file 10.
c
c                 = .false. if output is to the default output file
c                           only.
c
c     mou         integer
c                 unit number for line printer output
c
c     elmin       real
c                 L_min = minimum sum of squared residuals for the nonnegatively
c                         constrained least squares problem, 
c                         i.e.,  (b - A xmin)**T * (b - A xmin) .
c                 Computed herein if  xmin = 0-vector
c
c     jmin        integer(n)
c                 basic and nonbasic indices for solution to the nonnegativelu
c                         constrainted least squares problem.
c                 Computed herein if  xmin = 0-vector
c
c     kmin        integer
c                 no. vectors in the basis set for solution to the nonnegatively
c                         constrainted least squares problem.
c                 Computed herein if  xmin = 0-vector
c
c                 The column numbers of the basis corresponding to  xmin  are in
c                 the array   jmin(j), j=1,kmin.  The non-basis column numbers
c                 are in the array  jmin(j), j=kmin+1,n
c
c
c   Output parameters:
c   ------------------ 
c
c     philo        real(nw)
c                  the lower bound for each window.
c                  this array can share storage with  phlo  if
c                  the  phlo  values do not need to be saved.
c
c     phiup        real(nw)
c                  the upper bound for each window.
c                  this array can share storage with  phup  if
c                  the  phup  values do not need to be saved.
c
c     cratlo       real(nw)
c                  a convergence code for each window.
c                  in case of no unusual conditions, these are
c                  convergence ratios:  abs((philo - phihi) / phi ).
c                  (a negative sign means the interval (philo,phihi)
c                  failed to bracket the root  phi.)  this measures...
c
c               otherwise,
c              = 1111111  initial bound  phlo  or  phup  had
c                         L(phi) .le. amusq .
c              = 2222222  iteration stopped because interval became
c                         very small (no error, but probably
c                         expensive).
c              = 3333333  (min  L(phi) ) .ge. amusq.
c                         no confidence interval. 
c              = 6666666  LS  cycled (i.e., computed without terminating
c                         successfully) trying to evaluate  L(phimin),
c                         and no confidence interval was computed.
c              = 7777777  LS  cycled during  BRAKET  iteration,
c                         so the computed confidence interval is
c                         probably not optimal.
c
c     cratup       real(nw)
c                  a convergence code for each window.
c                  see cratlo.
c
c     phiest       real(nw)
c                  the point estimate for each window.
c                  should get  philo.le.phiest.le.phiup
c
c     nbctlo       integer(nw)
c                  the number of basis changes for each window in
c                  the quadratic programming iteration  LS. 
c                  this is a rough measure of the work in computing
c                  each lower bound.
c
c     nbctup       integer(nw)
c                  the number of basis changes for the upper bounds.
c                  see nbctlo.
c
c     ntotbc       integer
c                  sum of entries in  nbctlo  and nbctup.
c
c
c  Scratch storage:
c  ----------------
c
c     work       real(2*(mA+1)*n + 2*max(mA,n) + 4*mA + 5*n + 5)
c
c     iwork      integer(3*n) 
c
c
c  Subroutines used:
c  ----------------- 
c
c     This package:  BRCLS1, BASPR, BRAKET, EVALFN, LS,
c                    VECPR, GSDEL, GSERR, GSINS, GSDERR
c
c     Linpack:       STRSL
c
c     Linpack(BLAS): SCOPY, SDOT, SAXPY, SSCAL
c
c     General:       R1MACH (provides machine constants)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
       logical ipr
c
       dimension  A(ldA,n),   WIND(ldwind,nw),        b(mA),
     1            phlo(nw),   phup(nw),   
     1            philo(nw),  phiup(nw),  cratlo(nw), cratup(nw), 
     1            phiest(nw), nbctlo(nw), nbctup(nw),
     1            work(1),    iwork(1),   jmin(1)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  Divide scratch storage here
c
       iwt   = 1
       ix    = iwt    + n
       ixp   = ix  + n
       iz    = ixp    + n + 1 
       ig    = iz     + n + 1 
       id    = ig     + mA
       iad   = id     + max(mA,n)
       ires  = iad    + max(mA,n)
       ibmod = ires   + mA+1
       irr   = ibmod  + mA+1
       iQ    =  irr   + mA+1
       iR    = iQ     + (mA+1)*n
c
       ijb   = 1
       ijnb  = ijb    + n
c
c   ----------------------------------------------------------------------
      call BRCLS1 ( A,          ldA,         n,           mA,
     1              WIND,       ldwind,      1,           nw,
     2              amusq,      b,           phlo,        phup,
     3              epconv,     iprint,      ipr,         mou,
     4              philo,      phiup,       cratlo,      cratup,
     5              phiest,     nbctlo,      nbctup,      ntotbc,
     6              work(iwt),  work(ix),    xmin,        work(ixp),
     7              work(iz),   work(ig),    work(id),    work(iad),
     8              work(ires), work(ibmod), work(irr),   work(iQ),
     8                                                    work(iR),   
     9              elmin,      iwork(ijb),  iwork(ijnb), jmin,       
     9                                                    kmin  )
c   ----------------------------------------------------------------------
c
c     --------
       return
c     --------
c
       end
c
c
c*******************************************************************************
c
c
      subroutine BRCLS1 ( A,      ldA,    n,      mA,
c
     1                    WIND,   ldwind, iwlo,   iwup,
c
     2                    amusq,  b,      phlo,   phup,
c
     3                    epconv, iprint, ipr,    mou,
c
     4                    philo,  phiup,  cratlo, cratup,
c
     5                    phiest, nbctlo, nbctup, ntotbc,
c
     6                    wt,     x,      xmin,   xp,
c
     7                    z,      g,      d,      ad,
c
     8                    resid,  bmod,   rr,     Q,      R,      
c
     9                    lmin,   jb,     jnb,    jmin,   kmin  )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  See documentation in subroutine BRACLS.
c
c  Scratch storage:
c  ---------------- 
c
c     wt      real(n)
c             window vector.
c
c     x       real(n)
c             solution vector.
c
c     xp      real(n+1)
c             used in  LS .
c
c     z       real(n+1)
c             used in  LS .
c
c     g       real(mA)
c             used in BRAKET, LS.
c
c     d       real(max(mA,n)) 
c             used in  BRAKET.
c
c     ad      real(max(mA,n)) 
c             used in  BRAKET.
c
c     resid   real(mA+1)
c             residual for least squares problem. 
c
c     bmod    real(mA+1)
c             stores a modified  b .
c
c     rr      real(mA+1)
c             working storage used in subroutine GSINS .
c
c     Q       real(mA+1,n)
c             stores orthogonal factor of part of  A.
c
c     R       real(mA+1,n)
c             stores triangular factor of part of  A.
c
c     jb      integer(n)
c             indices of "basic" variables.
c
c     jnb     integer(n)
c             indices of "nonbasic" variables, the zero ones.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      logical  ipr, loflag, upflag, ip
c
      real    lphil0, lphir0, lphmin, lmin, lofphi, lofphl, lofphr
c
      dimension  A(ldA,n), b(ldA), wt(n)
c
      dimension  phlo(iwup), phup(iwup)
c
      dimension  WIND(ldwind,iwup), Q(mA+1,n),
     1           R(mA+1,n),        x(n),       bmod(mA+1), rr(mA+1),
     2           xp(n+1),          z(n+1),     g(n),       resid(ldA)
c
      dimension  jb(n), jnb(n), jmin(n) 
c
      dimension  philo(iwup),  phiup(iwup),  phiest(iwup),
     1           xmin(n),      d(mA),        nbctlo(iwup), ad(mA),
     2           nbctup(iwup), cratlo(iwup), cratup(iwup)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  Initialization of some parameters
c
      if ( ipr ) then
                      iw10 =  1
                      ip10 = -1
              else
                      iw10 = 0
                      ip10 = 0
      end if
c
       ip = ipr
c
       call SCOPY ( max( mA, n ), 0.0, 0, d, 1 ) 
c
      do 40 j=1,n
         jb(j) = 0
         jnb(j) = 0
   40 continue
c
      do 60 j=iwlo,iwup
         philo(j) = 0.0
         phiup(j) = 0.0
         cratlo(j) = 0.0
         cratup(j) = 0.0
         nbctlo(j) = 0
         nbctup(j) = 0
   60 continue
c
      if ( iprint .gt. 0 ) then
             print 80
   80        format ('1')
             if (ipr) write (mou,80) 
      end if
c
c  Set the  LS  convergence constants
c
c        epsmi = smallest interval convergence constant
c
c       epconv = relative precision for converged bounds  ( was previously 
c                obtained from subroutine INPUT  ) 
c
      epsmi = 100.0*R1MACH(1) 
c
      if ( iprint .gt. 0 ) then
             print 94, epsmi
   94    format (///' Smallest  LS  interval convergence const.  =',
     1                1pe11.4/)
             if (ipr) write (mou,94) epsmi
      end if
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  If it hasn't aready been done, then solve for
c
c         xmin = solution of the constrained least squares problem
c
c                         {                        |          }
c            L_min =  min { (b - A x)**T (b - A x) | x .ge. 0 }
c                      x  {                        |          }
c
c
      xmnorm = SNRM2( n, xmin, 1 )
c
c                            -----------
      if ( xmnorm .ne. 0.0 )  go to 100
c                            -----------
c
c  Use subr. LS to solve the non-negatively constrained least squares
c  problem by setting equality constraint to the last row of A and ad 
c  with weight alpha = 1.0 .
c
      k = 0
c
      do j = 1, n
         wt(j) = A(mA,j)
         jnb(j) = j
      end do
c
      alpha = 1.0
c
c   -----------------------------------------------------------
      call LS ( A,      ldA,    mA-1,   n,      b,      b(mA),
c
     1          wt,     alpha,  iprint, ipr,    mou,
c
     2          k,      jb,     jnb,    Q,      mA+1,   
c
     3          R,      mA+1,   xmin,   resid,  rnorm,  g, 
c
     4          nbch,   icycle, bmod,   rr,     xp,     z     )
c   ------------------------------------------------------------
c
      if ( icycle .eq. 1 ) then
                           if ( iprint .gt. 0 ) then
                                                print 101 
                                       if (ipr) write (mou,101)
                           end if
                 k = 0
      end if
c
  101 format (//' Cycling occurred in initial constrained point ',
     1              'estimation problem (calc. of xmin)'/)
c
c
c  Save column numbers of basis for  xmin  in array  jmin(j), j = 1, kmin.  
c            Save non-basis column numbers in array  jmin(j), j = kmin+1, n
c
      kmin = k
c
      call SCOPY ( k, jb, 1, jmin, 1 )
c
      if ( k .lt. n )  call SCOPY ( n-k, jnb, 1, jmin(k+1), 1 )
c
      sumsqr = rnorm*rnorm + resid(1)*resid(1)
c
      lmin = sumsqr 
c
      if ( iprint .gt. 0 ) then
                           print 104,sumsqr
                  if (ipr) write (mou,104) sumsqr
c                -------------------------------------------
                  call VECPR ( xmin, n, 'xmin', iw10, mou )
c                -------------------------------------------
      end if
c
  104 format (///5x,'(',33x,':',12x,')'/5x,'( ::    -1',19x,
     1        ':: 2 :',12x,')'/' min ( :: (S)  *( (K)*(x) - (y)',
     2        ' ) ::   : (x).ge.(0) )  =',1pe16.6/'  x  ( ::',25x,
     3        '::   :',12x,')'/5x,'(',33x,':',12x,')',
     4        ///' The constrained solution vector is' )
c
      if ( amusq .lt. sumsqr ) then
                               print 106, sumsqr
                     if (ipr)  write (mou,106) sumsqr 
      end if
c
  106 format (///' ********************************* Warning ***',
     1           '*************************'///' The problem is ',
     2           'inconsistent -- mu**2 .lt.',1pe15.6,'  !'///
     3           '**************************************************',
     4           '*********************'// )
c
      if (iprint .ge. 4)   then
c            ---------------------------------------------
              call VECPR ( resid, mA, 'rsid', ip10, mou )
c
              call VECPR ( g, n, 'grad', ip10, mou )
c            ---------------------------------------------
      end if
c
      if ( (iprint .ge. 4) .or. (.not. ipr) )
c                    --------------------------------------------------
     1                call BASPR ( jb, k, jnb, n-k, iprint, ipr, mou )
c                    --------------------------------------------------
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
  100 continue
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c         Now do the  BRAKET-LS  quadratic programming procedure
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  Set maximum no. of  LS  iterations to allow while searching for a given root.
c
      maxit = 2*mA
c
      if ( iprint .gt. 0 ) then
            print 80
            if (ipr) write (mou,80)
      end if
c
      ntotbc = 0
c
c
c+++++++++++++++++++++++++ Enter the window loop +++++++++++++++++++++++++++++++
c
c
      do 1000 iw = iwlo, iwup
c
c
c  Generate window no. iw  in array  wt(j), j=1,n.
c
       call SCOPY ( n, WIND(1,iw), 1, wt, 1 ) 
c
      if ( ipr ) print 108, iw
  108 format (///' Window number',i4,'  ----------------------------',
     1           '--------------------------------') 
c
      if ( iprint .gt. 0 ) then
c
            if (ipr) write (mou,108) iw
            if ( .not. ipr ) print 108, iw
c
            if ( iprint .ge. 2 ) call VECPR ( wt, n, 'w**T', iw10, mou )
c
      end if
c
c  Calculate point estimate  phiest = w**T xmin
c
      phiest(iw) = SDOT ( n, wt, 1, xmin, 1 ) 
c
      wmax = 0.0
      do 110 j=1,n
         if ( abs( wt(j) ) .gt. wmax )  wmax = abs( wt(j) )
  110 continue
c
c  Calculate weight alpha to make
c
c          alpha w**T x  = alpha phi 
c
c          an equality in the constrained least squares problem
c
      eta = R1MACH(4)
      gam = 0.0
      do 112 j=1,n
         anorm = abs( wt(j) )
         do 111 i=1,mA
            anorm = anorm + abs( A(i,j) )
  111    continue
         if (anorm.gt.gam) gam = anorm
  112 continue
      alpha = sqrt( sqrt( gam ) / ( 1.0e-4*eta ) )
c
      if ( ipr .and. (iprint .ge. 1) ) write (mou,114) alpha
      if (iprint .ge. 2 ) then
                     print 114, alpha
  114                format (//' The weight of the equality',
     1                         ' constraint in  LS  is ',
     1                         ' alpha =',1pe16.6/) 
      end if
c
c  Initialize  phil0 , L(phil0) , phir0 , L(phir0)
c
      phil0 = phlo(iw)
      phir0 = phup(iw)
c
      nbctlo(iw) = 0
        k = 0
c
c  Call subr. LS  to evaluate L(phil0)
c
c   ------------------------------------------------------------
       call LS ( A,      ldA,    mA,     n,      b,      phil0,
c
     1           wt,     alpha,  iprint, ipr,    mou, 
c
     2           k,      jb,     jnb,    Q,      mA+1,   
c
     3           R,      mA+1,   x,      resid,  rnorm,  g,
c
     4           nbch,   icycle, bmod,   rr,     xp,     z     )
c   -------------------------------------------------------------
c
       lphil0 = rnorm*rnorm
c
      if ( iprint .gt. 0 ) then
             if (iprint.ge.1)  print 130, phil0,lphil0
 130      format (//' Initial lower bound,   phil0  ,  L(phil0)   = ',
     1              1p2e16.6/)
             if(ipr) write(mou,130) phil0, lphil0
      end if
c
      if ( ( icycle .eq. 1 ) .and. ( iprint .gt. 0 ) ) then
             print 133, iw, phil0, lphil0 
  133    format (/' Cycling while evaluating init. L(phil0) for ',
     1            'window ',i3,/' Stopped with  phil0 ,  L(phil0) =',
     2            1p2e15.6/)
             if (ipr) write (mou,133) iw, phil0, lphil0
      end if
c
      nbctlo(iw) = nbctlo(iw) + nbch
c
      loflag = .false.
c
      nbctup(iw) = 0
      k = 0
c
c  Call subroutine  LS  to evaluate  lphir0 = L(phir0)
c
c   -----------------------------------------------------------
      call LS ( A,      ldA,    mA,     n,      b,      phir0,
c
     1          wt,     alpha,  iprint, ipr,    mou,
c
     2          k,      jb,     jnb,    Q,      mA+1,   
c
     3          R,      mA+1,   x,      resid,  rnorm,  g,
c
     4          nbch,   icycle, bmod,   rr,     xp,     z    )
c   ------------------------------------------------------------
c
      lphir0 = rnorm*rnorm
c
      if ( iprint .gt. 0 ) then
             if (iprint.ge.1) print 135, phir0, lphir0 
 135     format (/' Initial upper bound,   phir0  ,  L(phir0)   = ',
     1            1p2e16.6/)
             if (ipr) write (mou,135)  phir0, lphir0
      end if
c
      if ( ( icycle .eq. 1 ) .and. ( iprint .gt. 0 ) ) then
                        print 140, iw, lphir0
  140      format (/' Cycling while evaluating init. L(phir0) for',
     1              ' window ',i3/' Stopped with  L(phir0) =',1pe11.3/)
                        if (ipr) write(mou,140) iw, lphir0
      end if
c
      nbctup(iw) = nbctup(iw) + nbch
c
      upflag = .false.
c
c
c  Check whether  L(phil0) .le. mu**2
c
      if ( lphil0 .le. amusq ) then
c
                philo(iw) = phil0
c
                if (iprint.ge.1) print 150, phil0
 150            format (//' Lower bound set to initial value',1pe18.6/)
                if( ipr .and. ( iprint .gt. 0 ) ) write(mou,150) phil0
c
                cratlo(iw) = 1111111.0 
                loflag = .true.
      end if
c
c  Check whether  L(phir0) .le. mu**2
c
      if ( lphir0 .le. amusq ) then
c
           phiup(iw) = phir0
c
           if (iprint.ge.1) print 160, phir0
 160       format (//' Upper bound set to initial value',1pe18.6/)
           if ( ipr .and. ( iprint .gt. 0 ) ) write(mou,160) phir0
c
           cratup(iw) = 1111111.0 
           upflag = .true.
c
c        Has philo(iw) already been set ? 
c
           if ( loflag ) go to 510
c
           phil = phil0 
           lofphl = lphil0
           phir = phir0 
           lofphr = lphir0
           phi = phir
           lofphi = lofphr
           if ( iprint .ge. 2 ) print 260, k, phi, lofphi
           if ( ipr .and. ( iprint .ge. 1 ) )
     1                            write (mou,260) k,phi,lofphi
c
           go to 300
c                     i.e., go on to compute philo(iw)
      end if
c
c  Program gets to this point only if  L(phir0) .gt. mu**2
c
c  Test whether  philo(iw)  has already been set 
c
      if (loflag) then
                  phil = phil0
                  lofphl = lphil0
                  phir = phir0
                  lofphr = lphir0
                  phi = phir
                  lofphi = lofphr
                  if ( iprint .ge. 2 ) print 260, k, phi, lofphi
                  if ( ipr .and. ( iprint .ge. 1) )
     1                                 write (mou,260) k, phi, lofphi 
                  go to 500
c                           i.e., go on to compute phiup(iw)
      end if
c
c  Program gets here only if
c
c       ( L(phil0) .gt. mu**2 )  .and.  ( L(phir0) .gt. mu**2 )
c
c   In this case it is necessary to estimate  phimin,
c   i.e., to find that value  phi = phimin  such that
c
c                      {                        |                          }
c     L(phimin) = min  { (b - A x)**T (b - A x) | x .ge. 0 ,  w**T x = phi }
c                x,phi {                        |                          }
c
c  -- or, failing that, to find some  phimin = phimid  such that
c
c         (phil0 .lt. phimid .lt. phir0) .and. ( L(phimid) .lt. mu**2 )
c
c  Try to estimate  phimin  by     phimin = wt**T xmin
c
      k = kmin
c
      phimin = phiest(iw)
c
c
      if ( iprint .gt. 0 ) then
             if ( iprint .ge. 1 ) print 170, phimin
 170     format (/' xmin  gives    phimin  =  w**T xmin  =',1pe16.6/)
             if (ipr) write (mou,170) phimin 
      end if
c
c
c  Check whether phimin lies in the initial interval
c
      if ( ( phimin .le. phil0 ) .or. ( phimin .ge. phir0 )
     1                           .or. ( kmin .le. 0 )       )
     1       then
c
             if ( iprint .gt. 0 ) then
c
                   print 180, iw, kmin, phil0, phimin, phir0
  180              format (/' For window no.',i4,3x,'kmin ,philo ,',
     1                      ' phimin , phiro  =',i4,1p3e15.6//
     1                      ' Alternate method would be necessary',
     1                      ' to calculate phimin.'//
     1                      ' The program fails to improve the ',
     1                      ' initial bounds for this window.')
c
                   if (ipr) write (mou,180) iw,kmin,phil0,phimin,phir0 
c
             end if
c
      end if
c
c  Program reaches this point only if  w**T xmin  gave acceptable value for 
c  phimin.
c
c  Set up (m+1) x k matrix  Q  and k x k upper triangle  R  for basis set 
c
c  which gives   phi = phimin ,
c
c  i.e., perform the Gramm-Schmidt orthogonalization
c
c         (  alpha w1**T )
c         (              )                (     )
c         (     A1       )                (  Q  )
c         (              )      G-S       (     )
c         (--------------)  ----------->  (-----)
c         (              )                (     )
c         (      I       )                (  R  )
c         (              )                (     )
c
c     where
c            A1  = matrix of basis columns from kernel matrix  A
c
c             w1 = vector of corresponding elements of window vector
c
c              I = k x k identity matrix
c
      call SCOPY ( k, jmin, 1, jb, 1 )
c
      if ( k .lt. n )  call SCOPY ( n-k, jmin(k+1), 1, jnb, 1 )
c
      do 210 j=1,k
c
         jx = jb(j)
         Q(1,j) = alpha*wt(jx) 
c
         call SCOPY ( mA, A(1,jx), 1, Q(2,j), 1 )
c
c       -------------------------------------------------------------
         call GSINS ( Q, mA+1, mA+1, j-1, Q(1,j), R(1,j), rr, info )
c       -------------------------------------------------------------
c
         if ( info .ne. 0 )  call GSERR ( info, j, j )
c
  210 continue
c
c
c  Initialize   phi = phimin  and evaluate  L(phi)  by projecting (mA+1)-vector
c
c       ( alpha*phi , b**T )**T   onto the column space of matrix  Q  
c
      phi = phimin
c
      resid(1) = alpha*phi
c
      call SCOPY ( mA, b, 1, resid(2), 1 )
c
c    ------------------------------------------------------
      call GSINS ( Q, mA+1, mA+1, k, resid, xp, rr, info )
c    ------------------------------------------------------
c
       if ( info .ne. 0 )  call GSERR ( info, n+1, k+1 )
c
      lofphi = lmin 
      lphmin = lofphi
c
      if ( iprint .ge. 1 ) print 230, lphmin
 230  format (/14x,' L(phimin) =',1pe16.6/)
      if ( ipr .and. ( iprint .ge. 1 ) ) write (mou,230) lphmin
c
c
c  Check whether   L(phimin) .lt. mu**2 
c
c     If not, set  ( philo(iw),phiup(iw) )  to  ( phil0,phir0 ) ,
c
c     and go on to the next window 
c
c
      if ( lphmin .ge. amusq ) then
c
                   philo(iw) = phil0
                   phiup(iw) = phir0
c
                   if ( iprint .gt. 0 ) then
                         if (iprint.ge.1) print 150, phil0
                         if (ipr) write (mou,150) phil0
                         if (iprint.ge.1) print 160, phir0
                         if (ipr) write (mou,160) phir0
                   end if
c
                   cratlo(iw) = 3333333.0
                   cratup(iw) = 3333333.0
c
                   go to 510
c
      end if
c
c
       if ( iprint .ge. 2 ) print 260, k,phi,lofphi
  260 format (1h ///1h ,'initial values   k/phi/L(phi)',i4, 
     1        1p2e16.6/)
      if ( ipr .and. ( iprint .gt. 0  )) write(mou,260) k,phi,lofphi 
c
c  Initialize  phil , phir , L(phil) , and  L(phir)
c
      phil = phil0
      lofphl = lphil0
      phir = phimin 
      lofphr = lphmin
c
      if ( iprint .ge. 2 ) print 280
  280 format (1h /) 
      if ( ipr .and. ( iprint .ge. 1) ) write (mou,280)
c
c
  300 continue
c
c  Call subroutine BRAKET to do the boxing iteration for lower bound  philo(iw)
c
c   ---------------------------------------------------------------
      call BRAKET ( mA,     n,      A,      ldA,    b,      wt,
c
     1              alpha,  k,      jb,     jnb,    Q,      mA+1,
c
     2              R,      mA+1,   x,      resid,  bmod,   rr,
c
     3              xp,     z,      g,      epconv, epsmi,
c
     4              amusq,  maxit,  nbct,   crat,   phir,
c
     5              phil,   phi,    lofphi, lofphr, lofphl, iw,
c
     6              ad,     d,     'lower', iprint, ipr,    mou  )
c   ---------------------------------------------------------------
c
       nbctlo(iw) = nbctlo(iw) + nbct
       cratlo(iw) = crat
c
       philo(iw) = phi
c
         if ( crat .eq. 9999999.0 )  philo(iw) = phil 
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  Now compute the upper bound  phiup(iw)
c
      if ( upflag ) go to 510
c
c
c  Again set up (m+1) x k matrix  Q  and k x k upper triangle  R  for 
c
c     basis set which gives  phi = phimin
c
c    i.e., again perform the Gramm-Schmidt orthogonalization
c
c         ( alpha w1**T )
c         (             )                (     )
c         (     A1      )                (  Q  )
c         (             )      g-s       (     )
c         (-------------)  ----------->  (-----)
c         (             )                (     )
c         (      I      )                (  R  )
c         (             )                (     )
c
c     where
c             A1 = matrix of basis columns from kernel matrix  A
c
c             w1 = vector of corresponding elements of window vector
c
c              i = k x k identity matrix
c
c
      k = kmin
c
      call SCOPY ( k, jmin, 1, jb, 1 )
c
      if ( k .lt. n )  call SCOPY ( n-k, jmin(k+1), 1, jnb, 1 )
c
      do 410 j=1,k
c
         jx = jb(j)
         Q(1,j) = alpha*wt(jx) 
c
         call SCOPY ( mA, A(1,jx), 1, Q(2,j), 1 )
c
c       -------------------------------------------------------------
         call GSINS ( Q, mA+1, mA+1, j-1, Q(1,j), R(1,j), rr, info )
c       -------------------------------------------------------------
c
         if ( info .ne. 0 )  call GSERR ( info, j, j )
c
  410 continue
c
c  Initialize  phi , lofphi , phil , lofphl , and  phir
c
      phi = phimin
      lofphi = lphmin
c
      if ( iprint .ge. 2 ) print 260, k,phi,lofphi
      if ( ( iprint .ge. 1 ) .and. ipr ) write(mou,260) k,phi,lofphi 
c
c
      phil = phimin 
      lofphl = lphmin
      phir = phir0
      lofphr = lphir0
c
      if ( iprint .ge. 2 ) print 280
      if ( ipr .and. ( iprint .ge. 1 ) ) write (mou,280)
c
  500 continue
c
c  Call subroutine BRAKET to do boxing iteration for upper bound phiup(iw)
c
c   --------------------------------------------------------------
      call BRAKET ( mA,     n,      A,      ldA,    b,      wt,
c
     1              alpha,  k,      jb,     jnb,    Q,      mA+1,
c
     2              R,      mA+1,   x,      resid,  bmod,   rr,
c
     3              xp,     z,      g,      epconv, epsmi,
c
     4              amusq,  maxit,  nbct,   crat,   phil,
c
     5              phir,   phi,    lofphi, lofphl, lofphr, iw,
c
     6              ad,     d,     'upper', iprint, ipr,    mou  )
c   --------------------------------------------------------------
c
       nbctup(iw) = nbctup(iw) + nbct
       cratup(iw) = crat
c
       phiup(iw) = phi
c
         if ( crat .eq. 9999999.0 )  phiup(iw) = phir 
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
  510 continue
c
      ntotbc = ntotbc + nbctlo(iw) + nbctup(iw)
c
c
 1000 continue
c
c+++++++++++++++++ End of window loop ++++++++++++++++++++++++++++++++++++++++++
c
c     --------
       return
c     --------
c
       end
c
c
c*******************************************************************************
c
c
      subroutine BRAKET ( mA,     n,      A,      ldA,    b,      wt,
c
     1                    alpha,  k,      jb,     jnb,    Q,      ldQ,          
c
     2                    R,      ldR,    x,      resid,  bmod,   rr,
c
     3                    xp,     z,      g,      epconv, epsmi,        
c
     4                    amusq,  maxit,  nbct,   crat,   phineg,
c
     5                    phipos, phi,    lofphi, lofphn, lofphp, iw,
c
     6                    ad,     d,      lowup,  iprint, ipr,    mou  )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  A subroutine to solve the piecewise quadratic equation
c
c               L(phi)  -  amusq  =  0.0
c
c        where
c
c               L(phi) = piecewise quadratic function defined
c                        for any value of phi by the subroutine LS
c
c        and
c                amusq = constant picked to give desired confidence level
c
c  Subroutine LS computes
c
c                    {                        |                         }
c     L(phi)  =  min { (b - A x)**T (b - A x) | w**T x = phi , x .ge. 0 }
c                 x  {                        |                         }
c
c  for any given phi.   This subroutine uses quadratic extrapolation
c  coupled with bisection to solve 
c
c                     L(phi)  =  amusq   
c
c  for a root phi contained in the initial interval  (phineg,phipos).
c
c  On Input
c  --------
c
c      phineg = a  bound for phi
c
c      lofphn = L(phineg) < amusq
c
c      phipos = a  bound for phi
c
c      lofphp = L(phipos) > amusq
c
c         phi = initial estimate of desired root
c
c      lofphi = L(phi)
c
c  on output,
c
c        phi = converged value of desired root
c
c     lofphi = L(phi)
c
c       nbct = total number of basis changes used to get phi
c
c       crat = convergence ratio, or a code for special
c                      circumstances encountered in seeking phi
c
c
      external  LS
c
      logical  ipr
c
      character*5   lowup
      character*15  bisect, corrdn, corrup, chcode, reset, cycl
c
      real  lofphi, lofphn, lofphp,
     1      lpmimo, lpmimu
c
      dimension  A(ldA,n), Q(ldQ,n), R(ldR,n)
c
      dimension  b(ldA),   wt(n),  x(n), resid(ldA), bmod(ldA+1),
     1           xp(n+1),  z(n+1), g(n), ad(mA),     d(mA), 
     2           rr(ldA+1)
c
      integer  jb(n), jnb(n)
c
      data bisect /'bisection      '/
      data corrup /'correction up  '/
      data corrdn /'correction down'/
      data reset  /'reset to delmin'/
      data cycl   /'cycling  change'/
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      epsl = 10.0*R1MACH(4)
      phiold = phi
c
      nbct = 0
      delphi = 0.0
      lpmimu = lofphi - amusq 
c
c++++++++++++  Enter iteration for solving   L(phi) = amusq  +++++++++++++++++++
c
      do 400 itc=1,maxit
c
c
      phiold = phi
      dphold = delphi
c
c  On first pass through loop, or after cycling  (i.e. any time k=0),
c  get new  phi  by bisection.
c
      if ( (itc .eq. 1 ) .or. ( k .eq. 0 ) ) go to 263
c
c  Calculate candidate correction  delphi  by quadratic extrapolation.
c
c  First solve triangular system
c
c                              ( alpha )
c                   R d = Q**T (-------)
c                              (   0   )
c
c  to get the n-vector d
c
      call SCOPY ( k, Q(1,1), ldQ, d, 1 )
c
      call SSCAL ( k, alpha, d, 1 )
c
c    ---------------------------------------
      call STRSL ( R, ldR, k, d, 01, info )
c    ---------------------------------------
c
      if ( ( info .ne. 0 ) .and. ( iprint .gt. 0 ) ) then 
                 print 235,info
  235            format (1h ,'STRSL returned  info = ',i3/) 
                 if (ipr) write (mou,235) info
      end if
c
c  Calculate   daad = d**T A**T A d
c
      do 250 i=1,mA 
c
         ad(i) = 0.0
c
  250 continue
c
      do 260 j=1,k
c
         ja = jb(j)
c  
         call SAXPY ( mA, d(j), A(1,ja), 1, ad, 1 )
c
  260 continue
c
      daad = SDOT ( mA, ad, 1, ad, 1 )
c
c  If  daad = 0  get new phi by bisection
c
      if ( daad .eq. 0.0 ) go to 263
c
c  Calculate   dar0 = d**T A**T b - A x
c
      dar0 = SDOT ( mA, ad, 1, resid(2), 1 )
c
c  Solve the quadratic equation
c
c           daad*delphi**2 - 2*dar0*delphi + lofphi - amusq  =  0
c
c  for delphi
c
      discr = dar0**2 - daad*(lofphi - amusq)
c
      if ( discr .le. 0.0 ) then
             delphi = dar0/daad
      else
             discr = sqrt(discr)
             if ( dar0 .lt. 0.0 ) delph1 = ( dar0 - discr ) / daad
             if ( dar0 .ge. 0.0 ) delph1 = ( dar0 + discr ) / daad
c
             delph2 = ( lofphi - amusq) / ( delph1*daad)
c
             if ( phipos .lt. phineg )  then
                               delphi = amin1(delph1,delph2)
             else
                               delphi = amax1(delph1,delph2)
             end if
      end if
c
c  End quadratic extrapolation computation
c
c  Check whether new  delphi  is very similar to previous delphi. 
c  If so iteration is not converging quadratically, so get the new  phi  
c  by bisection.
c
      if ( delphi .ne. 0.0 ) then 
             if ( abs( 1.0 - abs( dphold/delphi ) ) .lt. 0.01 ) 
     1                                                      go to 263 
             if ( delphi .lt. 0.0 ) chcode = corrdn
             if ( delphi .gt. 0.0 ) chcode = corrup
       end if
c
c  Check whether new  delphi  larger than minimum acceptable value.  
c  If not reset  delphi  to that value (delmin).
c
      delmin = epconv*abs(phi)*0.9
c
      if ( abs( delphi ) .lt. delmin ) then
                    delphi = sign(1.0,delphi)*delmin
                    chcode = reset
      end if
c
c  Get new  phi  by correcting the previous phi
c
      phi = phi + delphi
c
c  If new phi lies in interval  (phineg,phipos)  use it.
c  Otherwise get new  phi  by bisection
c
      if ( ( (phi .ge. phineg) .and. (phi .le. phipos) ) .or.
     1     ( (phi .le. phineg) .and. (phi .ge. phipos) ) ) go to 267
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
 263  continue
c
c  Calculate new  phi  by bisection
c
      chcode = bisect
      phi = (phineg + phipos)/2.0
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
 267  continue
c
c  Print summary for new  phi , delphi , etc.
c
       if ( iprint .ge. 3 )
     1      print 268, iw,itc,phineg,phipos,phiold,delphi,phi,chcode
  268 format (' wind/it',2i4,3x,'phineg/phipos',1p2e15.6/
     1        11x,'phiold/delphi/phi',1p3e15.6/56x,a15)
       if ( ( iprint .ge. 2 ) .and. ipr )
     1              write(mou,268) iw,itc,phineg,phipos,phiold,delphi, 
     1                                 phi,chcode 
c
c  Call subr.  EVALFN  to evaluate  L(phi)  for new  phi
c
c  EVALFN calls  LS  to compute the function value at  phi.
c
c  If  LS fails,  other points in the interval (phipos,phineg)
c  are chosen as alternatives.
c
       phitry = phi
c
c   ---------------------------------------------------------------
       call EVALFN ( phipos, phineg, lofphp, lofphn, amusq, 
c
     1               ierror, mA,     n,      A,      ldA,    b,
c
     2               phi,    wt,     alpha,  k,      jb,     jnb,
c
     3               Q,      ldQ,    R,      ldR,    x,
c
     4               resid,  nbch,   bmod,   rr,     xp,     z,
c
     5               g,      rnorm,  infor,  iprint, ipr,    mou  )
c   ---------------------------------------------------------------
c
       lofphi = rnorm*rnorm
       lpmimo = lpmimu
       lpmimu = lofphi - amusq
       nbct = nbct + nbch
c
       if ( ierror .eq. 1 )  then 
               if ( iprint .gt. 0 ) then
                      print 270, phitry,iw,lowup,phi,lofphi
  270                 format (' Cycling at phi = ',e13.6,
     1                       ' while searching for window no.',i3, 
     1                 1x,a5,' bound'/' Changed phi  ,  L(phi)  = ',
     1                       1p2e13.6/)
                       if (ipr) 
     1                     write (mou,270) phitry,iw,lowup,phi,lofphi
                end if
                chcode = cycl 
       end if
c
       if ( ierror .eq. 2 )  then 
              if ( iprint .gt. 0 ) then 
                      print 280, iw,lowup,phi,lofphi
  280                format ( 
     1               ' Cycling while searching for window no. ',
     1                 i3, 1x,a5,' bound'/' phi  ,  L(phi)  =',
     2                 1p2e13.6/)
                      if (ipr)
     1                      write (mou,280) iw,lowup,phi,lofphi
                      print 282, phipos,lofphp, phineg,lofphn
  282                 format (' BRAKET terminated with phipos,',
     1                     ' L(phipos), phineg, L(phineg) ',
     1                     1p4e13.6/)
                      if (ipr) write (mou,282)
     1                        phipos,lofphp, phineg,lofphn
               end if
c
               k = 0
c
               phi = phipos
               lofphi = lofphp
               crat = 7777777.0
c
c             --------
               return
c             --------
      end if
c
c  Print summary for new  L(phi) 
c
      if ( iprint .ge. 3 ) print 285, nbch,k,phi,lofphi
  285 format (11x,'nbch/k',2i4,8x,'phi/L(phi)',
     1        1p2e15.6/)
      if ( ( iprint .ge. 2 ) .and. ipr ) 
     1                             write (mou,285) nbch,k,phi,lofphi
c
c
c  If new  delphi  was gotten by resetting to  delmin,
c  and if new L(phi) on same side of  amusq  as old  L(phi),
c  then go back and recalculate new  phi  by bisection
c
      if ( ( sign( 1.0, lpmimo ) .eq. sign( 1.0, lpmimu ) ) 
     1     .and.                   ( chcode .eq. reset ) ) go to 263
c
c
c  The new  phi  is acceptable.  
c  Replace either  phineg  or  phipos with the new  phi  so that the root 
c  is still bracketed.
c
      if  ( lpmimu .ge. 0.0 ) then
                                         phipos = phi
                                         lofphp = lofphi
                                    else
c
                                         phineg = phi
                                         lofphn = lofphi
      end if
c
c+++++++++++++++++++++++ Now check for convergence +++++++++++++++++++++++++++++
c
      phint = abs(phipos - phineg)
c
c  +++++++++++++++++++
c  Convergence test 1:  Absolute length of interval small
c  +++++++++++++++++++
c
      if ( phint .le. epsmi ) then
c
            if ( iprint .gt. 0 ) then
                      if( iprint .ge. 1 ) print 287
  287           format (/' Iteration stopped because interval width',
     1                  ' is very small'/) 
                      if (ipr) write (mou,287)
                      if ( iprint .ge. 2 ) print 330, lowup,
     1                                                phi,lofphi,nbct
                      if (ipr) 
     1                write (mou,330) lowup,phi,lofphi,nbct 
             end if
c
             crat = 2222222.0
c
             go to 310
c
      end if
c
c  +++++++++++++++++++
c  Convergence test 2:  Relative length of interval small
c  +++++++++++++++++++
c
      phites = amax1 ( abs(phi) , abs(phineg) , abs(phipos) )
c
      if ( abs( phint / phites ) .le. epconv ) then
c
              if ( iprint .gt. 0 ) then 
                     if ( iprint .ge. 1 ) print 289, phineg, phipos
  289           format (/1h ,'phi neg approx = phi pos',
     1                         4x,'phineg/phipos',1p2e15.6/)
                     if (ipr)
     1                   write (mou,289) phineg,phipos
              end if
c
              go to 310
c
      end if
c
c  +++++++++++++++++++
c  Convergence test 3:  Function value close to zero
c  +++++++++++++++++++
c
      if ( abs( (lofphi - amusq) / amusq ) .gt. epsl ) go to 400
c
      if ( iprint .gt. 0 ) then
            if ( iprint .ge. 1 ) print 300
  300       format (/1h ,'L(phi) approx = mu**2'/)
            if (ipr) write (mou,300)
      end if
c
  310 if ( phi .ne. 0.0 ) then
             crat = abs( (phipos - phineg) / phi )
             crat = - sign( crat, (lofphp - amusq)*(lofphn - amusq) )
      else
             crat = R1MACH(2) 
      end if
c
      if ( iprint .gt. 0 ) then
             if ( iprint .ge. 1 )
     1                    print 330, lowup, phi, lofphi, nbct
  330     format (//' ',a5,' Bound root at   phi , L(phi)  = ',
     1            1p2e16.6/' Total no. basis changes =',i6//)
             if (ipr) write(mou,330) lowup, phi, lofphi, nbct
      end if
c
      if (phi.ne.0.0) then
c
             crat = abs( (phipos - phineg) / phi )
             crat = - sign ( crat, (lofphp - amusq)*(lofphn - amusq) ) 
      else
             crat = R1MACH(2) 
c
      end if
c
c    --------
      return
c    --------
c
  400 continue
c
c+++++++++++++++ End of iteration solving  L(phi) = amusq  +++++++++++++++++++++
c
      if ( iprint .gt. 0 ) then
             print 450, maxit,nbct
  450        format (///' No root found in',i5,'  iterations'/,
     1                  ' Total no. basis changes =',i6//)
             if (ipr) write (mou,450) maxit,nbct
      end if
c
      crat = 9999999.0
c
      phi = phipos
      lofphi = lofphp
c
c    --------
      return
c    --------
c
      end 
c
c
c*******************************************************************************
c
c
       subroutine EVALFN ( phipos, phineg, lofphp, lofphn, amusq, 
c
     1                     ierror, mA,     n,      A,      ldA,   b, 
c
     2                     phi,    wt,     alpha,  k,      jb,    jnb,
c
     3                     Q,      ldQ,    R,      ldR,    x, 
c
     4                     resid,  ntbch,  bmod,   rr,     xp,    z, 
c
     5                     g,      rnorm,  infor,  iprint, ipr,   mou )       
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
       real  phipos, phineg, lofphp, lofphn, amusq, A(ldA,n),
c
     1       b(ldA), phi, wt(n), alpha, Q(ldQ,n), R(ldR,n), rr(ldA+1),
c
     1       x(n), resid(ldA), bmod(ldA+1), xp(n), z(n), g(n), rnorm
c
       integer  ierror, n, ldA, k, jb(n), jnb(n), ldQ, ldR, ntbch,
c
     1          infor, iprint 
c
       logical  ipr
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c   This subroutine uses LS to evaluate  L(phi).
c
c       If  LS  cycles with a value greater than  amusq, then
c       the subroutine tries  ktry  alternate points in the 
c       interval  (phizro,phipos) in order to get a successful
c       value from  LS  or a cycled value less than  amusq. 
c
c       On return
c       ---------
c                 ierror = 0  means  LS was successful at  phi
c                        = 1  means  LS was successful at an alternate
c                                       point
c                        = 2  means  LS failed at all trial points
c
c                 phi is the point at which the function was evaluated,
c                         or  phizro  if  LS  failed at all trial points.
c
c       phizro  is the zero of the linear interpolant between
c          (phipos,lofphp-amusq)  and (phineg,lofphn-amusq) .
c
c       Since the function is convex,  it is guaranteed to be
c       less than  amusq on the interval (phizro,phineg), and thus
c       the desired root is bracketed by (phipos,phizro) .
c
c       Note:  this routine can be used to search for a noncycling point
c       in any desired interval (a,b)  by setting 
c            phipos = a         lofphp = 1          amusq = 0
c            phineg = 2b-a      lofphn = -1 .
c       This gives phizro = b, and the search will start at  b .
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
       ktry = 25
c
       ntbch = 0
       ierror = 0
c
       phizro = ( (lofphp - amusq)*phineg - (lofphn - amusq)*phipos )
     1                  / (lofphp - lofphn)
c
       do 200 kount = 1, ktry 
c
c   ---------------------------------------------------------------
          call LS ( A,      ldA,    mA,     n,      b,      phi,
c
     1              wt,     alpha,  iprint, ipr,    mou,
c
     2              k,      jb,     jnb,    Q,      ldQ,    
c
     3              R,      ldR,    x,      resid,  rnorm,  g,
c
     4              nbch,   infor,  bmod,   rr,     xp,     z    )
c   ----------------------------------------------------------------
c
          ntbch = ntbch + nbch
c
c                             --------
          if ( infor .ne. 1 )  return
c                             --------
c
          ierror = 1
          phi = phizro + (kount)*(phipos - phizro) / (ktry+1)
          k = 0
c
  200  continue
c
       ierror = 2
c
c     --------
       return
c     --------
c
       end
