c--------------------------------------------------c
c     qmcbv.f computes vector of critical values   c
c       for constrained mcb inference              c
c       with 1 or 2 factors                        c
c     last revision: 12/10/90 by jch               c
c--------------------------------------------------c
      subroutine qmcbv(km1,fact1,cx,nu,nalpha,cc,d)
      implicit doubleprecision (a-h,o-z)
c6    implicit real*8 (a-h,o-z)
      external glv,ghb1,gh2d
      dimension cc(nalpha),d(nalpha)
      dimension fact1(km1),fact2(100)
      dimension a(100),b(100),c(100)
      dimension ph(48),wh(48),whln(48),zh(48)
      dimension pl(48),wl(48)
      dimension ah(7),bh(12)
      common/list2/a,b,c
      common/list5/dinfnu
      common/list6/ph,wh,zh,whln
      common/list7/pl,wl
      common/list8/zero,zsr,zlog,zln
      common/list9/sr2,pi,srpi,pib2
      common/list10/iin,iout,iouti
      common/list12/ah,bh
c     ----------------------------------
c     initialize commonly used constants
c     ----------------------------------
      call commgq
c     ----------------------------------------
c     initialize default tolerence (eps) on cc
c     ----------------------------------------
      eps=0.20d-04
c     --------------------------------------------------------------
c     initialize default tolerence (dstop) on half width of interval
c     --------------------------------------------------------------
      dstop=0.5d-03
c     -------------------------------------------------------
c     initialize default bound on number of iteration (itmax)
c     -------------------------------------------------------
      itmax=9
      imc=21
      ifact=1
      k=km1+1
      dnu=dble(float(nu))
c
c     -----------------------
c     initialize correlations
c     -----------------------
      do 10 j=1,km1
c       ----------------------
c       a = factor pattern
c       c = 1/sqrt(uniqueness)
c       ----------------------
        a(j)=fact1(j)
        if (ifact.eq.2) goto 5
        c(j)=1.0d0/dsqrt(1.0d0-a(j)*a(j))
        goto 10
    5   b(j)=fact2(j)
        c(j)=1.0d0/dsqrt(1.0d0-a(j)*a(j)-b(j)*b(j))
   10 continue
c
c     -------------------
c     find critical value
c     -------------------
      do 1000,i=1,nalpha
c         -----------------------------------------------
c         obtain upper and lower bounds on critical value
c         -----------------------------------------------
          call bdb11(km1,cc(i),nu,dnu,dhigh,dlow)
c         call bdmcb1(k,cc(i),nu,dnu,dhigh,dlow)
c
c         ------------------------------------------------------------
c         compute critical value d by the modified regula falsi method
c         ------------------------------------------------------------
 15       if (ifact.eq.2) goto 20 
c
          call mrgfls(ghb1,imc,k,glv,cc(i),dhigh,dlow,w,dnu,-1.0d0,
     *                dstop,eps,itmax,iflag,iter)
          if(iflag.lt.0) stop
          d(i)=w
c
          goto 30
   20     glhigh=glv(gh2d,imc,k,dhigh,dnu,-1.0d0)
          gllow=glv(gh2d,imc,k,dlow,dnu,-1.0d0)
c
          call mrgfls(gh2d,imc,k,glv,cc(i),dhigh,dlow,w,dnu,-1.0d0,
     *                dstop,eps,itmax,iflag,iter)
          if(iflag.lt.0) stop
          d(i)=w
c
   30 continue
 1000 continue
      return
      end
