c
c*******************************************************************************
c
      subroutine CFRDIT ( A,      mr,     mc,     AQWRK,  nr,    nc,  
c
     1                    b,      WIND,           mw,     xlo,   xup,
c
     2                    philo,          phiup,          phiest,
c
     3                    x,      qdpwk,  U,      maxmn,  VT,    sig, 
c
     4                    wt,     SigUTb, VTQm1d, VTQwsq, sigsq, Qinvd,
  
     5                    ajnorm, sg2ptn, amusq,  itt,    nw,    svdwrk,
c
     6                    lsvwrk, p,      q,      ipr,    iprnt, mou,
c
     7                    UTb,    EVwork  )        
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      logical  ipr
c
      external  FERDIT, VECPR, VECPRD
c
c**** double precision  Qinvd(maxmn), sg2ptn(mc),
      dimension         Qinvd(maxmn), sg2ptn(mc),
     1                  xlo(mc),      xup(mc),      qdpwk(mc), 
     2                  VTQwsq(mc),   SigUTb(mc),   VTQm1d(mc),
     3                  sigsq(mc),    UTb(maxmn),   EVwork(mc)
c
c**** double precision  dble
c
      dimension  A(mr,mc),       AQWRK(mr,mc), VT(maxmn,maxmn),
     1           U(maxmn,maxmn), WIND(mc,mw),  svdwrk(lsvwrk)
c
      dimension  b(mr),       philo(mw), phiup(mw),  phiest(mw), x(mc),
     1           sig(maxmn ), wt(mc),    ajnorm(mc), p(mc),      q(mc) 
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      if ( ipr ) then
                      iw10 = 1
                 else
                      iw10 = 0
      end if
c
      do 400 j=1,nc 
         xlo(j) = p(j)
         xup(j) = q(j) 
  400 continue
c
      if ( iprnt .gt. 0 ) then
                           print 430
                  if (ipr) write (mou,430)
c                --------------------------------------------
c****             call VECPRD ( xlo, nc, 'p0  ', iw10, mou )
                  call VECPR  ( xlo, nc, 'p0  ', iw10, mou )
c                --------------------------------------------
c****             call VECPRD ( xup, nc, 'q0  ', iw10, mou )
                  call VECPR  ( xup, nc, 'q0  ', iw10, mou )
c                --------------------------------------------
      end if
c
  430 format (' '/' Initial FERDIT bounds on x')
c
c   ----------------------------------------------------------
      call FERDIT ( itt,   A,      AQWRK,  WIND,   b,      
c
     1              mr,    mc,     nr,     nc,     nw,
c
     2              maxmn, mw,     amusq,  xlo,    xup,      
c
     3              philo, phiup,  phiest, wt, 
c
     4              qdpwk, x,      U,      VT,     sig,
c
     5              sigsq, SigUTb, VTQm1d, VTQwsq, EVwork, 
c
     6              Qinvd, ajnorm, sg2ptn, svdwrk, lsvwrk,
c
     7              ipr,   iprnt,  mou,    UTb   ) 
c   -----------------------------------------------------------
c
      if ( iprnt .gt. 0 ) then
                           print 120
  120               format (' '/' '/' Final FERDIT bounds on the phi')
                           if (ipr) write (mou,120)
c                         ----------------------------------------------
                           call VECPR ( philo,  nw, 'phlo', iw10, mou )
c                         ----------------------------------------------
                           call VECPR ( phiest, nw, 'phes', iw10, mou )  
c                         ----------------------------------------------
                           call VECPR ( phiup,  nw, 'phhi', iw10, mou )
c                         ----------------------------------------------
      end if
c
      return
c
      end 
c
c
c*******************************************************************************
c
      subroutine FERDIT ( itt,    A,      AQWRK,  WIND,   b,      
c
     1                    mr,     mc,     m,      n,      nwnd,
c
     2                    maxmn,  mw,     musq,   xlo,    xup,    
c
     3                    philo,  phiup,  phiest, w,
c
     4                    qdp,    x,      U,      VT,     sig,
c
     5                    sigsq,  SigUTb, VTQm1d, VTQwsq, EVwork,  
c
     6                    Qinvd,  ajnorm, sg2ptn, svdwrk, lsvwrk,
c
     7                    ipr,    iprint, mou,    UTb     )
c
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c     AQWRK is the A-matrix (already regularized and scaled with the errors)
c         b is the rhs
c       qdp is for (hopefully) shrinking upper bounds
c         x is estimated solution
c         w is window vector
c         m is number of rows of AQWRK
c         n is number of cols of AQWRK
c        mr is first dimension of AQWRK, m is # rows of AQWRK
c        mc is second dim of AQWRK, n is number of columns
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      external  ARRAY,  EVAL,   VECPR,  VECPRD, WDTH,   ZERO
c
      logical  ipr
c
c**** double precision  qdp(mc),      xlo(mc),    xup(mc),
      dimension         qdp(mc),      xlo(mc),    xup(mc),
     1                  Qinvd(maxmn), VTQwsq(mc), VTQm1d(mc), 
     2                  SigUTb(mc),   sg2ptn(mc), sigsq(mc), 
     3                  UTb(maxmn),   EVwork(mc)
c
c**** double precision  qtemp, Qdtmp, pnew,  qnew,   Aij,   sumz, vi,
      real              qtemp, Qdtmp, pnew,  qnew,   Aij,   sumz, vi,
c****1                  dble,  dmax1, dmin1, VTQdsq, bTU2b, sumr,
     1                                       VTQdsq, bTU2b, sumr,
     2                  xj,    VTij,  zi,    ri,     sgi,   tnd
c
      real  musq 
c
      dimension  A(mr,mc),    AQWRK(mr,mc),   svdwrk(lsvwrk), 
     1           WIND(mc,mw), U(maxmn,maxmn), VT(maxmn,maxmn)
c 
      dimension  x(mc),     b(mr),     sig(maxmn),  w(mc),  
     1           phiup(mc), philo(mc), phiest(mc),  ajnorm(mc)
