c
c*******************************************************************************
c
      subroutine QCALC ( A,      mr,     nr,     nc,     b, 
c
     1                   p,      q,      WIND,           mc,
c
     2                   nwtot,  phlo,   phup,   ipr,
c
     3                   iprint, mou  )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  Subroutine to compute initial bounds  phlo(k)  and  phup(k)  for the averages
c
c                 phi_k = w_k**T x  ,  k = 1, 2, ... , nwtot,
c
c       where  w_k  is the kth column of the nc x nwtot matrix  WIND , 
c 
c                   and  x  is assumed to be non-negative,
c
c  If  A , b , and  x  are all non-negative, It also computes a vector  q  of
c
c                      upper bounds for the elements of  x
c 
c-------------------------------------------------------------------------------
c
      external  R1MACH, VECPR
c
      real  mu
c
      logical  ipr
c
      dimension  A(mr,nc),    WIND(mc,nwtot), b(nr), p(nc), q(nc), 
     1           phlo(nwtot), phup(nwtot)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      if (ipr) then 
                    iw10 = 1
               else 
                    iw10 = 0
      end if
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  Check whether q has already been set.  If so, go to 100 to calc phlo & phup.
c                                         If not, compute q below assuming that
c                                         A , x , b , and  WIND non-negative.
c
      do 310 j=1,nc
c                                    -----------
         if ( q(j) .lt. r1mach( 2 ) ) go to 100
c                                    -----------
  310 continue
c
c  Check whether  A  and  b  are both non-negative.
c  If so  compute  q ,  phlo ,  and  phup , assuming  x  is nonnegative  
c  If not, do not change input  q , but use it to compute  phlo  and  phup
c
      do 5 i=1,nr
c                             -----------
         if ( b(i) .lt. 0.0 )  go to 100
c                             -----------
         do 5 j=1,nc
c                                    -----------
              if ( A(i,j) .lt. 0.0 )  go to 100
c                                    -----------
    5 continue
c

c
c  Set  mu  for 99.95% (~100%) simultaneous confidence in q-bounds
c  (This is not the value of mu to be used in calculating the final bounds). 
c
      twonm1 = 2 * max0( nr, nc ) - 1
c
      mu = ( 1.0 / sqrt( 2.0 ) ) * ( 3.485 + sqrt( twonm1 ) )
c
      print 14, mu
      if (ipr) write(mou,14) mu
   14 format('1',' ***   For q-box bounds in subr. QCALC, mu  =', 
     1            1pe11.3,'   ***'/)
c
c  Calculate q-vector (i.e. diagonal elements of diagonal matrix Q)
c
      do 21 j=1,nc
c
         q(j) = R1MACH( 2 )
c
         do 20 i=1,nr
c
            if ( A(i,j) .gt. 0.0 ) then
c
                   qmin = ( amax1( b(i), 0.0 ) + mu ) / A(i,j)
c
                   if ( qmin .lt. q(j) )  q(j) = qmin
c
            end if
c
   20    continue
c
   21 continue
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c---------------
 100  continue
c---------------
c
c                     Compute  phlo  and  phup
c
      if ( iprint .gt. 0 ) then
             if ( .not. ipr ) print *, '1'
             if ( ipr )   write(mou,*) '1'
             print 22
             if ( ipr ) write (mou,22)
             call VECPR ( q, nc, 'q-vector', iw10, mou )
      end if
   22 format ( '0'/' q-vector and initial  phlo , phup  bounds' )
c
      do 200 iw=1,nwtot
c
         phlo(iw) = 0.0
         phup(iw) = 0.0
c
         do 150 j=1,nc
c 
            if ( WIND(j,iw) .ge. 0.0 ) then
c
                          phlo(iw) = phlo(iw) + WIND(j,iw)*p(j)
                          phup(iw) = phup(iw) + WIND(j,iw)*q(j)
                   else
                          phlo(iw) = phlo(iw) + WIND(j,iw)*q(j)
                          phup(iw) = phup(iw) + WIND(j,iw)*p(j)
            end if
c
  150    continue
c
  200 continue
c
      if ( iprint .gt. 0 ) then
c
               call VECPR ( phlo, nwtot, 'init. phi_lo', iw10, mou )
               call VECPR ( phup, nwtot, 'init. phi_up', iw10, mou )
      end if
c
c    --------
      return
c    --------
c
      end 
