shuffle.f90 2.04 KB
Newer Older
1

2 3
MODULE shuffle
  use globals
Luca's avatar
Luca committed
4
  use utils
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
Luca's avatar
Luca committed
11 12 13
    complex(kind=prec) :: a, m(:,:)
    integer :: i
    complex(kind=prec) :: 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
Luca's avatar
Luca committed
21 22
    complex(kind=prec) :: m1(:,:), m2(:,:)
    complex(kind=prec) :: 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)
Luca's avatar
Luca committed
28
    complex(kind=prec) :: v1(:), v2(:)
29
    integer :: i
Luca's avatar
Luca committed
30
    complex(kind=prec) :: 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))
Luca's avatar
Luca committed
33
    complex(kind=prec) :: p1(product((/(i,i=1,size(v1)+size(v2)-1)/))/  & 
34 35
      (product((/(i,i=1,size(v1)-1)/))*product((/(i,i=1,size(v2))/))), & 
      size(v1) + size(v2) - 1)
Luca's avatar
Luca committed
36
    complex(kind=prec) :: p2(product((/(i,i=1,size(v1)+size(v2)-1)/))/  & 
37 38
      (product((/(i,i=1,size(v1))/))*product((/(i,i=1,size(v2)-1)/))), & 
      size(v1) + size(v2) - 1)
Luca's avatar
Luca committed
39
    complex(kind=prec) :: alpha, beta, w1(size(v1)-1), w2(size(v2)-1)
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

    res = 0
    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
58

59 60
  END FUNCTION shuffle_product

61 62
END MODULE shuffle

Luca's avatar
Luca committed
63 64 65 66
PROGRAM test
  use utils
  use shuffle
  implicit none
67

Luca's avatar
Luca committed
68 69
  complex(kind=prec) :: v1(3), v2(2)
  integer :: amount_shuffles
70

Luca's avatar
Luca committed
71 72
  v1 = cmplx((/1,2,3/))
  v2 = cmplx((/4,5/))
73

Luca's avatar
Luca committed
74
  call print_matrix(shuffle_product(v1,v2))
75

Luca's avatar
Luca committed
76
END PROGRAM test
77