c
c-------------------------------------------------------------------------------
c
      if ( ipr ) then
                      iw10 = 1
                 else
                      iw10 = 0
      end if
c
      eps = r1mach( 4 )
c
      tol = 5.0e-2
c
      minmn = min0( m , n )
c
c  Find initial center of the Q-box and compute   Q**-1 d  
c
      do 50 j=1,n
c
c****    qtemp = ( xup(j) - xlo(j) ) / 2.0d0
         qtemp = ( xup(j) - xlo(j) ) / 2.0
         qdp(j) = qtemp
c                      -- qdp is half-lengths of sides of Q-box
c
c****    Qdtmp = ( xup(j) + xlo(j) ) / 2.0d0
         Qdtmp = ( xup(j) + xlo(j) ) / 2.0
         Qdtmp = Qdtmp / qtemp
c
         Qinvd(j) = Qdtmp
c                        -- Qinvd is the vector  Q**-1 d
   50 continue
c
      if ( iprint .gt. 1 ) then
                           write(mou,51) eps, tol
   51         format(' '/' In FERDIT,    eps =', 1pe11.4 /
     1                   '               tol =', 1pe11.4  )
c****         call VECPRD ( qdp,   n, 'Q   ', iw10, mou )
c****         call VECPRD ( Qinvd, n, 'Q-1d', iw10, mou )
              call VECPR  ( qdp,   n, 'Q   ', iw10, mou )
              call VECPR  ( Qinvd, n, 'Q-1d', iw10, mou )
      end if
c
c  Form  AQWRK = A * diag( q_1 , q_2 , ... , q_n )
c
      do 60 i=1,m
      do 60 j=1,n
c****  Aij = dble( A(i,j) )
       Aij =       A(i,j) 
c****  AQWRK(i,j) = sngl( Aij * qdp(j) )
       AQWRK(i,j) =       Aij * qdp(j) 
   60 continue
c
c  ajnorm(j) is the Euclid norm of j'th column of  AQWRK 
c
      do 65 j=1,n
       ajnorm(j) = snrm2 ( m, AQWRK(1,j), 1 )
   65 continue
      if( iprint .gt. 1 ) 
     1                    call VECPR ( ajnorm, n, 'anrm', iw10, mou )
c
c
c++++++++++++++++++++ Begin FERDIT iteration loop ++++++++++++++++++++++++++++++
c
c
      do 1000 itr=1,itt
c
c
         if ( iprint .gt. 0 ) then
               print 1002, itr
 1002          format(' '/' '/'0******************** Beginning FERDIT',
     1                       ' iteration ',i3,'   ***************' )
               if (ipr) write(mou,1002) itr
         end if
c                         
         if ( iprint .gt. 3 ) 
     1                call ARRAY ( AQWRK, m, n, mr, 'A*Q ', iw10, mou )
c                       
c
c  Calculate singular value decomp. of AQWRK using LAPACK subrout.  SGESVD
c
         if ( ipr ) print 19
   19    format(' '/' Computing singular value decomp. in  FERDIT' )
c
c   ----------------------------------------------------------
         call SGESVD ( 'A',    'A',    m,     n,  AQWRK,  mr,
c
     1                 sig,    U,      maxmn, VT, maxmn, 
c
     2                 svdwrk, lsvwrk, info   )
c   ----------------------------------------------------------
c
      if ( iprint .gt. 0 ) then
               print 12
               if (ipr) write (mou,12) 
               call ARRAY ( sig, minmn, 1, minmn, 'sig ', iw10, mou )
      end if
c
   12 format (' '/' '/' Singular values computed in  FERDIT' ) 
c
      if ( info .ne. 0 ) then
              print 9,info 
    9         format (' ' /' '/' SVD (SGESVD) failed in FERDIT ',
     1                ' **********  info  =',i5, '  ************'/)
              if (ipr) write(mou,9) info
              if (info .gt. 0)
     1                call VECPR ( svdwrk(2), minmn, 'e   ', iw10, mou )
c            --------
              return
c            --------
      end if
c
      if ( m .lt. n ) then
c
         do 14 i=m+1,n
            sig(i) = 0.0
   14    continue
c
      end if
c
c
c  Transpose the matrix in the array  U ,  i.e., set  U = U**T
c
      mm1 = m - 1
      do 11 i=1,mm1 
         ip1 = i + 1
         do 10 j=ip1,m
            Uij = U(i,j)
            U(i,j) = U(j,i)
            U(j,i) = Uij
   10    continue
   11 continue
c
c  Form single precision  AQWRK = A * diag( q_1 , q_2 , ... , q_n )
c
      do 13 i=1,m
      do 13 j=1,n
c****   Aij = dble( A(i,j) )
        Aij =       A(i,j) 
        qtemp = qdp(j)
        Aij = Aij*qtemp
c****   AQWRK(i,j) = sngl( Aij )
        AQWRK(i,j) =       Aij 
   13 continue
c
c   Form the vectors  UTb = U**T * b  and  SigUTb = Sig * U**T * b 
c
c                  ( Array U contains matrix  U**T )
c
      do 18 i=1,m
c****    sumz = 0.0d0
         sumz = 0.0
         do 17 j=1,m
c****        sumz = sumz + dble( U(i,j) ) * dble( b(j) )
             sumz = sumz +       U(i,j)   *       b(j) 
   17    continue
         UTb(i) = sumz
c****    if ( i .le. minmn ) SigUTb(i) = sumz * dble( sig(i) )
         if ( i .le. minmn ) SigUTb(i) = sumz *       sig(i) 
   18 continue
c
      if ( m .lt. n ) then
            do 23 j=m+1,n
c****          SigUTb(j) = 0.0d0
               SigUTb(j) = 0.0
   23       continue
      end if
c
c****  if ( iprint .gt. 1 )  call VECPRD ( SigUTb, n, 'SUTb', 1, mou )
       if ( iprint .gt. 1 )  call VECPR  ( SigUTb, n, 'SUTb', 1, mou )
