2

I'm writing a spares grid code and need to combine N 1-dimensional grid points (written in vector form) into the an array of all possible points. For example one can mix two vectors (a,b) with (c,d,e) giving the following points:

(a,c) (a,d) (a,e) (b,c) (b,d) (b,e)

Matlab has a function called combvec:

http://www.mathworks.co.uk/help/nnet/ref/combvec.html

I'm writing this code in FORTRAN however I can't find the underlying algorithm. The code needs to take in N (N>1) vectors (i.e 2,3...N) and each can be a different length. Does anyone know of an algorithm?

Alexander Vogt
  • 17,879
  • 13
  • 52
  • 68
user2350366
  • 491
  • 1
  • 4
  • 20
  • @bdecaf, I guess it's vectors that can be used to create a sparse matrix. Similar to `sparse(vector1, vector2, values, size1, size2)` in MATLAB (I see you're a MATLAB-guy) – Stewie Griffin Jun 23 '14 at 11:30
  • yeah was thinking about that - that's why I erased the comment. – bdecaf Jun 23 '14 at 11:34
  • But to get back to the question what you ask is called [cartesian product](http://en.wikipedia.org/wiki/Cartesian_product). Searching for it leads to many results including SO:[generate a matrix of possible combinations using fortran](http://stackoverflow.com/questions/14572966/generate-a-matrix-of-possible-combinations-using-fortran) – bdecaf Jun 23 '14 at 11:36
  • Would the question [Complete set of combinations combining 3 set](http://stackoverflow.com/q/24259125/2545927) with answers be of help? That question relates to Matlab and I am assuming it is relevant as you have added the tag `matlab`. – kkuilla Jun 23 '14 at 12:05
  • [Here](http://stackoverflow.com/q/21895335/2586922)'s a solution in Matlab. It shouldn't be hard to translate it into Fortran – Luis Mendo Jun 23 '14 at 16:34

2 Answers2

1

I don't know Fortran, but since you say you can't find the underlying algorithm I'm assuming you will be able to write this yourself once you know the algorithm. It's quite easy actually. The pseudocode would be something like this (assuming there are no duplicates):

index = 0   ! or 1
for each element in first vector
    for each element in second vector
        matrix(index,1) = current element of first vector
        matrix(index,2) = current element of second vector
        index = index + 1
    end for
end for

This should give you a matrix similar to the one you would get using combvec. There are probably more efficient ways to do this, but as I don't know the details of Fortran I can't help you there unfortunately. In Matlab you would of course vectorize this.

Good luck =)

Stewie Griffin
  • 14,889
  • 11
  • 39
  • 70
0

The following function should do what you want, I think. It is as simple as possible, taking as input a rank 1 array containing the cardinalities of the data sets and returning a rank 2 array, one column for each set, containing the indices for that set. The expression 1 + mod((i-1)/rep, N) represents the ith element of a sequence of integers 1,2,...,N, with each element repeated rep times.

! (requires explicit interface)
pure function cartprod(v) result(cp)
  integer, intent(in) :: v(1:)
  integer :: cp(product(v), size(v,1))

  integer :: i, j, p, rep

  p = product(v)

  do j = 1, size(v,1) 
    rep = p / product(v(1:j))
    do i = 1, p
      cp(i,j) = 1 + mod((i-1)/rep, v(j))
    enddo
  enddo
end function

Suppose you have defined a dynamic length vector as follows, you could directly obtain the matrix of combinations:

module dynamic_vector
  implicit none

  type :: d_vector
    integer, allocatable :: val(:)
  end type

contains

  pure function combvec(v) result(cv)
    type(d_vector), intent(in) :: v(1:)
    integer, allocatable :: cv(:,:)

    integer :: i, j, prod, rep, len, sizes(size(v,1))

    len = size(v,1)

    ! Determine sizes of the vectors, loop is necessary because we may not
    ! reference v%val if v is an array and val is allocatable.
    do i = 1, len
      sizes(i) = size(v(i)%val,1)
    enddo

    prod = product(sizes)

    ! Allocate and fill the output matrix
    allocate(cv(prod, len))

    do j = 1, len
      rep = prod / product(sizes(1:j))
      do i = 1, prod
        cv(i,j) = v(j)%val(1 + mod((i-1)/rep, sizes(j)))
      enddo
    enddo
  end function
end module

A short test program:

program test
  use dynamic_vector
  implicit none

  type(d_vector) :: foo(2)
  integer :: i, bar(:,:)
  allocatable :: bar

  allocate(foo(1)%val, source = [1,2])
  allocate(foo(2)%val, source = [3,4,5])

  bar = combvec(foo)
  write(*,'(2(I0,X))') (bar(i,:), i = 1, 6)
end program

Result:

1 3
1 4
1 5
2 3
2 4
2 5
sigma
  • 2,758
  • 1
  • 14
  • 18
  • By the way, gfortran 4.8.2 gives me an internal compiler error upon compiling the module and program, while the Intel compiler version 14 is fine. Also, I have not bothered guarding against unallocated arrays in `combvec` ;). – sigma Jun 24 '14 at 21:30
  • Verified that gfortran has an [existing bug](https://gcc.gnu.org/bugzilla/show_bug.cgi?id=45440) with statements like `allocate(x, source = [1,2])`, where the extent of `x` is left out. This does work: `allocate(x(2), source = [1,2])`. – sigma Jun 25 '14 at 10:02