Worksharing Constructs (Others)¶
Most of the time, we end up having more than one loop — a nested loop — where two or three loops will be next to each other. OpenMP provides a clause for handling this kind of situation with collapse. To understand this, we will now study matrix multiplication, which involves a nested loop.
Collapse¶
The collapse clause can be used for nested loops — an entire part of the iteration will be divided by the available number of threads. The figure below shows an example of not using the collapse clause — only the outer loop is parallelised.
This is not what we want. With the collapse clause, we notice that all available threads execute every single iteration, as seen in the figure below.
Collapse syntax
#pragma omp parallel
#pragma omp for collapse(2)
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
{ cout << " Thread id" << " " << omp_get_thread_num() << endl; }
// Or
#pragma omp parallel for collapse(2)
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
{ cout << " Thread id" << " " << omp_get_thread_num() << endl; }
!$omp parallel
!$omp do collapse(2)
do i = 1, n
do j = 1, n
print*, 'Thread id', omp_get_thread_num()
end do
end do
!$omp end do
!$omp end parallel
!! Or
!$omp parallel do collapse(2)
do i = 1, n
do j = 1, n
print*, 'Thread id', omp_get_thread_num()
end do
end do
!$omp end parallel do
Examples and Questions: Collapse
#include <iostream>
#include <omp.h>
using namespace std;
int main()
{
int N=5;
#pragma omp parallel
#pragma omp for collapse(2)
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
cout << "Outer loop id " << i << " Inner loop id "<< j
<< " Thread id " << omp_get_thread_num() << endl;
return 0;
}
program main
use omp_lib
implicit none
integer :: n, i, j
n=5
!$omp parallel
!$omp do collapse(2)
do i = 1, n
do j = 1, n
print*, 'Outer loop id ', i, 'Inner loop id ', j, 'Thread id', omp_get_thread_num()
end do
end do
!$omp end do
!$omp end parallel
end program main
Outer loop id 4 Inner loop id 2 Thread id 16
Outer loop id 1 Inner loop id 4 Thread id 3
Outer loop id 5 Inner loop id 1 Thread id 20
Outer loop id 4 Inner loop id 1 Thread id 15
Outer loop id 2 Inner loop id 1 Thread id 5
Outer loop id 3 Inner loop id 1 Thread id 10
Outer loop id 1 Inner loop id 2 Thread id 1
Outer loop id 2 Inner loop id 3 Thread id 7
Questions
Can you add any of the scheduling clauses here, such as
static,dynamic, etc.?Is it really necessary to use them when you use
collapse, or is it dependent on other factors such as the nature of the computation and available threads?
Reduction¶
The reduction clauses are data-sharing attribute clauses that can be used to perform some forms of repetition calculations in the parallel region.
It can be used for arithmetic reductions:
+,*,-,max,minAnd also with logical operator reductions in C:
&,&&,|,||,^
Reduction syntax
#pragma omp parallel
#pragma omp for reduction(+:sum)
for(int i = 0; i < N; i++)
{ sum +=a[i]; }
// Or
#pragma omp parallel for reduction(+:sum)
for(int i = 0; i < N; i++)
{ sum += a[i]; }
!$omp parallel
!$omp do reduction(+:sum)
do i = 1, n
sum = sum + a(i)
end do
!$omp end do
!$omp end parallel
!! Or
!$omp parallel do reduction(+:sum)
do i = 1, n
sum = sum + a(i)
end do
!$omp end parallel do
Examples and Question: Reduction
#include <iostream>
#include <omp.h>
using namespace std;
int main()
{
int sum, N = 10;
float *a = (float*)malloc(sizeof(float) * N);
#pragma omp parallel for reduction(+:sum)
for(int i = 0; i < N; i++)
{
a[i] = i;
sum += a[i];
}
cout << "Sum is "<< sum << endl;
return 0;
}
program main
use omp_lib
implicit none
real(8), dimension(:), allocatable :: a
integer :: n, i, sum
n=10
allocate(a(n))
!$omp parallel do reduction(+:sum)
do i = 1, n
a(i) = i
sum = sum + a(i)
end do
!$omp end parallel do
print *, 'Sum is ', sum
end program main
Questions
What happens if you do not use the reduction clause? Do we still get the correct answer?
Matrix Multiplication¶
In this example, we consider a square matrix where M=N is equal for both A and B matrices. We create a 1D array to represent a 2D matrix. We must use the collapse clause since matrix multiplication deals with 3 loops.
Matrix multiplication function call
for(int row = 0; row < width ; ++row)
for(int col = 0; col < width ; ++col)
{
float sum=0;
for(int i = 0; i < width ; ++i)
{ sum += a[row*width+i] * b[i*width+col]; }
c[row*width+col] = sum;
}
do row = 0, width-1
do col = 0, width-1
real(8) :: sum
sum=0
do i = 0, width-1
sum = sum + (a((row*width)+i+1) * b((i*width)+col+1))
enddo
c(row*width+col+1) = sum
enddo
enddo
Questions and Solutions¶
Examples: Matrix Multiplication
#include<stdio.h>
#include<stdlib.h>
#include<omp.h>
void Matrix_Multiplication(float *a, float *b, float *c, int width)
{
for(int row = 0; row < width ; ++row)
for(int col = 0; col < width ; ++col)
{
float sum=0;
for(int i = 0; i < width ; ++i)
{ sum += a[row*width+i] * b[i*width+col]; }
c[row*width+col] = sum;
}
}
int main()
{
printf("Programme assumes that matrix size is N*N \n");
printf("Please enter the N size number \n");
int N = 0; scanf("%d", &N);
float *a, *b, *c;
a = (float*)malloc(sizeof(float) * (N*N));
b = (float*)malloc(sizeof(float) * (N*N));
c = (float*)malloc(sizeof(float) * (N*N));
for(int i = 0; i < (N*N); i++)
{ a[i] = 1.0f; b[i] = 2.0f; }
clock_t start = clock();
Matrix_Multiplication(a, b, c, N);
clock_t end = clock();
printf("Time measured: %.3f seconds.\n", (double)(end - start)/CLOCKS_PER_SEC);
for(int i = 0; i < N; i++)
{ for(int j = 0; j < N; j++) printf("%f ", c[j]); printf("\n"); }
free(a); free(b); free(c);
return 0;
}
#include<stdio.h>
#include<stdlib.h>
#include<omp.h>
void Matrix_Multiplication(float *a, float *b, float *c, int width)
{
for(int row = 0; row < width ; ++row)
for(int col = 0; col < width ; ++col)
{
float sum=0;
for(int i = 0; i < width ; ++i)
{ sum += a[row*width+i] * b[i*width+col]; }
c[row*width+col] = sum;
}
}
int main()
{
printf("Programme assumes that matrix size is N*N \n");
printf("Please enter the N size number \n");
int N = 0; scanf("%d", &N);
float *a, *b, *c;
a = (float*)malloc(sizeof(float) * (N*N));
b = (float*)malloc(sizeof(float) * (N*N));
c = (float*)malloc(sizeof(float) * (N*N));
for(int i = 0; i < (N*N); i++)
{ a[i] = 1.0f; b[i] = 2.0f; }
// ADD YOUR PARALLEL REGION
Matrix_Multiplication(a, b, c, N);
for(int i = 0; i < N; i++)
{ for(int j = 0; j < N; j++) printf("%f ", c[j]); printf("\n"); }
free(a); free(b); free(c);
return 0;
}
#include<stdio.h>
#include<stdlib.h>
#include<omp.h>
void Matrix_Multiplication(float *a, float *b, float *c, int width)
{
#pragma omp for collapse(2)
for(int row = 0; row < width ; ++row)
for(int col = 0; col < width ; ++col)
{
float sum=0;
for(int i = 0; i < width ; ++i)
{ sum += a[row*width+i] * b[i*width+col]; }
c[row*width+col] = sum;
}
}
int main()
{
printf("Programme assumes that matrix size is N*N \n");
printf("Please enter the N size number \n");
int N = 0; scanf("%d", &N);
float *a, *b, *c;
a = (float*)malloc(sizeof(float) * (N*N));
b = (float*)malloc(sizeof(float) * (N*N));
c = (float*)malloc(sizeof(float) * (N*N));
for(int i = 0; i < (N*N); i++)
{ a[i] = 1.0f; b[i] = 2.0f; }
omp_set_num_threads(omp_get_max_threads());
double start = omp_get_wtime();
#pragma omp parallel
Matrix_Multiplication(a, b, c, N);
double end = omp_get_wtime();
printf("Time measured: %.3f seconds.\n", end - start);
for(int i = 0; i < N; i++)
{ for(int j = 0; j < N; j++) printf("%f ", c[j]); printf("\n"); }
free(a); free(b); free(c);
return 0;
}
Compilation and Output
// compilation
$ gcc Matrix-multiplication.c -o Matrix-Multiplication-Serial
// execution
$ ./Matrix-Multiplication-Serial
Programme assumes that the matrix (square matrix) size is N*N
Please enter the N size number
4
8 8 8 8
8 8 8 8
8 8 8 8
8 8 8 8
// compilation
$ gcc -fopenmp Matrix-multiplication-Solution.c -o Matrix-Multiplication-Solution
// execution
$ ./Matrix-Multiplication-Solution
Programme assumes that the matrix (square matrix) size is N*N
Please enter the N size number
4
8 8 8 8
8 8 8 8
8 8 8 8
8 8 8 8
// compilation
$ gfortran -fopenmp Matrix-multiplication-Solution.f90 -o Matrix-Multiplication-Solution
// execution
$ ./Matrix-Multiplication-Solution
Programme assumes that the matrix (square matrix) size is N*N
Please enter the N size number
4
8 8 8 8
8 8 8 8
8 8 8 8
8 8 8 8
Questions
Right now, we are dealing with square matrices. Could you write a code for a different matrix size while still fulfilling the matrix multiplication condition?
Could you use any one of the loop scheduling methods, such as
dynamicorstatic? Do you see any performance gain?