c
c   Form  bTU2b = b**T * U_2 * U_2**T * b   where  U = ( U_1 , U_2 )
c
c**** bTU2b = 0.0d0
      bTU2b = 0.0
      if ( m .gt. n ) then
             do 21 i=n+1,m
                bTU2b = bTU2b + UTb(i)**2
   21        continue
      end if
      if ( iprint .gt. 1 ) then
                           print 22, bTU2b
   22        format (' '/' '/' || U_2 * b ||**2  =  ',
     1                       'b**T * U_2 * U_2**T * b  =', 1pd16.8 )
             if ( ipr ) write(mou,22) bTU2b
      end if
c                                              minmn        2
c   Form  VTQm1d = V**T * Q**-1 * d ,  VTQdsq = SUM (VTQm1d)
c                                               i=1         i
c                                     and  
c                                     vector sigsq of squared singular values
c
c**** VTQdsq = 0.0d0 
      VTQdsq = 0.0 
      do 30 i=1,n
c****    sumr = 0.0d0
         sumr = 0.0
         do 25 j=1,n
c****       sumr = sumr + dble( VT(i,j) ) * Qinvd(j)
            sumr = sumr +       VT(i,j)   * Qinvd(j)
   25    continue
         VTQm1d(i) = sumr
         if ( i .le. minmn )  VTQdsq = VTQdsq + sumr**2
c****    sigsq(i) = ( dble( sig(i) ) )**2
         sigsq(i) = (       sig(i)   )**2
   30 continue
c
c****   if ( iprint .gt. 1 )   call VECPRD ( VTQm1d, n, 'VTQd', 1, mou )
        if ( iprint .gt. 1 )   call VECPR  ( VTQm1d, n, 'VTQd', 1, mou )
c
c   Now loop thru the nwnd window vectors, finding best tau for each
c
c   nfail is a counter telling how many times in an
c     iteration we failed to find tau.  If we fail for all
c       of the window vectors we quit.
c
      nfail = 0
c
c
c+++++++++++++++ Enter window vector loop ++++++++++++++++++++++++++++++++++++++
c
      tauav = 0.0
c
c
      do 950 nw=1,nwnd
c
c
        icnt = 0
c
      do 400 i=1,n
       w(i) = WIND(i,nw)
  400 continue
c
      if ( ipr ) then
                 print 410, nw
  410         format (' '/' '/' beginning window', i4,
     1                '  ---------------------------------' )
                 if ( iprint .gt. 0 ) write (mou,410) nw
      end if
c
c  Choose upper bound for tau  
c
      sum = 0.0
      do 300 j=1,n
         sum = sum + ajnorm(j) * abs( w(j) )
  300 continue
c
      sighi = sum
c
c  Set lower bound for tau to  10.0 * eps_mach * (largest singular value).
c
      siglo = 10.0 * eps * sig(1)
c
       savelo = siglo
       savehi = sighi
c
c  Form  (v)_i**2 = (VT * q * w)_i**2
c
      do 40 i=1,n
c****    vi = 0.0d0
         vi = 0.0
         do 35 j=1,n
c****       vi = vi + dble( VT(i,j) ) * qdp(j) * dble( w(j) )
            vi = vi +       VT(i,j)   * qdp(j) *       w(j) 
   35    continue
         VTQwsq(i) = vi**2
   40 continue
c**** if ( iprint .gt. 1 ) call VECPRD ( VTQwsq, n, 'VTQw', iw10, mou )
      if ( iprint .gt. 1 ) call VECPR  ( VTQwsq, n, 'VTQw', iw10, mou )
c
c
   42 continue
c
      if ( siglo .lt. savelo ) savelo = siglo
      if ( sighi .gt. savehi ) savehi = sighi
c
c   --------------------------------------------------------------
      call ZERO ( siglo,  sighi,  eps,    tol,    tau,
c
     1            fval,   iflag,  AQWRK,  
c
     2            qdp,            VT,     maxmn,  sig,
c
     3            SigUTb, VTQm1d, UTb,    musq,   mr,
c
     4            mc,     m,      n,      nw,     VTQwsq,
c
     5            sigsq,  elsq,   bTU2b,  ajnorm, EVwork,  
c
     6            VTQdsq, sg2ptn, ipr,    iprint, mou     ) 
c   ---------------------------------------------------------------
c
      if ( iprint .gt. 1 )  then
                                  print 39, iflag, tau, fval
   39      format (' '/' On return from  ZERO , iflag =',i2 /
     1                 '                          tau =', 1pe13.5/
     2                 '                        dl/dt =', 1pe13.5  )
             if (ipr) write (mou,39) iflag, tau, fval
      end if
c
c
c------  iflag=0  Indicates a normal return from ZERO.
c                 In this case the zero is returned in  tau  and the
c                 corresponding function value, which should be very small,
c                 is returned in  fval .
c
      if ( iflag .eq. 0 )  go to 43
c
c------  iflag=1  Indicates that the two function values had same sign.
c                 The subroutine returns with  tau  set to the end point
c                 ( siglo  or  sighi ) which give the smallest absolute
c                 value for the function.  The corresponding function value
c                 is returned in  fval .
c
c                 This indicates one of the following possibilities:  
c
c                 1. a zero exists, but the initial interval (siglo,sighi)
c                    did not enclose it.  (both function values may be either
c                    positive or negative)
c
c                 2. the problem is well-conditioned, so tau=0 is best value.
c                    (both function values must be positive)
c
c                 3. the Q-box bounds are the best obtainable, so tau should
c                    be very large.  (both function values must be negative)
c
c           if the devivative evaluated at  siglo  is: 
c
c                 a. positive, then 1. or 2. is true,
c
c                 b. negative, then 1. or 3. is true.
c
      if ( iprint .gt. 1 ) then
                                print 41, siglo, sighi
   41          format (' '/' Initial Interval  ( a , b ) = (',
     1               1pe13.5,' ,',1pe13.5,' )  did not bracket a zero')
               if (ipr) write (mou,41) siglo, sighi
      end if
c
      if ( iprint .gt. 2 ) write(mou,46)
   46    format (' '/' In FERDIT, No. 1 call to EVAL')
