      subroutine su3reunit(am) ! Bernd Berg Mar 04 2008.
c Alexei Bazavov Oct 20 2005. Reunitarization of SU(3) matrix.
      include '../../Libs/Fortran/implicit.sta'
      include '../../Libs/Fortran/constants.par'
      dimension am(1:18)

c normalize first column
      det=am(1)*am(1)+am(10)*am(10)+am(4)*am(4)+am(13)*am(13)+
     &    am(7)*am(7)+am(16)*am(16)
      det=dsqrt(det)
      am(1)=am(1)/det
      am(10)=am(10)/det
      am(4)=am(4)/det
      am(13)=am(13)/det
      am(7)=am(7)/det
      am(16)=am(16)/det

c orthogonalize second column
      scalreal=am(2)*am(1)+am(11)*am(10)+am(5)*am(4)+am(14)*am(13)+
     &         am(8)*am(7)+am(17)*am(16)
      scalimag=am(11)*am(1)-am(2)*am(10)+am(14)*am(4)-am(13)*am(5)+
     &         am(17)*am(7)-am(16)*am(8)
      am12r=am(2)-(scalreal*am(1)-scalimag*am(10))
      am12i=am(11)-(scalreal*am(10)+scalimag*am(1))
      am22r=am(5)-(scalreal*am(4)-scalimag*am(13))
      am22i=am(14)-(scalreal*am(13)+scalimag*am(4))
      am32r=am(8)-(scalreal*am(7)-scalimag*am(16))
      am32i=am(17)-(scalreal*am(16)+scalimag*am(7))

c normalize second column
      det=am12r*am12r+am12i*am12i+am22r*am22r+am22i*am22i+
     &    am32r*am32r+am32i*am32i
      det=dsqrt(det)
      am(2)=am12r/det
      am(11)=am12i/det
      am(5)=am22r/det
      am(14)=am22i/det
      am(8)=am32r/det
      am(17)=am32i/det

c build third column as cross-product of first and second
      am(3)=am(4)*am(8)-am(13)*am(17)-am(5)*am(7)+am(14)*am(16)
      am(12)=-am(4)*am(17)-am(8)*am(13)+am(5)*am(16)+am(7)*am(14)
      am(6)=am(7)*am(2)-am(16)*am(11)-am(1)*am(8)+am(10)*am(17)
      am(15)=-am(7)*am(11)-am(2)*am(16)+am(1)*am(17)+am(8)*am(10)
      am(9)=am(1)*am(5)-am(10)*am(14)-am(2)*am(4)+am(11)*am(13)
      am(18)=-am(1)*am(14)-am(5)*am(10)+am(2)*am(13)+am(4)*am(11)


      return
      end
