Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
! =================================================================================================
! This file is part of the CodeVault project. The project is licensed under Apache Version 2.0.
! CodeVault is part of the EU-project PRACE-4IP (WP7.3.C).
!
! Author(s):
! Mariusz Uchronski <mariusz.uchronski@pwr.edu.pl>
!
! This example demonstrates the use of Fortran 90 with OpenMP for matrix-matrix multiplication.
! The example is set-up to perform single precision matrix-matrix multiplication.
!
! See [http://openmp.org] for the full OpenMP documentation.
!
! =================================================================================================
program gemm_openmp
implicit none
integer,parameter :: seed = 86456
integer, parameter :: n = 400
integer, parameter :: m = 600
integer, parameter :: k = 800
real, allocatable :: a(:)
real, allocatable :: b(:)
real, allocatable :: c(:)
allocate(a(n*m))
allocate(b(m*k))
allocate(c(n*k))
call fill_random(a, n, m)
call fill_random(b, m, k)
call run_gemm_openmp(a, b, c, n, m, k)
deallocate(a)
deallocate(b)
deallocate(c)
end program gemm_openmp
subroutine fill_random(a, n, m)
implicit none
integer, intent(in) :: n, m
real, dimension(n*m), intent(inout) :: a
integer :: i, j
do i=1, n
do j=1, m
a((i-1)*m+j) = rand()
end do
end do
end subroutine fill_random
subroutine print_array(a, n, m)
implicit none
integer, intent(in) :: n, m
real, dimension(n*m), intent(in) :: a
integer :: i, j
do i=1, n
do j=1, m
print '(F8.4,$)', a((i-1)*m+j)
end do
print *, ''
end do
end subroutine print_array
subroutine run_gemm_openmp(a, b, c, n, m, k)
integer, intent(in) :: n, m , k
real, dimension(n*m), intent(in) :: a
real, dimension(m*k), intent(in) :: b
real, dimension(n*k), intent(inout) :: c
integer :: i, j, l
real :: summ
!$OMP PARALLEL DO DEFAULT (NONE) &
!$OMP SHARED(n,m,k,a,b,c) PRIVATE(i,j,l,summ)
do i=1, n
do j=1, k
summ = 0
do l=1, m
summ = summ + a((i-1)*m+l) * b((l-1)*k+j)
end do
c((i-1)*k+j) = summ
end do
end do
!$OMP END PARALLEL DO
end subroutine