!--------------------------------------------------------------------------------
  subroutine rotate_epmat ( cz1, cz2, xq, iq)
  !--------------------------------------------------------------------------------
  !
  ! 1) rotate the electron-phonon matrix from the cartesian representation
  !    of the first qpoint of the star to the eigenmode representation (using cz1).
  ! 
  ! 2) rotate the electron-phonon matrix from the eigenmode representation
  !    to the cartesian representation of the qpoint iq (with cz2).
  !
  ! Feliciano Giustino, UCB Feb 2007
  !
  !--------------------------------------------------------------------------
  !
#include "f_defs.h"
  USE kinds, only : DP
  use io_global, only : stdout
  use el_phon, only : epmatq
  use phcom, only : nmodes
  use control_flags, only : iverbosity
  use pwcom, only : nbnd, nks, nr1, nr2, nr3, at, bg
  USE ions_base, ONLY : amass, tau, nat, ntyp => nsp, ityp
  implicit none
  !
  integer :: iq_first, iq
  !  originating q point of the star
  !  teh current qpoint
  !  the sym op leading from the first q in the star to this q point
  real(kind=DP) :: xq(3)
  !  the rotated q vector
  integer :: s(3,3,48), g0vec(3)
  !  the symmetry operation for the eigenmodes
  !  the rotated of each atom
  !  the folding G-vector
  !  the index of the inverse operation
  complex(kind=DP) :: cz1( nmodes, nmodes), cz2(nmodes, nmodes)
  !  the eigenvectors for the first q in the star
  !  the rotated eigenvectors, for the current q in the star
  !
  complex(kind=DP), parameter :: cone = (1.d0, 0.d0), czero = (0.d0, 0.d0)
  real(kind=DP), parameter :: zero = 0.d0, twopi = 6.28318530717959, &
      rydcm1 = 13.6058d0 * 8065.5d0
  !
  ! work variables 
  !
  complex(kind=DP) :: eptmp( nmodes)
  integer :: imode, jmode, mu, nu, ism1, na, ik, ibnd, jbnd
  real(kind=DP) :: massfac
  complex(kind=DP) :: cz_tmp(nmodes,nmodes)

  !
  ! the mass factors: 
  !  1/sqrt(M) for the  direct transform
  !  sqrt(M)   for the inverse transform 
  !
  ! if we set cz1 = cz2 here and we calculate below
  ! cz1 * cz2 we find the identity
  !
  do mu = 1, nmodes
    na = (mu - 1) / 3 + 1
    massfac = sqrt(amass(ityp(na)))
    cz1 (mu, :) = cz1 (mu, :) / massfac
    cz2 (mu, :) = cz2 (mu, :) * massfac
  enddo
  !
  ! the inverse transform also requires the hermitian conjugate
  !
!
! unbelievable: xlf90 does not like the assignment:
!  cz2 = conjg ( transpose ( cz2 ) )
!
! I need to use an auxiliary array! On pgi everything was fine
! Should we grep "transpose" trhouhout the code ?
!
  cz_tmp = conjg ( transpose ( cz2 ) )
  cz2 = cz_tmp
  ! 
  !  ep_mode (j) = cfac * sum_i ep_cart(i) * u(i,j)
  !
  do ibnd = 1, nbnd
   do jbnd = 1, nbnd
    do ik = 1, nks
       !
       !  bring e-p matrix from the cartesian representation of the
       !  first q in the star to the corresponding eigenmode representation
       !
!
! the following does not work, apparently one cannot overwrite with zgemm...
! (epmatq is left unchanged at the exit)
!
!       call zgemv ('t', nmodes, nmodes, cone, cz1, nmodes,  &
!          epmatq (ibnd, jbnd, ik, :, iq), 1, czero, epmatq (ibnd, jbnd, ik, :, iq), 1 )
!
       call zgemv ('t', nmodes, nmodes, cone, cz1, nmodes,  &
          epmatq (ibnd, jbnd, ik, :, iq), 1, czero, eptmp, 1 )
       !
       ! rotate epmat in the cartesian representation for this q in the star
       !
       call zgemv ('t', nmodes, nmodes, cone, cz2, nmodes, &
          eptmp, 1, czero, epmatq (ibnd, jbnd, ik, :, iq), 1 )
       !
    enddo
   enddo
  enddo
  !
  return
  end subroutine rotate_epmat
  !