c
c
c*********************************************************************
c
c  subroutine  f e r d  for single dimensional confidence intervals
c
c
c-----       subroutine  f e r d  with call loop
c
c
      subroutine ferd ( alo, aup,  b,     s,     w,    ht,    q,
c
     1                  ut,  utmp, atemp, istar, plab, blab,  elab,
c
     2                  mr,  mc,   nr,    nc,    nw,   nopt5, nbits,
c
     3                  taw, kr,   mu,    plo,   pup,  nscr,  nrcw,
c
     4                  id,  pest  )
c
c-----------------------------------------------------------------------
c
      real  mu
c
      character*2 isym(3)
c
c
      dimension  alo(mr,nc),aup(mr,nc),b(nr),s(nr),w(nc),ht(mr,nc),
c
     1         q(nc),ut(nr),utmp(nrcw,10),atemp(nc),istar(nc),plab(nw),
c
     2         blab(nr),elab(nc),plo(kr),pup(kr),id(20), pest(kr)
c
      data isym/ 'l ', 'u ', '= ' /
c
c----------------------------------------------------------------------
c
c  Calculate q vector
c
      do 20 j=1,nc
         istar(j) = 0
c                     ** Addition to allow for zero columns in matrix  5/7/80
         q(j) = r1mach( 2 )
         do 20 i=1,nr
            b(i) = amax1( b(i), 0.0 )
            if ( alo(i,j) ) 10, 20, 10
   10       qmin = ( b(i) + mu*s(i) ) / alo(i,j)
            if ( q(j) - qmin ) 20, 20, 15
   15       q(j) = qmin
            istar(j) = i
   20 continue
c
      do 22 j=1,nc
         if ( istar(j) .eq. 0.0 ) q(j) = 0.0
   22 continue
c              ** Addition for zero columns in matrix   5/7/80
c
c
      do 30 i=1,nr
c
c     calculate effective matrix errors (using ut vector)
c
      ut(i)=0.0
      do 30 j=1,nc
      ut(i)=ut(i)+(aup(i,j)-alo(i,j))*q(j)
   30 continue
c
c     normalize matrix (using average matrix q, s, and ut from above) 
c
      enorm=0.0
      do 40 i=1,nr
      do 40 j=1,nc
      ht(i,j)=(q(j)*((alo(i,j)+aup(i,j))/2.0))/(s(i)+ut(i)) 
      enorm=enorm+ht(i,j)*ht(i,j)
   40 continue
      enorm=sqrt(enorm)
      conind=enorm*10.0*0.5**nbits
c
c     increase taw for large condition number
c
      taw=amax1(taw,conind)
c
      write (6,905) conind
      print *, ' '
      print *, ' taw = ', taw
      print *, ' '
c
c  Calculate transformation matrix
c    ---------------------------------------------------------------
      call ginv ( ht, mr, nr, nc, taw/sqrt(float(nc)), alo, atemp )
c    ---------------------------------------------------------------
c
      read(nscr) alo
      rewind nscr
c
c  Unnormalize transformed matrix
c
      do 50 i=1,nr
         do 50 j=1,nc
            ht(i,j) = ( q(j)*ht(i,j) ) / ( s(i) + ut(i) )
   50 continue
c
c  Compute the estimated solution vector x (in array w)
c
      do 55 j=1,nc
         w(j) = 0.0
         do 55 i=1,nr
            w(j) = w(j) + HT(i,j)*b(i)
   55 continue
      call VECPR( w, nc, '   x', 1 )
c
c++++  Begin window loop  ++++
c
      iter=19
c
      do 200 k=1,nw 
c
      read(nscr) w
c
      iter=iter+1
c
c     calculate ut=wth
c
      do 60 i=1,nr
      ut(i)=0.0
   60 continue
      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 point estimate for current window
c
      pest(k) = 0.0 
      do 110 i=1,nr 
         pest(k) = pest(k) + ut(i)*b(i) 
  110 continue
c
c---- lower estimate for window
c
      call loop(w,b,s,q,alo,aup,mr,nr,nc,philo,err,istar,ut,utmp(1,1),
     1utmp(1,2),utmp(1,3),que)
      utmp(k,6)=err 
      utmp(k,7)=philo
      plo(k) = philo
      qlo=que