c
c   --------------------------------------------------------------
      call EVAL ( siglo,  dsiglo, AQWRK,  
c
     1            qdp,            VT,     maxmn,  sig,
c
     2            SigUTb, VTQm1d,         musq,   mr,
c
     3            mc,     m,      n,      nw,     VTQwsq,
c
     4            sigsq,  elsq,   bTU2b,  ajnorm, UTb,   EVwork,
c
     5            VTQdsq, sg2ptn, ipr,    iprint, mou   )
c   ---------------------------------------------------------------
c
c  If  dsiglo  is positive, reset  siglo = siglo / 2.0
c           and try again ( up to 10 times ).  If no bracket is found in
c             10 tries, problem is probably well conditioned.  in this 
c                case, take  tau = siglo
c
      if ( dsiglo .gt. 0.0 ) then
                                  siglo = siglo / 2.0
                                  icnt = icnt + 1 
                                  if ( icnt .le. 10 ) go to 42
                                  tau = siglo
                                  go to 44
      end if
c
      if ( dsiglo .le. 0.0 ) then
c
c  Otherwise, increase  sighi  and try again ( up to 5 times )
c          after 5 unsuccessful tries, set  tau = sighi  because the
c          the  Q-box bounds are probably the best available
c
                     sighi = 10.0 * sighi
                     icnt = icnt + 1
                     if ( icnt .le. 5 ) go to 42
                     tau = sighi
c
      end if
c
   44 continue
c
      if ( iprint .gt. 1 ) then
                                print 2001, nw, savelo, savehi
 2001    format(' For window ',i3,' No Zero in the Interval  (',
     1                                  1pe13.5, ' ,', 1pe13.5, ' )')
             if (ipr) write (mou,2001) nw, savelo, savehi
             print 2002, tau
 2002        format(' '/'    tau set to   tau =', 1pe13.5)
             if (ipr) write (mou,2002) tau
      end if
c
      if ( iprint .gt. 1 )  write(mou,48)
   48    format (' '/' In FERDIT, No. 2 call to EVAL')
c
      nfail=nfail+1 
c
c   ---------------------------------------------------------------
      call EVAL ( tau,    val,    AQWRK,  
c
     1            qdp,            VT,     maxmn,  sig,
c
     2            SigUTb, VTQm1d,         musq,   mr,
c
     3            mc,     m,      n,      nw,     VTQwsq,
c
     4            sigsq,  elsq,   bTU2b,  ajnorm, UTb,   EVwork,
c
     5            VTQdsq, sg2ptn, ipr,    iprint, mou   )
c   ---------------------------------------------------------------
c
c
c+++++  Transfer to Here on Normal Return from ZERO  +++++ 
c
  43  continue
c
c
c++++  Now compute biased estimate of  x  and bounds on   phi = w**T x
c
c------  x' = V ( Sig**2 + (tau/n)I )**-1 ( Sig U**T b + (tau/n) V**T Q**-1 d )
c
c------  x = Q x'
c
c**** tnd = dble( tau ) / dble( n )
      tnd =       tau   /       n 
      do 20 j=1,n
c****    xj = 0.0d0 
         xj = 0.0
        do 15 i=1,n 
c****       VTij = dble( VT(i,j) )
            VTij =       VT(i,j) 
            zi = SigUTb(i)
            ri = VTQm1d(i)
            sgi = sg2ptn(i) 
            xj = xj + ( VTij*(zi + tnd*ri) ) / sgi
   15   continue
c****    x(j) = sngl( xj * qdp(j) )
         x(j) =       xj * qdp(j) 
   20 continue
c
      if ( iprint .gt. 1 ) call ARRAY ( x, n, 1, mc, 'x   ', iw10, mou )
c
c------ phi = w**T x
c
      sum = 0.0
      do 47 j=1,n
         sum = sum + w(j)*x(j)
   47 continue
      phiest(nw) = sum
c
c------ Call subr. WDTH to compute width of confidenced interval for phi
c
      if ( iprint .gt. 1 ) write(mou,52)
   52   format (' '/' In FERDIT, No. 1 call to WDTH')
c
c   ---------------------------------------------------------
      call WDTH ( wid,    nw,  VTQwsq, mc,      n,
c
     1            tau,                 VTQm1d,  SigUTb,
c
     2            sg2ptn, ipr, iprint, mou,     m  )
c   ----------------------------------------------------------
c
      phunew = sum + wid
c
      phlnew = sum - wid
c
      if ( iprint .gt. 1 ) then
                                print 1001, wid
 1001         format (' '/' In FERDIT,   L**2 =', 1pe11.4 )
              if (ipr) write (mou,1001) wid
      end if
c
      phiup(nw) = amin1 ( phiup(nw), phunew )
c
      philo(nw) = amax1 ( philo(nw), phlnew )
c
      if ( iprint .gt. 1 ) then
                       print 1004, philo(nw), phiest(nw), phiup(nw)
 1004        format(' '/' In FERDIT,   philo =', 1pe12.4 /
     1                  '             phiest =', 1pe12.4 /
     2                  '              phiup =', 1pe12.4   ) 
            if (ipr) write(mou,1004) philo(nw), phiest(nw), phiup(nw)
      end if
c
c++++  Update xup, xlo, bounds on x
c
      do 80 j=1,n
c 
        if ( iprint .gt. 1 )  write(mou,74)
   74   format(' '/' In FERDIT, No. 2 call to WDTH (repeat n times)') 
c
        do 75 jj=1,n
c****      VTQwsq(jj) = ( dble( VT(jj,j) ) * qdp(j) )**2
           VTQwsq(jj) = (       VT(jj,j)   * qdp(j) )**2
   75   continue
c
c
c   -----------------------------------------------------------
        call WDTH ( wid,    nw,   VTQwsq,  mc,     n,
c
     1              tau,                   VTQm1d, SigUTb,
c
     2              sg2ptn, ipr,  iprint,  mou,    m  )
c   ------------------------------------------------------------
c
c****   pnew = dble( x(j) ) - dble( wid )
        pnew =       x(j)   -       wid 
c****   qnew = dble( x(j) ) + dble( wid )
        qnew =       x(j)   +       wid 
