      subroutine su3_b2nhitupdt(nhit,i,j,iacpt,alsu3,bstaple)
C Bernd Berg Mar 6 2009; Alexei Bazavov Oct 20 2005.
C One update attempt. Accepted or rejected.
C Uses bstaple=beta*staple (beta differs with plaquettes).
C Cabbibo-Marinari update of SU(2) subgroup in SU(3)
c i,j:    col/raw to update (i<j)
c iacpt:  number of accepted proposals.
      include '../../Libs/Fortran/implicit.sta'
      include '../../Libs/Fortran/constants.par'
      parameter(eps=one/ten**16) ! less than eps means =0
      dimension alsu3(18),bstaple(18),rp(18),ru1(18),rr(18),rl(18)
      ltest=.true.
      ltest=.false.

c Copy link into matrix
      call su3copy_m_m(alsu3,rl)
c Multiply link by bstaple
      call su3mult_m_m_m(rl,bstaple,rp) 
c Extract i,j part and project onto SU(2)
      i0=4*i-3
      i1=3*(i-1)+j
      i2=3*(j-1)+i
      i3=4*j-3
      r0=(rp(i0)+rp(i3))/two
      r1=(rp(i1+9)+rp(i2+9))/two
      r2=(rp(i1)-rp(i2))/two
      r3=(rp(i0+9)-rp(i3+9))/two

c Calculate determinant
      det=r0*r0+r1*r1+r2*r2+r3*r3
      det=sqrt(det)

c Generate a0 with FHKP heatbath
      alpha=(two/three)*det ! effective beta is 2/3 * beta
      call su2_a0nhit(nhit,alpha,a0,lacpt)
      if(.not.lacpt) return ! Reject a0.

      iacpt=iacpt+1 ! Accept.
c Normalize staple and take inverse (herm. conj. matrix)
      r0=r0/det
      r1=r1/det
      r2=r2/det
      r3=r3/det
c Create SU(3) matrix out of this subgroup
      call su3init1(ru1)
      ru1(i0)=r0
      ru1(i0+9)=-r3
      ru1(i1)=-r2
      ru1(i1+9)=-r1
      ru1(i2)=r2
      ru1(i2+9)=-r1
      ru1(i3)=r0
      ru1(i3+9)=r3

c Generate vector on 3-D sphere
      call ranmar(xr)
      cos_theta=two*xr-one
      call ranmar(xr)
      phi=tpi*xr
      sin_theta=sqrt(one-cos_theta*cos_theta)
      aaa=sqrt(abs(one-a0*a0)) ! Avoids -1.D-300 or so in sqrt.
      a1=aaa*sin_theta*cos(phi)
      a2=aaa*sin_theta*sin(phi)
      a3=aaa*cos_theta

c Create (a0,a) as SU(3) matrix rp
      call su3init1(rp)
      rp(i0)=a0
      rp(i0+9)=a3
      rp(i1)=a2
      rp(i1+9)=a1
      rp(i2)=-a2
      rp(i2+9)=a1
      rp(i3)=a0
      rp(i3+9)=-a3
c Multiply rp matrix by ru1 matrix on the right
      call su3mult_m_m_m(rp,ru1,rr)
c Update link by multiplying by resulting matrix (rr) on the left
      call su3mult_m_m_m(rr,rl,alsu3)

      if(ltest) stop "su3_b2update: ltest."
      return
      end
