orthogonalize_u Subroutine

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

    implicit none

    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
    integer :: LWORK

    if (n .lt. m) call io_error('n<m')
    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 io_error(' ERROR: IN ZGESVD IN orthogonalize_u')
    endif
    deallocate (smat, eig, WORK, RWORK)
    ! 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)
      enddo
    enddo
    enddo
    deallocate (evecl, evecr)

    return
  end subroutine orthogonalize_u