c
c****   xlo(j) = dmax1( xlo(j), pnew )
        xlo(j) = amax1( xlo(j), pnew )
c****   xup(j) = dmin1( xup(j), qnew )
        xup(j) = amin1( xup(j), qnew )
c
   80 continue
c
c
      if ( iprint .gt. 1 ) then
                               print 1005
 1005        format ( ' '/ 'In FERDIT,  Q - box  Updated' )
                            if ( ipr ) write (mou,1005)
      end if
c
c++++  Check whether new bounds on solution vector  x  give better bounds
c      for  phi = w**T x.  If so, Update  philo,  phiup 
c
      sumhi = 0.0
      sumlo = 0.0
c
      do 81 j=1,n
c
         if( w(j) .le. 0.0 ) go to 83
c
         sumhi = sumhi + w(j)*xup(j)
         sumlo = sumlo + w(j)*xlo(j)
         go to 81
c
   83    continue
         sumhi = sumhi + w(j)*xlo(j)
         sumlo = sumlo + w(j)*xup(j)
c
   81 continue
c
      phiup(nw) = amin1 ( phiup(nw), sumhi )
      philo(nw) = amax1 ( philo(nw), sumlo )
c
      phiest(nw) = amin1 ( phiest(nw), sumhi )
      phiest(nw) = amax1 ( phiest(nw), sumlo )
c
      if ( iprint .gt. 1 ) then
                         print 1004, philo(nw), phiest(nw), phiup(nw)
                         if (ipr) write(mou,1004) philo(nw),
     1                                          phiest(nw), phiup(nw)
      end if
c
c++++ Update average value of tau
c
      tauav = ( float( nw - 1 )*tauav + tau ) / float( nw )
c
      if ( iprint .gt. 0 ) then
                               print 1003, tau, tauav
 1003     format(' '/'    tau=', 1pe11.4, '   tau av =', 1pe11.4 )
              if (ipr) write(mou,1003) tau, tauav
      end if
c
      if ( (nw .eq. nwnd) .and. (iprint .gt. 1) ) then
            print 2003
 2003       format(' '/' '/' '/' '/' ',8x,'w',6x,'Solution x',5x,
     1                             'Lower bd xlo',4x,'Upper bd xup')
            if (ipr) write(mou,2003)
            do 82 j=1,n
                print 2004, j, w(j), x(j), xlo(j), xup(j)
 2004           format(' ',i3,f8.2,1pe16.7,1p2d16.7)
                if (ipr) write(mou,2004) j, w(j), x(j), xlo(j), xup(j)
   82       continue
      end if
c
c
  950 continue
c
c
c+++++++++++++  End of window loop  ++++++++++++++++++++++++++++++++++++++++++++
c
c++++ Compute a final point estimate of the biased solution vector using 
c        the last obtained value of tau (value for the final window)
c
cccccc      tau = tauav
cccccc------------------activate above statement to use average value of tau
c
c------  x' = V ( Sig**2 + (tau/n)I )**-1 ( Sig U**T b + (tau/n) V**T Q**-1 d )
c
c------  x = Q x'
c
c**** tnd = dble( tau ) / dble( n )
      tnd =       tau   /       n 
      do 1020 j=1,n
c****    xj = 0.0d0 
         xj = 0.0 
        do 1015 i=1,n 
c****       VTij = dble( VT(i,j) )
            VTij =       VT(i,j) 
            zi = SigUTb(i)
            ri = VTQm1d(i)
            sgi = sg2ptn(i) 
            xj = xj + ( VTij*(zi + tnd*ri) ) / sgi
 1015   continue
c****    x(j) = sngl( xj * qdp(j) )
         x(j) =       xj * qdp(j) 
 1020 continue
c
      if ( iprint .gt. 1 ) call ARRAY ( x, n, 1, mc, 'x   ', iw10, mou )
c
      if ( iprint .gt. 0 ) then
             print 2006, itr
 2006        format(///' ','Window',5x,'philo',11x,'phiest',10x,'phiup',
     1                                            8x,'for',' iter',i4 )
             if (ipr) write(mou,2006) itr
             do 955 nw=1,nwnd
                print 2005, nw, philo(nw), phiest(nw), phiup(nw)
 2005           format(' ',2x,i4,5x,3(1pe11.4,5x))
                if (ipr) write(mou,2005) nw,philo(nw),
     1                                            phiest(nw),phiup(nw)
  955        continue
      end if
c
c++++ Update  Q-values ,  Q**-1 * d ,   A * Q
c
      do 90 j=1,n
c
c****    qtemp = ( xup(j) - xlo(j) ) / 2.0d0
         qtemp = ( xup(j) - xlo(j) ) / 2.0
         qdp(j) = qtemp 
c
c****    Qdtmp = ( xup(j) + xlo(j) ) / 2.0d0
         Qdtmp = ( xup(j) + xlo(j) ) / 2.0
         Qdtmp = Qdtmp / qtemp
c
         Qinvd(j) = Qdtmp
c
         do 85 i=1,m
c****       Aij = dble( A(i,j) )
            Aij =       A(i,j) 
            Aij = Aij * qtemp 
c****       AQWRK(i,j) = sngl( Aij )
            AQWRK(i,j) =       Aij 
   85    continue
c
   90 continue
c
      do 92 j=1,n
         ajnorm(j) = snrm2 ( m, AQWRK(1,j), 1 ) 
   92 continue
c
      if ( iprint .gt. 1 ) then
              write (mou,93) 
   93        format(' '/' '/'Updating Q-values in FERDIT,  ')
c****         call VECPRD ( qdp,    n, 'Q   ', iw10, mou )
c****         call VECPRD ( Qinvd,  n, 'Q-1d', iw10, mou )
              call VECPR  ( qdp,    n, 'Q   ', iw10, mou )
              call VECPR  ( Qinvd,  n, 'Q-1d', iw10, mou )
c
              call VECPR  ( ajnorm, n, 'anrm', iw10, mou )
      end if
c
c                            --------
      if ( nfail .eq. nwnd )  return
