      subroutine su3_update_su2(idl,isl,st,i,j,itry)
C Alexei Bazavov / Bernd Berg Mar 6 2009.
C Cabbibo-Marinari update of SU(2) subgroup in SU(3)
c idl,isl:   links.
c st:        staple.
c i,j:       col/raw to update (i<j).
c itry:      number of RUA steps.
      include '../../Libs/Fortran/implicit.sta'
      include '../../Libs/Fortran/constants.par'
      parameter(eps=one/ten**16) ! less than eps means =0
      include '../../Libs/Fortran_par/lat.par'
      include '../../Libs/Fortran_par/mc.par'
      include '../../Libs/Fortran_par/common_lat.f'
      include '../../Libs/SU3_par/common_su3.f'

      dimension st(1:18),rp(1:18),ru1(1:18),rr(1:18)
      dimension rl(1:18)

c Copy link into matrix
      call su3copy_m_m(alsu3(1,idl,isl),rl)
c Multiply link by staple
      call su3mult_m_m_m(rl,st,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 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 a0 with FHKP heatbath
      alpha=(two*beta*det)/three ! effective beta is 2/3 * beta
      a0tr=su2_a0(alpha,itry)

c Generate vector on 3-D sphere
      call ranmar(xr)
      cos_theta=two*xr-one
      call ranmar(xr)
      phi=2*PI*xr
      sin_theta=sqrt(1-cos_theta*cos_theta)
      aaa=sqrt(abs(one-a0tr*a0tr)) ! Avoids -1.D00 on 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)=a0tr
      rp(i0+9)=a3
      rp(i1)=a2
      rp(i1+9)=a1
      rp(i2)=-a2
      rp(i2+9)=a1
      rp(i3)=a0tr
      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(1,idl,isl))

      return
      end
