c
c**********************************************************************
c
      subroutine ferdo ( A, b, s, w, HT, q, ut, atemp, istar, plab, mr,
c
     1                   mc, nr, nc, nw, nbits, tau, kr, mu, plo, pup,
c
     2                   nscr, pest )
c
c----------------------------------------------------------------------
c
c  computational guts of ferdo - calcs plo-pup solution
c
c
c---------------------------------------------------------------------
c
      real  mu, kappa
c
      dimension  A(mr,mc), b(mr), s(mr), w(mc), HT(mr,mc), q(mc),
c
     1           ut(mr), atemp(mc), plo(kr), pup(kr), pest(kr)
c
      dimension  istar(mc), plab(kr)
c
c----------------------------------------------------------------------
c
c**  Set  kappa = 1.0  for  alpha = 0.6826
c**                    i.e., for 68.26% 1-at-a-time confidencde intervals 
      kappa = 1.0
c
c** Calculate q vector
c
      do 21 j=1,nc
         istar(j) = 0
         q(j) = r1mach( 2 )
         do 20 i=1,nr
            if ( A(i,j) ) 10, 20, 10
   10       qmin = ( amax1( b(i), 0.0 ) + mu*s(i) ) / A(i,j)
            if ( q(j) - qmin ) 20, 20, 15
   15       q(j) = qmin
            istar (j)= i
   20    continue
   21 continue
c
c** Calculate  HT     HT = S**(-1) A Q
c
      enorm = 0.0
      do 25 i=1,nr
         do 25 j=1,nc
            HT(i,j) = ( q(j)*A(i,j) ) / s(i)
            enorm = enorm + HT(i,j)**2
   25 continue
c
c**  GINV replaces HT by (HT) * ( (HT)**T*(HT) + (tau**2/nc) I )**-1
c
c    -------------------------------------------------------------
      call GINV ( HT, mr, nr, nc, tau/sqrt(float(nc)), A, atemp )
c    -------------------------------------------------------------
c
      read (nscr) A 
      rewind nscr
c
c** Make HT the transpose of solution matrix
c**    Un-normalize, i.e., replace  HT  by  S**-1 HT Q
c
      do 50 i=1,nr
         do 50 j=1,nc
            HT(i,j) = ( q(j)*HT(i,j) ) / s(i)
   50 continue
c
c** Now have  (HT) = transpose of
c                           ( A**T S**-2 A + (tau**2/nc) Q**-2 )**-1 A**T S**-2
      write (6,907) 
      runsum=0.0
c
c** Window loop begins here 
c
      do 500 k=1,nw 
c
         read(nscr) w
c
         do 60 i=1,nr
         ut(i) = 0.0
   60    continue
c
c**      u**T = w**T H
c
         do 100 j=1,nc 
            if( w(j) ) 70, 100, 70
   70       do 80 i=1,nr
               ut(i) = ut(i) + w(j)*HT(i,j)
   80       continue
  100    continue
c
c**      Compute  u**T b , the point estimate,
c            and  u**T S**2 u , the variance of the estimate
c
         ussum = 0.0
           utb = 0.0
         do 150 i=1,nr 
              utb = utb + ut(i)*b(i)
            ussum = ussum + ( ut(i)*s(i) )**2
  150    continue
c
c**      Save the point estimate in array pest
c
         pest(k) = utb 
c
c**      Compute the bias term
c
         wuaq = 0.0
         do 200 j=1,nc
c
            sum = 0.0
            do 175 i=1,nr 
               sum = sum + ut(i)*A(i,j)
  175       continue
c
            wuaq = wuaq + q(j)*abs( w(j) - sum )
c                                                *** wuaq estimates bias error
  200    continue
c
c**      Compute epslon = bias + kappa*(standard deviation) 
c
         epslon = kappa*sqrt( ussum ) + wuaq
c
c**      Compute the confidence bounds
c
         phiup = utb + epslon
         philo = utb - epslon
c
         err1 = sqrt( ussum )
         err2 = wuaq
         sp = 100. * epslon / utb
c
c**      The following section computes the running integral of the spectrum
c**          from 0 to E and approximate error ( sqrt( 2.* pi ) = 2.506....)
c
         if ( k .gt. 5 )  runsum = runsum + utb*( plab(k) - plab(k-1) )
         write (6,904) plab(k),k,philo,phiup,utb,sp,runsum,err1,err2
c
         pup(k) = phiup
         plo(k) = philo
c
  500 continue
c
c** Window loop ends here
c
c    --------
      return
c    --------
c
  904 format(1x,1pe13.3,i5,3e13.3,0pf8.1,1pe18.8,2e13.3)
  907 format (1h0,5x,6hEnergy,14x,3hplo,10x,3hpup,9x,4hpave,6x,5hsp((),
     1 8x,6hrunsum,10x,4herr1,9x,4herr2/)
c
c
      end 
