canonicalize Subroutine

    pure subroutine canonicalize(V, lambda)
        !! This routine canonicalizes the eigenvector matrix V and the
        !! corresponding eigenvalues lambda
        !!
        !! This procedure is particularly useful for tests,
        !! where degenerate Eigenspaces might result in different
        !! Eigenvectors.
        real(kind=dp), intent(inout) :: V(:, :), lambda(:)
        real(kind=dp), allocatable :: norms_projections(:)
            !! The norms of the projections
        integer :: low, i, j, d
        integer, allocatable :: idx(:), dimensions(:)
        real(kind=dp), allocatable :: projections(:, :)
        debug_function_name('canonicalize')

        ASSERT(size(V, 1) == size(V, 2) .and. size(V, 1) == size(lambda))

        allocate(norms_projections(size(V, 2)))

        idx = [(i, i = 1, size(lambda, 1))]
        call sort(lambda, idx)
        V(:, :) = V(:, idx)

        call determine_eigenspaces(lambda, dimensions)

        low = 1
        do j = 1, size(dimensions)
            d = dimensions(j)
            if (d == 1) then
                V(:, low) = V(:, low) / norm(V(:, low))
                low = low + d
                cycle
            end if
            ! Ensure that the basis of each subspace is orthonormal
            V(:, low : low + d - 1) = gram_schmidt(V(:, low : low + d - 1))

            projections = project_canonical_unit_vectors(V(:, low : low + d - 1))

            do i = 1, size(norms_projections)
                norms_projections(i) = norm(projections(:, i))
            end do

            idx = [(i, i = 1, size(norms_projections))]
    

        ! merge_sort
        block
            integer, dimension(:), allocatable :: tmp
            integer :: current_size, left, mid, right
            integer, parameter :: along = 1, &
                bubblesort_size = 20

            associate(X => idx)
                ! Determine good number of starting splits
                block
                    integer :: n_splits
                    n_splits = 1
                    do while (size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0) > bubblesort_size)
                        n_splits = n_splits + 1
                    end do
                    current_size = size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0)
                end block

                ! Reuse this variable as tmp for swap in bubble_sort
                ! and for merge in merge_sort.
                block
                    integer :: max_buffer_size, n_merges
                    n_merges = ceiling(log(real(size(X, along)) / real(current_size)) / log(2.0))
                    max_buffer_size = current_size * merge(2**(n_merges - 1), 1, n_merges >= 1)
                    allocate(tmp(max_buffer_size))
                end block

                ! Bubble sort bins of size `bubblesort_size`.
                do left = lbound(X, along), ubound(X, along), current_size
                    right = min(left + bubblesort_size - 1, ubound(X, along))
                        ! bubblesort
    block
      integer :: n, i
      associate(V => X(left : right))
        do n = ubound(V, 1), lbound(V, 1) + 1, -1
          do i = lbound(V, 1), ubound(V, 1) - 1
            if (.not.   comp_idx(V(i), V(i + 1))) then
                  ! swap
    block
        associate(tmp => tmp(1))
            tmp = V(i)
            V(i) = V(i + 1)
            V(i + 1) = tmp
        end associate
    end block
            end if
          end do
        end do
      end associate
    end block
                end do

                do while (current_size < size(X, along))
                    do left = lbound(X, along), ubound(X, along), 2 * current_size
                        right = min(left + 2 * current_size - 1, ubound(X, along))
                        mid = min(left + current_size, right) - 1
                        tmp(: mid - left + 1) = X(left : mid)
                            ! merge
    block
        integer :: i, j, k

        associate(A => tmp(: mid - left + 1), B => X(mid + 1 : right), C => X(left : right))

            if (size(A) + size(B) > size(C)) then
                error stop
            end if

            i = lbound(A, 1)
            j = lbound(B, 1)
            do k = lbound(C, 1), ubound(C, 1)
                if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then
                    if (  comp_idx(A(i), B(j))) then
                        C(k) = A(i)
                        i = i + 1
                    else
                        C(k) = B(j)
                        j = j + 1
                    end if
                else if (i <= ubound(A, 1)) then
                    C(k) = A(i)
                    i = i + 1
                else if (j <= ubound(B, 1)) then
                    C(k) = B(j)
                    j = j + 1
                end if
            end do
        end associate
    end block
                    end do
                    current_size = 2 * current_size
                end do
            end associate
        end block
    

        ! merge_sort
        block
            integer, dimension(:), allocatable :: tmp
            integer :: current_size, left, mid, right
            integer, parameter :: along = 1, &
                bubblesort_size = 20

            associate(X => idx(: d))
                ! Determine good number of starting splits
                block
                    integer :: n_splits
                    n_splits = 1
                    do while (size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0) > bubblesort_size)
                        n_splits = n_splits + 1
                    end do
                    current_size = size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0)
                end block

                ! Reuse this variable as tmp for swap in bubble_sort
                ! and for merge in merge_sort.
                block
                    integer :: max_buffer_size, n_merges
                    n_merges = ceiling(log(real(size(X, along)) / real(current_size)) / log(2.0))
                    max_buffer_size = current_size * merge(2**(n_merges - 1), 1, n_merges >= 1)
                    allocate(tmp(max_buffer_size))
                end block

                ! Bubble sort bins of size `bubblesort_size`.
                do left = lbound(X, along), ubound(X, along), current_size
                    right = min(left + bubblesort_size - 1, ubound(X, along))
                        ! bubblesort
    block
      integer :: n, i
      associate(V => X(left : right))
        do n = ubound(V, 1), lbound(V, 1) + 1, -1
          do i = lbound(V, 1), ubound(V, 1) - 1
            if (.not.   V(i) <= V(i + 1)) then
                  ! swap
    block
        associate(tmp => tmp(1))
            tmp = V(i)
            V(i) = V(i + 1)
            V(i + 1) = tmp
        end associate
    end block
            end if
          end do
        end do
      end associate
    end block
                end do

                do while (current_size < size(X, along))
                    do left = lbound(X, along), ubound(X, along), 2 * current_size
                        right = min(left + 2 * current_size - 1, ubound(X, along))
                        mid = min(left + current_size, right) - 1
                        tmp(: mid - left + 1) = X(left : mid)
                            ! merge
    block
        integer :: i, j, k

        associate(A => tmp(: mid - left + 1), B => X(mid + 1 : right), C => X(left : right))

            if (size(A) + size(B) > size(C)) then
                error stop
            end if

            i = lbound(A, 1)
            j = lbound(B, 1)
            do k = lbound(C, 1), ubound(C, 1)
                if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then
                    if (  A(i) <= B(j)) then
                        C(k) = A(i)
                        i = i + 1
                    else
                        C(k) = B(j)
                        j = j + 1
                    end if
                else if (i <= ubound(A, 1)) then
                    C(k) = A(i)
                    i = i + 1
                else if (j <= ubound(B, 1)) then
                    C(k) = B(j)
                    j = j + 1
                end if
            end do
        end associate
    end block
                    end do
                    current_size = 2 * current_size
                end do
            end associate
        end block

            ! reorthogonalize
            V(:, low : low + d - 1) = gram_schmidt(projections(:, idx(: d)))

            where (near_zero(V(:, low : low + d - 1)))
                V(:, low : low + d - 1) = 0._dp
            end where
            low = low + d
        end do

        contains

            pure logical function comp_idx(i, j)
                integer, intent(in) :: i, j
                comp_idx = norms_projections(i) >= norms_projections(j)
            end function
    end subroutine canonicalize