symmetrize_ukirr Subroutine

  subroutine symmetrize_ukirr(num_wann, num_bands, ir, ndim, umat, &
                              sitesym, stdout, error, comm, n)
    !================================================!
    !
    !  calculate u~(k)=1/N_{R'} \sum_{R'} d^{+}(R',k) u(k) D(R',k)
    !  where R'k = k
    !  and orthonormalize it
    !
    !================================================!

    use w90_wannier90_types, only: sitesym_type
    use w90_error, only: w90_error_type, set_error_fatal, set_error_unconv
    implicit none

    ! arguments
    type(sitesym_type), intent(in) :: sitesym
    type(w90_comm_type), intent(in) :: comm
    type(w90_error_type), allocatable, intent(out) :: error

    integer, intent(in) :: ir, ndim
    integer, intent(in) :: num_bands
    integer, intent(in) :: num_wann
    integer, intent(in) :: stdout
    integer, optional, intent(in) :: n

    complex(kind=dp), intent(inout) :: umat(ndim, num_wann)

    ! local variables
    integer :: isym, ngk, i, iter, ntmp
    integer, parameter :: niter = 100
    real(kind=dp)    :: diff
    complex(kind=dp) :: usum(ndim, num_wann)
    complex(kind=dp) :: cmat_sub(ndim, num_wann)
    complex(kind=dp) :: cmat(ndim, num_wann)
    complex(kind=dp) :: cmat2(num_wann, num_wann)

    !write(stdout,"(a)") '-- symmetrize_ukirr --'
    if (present(n)) then
      if (ndim .ne. num_bands) then
        call set_error_fatal(error, 'ndim!=num_bands', comm)
        return
      end if
      ntmp = n
    else
      if (ndim .ne. num_wann) then
        call set_error_fatal(error, 'ndim!=num_wann', comm)
        return
      end if
      ntmp = ndim
    end if

    ngk = count(sitesym%kptsym(:, ir) .eq. sitesym%ir2ik(ir))
    if (ngk .eq. 1) then
      call orthogonalize_u(ndim, num_wann, umat, ntmp, error, comm)
      return
    end if

    do iter = 1, niter
      usum(:, :) = 0
      cmat2(:, :) = 0
      do i = 1, num_wann
        cmat2(i, i) = cmat2(i, i) + ngk
      end do
      do isym = 1, sitesym%nsymmetry
        if (sitesym%kptsym(isym, ir) .ne. sitesym%ir2ik(ir)) cycle
        !
        ! cmat = d^{+}(R,k) U(k) D(R,k)
        ! size of umat: umat(ndim,num_wann)
        !
        ! cmat_sub = U(k) D(R,k)
        call zgemm('N', 'N', ntmp, num_wann, num_wann, cmplx_1, umat, ndim, &
                   sitesym%d_matrix_wann(:, :, isym, ir), num_wann, cmplx_0, cmat_sub, ndim)
        ! cmat = d^{+}(R,k) * cmat_sub
        call zgemm('C', 'N', ntmp, num_wann, ntmp, cmplx_1, sitesym%d_matrix_band(:, :, isym, ir), &
                   ndim, cmat_sub, ndim, cmplx_0, cmat, ndim)

        usum(:, :) = usum(:, :) + cmat(:, :)
        ! check
        cmat2(:, :) = cmat2(:, :) - &
                      matmul(conjg(transpose(umat(:ntmp, :))), cmat(:ntmp, :))
      end do ! isym
      diff = sum(abs(cmat2))
      if (diff .lt. sitesym%symmetrize_eps) exit
      if (iter .eq. niter) then
        write (stdout, "(a)") 'Error in symmetrize_u: not converged'
        write (stdout, "(a)") 'Either eps is too small or specified irreps is not'
        write (stdout, "(a)") '  compatible with the bands'
        write (stdout, "(a,2e20.10)") 'diff,eps=', diff, sitesym%symmetrize_eps
        call set_error_unconv(error, 'symmetrize_ukirr: not converged', comm)
        return
      end if
      usum = usum/ngk
      call orthogonalize_u(ndim, num_wann, usum, ntmp, error, comm)
      if (allocated(error)) return

      umat(:, :) = usum
    end do ! iter

    return
  end subroutine symmetrize_ukirr