orthogonalize_u Subroutine

private subroutine orthogonalize_u(ndim, m, u, n, error, comm)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: ndim
integer, intent(in) :: m
complex(kind=dp), intent(inout) :: u(ndim,m)
integer, intent(in) :: n
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm

Calls

proc~~orthogonalize_u~~CallsGraph proc~orthogonalize_u orthogonalize_u proc~set_error_dealloc set_error_dealloc proc~orthogonalize_u->proc~set_error_dealloc proc~set_error_fatal set_error_fatal proc~orthogonalize_u->proc~set_error_fatal zgesvd zgesvd proc~orthogonalize_u->zgesvd proc~comms_sync_error comms_sync_error proc~set_error_dealloc->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_dealloc->proc~set_base_error proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error

Called by

proc~~orthogonalize_u~~CalledByGraph proc~orthogonalize_u orthogonalize_u proc~symmetrize_ukirr symmetrize_ukirr proc~symmetrize_ukirr->proc~orthogonalize_u proc~sitesym_dis_extract_symmetry sitesym_dis_extract_symmetry proc~sitesym_dis_extract_symmetry->proc~symmetrize_ukirr proc~sitesym_symmetrize_u_matrix sitesym_symmetrize_u_matrix proc~sitesym_symmetrize_u_matrix->proc~symmetrize_ukirr proc~dis_extract dis_extract proc~dis_extract->proc~sitesym_dis_extract_symmetry proc~dis_extract->proc~sitesym_symmetrize_u_matrix proc~dis_main dis_main proc~dis_main->proc~sitesym_symmetrize_u_matrix proc~dis_main->proc~dis_extract proc~internal_find_u internal_find_u proc~dis_main->proc~internal_find_u proc~internal_find_u->proc~sitesym_symmetrize_u_matrix proc~overlap_project overlap_project proc~overlap_project->proc~sitesym_symmetrize_u_matrix proc~w90_disentangle~2 w90_disentangle proc~w90_disentangle~2->proc~dis_main proc~w90_project_overlap~2 w90_project_overlap proc~w90_project_overlap~2->proc~overlap_project proc~w90_disentangle w90_disentangle proc~w90_disentangle->proc~w90_disentangle~2 proc~w90_project_overlap w90_project_overlap proc~w90_project_overlap->proc~w90_project_overlap~2 program~wannier wannier program~wannier->proc~w90_disentangle~2 program~wannier->proc~w90_project_overlap~2

Source Code

  subroutine orthogonalize_u(ndim, m, u, n, error, comm)
    !================================================!

    implicit none

    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm
    integer, intent(in) :: ndim, m
    complex(kind=dp), intent(inout) :: u(ndim, m)
    integer, intent(in) :: n

    complex(kind=dp), allocatable :: smat(:, :), evecl(:, :), evecr(:, :)
    complex(kind=dp), allocatable :: WORK(:)
    real(kind=dp), allocatable :: eig(:), RWORK(:)
    integer :: INFO, i, j, l, LWORK, ierr

    if (n .lt. m) then
      call set_error_fatal(error, 'n<m', comm)
      return
    end if
    allocate (smat(n, m)); smat(1:n, 1:m) = u(1:n, 1:m)
    allocate (evecl(n, n), evecr(m, m))
    allocate (eig(min(m, n)))

    allocate (RWORK(5*min(n, m)))
    LWORK = 2*min(m, n) + max(m, n)
    allocate (WORK(LWORK))

    ! Singular-value decomposition
    call zgesvd('A', 'A', n, m, smat, n, eig, evecl, n, evecr, m, WORK, LWORK, RWORK, INFO)
    if (info .ne. 0) then
      call set_error_fatal(error, ' ERROR: IN ZGESVD IN orthogonalize_u', comm)
      return
    end if
    deallocate (smat, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating smat in orthogonalize_u', comm)
      return
    end if
    deallocate (eig, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating eig in orthogonalize_u', comm)
      return
    end if
    deallocate (WORK, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating WORK in orthogonalize_u', comm)
      return
    end if
    deallocate (RWORK, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating RWORK in orthogonalize_u', comm)
      return
    end if

    ! u_matrix is the initial guess for the unitary rotation of the
    ! basis states given by the subroutine extract
    u = 0
    do j = 1, m
    do l = 1, m
      do i = 1, n
        u(i, j) = u(i, j) + evecl(i, l)*evecr(l, j)
      end do
    end do
    end do

    deallocate (evecl, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating evecl in orthogonalize_u', comm)
      return
    end if
    deallocate (evecr, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating evecr in orthogonalize_u', comm)
      return
    end if
  end subroutine orthogonalize_u