c
c     upper estimate for window
c
  125 do 150 j=1,nc 
      w(j)=-w(j)
  150 continue
      do 160 i=1,nr 
      ut(i)=-ut(i)
  160 continue
      call loop(w,b,s,q,alo,aup,mr,nr,nc,phiup,err,istar,ut,utmp(1,1),
     1utmp(1,2),utmp(1,4),que)
      phiup=-phiup
      utmp(k,6)=utmp(k,6)+err 
      utmp(k,8)=phiup
      pup(k) = phiup
      qup=que
      utmp(k,9)=(utmp(k,6)/(utmp(k,8)-utmp(k,7)))*100.0
      utmp(k,10)=100.0*(qup+qlo)/(phiup-philo)
c
c     optional plot of upper and lower synthetic windows
c
      if(nopt5-3)200,170,190
  170 if(iter-20)200,180,200
  180 iter=0
  190 do 195 j=1,nc 
      utmp(j,5) = -w(j)
      utmp(j,4)=-utmp(j,4)
  195 continue
c
      write(6,907) k,plab(k)
      call pplot3(utmp(1,3),nrcw,1,nc,0.004,3,4hener,elab,isym,6)
c
  200 continue
c
      if(nopt5-1)700,700,220
c
c     start b adjusted loop
c
  220 write (6,908) id
      write (6,902) 
      factor=0.0
c                   t
c calculate aave(.j) h
c
      do 400 l=1,nr 
      do 250 i=1,nr 
      sum=0.0
      do 225 j=1,nc 
      sum=sum+((alo(l,j)+aup(l,j))/2.0)*ht(i,j)
  225 continue
      ut(i)=sum
  250 continue
c
c     lower estimate for b-adj
c
      do 260 j=1,nc 
      atemp(j)=alo(l,j)
  260 continue
      call loop(atemp,b,s,q,alo,aup,mr,nr,nc,philo,err,istar,ut,
     1utmp(1,1),utmp(1,2),utmp(1,3),que)
      utmp(l,4)=philo
c
c     upper estimate for b-adj
c
      do 270 i=1,nr 
      ut(i)=-ut(i)
  270 continue
      do 280 j=1,nc 
      atemp(j)=-aup(l,j)
  280 continue
      call loop(atemp,b,s,q,alo,aup,mr,nr,nc,phiup,err,istar,ut,
     1utmp(1,1),utmp(1,2),utmp(1,3),que)
      phiup=-phiup
      utmp(l,5)=phiup
      badj=0.5*(philo+phiup)
      delb=b(l)-badj
      delbsl=delb/s(l)
      pctdev=delbsl*100.0
      factor=factor+delbsl**2 
      write (6,903) l,blab(l),b(l),s(l),philo,phiup,delb,pctdev
  400 continue
      factor=sqrt(factor/float(nr))
c
      do 450 l=1,nr 
      utmp(l,3) = utmp(l,4)
      utmp(l,4) = utmp(l,5)
  450 utmp(l,5) = b(l)
      write(6,906)
      call pplot3(utmp(1,3),nrcw,1,nr,0.004,3,4hp.h.,blab,isym,6)
      write (6,904) factor
c
  700 if(nopt5.gt.1)  write(6,908) id
      write (6,900) 
      do 750 i=1,nw 
      pavg = ( utmp(i,7) + utmp(i,8) ) / 2.0
      sp = 100. * (utmp(i,8)-pavg) / pavg
  750 write(6,901) i,plab(i),plo(i),pup(i),pavg,sp,utmp(i,6),utmp(i,9),
     1 utmp(i,10)
c
  800 return