c                            --------
c
c                             - Failed to find tau for all the window vectors
c
 1000 continue
c
c
c+++++++++++++++ End of FERDIT iteration loop ++++++++++++++++++++++++++++++++++
c
c
      if ( iprint .gt. 0 ) then
                                print 1010
 1010           format (//' End of FERDIT Iterations '/)
                               if (ipr) write (mou,1010)
      end if
c
c
      return
c
      end 
c
c
c*******************************************************************************
c
c
      subroutine EVAL ( tau,    val,    AQWRK,  
c
     1                  qdp,            VT,     maxmn,  sig,
c
     2                  SigUTb, VTQm1d,         musq,   mr, 
c
     3                  mc,     m,      n,      nw,     VTQwsq,
c
     4                  sigsq,  elsq,   bTU2b,  ajnorm, UTb,   EVwork,
c
     5                  VTQdsq, sg2ptn, ipr,    iprint, mou   )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      common /evalbl/ taud, musqd, rho
c
      logical  ipr
c
      real  musq
c
c**** double precision  rho,   taud,  tnd,   musqd,  tnosg2, sgosg2,
      real              rho,   taud,  tnd,   musqd,  tnosg2, sgosg2,
     1                  sum1,  sum2,  sum3,  sum4,   sum5,   sum6,
     2                  sum7,  zptnr, vsq,   zptnr2, sg2tn2,
c****3                  vald,  dble,  bTU2b, VTQdsq 
     3                  vald,         bTU2b, VTQdsq 
c
c**** double precision  VTQwsq(mc), SigUTb(mc),  
      dimension         VTQwsq(mc), SigUTb(mc),  
     1                  VTQm1d(mc), UTb(maxmn),   EVwork(mc),
     2                  qdp(mc),    sigsq(maxmn), sg2ptn(mc)
c
      dimension  AQWRK(mr,mc), ajnorm(mc), VT(maxmn,maxmn), sig(maxmn)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      save
c
c**** taud = dble( tau )
      taud =       tau 
c**** musqd = dble( musq )
      musqd =       musq 
c**** tnd = taud / dble( n )
      tnd = taud /       n 
c
c**** sum1 = 0.0d0
c**** sum2 = 0.0d0
c**** sum3 = 0.0d0
c**** sum4 = 0.0d0
c**** sum5 = 0.0d0
      sum1 = 0.0
      sum2 = 0.0
      sum3 = 0.0
      sum4 = 0.0
      sum5 = 0.0
c
      minmn = min0( m, n )
c
      do 10 i=1,n
c
         zptnr  = SigUTb(i) + tnd * VTQm1d(i)
         zptnr2 = zptnr**2 
c
         vsq = VTQwsq(i)
c
         sg2ptn(i) = sigsq(i) + tnd
         sg2tn2 = sg2ptn(i)**2
c
         sum2 = sum2 + vsq / sg2tn2
c
         if ( i .le. minmn ) then
c
               sum3 = sum3 + zptnr * VTQm1d(i) / sg2ptn(i)
c
               sum4 = sum4 + zptnr2 / sg2tn2
c
         end if
c
         sum5 = sum5 + vsq / sg2ptn(i)
c
   10 continue
c
c++++  Now compute rho
c
c**** sum6 = 0.0d0
c**** sum7 = 0.0d0
      sum6 = 0.0
      sum7 = 0.0
c
      do 13 i=1,minmn
c
         tnosg2 = tnd / sg2ptn(i)
c****    sgosg2 = dble( sig(i) ) / sg2ptn(i)
         sgosg2 =       sig(i)   / sg2ptn(i)
c
         EVwork(i) = tnd * sig(i) * VTQm1d(i) / sg2ptn(i)
c
         sum6 = sum6 + ( EVwork(i) - tnosg2*UTb(i) )**2
c 
         sum7 = sum7 + ( sgosg2*UTb(i) 
c
c****1                     - dble( sig(i) )*sgosg2*VTQm1d(i) )**2 
     1                     -       sig(i)  *sgosg2*VTQm1d(i) )**2 
c
   13 continue
c
      rho = sum6 + tnd*sum7
c
      if ( m .gt. n ) rho = rho + bTU2b 
c
c--- Compute   sum1 = ( confidence interval width scaling constant )
c
      sum1 = taud + musqd - rho
c
      if ( iprint .gt. 1 ) then
                           print 991, taud, sum6, sum7, rho, 
     1                                sum1, sum2, sum3, sum4, sum5
  991         format(' '/ 10x, ' in EVAL,  taud =', 1pd16.8,
     1                  / 10x, '           sum6 =', 1pd16.8,
     1                  / 10x, '           sum7 =', 1pd16.8,
     2                  / 10x, '            rho =', 1pd16.8,
     3                  / 10x, '           sum1 =', 1pd16.8 /
     4                    10x, '           sum2 =', 1pd16.8 /
     5                    10x, '           sum3 =', 1pd16.8 /
     6                    10x, '           sum4 =', 1pd16.8 /
     7                    10x, '           sum5 =', 1pd16.8  )
              if (ipr) write(mou,991) taud, sum6, sum7, rho, 
     1                                sum1, sum2, sum3, sum4, sum5
      end if 
c
      sum4 = sum4 + VTQdsq 
c
c  Compute  elsq ,  the interval width L**2 corresponding to tau
c
c**** elsq = sngl( sum1 * sum5 )
      elsq =       sum1 * sum5 
c
c**** sum1 = - sum1 / dble( n )
      sum1 = - sum1 /       n 
      sum1 = sum1 * sum2
c
c**** sum4 = - sum4 / dble( n )
      sum4 = - sum4 /       n 
c
c**** sum3 = ( sum3 * 2.0d0 ) / dble( n )
      sum3 = ( sum3 * 2.0   ) /       n 
c
      if ( iprint .gt. 1 ) then
                               write(mou,1002) sum1, sum4, sum3
 1002       format(' '/ 10x, 'In EVAL, sum1=', 1pd16.8 /
     1                  10x, '         sum4=', 1pd16.8 /
     2                  10x, '         sum3=', 1pd16.8  )
      end if
