tree_rib_node_bisect Subroutine

    subroutine tree_rib_node_bisect(c_particle, order, div)
        !> Divide given cluster of spheres into two subclusters by inertial bisection
        !!
        !! @param[in] nsph: Number of all input spheres
        !! @param[in] csph: Centers of all input spheres
        !! @param[in] n: Number of spheres in a given cluster
        !! @param[inout] order: Indexes of spheres in a given cluster. On exit, indexes
        !!      `order(1:div)` correspond to the first subcluster and indexes
        !!      `order(div+1:n)` correspond to the second subcluster.
        !! @param[out] div: Break point of `order` array between two clusters.
        !! @param[inout] ddx_error: ddX error
        implicit none

        real(rp), intent(in) :: c_particle(:, :)
        integer(ip), intent(inout) :: order(:)
        integer(ip), intent(out) :: div
        
        real(rp) :: c(3),  s(3)
        real(rp), allocatable :: tmp_c(:, :), work(:)
        external :: dgesvd
        integer(ip) :: i, l, r, lwork, info, n, n_particle
        integer(ip), allocatable :: tmp_order(:)


        n = size(order)
        n_particle = size(c_particle, 2)

        allocate(tmp_c(3, n))
        allocate(tmp_order(n))

        ! Get coordinates in a contiguous array
        do i = 1, n
            tmp_c(:, i) = c_particle(:, order(i))
        end do
        ! Remove the geometrical center for the coordinates
        c = sum(tmp_c, 2) / n
        do i=1, n
            tmp_c(:,i) = tmp_c(:,i) - c
        end do
        
        !! Find right singular vectors
        ! Get proper size of temporary workspace
        lwork = -1
        call dgesvd('N', 'O', 3, n, tmp_c, 3, s, tmp_c, 3, tmp_c, 3, &
                    s, lwork, info)
        lwork = int(s(1))
        allocate(work(lwork))
        ! Get right singular vectors
        call dgesvd('N', 'O', 3, n, tmp_c, 3, s, tmp_c, 3, tmp_c, 3, &
                    work, lwork, info)
        if (info.ne.0) then
            call fmm_error('DGESVD failed in tree_node_bisect')
            return
        end if
        deallocate(work)
        !! Sort spheres by sign of the above scalar product, which is equal to
        !! the leading right singular vector scaled by the leading singular value.
        !! However, we only care about sign, so we take into account only the
        !! leading right singular vector.
        ! First empty index from the left
        l = 1
        ! First empty index from the right
        r = n
        ! Cycle over values of the singular vector
        do i = 1, n
            ! Positive scalar products are moved to the beginning of `order`
            if (tmp_c(1, i) > 0.0) then
                tmp_order(l) = order(i)
                l = l + 1
            ! Negative scalar products are moved to the end of `order`
            else
                tmp_order(r) = order(i)
                r = r - 1
            end if
        end do
        ! Set divider and update order
        div = r
        order = tmp_order
        deallocate(tmp_c)
        deallocate(tmp_order)
    end subroutine tree_rib_node_bisect