c
c
  900 format(1h0,13x,6henergy, 8x,3hplo,12x,3hpup,12x,4hpavg,7x,5hsp((),
     1 7x,3herr,7x,5he pct,2x,5hq pct/) 
  901 format(1h ,i5,5x,1pe10.3,3e15.4,0pf8.1,1pe15.4,0pf6.1,f7.1)
  902 format(1h0,2x,3hbin,7x,5hp.hgt, 9x,1hb,14x,2hsb,13x,3hblo,12x,
     13hbup,12x,4hdelb,8x,5hp dev/)
  903 format(1h ,i5,5x,1pe10.3,5e15.4,0pf7.1)
  904 format(1h0,22hinconsistancy factor =,e16.8) 
  905 format(1h0,4x,3h***,3x,21hcondition indicator =,e15.7,3x,3h***)   ju 21 76
  906 format(1h1,'linear plot of b and  b-adjusted bounds',/)
  907 format(////1h0,'window no.',i4,' and its bounds -- window',
     1 ' energy =',1pe12.4/)
  908 format(1h1,20a4)
c
c
      end 
c
c
c***********************************************************************
c
c
      subroutine loop ( w,b,s,q,alo,aup,mr,nr,nc,philo,err,istar,ut,
c
     1                  ulo,elo,synwlo,que)
c
c-----------------------------------------------------------------------
c
      dimension w(nc),b(nr),s(nr),q(nc),alo(mr,nc),aup(mr,nc),istar(nc),
c
     1          ut(nc),ulo(nc),elo(nc),synwlo(nc) 
c
c-----------------------------------------------------------------------
c
c     loop for subroutine ferd
c
c     ulo = u(transpose)
c
      do 100 i=1,nr 
      ulo(i)=ut(i)
  100 continue
c
c     uta = ulo * alo if -
c     uta = ulo * aup if +
c
c     elo = w - u(trans) * a
      do 200 j=1,nc 
      uta=0.0
      do 175 i=1,nr 
      if(ulo(i))125,175,150
  125 uta=uta+ulo(i)*alo(i,j) 
      go to 175
  150 uta=uta+ulo(i)*aup(i,j) 
  175 continue
      elo(j)=w(j)-uta
  200 continue
c
c     find largest negative component, ejmin, of elo * q, and its index jmin
c
  210 ejmin=0.0
      do 240 j=1,nc 
      ejqj=elo(j)*q(j)
      if(ejqj)220,240,240
  220 if(ejqj-ejmin)230,240,240
  230 ejmin=ejqj
      jmin=j
  240 continue
      is=istar(jmin)
      if (ejmin) 255,340,340
c
c     remove istar row from window error, elo
c
  255 do 270 j = 1,nc
      if (ulo(is)) 260,265,265
  260  elo(j) = elo(j) + ulo(is)*alo(is,j)
       go to 270
  265  elo(j) = elo(j) + ulo(is)*aup(is,j)
  270 continue
c
c     compute new ulo(is) to make elo(jmin) = 0.0 
c
      if (elo(jmin)) 275,280,280
  275  ulo(is) =(elo(jmin)/alo(is,jmin))*1.001                          ma 10 76
       go to 290
  280  ulo(is) =(elo(jmin)/aup(is,jmin))*0.999                          ma 10 76
  290 continue
c
c     put new istar row back into elo
c
      do 320 j = 1,nc
      if (ulo(is)) 305,310,310
  305  elo(j) = elo(j) - ulo(is)*alo(is,j)
       go to 320
  310  elo(j) = elo(j) - ulo(is)*aup(is,j)
  320 continue
c
c     continue this process until the window inequality is satisfied    ma 10 76
c
      go to 210
c
c     calculate new synthethetic window l(e)
c
  340 do 380 j=1,nc 
      sum1=0.0
      do 370 i=1,nr 
      if(ulo(i))350,370,360
  350 sum1=sum1+ulo(i)*alo(i,j)
      go to 370
  360 sum1=sum1+ulo(i)*aup(i,j)
  370 continue
      synwlo(j)=sum1
  380 continue
c
c     calculate phi and (u*s)**2                                        ma 10 76
c
      ulosq=0.0
      philo=0.0
      do 400 i=1,nr 
      philo=philo+ulo(i)*b(i) 
      ulosq=ulosq+(ulo(i)*s(i))**2
  400 continue
c
c  Set value of  mu = 1.0  (amu = 1.0)  for 68.26 percent 1-at-a-time intervals
c
       amu = 1.0
c
      err = amu*sqrt(ulosq)
      philo=philo-err
      que=qshoot(elo,q,nc)
      return
      end 
c
c
c************************************************************************
c
c
      function qshoot ( d, q, nc )
c
c------------------------------------------------------------------------
c
      dimension  d(nc), q(nc) 
c
c     function qshoot gives an estimate of error due to
c     quadratic undershoot of the synthetic window between
c     the comparison points . 
c
c-----------------------------------------------------------------------
c
      qshoot=0.0
      nn=nc-1
      do 40 i=2,nn
      der1=0.5*(d(i+1)-d(i-1))
      der2=d(i-1)-2.0*d(i)+d(i+1)
c     second derivative must be positive for a minimum to occur
      if(der2)40,40,10
   10 xmin=-der1/der2
c     does minimum occur within range
      if(abs(xmin)-1.0)20,40,40
   20 value=d(i)-(der1**2)/(2.0*der2)
c     has an undershoot occurred i. e., is minimum negative 
      if(value)30,40,40
   30 qshoot=qshoot-value*q(i)
   40 continue
      qshoot=4.0/9.0*qshoot
c
      return
c
c
      end 
