subroutine generate_k_way_intersections_from_lists(atom_lists, n_atoms_list, n_sets, k, combination, max_atoms, &
temp_intersections, temp_sets, temp_levels, intersection_count)
!! Generate all k-way intersections from atom lists
integer, intent(in) :: atom_lists(:, :), n_atoms_list(:), n_sets, k, max_atoms
integer, intent(inout) :: combination(:)
integer, intent(inout) :: temp_intersections(:, :), temp_sets(:, :), temp_levels(:)
integer, intent(inout) :: intersection_count
integer, allocatable :: current_intersection(:), temp_intersection(:)
integer :: current_n_intersect, temp_n_intersect
integer :: i, j, set_idx
logical :: has_intersection
call next_combination_init(combination, k)
do
! Compute intersection of all sets in this combination
has_intersection = .false.
! Start with the first set in the combination
set_idx = combination(1)
allocate (current_intersection(n_atoms_list(set_idx)))
current_intersection = atom_lists(1:n_atoms_list(set_idx), set_idx)
current_n_intersect = n_atoms_list(set_idx)
! Intersect with each subsequent set
do i = 2, k
set_idx = combination(i)
allocate (temp_intersection(current_n_intersect))
! Find intersection
temp_n_intersect = 0
do j = 1, current_n_intersect
if (any(atom_lists(1:n_atoms_list(set_idx), set_idx) == current_intersection(j))) then
temp_n_intersect = temp_n_intersect + 1
temp_intersection(temp_n_intersect) = current_intersection(j)
end if
end do
! Update current intersection
deallocate (current_intersection)
if (temp_n_intersect > 0) then
allocate (current_intersection(temp_n_intersect))
current_intersection = temp_intersection(1:temp_n_intersect)
current_n_intersect = temp_n_intersect
has_intersection = .true.
else
has_intersection = .false.
end if
deallocate (temp_intersection)
if (.not. has_intersection) exit
end do
! Store if we found an intersection
if (has_intersection .and. current_n_intersect > 0) then
intersection_count = intersection_count + 1
temp_intersections(1:current_n_intersect, intersection_count) = current_intersection(1:current_n_intersect)
temp_sets(1:k, intersection_count) = combination(1:k)
temp_levels(intersection_count) = k
end if
if (allocated(current_intersection)) deallocate (current_intersection)
! Get next combination
if (.not. next_combination(combination, k, n_sets)) exit
end do
end subroutine generate_k_way_intersections_from_lists