c
c**** sum3 = sum3 + sum4 + 1.0d0
      sum3 = sum3 + sum4 + 1.0
      sum3 = sum3 * sum5
c
      vald = sum1 + sum3
c**** val = sngl( vald )
      val =       vald 

c
      if ( iprint .gt. 1 ) then
                                print 1001, tau, val, elsq
 1001     format(' '/ 10x, 'In EVAL,  tau =', 1pe12.4 /
     1                10x, '        dL/dt =', 1pe12.4 /
     2                10x, '         L**2 =', 1pe12.4  )
             if (ipr) write(mou,1001) tau, val, elsq
      end if
c
c
      return
c
      end 
c
c
c*******************************************************************************
c
c
      subroutine WDTH ( wid,    nw,   VTQwsq, mc,      n, 
c
     1                  tau,                  VTQm1d,  SigUTb,
c
     2                  sg2ptn, ipr,  iprint, mou,     m  ) 
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      external  VECPRD
c
      common /evalbl/ taud, musqd, rho
c
      logical  ipr
c
c**** double precision  taud, musqd, rho,  widd, sum5, dsqrt
      real              taud, musqd, rho,  widd, sum5
c
c**** double precision  sg2ptn(mc), VTQwsq(mc), SigUTb(mc), VTQm1d(mc)
      dimension         sg2ptn(mc), VTQwsq(mc), SigUTb(mc), VTQm1d(mc)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      widd = taud + musqd - rho 
c
      if ( iprint .gt. 1 ) then
                           print 991, taud, musqd, rho, widd
  991         format(' '/ 10x, ' in WDTH,  taud =', 1pd16.8,
     1                  / 21x,            'musq =', 1pd16.8,
     2                  / 21x,            ' rho =', 1pd16.8,
     3                  / 21x,            'widd =', 1pd16.8  )
              if (ipr) write(mou,991) taud, musqd, rho, widd
      end if 
c
c**** if ( widd .lt. 0.0d0 ) then
      if ( widd .lt. 0.0   ) then
                                 print 1000, taud, musqd, rho, widd
 1000         format(' '/' Stoppimg in Subroutine WDTH' /' '/
     1                   '  taud = ', 1pd11.4 /
     2                   '  musq = ', 1pd11.4 /
     3                   '   rho = ', 1pd11.4 /
     4                   '  widd = ',1pd11.4 /' '/
     5                   ' Stopping because  sqrt(widd) =',
     6                   ' sqrt(taud + musq - rho)  undefined.'/' ' )
                      if (ipr)  write(mou,1000) taud, musqd, rho, widd
                      if (ipr) then
                                   ip10 = -1
                               else
                                   ip10 = 0
                      end if
c****                 call VECPRD ( VTQm1d, n, 'VTQd', ip10, mou )
c****                 call VECPRD ( VTQwsq, n, 'VQw2', ip10, mou )
c****                 call VECPRD ( SigUTb, n, 'SUTb', ip10, mou )
                      call VECPR  ( VTQm1d, n, 'VTQd', ip10, mou )
                      call VECPR  ( VTQwsq, n, 'VQw2', ip10, mou )
                      call VECPR  ( SigUTb, n, 'SUTb', ip10, mou )
c            ******
              stop
c            ******
      end if
c
c**** sum5 = 0.0d0
      sum5 = 0.0
c
      do 50 i=1,n
c
         sum5 = sum5 + VTQwsq(i) / sg2ptn(i)
c
   50 continue
c
c**** if ( sum5 .lt. 0.0d0 ) then
      if ( sum5 .lt. 0.0   ) then
                                 print 1001, sum5, widd
 1001         format(' '/' Stoppimg in Subroutine WDTH' /' '/
     1                   '   sum5 = ', 1pd11.4 /
     2                   '   widd = ', 1pd11.4 /' '/ 
     1                ' Stopping because', 
     3                '  sqrt( widd*sum5 )  is undefined.'/' '/)
                      if (ipr)  write(mou,1001) sum5, widd
                      if (ipr) then
                                   ip10 = -1
                               else
                                   ip10 = 0
                      end if
c****                 call VECPRD ( VTQm1d, n, 'VTQd', ip10, mou )
c****                 call VECPRD ( VTQwsq, n, 'VQw2', ip10, mou )
c****                 call VECPRD ( SigUTb, n, 'SUTb', ip10, mou )
                      call VECPR  ( VTQm1d, n, 'VTQd', ip10, mou )
                      call VECPR  ( VTQwsq, n, 'VQw2', ip10, mou )
                      call VECPR  ( SigUTb, n, 'SUTb', ip10, mou )
c            ******
              stop
c            ******
      end if
c
c
      widd = widd * sum5
c
c**** widd = dsqrt( widd )
      widd = sqrt( widd )
c
c**** wid = sngl( widd )
      wid =       widd 
c
      if ( iprint .gt. 1 ) then
                           print 1002, sum5, wid
 1002 format ( ' '/ 10x, 'In WDTH, sum5 =', 1pd16.8 /
     1              10x, '          wid =', 1pe13.5  )
              write(mou,1002) sum5, wid
      end if
c
      return
c
      end 
c
c
c*******************************************************************************
c
c
      subroutine ZERO ( a,      b,      machep, t,      f,
c
     1                  fb,     iflag,  AQWRK,   
c
     2                  qq,             VT,     maxmn,  sig,
c
     3                  SigUTb, VTQm1d, UTb,    musq,   mr, 
c
     4                  mc,     nr,     n,      nw,     VTQwsq,
c
     5                  sigsq,  elsqb,  bTU2b,  ajnorm, EVwork,  
c
     6                  VTQdsq, sg2ptn, ipr,    iprint, mou   )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      external  EVAL
c
      logical  ipr
c
      real  machep, musq, m
c
c**** double precision  qq(mc),    VTQwsq(mc), sg2ptn(mc),
      dimension         qq(mc),    VTQwsq(mc), sg2ptn(mc),
     1                  sigsq(mc), VTQm1d(mc), SigUTb(mc),
     2                             UTb(maxmn), EVwork(mc)
