shuffle.f90 2.06 KB
Newer Older
1

2 3
MODULE shuffle
  use globals
ulrich_y's avatar
ulrich_y committed
4
  use ieps
5 6 7 8 9 10
  implicit none

CONTAINS
  
  FUNCTION append_to_each_row(a, m) result(res)
    ! appends element a to each row of m
11
    type(inum) :: a, m(:,:)
Luca's avatar
Luca committed
12
    integer :: i
13
    type(inum) :: res(size(m,1),size(m,2)+1)
14 15 16 17
    do i=1,size(m,1)
      res(i,:) = [a,m(i,:)]
    end do
  END FUNCTION append_to_each_row
18

19 20
  FUNCTION stack_matrices_vertically(m1, m2) result(res)
    ! appends to matrix m1 the rows of matrix m2
21 22
    type(inum) :: m1(:,:), m2(:,:)
    type(inum) :: res(size(m1,1)+size(m2,1), size(m1,2))
23 24 25 26 27
    res(1:size(m1,1), :) = m1
    res(size(m1,1)+1:size(res,1),:) = m2 
  END FUNCTION stack_matrices_vertically

  RECURSIVE FUNCTION shuffle_product(v1, v2) result(res)
28
    type(inum) :: v1(:), v2(:)
29
    integer :: i
30
    type(inum) :: res(product((/(i,i=1,size(v1)+size(v2))/))/  &
31 32
      (product((/(i,i=1,size(v1))/))*product((/(i,i=1,size(v2))/))), & 
      size(v1) + size(v2))
33
    type(inum) :: alpha, beta, w1(size(v1)-1), w2(size(v2)-1)
34

35
    res = izero
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
    if(size(v1) == 0) then 
      res(1,:) = v2
      return
    else if(size(v2) == 0) then 
      res(1,:) = v1
      return
    end if

    alpha = v1(1)
    beta = v2(1)
    w1 = v1(2:size(v1))
    w2 = v2(2:size(v2))

    res = stack_matrices_vertically( &
      append_to_each_row(alpha, shuffle_product(w1, v2)), & 
      append_to_each_row(beta, shuffle_product(v1, w2)) )
Luca's avatar
Luca committed
52

53 54
  END FUNCTION shuffle_product

ulrich_y's avatar
ulrich_y committed
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72

  FUNCTION shuffle_with_zero(a) result(res)
    ! rows of result are shuffles of a with 0
    type(inum) :: a(:)
    type(inum) :: res(size(a)+1,size(a)+1)
    integer :: i,j, N
    N = size(a)+1
    do i = 1,N
      ! i is the index of the row
      ! j is the index of the zero
      j  = N+1-i
      res(i,j) = izero
      res(i,1:j-1) = a(1:j-1)
      res(i,j+1:N) = a(j:)
    end do
  END FUNCTION shuffle_with_zero


73 74
END MODULE shuffle

Luca's avatar
minor  
Luca committed
75 76 77 78
! PROGRAM test
!   use utils
!   use shuffle
!   implicit none
79

80
!   type(inum) :: v1(3), v2(2)
Luca's avatar
minor  
Luca committed
81
!   integer :: amount_shuffles
82

Luca's avatar
minor  
Luca committed
83 84
!   v1 = cmplx((/1,2,3/))
!   v2 = cmplx((/4,5/))
85

Luca's avatar
minor  
Luca committed
86
!   call print_matrix(shuffle_product(v1,v2))
87

Luca's avatar
minor  
Luca committed
88
! END PROGRAM test
89