      subroutine su3_getstaple(idl,isl,res_matr) 
c Alexei Bazavov / Bernd Berg Mar 04 2008.
c Calculate sum of staples from plqs containing link idl.
c idl, isl -- link coords.
c res_matr -- resulting matrix.
      include '../../Libs/Fortran/implicit.sta'
      include '../../Libs/Fortran/constants.par'
      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 res_matr(1:18),r1(1:18),r2(1:18)

      call su3init0(res_matr) ! Initialize matrix to 0.

c Loop on smaller dimensions.
      do idp=1,idl-1
c       building staple on forward plq
        call su3mult_m_mh_m(alsu3(1,idp,ipf(idl,isl)),
     &                      alsu3(1,idl,ipf(idp,isl)),r1)
        call su3mult_m_mh_m(r1,alsu3(1,idp,isl),r2)
c       add to resulting staple
        call su3add_m_m(res_matr,r2)
c       building staple on backward plq
        ii=ipb(idp,isl)
        call su3mult_mh_mh_m(alsu3(1,idp,ipf(idl,ii)),
     &                       alsu3(1,idl,ii),r1)
        call su3mult_m_m_m(r1,alsu3(1,idp,ii),r2)
c       add to resulting staple
        call su3add_m_m(res_matr,r2)
      end do

c Loop on bigger dimensions (coded to avoid if(idp.eq.idl) operation).
      do idp=idl+1,nd
c       building staple on forward plq
        call su3mult_m_mh_m(alsu3(1,idp,ipf(idl,isl)),
     &                      alsu3(1,idl,ipf(idp,isl)),r1)
        call su3mult_m_mh_m(r1,alsu3(1,idp,isl),r2)
c       add to resulting staple
        call su3add_m_m(res_matr,r2)
c       building staple on backward plq
        ii=ipb(idp,isl)
        call su3mult_mh_mh_m(alsu3(1,idp,ipf(idl,ii)),
     &                       alsu3(1,idl,ii),r1)
        call su3mult_m_m_m(r1,alsu3(1,idp,ii),r2)
c       add to resulting staple
        call su3add_m_m(res_matr,r2)
      end do

      return
      end