c
c**** double precision  VTQdsq, bTU2b
      real              VTQdsq, bTU2b
c
      dimension  AQWRK(mr,mc), VT(maxmn,maxmn), sig(maxmn), ajnorm(mc)
c
c++++++++ fb is the function value where the 'zero' f is found +++++++++++++++++
c
      iflag = 0
c
      sa = a
      sb = b
c
      if ( iprint .gt. 1 )  write(mou,970)
  970   format (' '/ 5x, 'In ZERO, Nos. 1 and 2 calls to EVAL' )
c
c   --------------------------------------------------------------
      call EVAL ( sa,     fa,     AQWRK,  
c
     1            qq,             VT,     maxmn,  sig,
c
     2            SigUTb, VTQm1d,         musq,   mr,
c
     3            mc,     nr,     n,      nw,     VTQwsq,
c
     4            sigsq,  elsqa,  bTU2b,  ajnorm, UTb,   EVwork, 
c
     5            VTQdsq, sg2ptn, ipr,    iprint, mou   )
c   ---------------------------------------------------------------
      call EVAL ( sb,     fb,     AQWRK,  
c
     1            qq,             VT,     maxmn,  sig,
c
     2            SigUTb, VTQm1d,         musq,   mr,
c
     3            mc,     nr,     n,      nw,     VTQwsq,
c
     4            sigsq,  elsqb,  bTU2b,  ajnorm, UTb,   EVwork,
c
     5            VTQdsq, sg2ptn, ipr,    iprint, mou   )
c   ----------------------------------------------------------------
c
       if ( iprint .gt. 1 ) then
                            print 990, elsqa, elsqb, fa, fb, sa, sb
  990     format(' '/,5x, 'In ZERO, L_1**2  =', 1pe13.5, 2x,
     1              ' L_2**2  =', 1pe13.5 / ' ' / 14x, 'dL_1/dt =',
     2              1pe13.5, 2x,  ' dL_2/dt =', 1pe13.5 / ' ' / 16x,
     3              'tau1  =', 1pe13.5, 2x, '   tau2  =', 1pe13.5 )
            if ( ipr ) write (mou,990) elsqa, elsqb, fa, fb, sa, sb
      end if
c
c
      if ( sign( 1.0, fa ) * sign( 1.0, fb ) .le. 0.0 )  go to 10
c
c
      iflag = 1
c
      if ( iprint .gt. 1 ) then
                                 print 1000
                                 if (ipr) write(mou,1000)
 1000         format(' '/ 5x, 'No Axis Crossing' )
      end if
c
      f = sb
c                                     --------
      if ( abs( fb ) .le. abs( fa ) )  return
c                                     --------
      fc = fb
      fb = fa
      fa = fc
      f  = sa
c
c  Return with the smallest (in abs value) function value fb and
c               and the value f which yields fb 
c    --------
      return
c    --------
c
c-------------
   10 continue
c-------------
c
      if ( iprint .gt. 1 ) then
                               print 1010
 1010      format ( ' '/ 5x, 'Root Bracketed' )
                               if ( ipr ) write (mou,1010)
      end if
c
c
  15  continue
c
      c = sa
      fc = fa
      e = sb - sa
      d = e
c
   20 continue
c
      if ( abs( fc ) .ge. abs( fb ) )  go to 30
c
      sa = sb
      sb = c
      c = sa
      fa = fb
      fb = fc
      fc = fa
c
   30 continue
c
      tol = ( 2.0 * machep + t ) * abs( sb )
      m = 0.5*( c - sb )
c
      if ( (abs(m) .le. tol) .or. (fb .eq. 0.0) )  go to 140
c 
      if ( (abs(e) .ge. tol) .and. (abs(fa) .gt. abs(fb)) )  go to 40
c
      e = m
      d = e
      go to 100
c
   40 continue
c
      s = fb / fa
c
      if ( sa .ne. 0.0 )  go to 50
c
      p = 2.0 * m * s
      q = 1.0 - s
      go to 60
c
   50 continue
c
      q = fa / fc
      r = fb / fc
      p = s*( 2.0*m*q*(q - r) - (sb - sa)*(r - 1.0) )
      q = (q - 1.0)*(r - 1.0)*(s - 1.0)
c
   60 continue
c
      if ( p .le. 0.0 )  go to 70
c
      q = -q
      go to 80
c
   70 continue
c
      p = -p
c
   80 continue
c
      s = e
      e = d
c
      if ( ( (2.0 * p) .ge. (3.0 * m * q - abs( tol*q )) )
c
     1            .or.   ( p .ge. abs( 0.5 * s * q ) )   )  go to 90
c
      d = p/q
      go to 100
c
   90 continue
c
      e = m
      d = e
c
  100 continue
c
      sa = sb
      fa = fb
c
      if ( abs( d ) .le. tol )  go to 110
c
      sb = sb + d
      go to 130
c
  110 continue
c
      if ( m .le. 0.0 )  go to 120
c
      sb = sb + tol 
      go to 130
c
  120 continue
c
      sb = sb - tol 
c
  130 continue
c
      if ( iprint .gt. 1 )  write(mou,980)
  980   format(' '/ 5x, 'In ZERO, No. 3 call to EVAL' )
c
c   -------------------------------------------------------------
      call EVAL ( sb,     fb,     AQWRK,  
c
     1            qq,             VT,     maxmn,  sig,
c
     2            SigUTb, VTQm1d,         musq,   mr,
c
     3            mc,     nr,     n,      nw,     VTQwsq,
c
     4            sigsq,  elsqb,  bTU2b,  ajnorm, UTb,   EVwork,
c
     5            VTQdsq, sg2ptn, ipr,    iprint, mou   )
c   --------------------------------------------------------------
c
      if ( (fb .gt. 0.0) .and. (fc .gt. 0.0) )  go to 15
c
      if ( (fb .le. 0.0) .and. (fc .le. 0.0) )  go to 15 
c
      go to 20
c
  140 continue
c
      f = sb
c
c
      return
c
      end 
