Store all r combinations of a list-CodePudding

I need a general program in fortran to obtain all possible combinations of r elements in a list of n elements. I have found this code that prints all the combinations (r=3, n =5) but I need them to be stored in an array.

I tried to record them as rows near the write statement but it does not work. Turning the recursive subprogram into a recursive function is also not working.

program combinations
  implicit none
  
  integer, parameter :: m_max = 3
  integer, parameter :: n_max = 5
  integer, dimension (m_max) :: comb
  character (*), parameter :: fmt = '(i0' // repeat (', 1x, i0', m_max - 1) // ')'
  
  call gen (1)
contains
  recursive subroutine gen (m)
    implicit none
    
    integer, intent (in) :: m
    integer :: n
    
    if (m > m_max) then
      write (*, fmt) comb
    else
      do n = 1, n_max
        if ((m == 1) .or. (n > comb (m - 1))) then
          comb (m) = n
          call gen (m   1)
        end if
      end do
    end if     
  end subroutine gen
end program combinations

CodePudding user response：

Firstly, mixing global variables and recursive procedures is a good way to cause a lot of unnecessary confusion and debugging, so let's turn comb and n_max into procedure arguments, use size(comb) to give m_max, and for now replace fmt with *:

program combinations
  implicit none

  integer :: comb(3)

  call gen(comb, 1, 5)
contains
  recursive subroutine gen(comb, m, n_max)
    integer, intent(inout) :: comb(:)
    integer, intent(in) :: m
    integer, intent(in) :: n_max

    integer :: n

    if (m > size(comb)) then
      write (*, *) comb
    else
      do n = 1, n_max
        if ((m == 1) .or. (n > comb(m - 1))) then
          comb(m) = n
          call gen(comb, m 1, n_max)
        end if
      end do
    end if
  end subroutine gen
end program combinations

The next thing to note is there's a subtle bug in your code. The line

if ((m == 1) .or. (n > comb (m - 1))) then

isn't guaranteed to work if m=1. Fortran does not guarantee short-circuiting of logical operators, so even if (m == 1) evaluates to .true., the (n > comb (m - 1)) could be evaluated, causing a segfault. Let's get around this by introducing a variable n_min, and calculating it correctly:

  recursive subroutine gen(comb, m, n_max)
    integer, intent(inout) :: comb(:)
    integer, intent(in) :: m
    integer, intent(in) :: n_max

    integer :: n
    integer :: n_min

    if (m > size(comb)) then
      write (*, *) comb
    else
      if (m == 1) then
        n_min = 1
      else
        n_min = comb(m-1)   1
      endif

      do n = n_min, n_max
        comb(m) = n
        call gen (comb, m 1, n_max)
      end do
    end if
  end subroutine gen

Okay, now we can start thinking about returning the combinations from gen. To do this, let's change gen from a subroutine into a function, and have it return a 2-D array. We're going to need to append one 2-D array onto another, so let's write a function to do that now:

  function append_combinations(input, new_combinations) result(output)
    integer, intent(in) :: input(:,:)
    integer, intent(in) :: new_combinations(:,:)
    integer, allocatable :: output(:,:)

    allocate(output(size(input,1), size(input,2) size(new_combinations,2)))
    output(:, :size(input,2)) = input
    output(:, size(input,2) 1:) = new_combinations
  end function

and now the whole program looks like

program combinations
  implicit none

  integer :: comb(3)
  integer, allocatable :: combs(:,:)
  integer :: i

  combs = gen(comb, 1, 5)

  write(*, *) ""
  do i=1,size(combs,2)
    write(*, *) combs(:,i)
  enddo
contains
  recursive function gen(comb, m, n_max) result(combs)
    integer, intent(inout) :: comb(:)
    integer, intent(in) :: m
    integer, intent(in) :: n_max
    integer, allocatable :: combs(:,:)

    integer :: n
    integer :: n_min
    integer, allocatable :: new_combs(:,:)

    if (m > size(comb)) then
      write (*, *) comb
      combs = reshape(comb, [size(comb),1])
    else
      if (m == 1) then
        n_min = 1
      else
        n_min = comb(m-1)   1
      endif

      allocate(combs(size(comb), 0))
      do n = n_min, n_max
        comb(m) = n
        new_combs = gen(comb, m 1, n_max)
        combs = append_combinations(combs, new_combs)
      end do
    end if
  end function gen

  function append_combinations(input, new_combinations) result(output)
    integer, intent(in) :: input(:,:)
    integer, intent(in) :: new_combinations(:,:)
    integer, allocatable :: output(:,:)

    allocate(output(size(input,1), size(input,2) size(new_combinations,2)))
    output(:, :size(input,2)) = input
    output(:, size(input,2) 1:) = new_combinations
  end function
end program combinations