unimore/hpc-2022-g3
unimore/hpc-2022-g3
1
Fork 0
mirror of https://github.com/Steffo99/unimore-hpc-assignments.git synced 2024-11-21 23:54:25 +00:00
Code Activity

HPC OpenMP Lab 2

Browse source
This commit is contained in:
Alessandro Capotondi 2021-04-15 18:47:41 +02:00
parent b3de6ba70d
commit 73b261c471
19 changed files with 2688 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
README.md
View file

@ -11,4 +11,5 @@ This repo contains the exercises and the tutorials used for Unimore's HPC class
### OpenMP Exercises
The exercises related to OpenMP programming model can be found in the folder `openmp`. Here the list of currectly available classes:
- `openmp\lab1`: OpenMP basics: *parallel*, *for-loop*, *sections*, and *tasking*.
- `openmp\lab2`: OpenMP Advanced: *reduction*, *tasking*, *optimizations*.

182
openmp/lab2/.solutions/fib-omp1.c Normal file
View file

@ -0,0 +1,182 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file fibonacci.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Recursive computation of Fibonacci
*
* @see https://en.wikipedia.org/wiki/Fibonacci_number
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "utils.h"
#define F_30 832040LL
#define F_40 102334155LL
#define F_50 12586269025LL
#define F_60 1548008755920LL
static int N;
static int CUTOFF;
#define SEPARATOR "------------------------------------\n"
// Parse command line arguments to set solver parameters
void parse_arguments(int argc, char *argv[]);
// Fibonacci Golden Model - DO NOT CHANGE!
unsigned long long fibonacci_g(unsigned long long n)
{
if (n < 2) return n;
return fibonacci_g(n - 2) + fibonacci_g(n - 1);
}
// Run the Fibonacci
unsigned long long fib(unsigned long long n)
{
if (n < 2)
return n;
unsigned long long x,y;
#pragma omp task shared(x)
x = fib(n - 2);
#pragma omp task shared(y)
y = fib(n - 1);
#pragma omp taskwait
return x+y;
}
int main(int argc, char *argv[])
{
parse_arguments(argc, argv);
printf(SEPARATOR);
printf("Number: %d\n", N);
printf("Cutoff: %d\n", CUTOFF);
printf(SEPARATOR);
// Run Jacobi solver
start_timer();
unsigned long long f_n;
#pragma omp parallel shared(f_n) num_threads(NTHREADS)
{
#pragma omp single nowait
{
f_n = fib(N);
}
}
stop_timer();
// Check error of final solution
unsigned long long g_n;
if(N==30)
g_n = F_30;
else if (N==40)
g_n = F_40;
else if (N==50)
g_n = F_50;
else if (N==60)
g_n = F_60;
else
g_n = fibonacci_g(N);
unsigned long long err = f_n - g_n;
printf(SEPARATOR);
printf("F(%d) = %llu\n", N, f_n);
printf("Error = %llu\n", err);
printf("Runtime = %lf ms\n", elapsed_ns() / 1E6);
printf(SEPARATOR);
return 0;
}
int parse_int(const char *str)
{
char *next;
int value = strtoul(str, &next, 10);
return strlen(next) ? -1 : value;
}
double parse_double(const char *str)
{
char *next;
double value = strtod(str, &next);
return strlen(next) ? -1 : value;
}
void parse_arguments(int argc, char *argv[])
{
// Set default values
N = 40;
CUTOFF = 20;
for (int i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "--number") || !strcmp(argv[i], "-n"))
{
if (++i >= argc || (N = parse_int(argv[i])) < 0)
{
printf("Invalid matrix order\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--cutoff") || !strcmp(argv[i], "-c"))
{
if (++i >= argc || (CUTOFF = parse_int(argv[i])) < 0)
{
printf("Invalid seed\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--help") || !strcmp(argv[i], "-h"))
{
printf("\n");
printf("Usage: ./jacobi [OPTIONS]\n\n");
printf("Options:\n");
printf(" -h --help Print this message\n");
printf(" -c --cutoff C Set task cutoff\n");
printf(" -n --number N Set the Fibonacci number\n");
printf("\n");
exit(0);
}
else
{
printf("Unrecognized argument '%s' (try '--help')\n", argv[i]);
exit(1);
}
}
}

186
openmp/lab2/.solutions/fib-omp2.c Normal file
View file

@ -0,0 +1,186 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file fibonacci.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Recursive computation of Fibonacci
*
* @see https://en.wikipedia.org/wiki/Fibonacci_number
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "utils.h"
#define F_30 832040LL
#define F_40 102334155LL
#define F_50 12586269025LL
#define F_60 1548008755920LL
#ifndef CUTOFF_DEF
#define CUTOFF_DEF 30
#endif
static int N;
static int CUTOFF;
#define SEPARATOR "------------------------------------\n"
// Parse command line arguments to set solver parameters
void parse_arguments(int argc, char *argv[]);
// Fibonacci Golden Model - DO NOT CHANGE!
unsigned long long fibonacci_g(unsigned long long n)
{
if (n < 2) return n;
return fibonacci_g(n - 2) + fibonacci_g(n - 1);
}
// Run the Fibonacci
unsigned long long fib(unsigned long long n)
{
if (n < 2)
return n;
unsigned long long x,y;
#pragma omp task if(n>CUTOFF) shared(x)
x = fib(n - 2);
#pragma omp task if(n>CUTOFF) shared(y)
y = fib(n - 1);
#pragma omp taskwait
return x+y;
}
int main(int argc, char *argv[])
{
parse_arguments(argc, argv);
printf(SEPARATOR);
printf("Number: %d\n", N);
printf("Cutoff: %d\n", CUTOFF);
printf(SEPARATOR);
// Run Jacobi solver
start_timer();
unsigned long long f_n;
#pragma omp parallel shared(f_n) num_threads(NTHREADS)
{
#pragma omp single nowait
{
f_n = fib(N);
}
}
stop_timer();
// Check error of final solution
unsigned long long g_n;
if(N==30)
g_n = F_30;
else if (N==40)
g_n = F_40;
else if (N==50)
g_n = F_50;
else if (N==60)
g_n = F_60;
else
g_n = fibonacci_g(N);
unsigned long long err = f_n - g_n;
printf(SEPARATOR);
printf("F(%d) = %llu\n", N, f_n);
printf("Error = %llu\n", err);
printf("Runtime = %lf ms\n", elapsed_ns() / 1E6);
printf(SEPARATOR);
return 0;
}
int parse_int(const char *str)
{
char *next;
int value = strtoul(str, &next, 10);
return strlen(next) ? -1 : value;
}
double parse_double(const char *str)
{
char *next;
double value = strtod(str, &next);
return strlen(next) ? -1 : value;
}
void parse_arguments(int argc, char *argv[])
{
// Set default values
N = 40;
CUTOFF = CUTOFF_DEF;
for (int i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "--number") || !strcmp(argv[i], "-n"))
{
if (++i >= argc || (N = parse_int(argv[i])) < 0)
{
printf("Invalid matrix order\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--cutoff") || !strcmp(argv[i], "-c"))
{
if (++i >= argc || (CUTOFF = parse_int(argv[i])) < 0)
{
printf("Invalid seed\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--help") || !strcmp(argv[i], "-h"))
{
printf("\n");
printf("Usage: ./jacobi [OPTIONS]\n\n");
printf("Options:\n");
printf(" -h --help Print this message\n");
printf(" -c --cutoff C Set task cutoff\n");
printf(" -n --number N Set the Fibonacci number\n");
printf("\n");
exit(0);
}
else
{
printf("Unrecognized argument '%s' (try '--help')\n", argv[i]);
exit(1);
}
}
}

188
openmp/lab2/.solutions/fib-omp3.c Normal file
View file

@ -0,0 +1,188 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file fibonacci.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Recursive computation of Fibonacci
*
* @see https://en.wikipedia.org/wiki/Fibonacci_number
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "utils.h"
#define F_30 832040LL
#define F_40 102334155LL
#define F_50 12586269025LL
#define F_60 1548008755920LL
#ifndef CUTOFF_DEF
#define CUTOFF_DEF 30
#endif
static int N;
static int CUTOFF;
#define SEPARATOR "------------------------------------\n"
// Parse command line arguments to set solver parameters
void parse_arguments(int argc, char *argv[]);
// Fibonacci Golden Model - DO NOT CHANGE!
unsigned long long fibonacci_g(unsigned long long n)
{
if (n < 2) return n;
return fibonacci_g(n - 2) + fibonacci_g(n - 1);
}
// Run the Fibonacci
unsigned long long fib(unsigned long long n)
{
if (n < 2)
return n;
if (n <= CUTOFF)
return fibonacci_g(n);
unsigned long long x,y;
#pragma omp task shared(x)
x = fib(n - 2);
#pragma omp task shared(y)
y = fib(n - 1);
#pragma omp taskwait
return x+y;
}
int main(int argc, char *argv[])
{
parse_arguments(argc, argv);
printf(SEPARATOR);
printf("Number: %d\n", N);
printf("Cutoff: %d\n", CUTOFF);
printf(SEPARATOR);
// Run Jacobi solver
start_timer();
unsigned long long f_n;
#pragma omp parallel shared(f_n) num_threads(NTHREADS)
{
#pragma omp single nowait
{
f_n = fib(N);
}
}
stop_timer();
// Check error of final solution
unsigned long long g_n;
if(N==30)
g_n = F_30;
else if (N==40)
g_n = F_40;
else if (N==50)
g_n = F_50;
else if (N==60)
g_n = F_60;
else
g_n = fibonacci_g(N);
unsigned long long err = f_n - g_n;
printf(SEPARATOR);
printf("F(%d) = %llu\n", N, f_n);
printf("Error = %llu\n", err);
printf("Runtime = %lf ms\n", elapsed_ns() / 1E6);
printf(SEPARATOR);
return 0;
}
int parse_int(const char *str)
{
char *next;
int value = strtoul(str, &next, 10);
return strlen(next) ? -1 : value;
}
double parse_double(const char *str)
{
char *next;
double value = strtod(str, &next);
return strlen(next) ? -1 : value;
}
void parse_arguments(int argc, char *argv[])
{
// Set default values
N = 40;
CUTOFF = CUTOFF_DEF;
for (int i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "--number") || !strcmp(argv[i], "-n"))
{
if (++i >= argc || (N = parse_int(argv[i])) < 0)
{
printf("Invalid matrix order\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--cutoff") || !strcmp(argv[i], "-c"))
{
if (++i >= argc || (CUTOFF = parse_int(argv[i])) < 0)
{
printf("Invalid seed\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--help") || !strcmp(argv[i], "-h"))
{
printf("\n");
printf("Usage: ./jacobi [OPTIONS]\n\n");
printf("Options:\n");
printf(" -h --help Print this message\n");
printf(" -c --cutoff C Set task cutoff\n");
printf(" -n --number N Set the Fibonacci number\n");
printf("\n");
exit(0);
}
else
{
printf("Unrecognized argument '%s' (try '--help')\n", argv[i]);
exit(1);
}
}
}

279
openmp/lab2/.solutions/jacobi-omp1.c Normal file
View file

@ -0,0 +1,279 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file jacobi.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief This code solves the steady state heat equation on a rectangular region.
* This code solves the steady state heat equation on a rectangular region.
* The sequential version of this program needs approximately
* 18/epsilon iterations to complete.
* The physical region, and the boundary conditions, are suggested
* by this diagram;
* W = 0
* +------------------+
* | |
* W = 100 | | W = 100
* | |
* +------------------+
* W = 100
* The region is covered with a grid of M by N nodes, and an N by N
* array W is used to record the temperature. The correspondence between
* array indices and locations in the region is suggested by giving the
* indices of the four corners:
* I = 0
* [0][0]-------------[0][N-1]
* | |
* J = 0 | | J = N-1
* | |
* [M-1][0]-----------[M-1][N-1]
* I = M-1
* The steady state solution to the discrete heat equation satisfies the
* following condition at an interior grid point:
* W[Central] = (1/4) * ( W[North] + W[South] + W[East] + W[West] )
* where "Central" is the index of the grid point, "North" is the index
* of its immediate neighbor to the "north", and so on.
*
* Given an approximate solution of the steady state heat equation, a
* "better" solution is given by replacing each interior point by the
* average of its 4 neighbors - in other words, by using the condition
* as an ASSIGNMENT statement:
* W[Central] <= (1/4) * ( W[North] + W[South] + W[East] + W[West] )
* If this process is repeated often enough, the difference between successive
* estimates of the solution will go to zero.
* This program carries out such an iteration, using a tolerance specified by
* the user, and writes the final estimate of the solution to a file that can
* be used for graphic processing.
* icensing:
* This code is distributed under the GNU LGPL license.
* odified:
* 18 October 2011
* uthor:
* Original C version by Michael Quinn.
* This C version by John Burkardt.
* eference:
* Michael Quinn,
* Parallel Programming in C with MPI and OpenMP,
* McGraw-Hill, 2004,
* ISBN13: 978-0071232654,
* LC: QA76.73.C15.Q55.
* ocal parameters:
* Local, double DIFF, the norm of the change in the solution from one iteration
* to the next.
* Local, double MEAN, the average of the boundary values, used to initialize
* the values of the solution in the interior.
* Local, double U[M][N], the solution at the previous iteration.
* Local, double W[M][N], the solution computed at the latest iteration.
*
*
* @see https://en.wikipedia.org/wiki/Jacobi_method
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include "utils.h"
static int N;
static int MAX_ITERATIONS;
static int SEED;
static double CONVERGENCE_THRESHOLD;
static FILE *data;
#define SEPARATOR "------------------------------------\n"
// Return the current time in seconds since the Epoch
double get_timestamp();
// Parse command line arguments to set solver parameters
void parse_arguments(int argc, char *argv[]);
// Run the Jacobi solver
// Returns the number of iterations performed
int run(double *A, double *xtmp)
{
int iter = 0, iterations_print = 1;
double err = 0.0;
do
{
err = 0.0;
#pragma omp parallel for reduction(max \
: err) num_threads(NTHREADS)
for (int i = 1; i < N - 1; i++)
{
for (int j = 1; j < N - 1; j++)
{
xtmp[i * N + j] = 0.25 * (A[(i - 1) * N + j] + A[(i + 1) * N + j] + A[i * N + j - 1] + A[i * N + j + 1]);
err = fmax(err, fabs(xtmp[i * N + j] - A[i * N + j]));
}
}
#pragma omp parallel for num_threads(NTHREADS)
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
A[i * N + j] = xtmp[i * N + j];
}
}
iter++;
#ifdef DEBUG
if (iter == iterations_print)
{
printf(" %8d %f\n", iter, err);
iterations_print = 2 * iterations_print;
}
#endif
} while (err > CONVERGENCE_THRESHOLD && iter < MAX_ITERATIONS);
return iter;
}
int main(int argc, char *argv[])
{
parse_arguments(argc, argv);
double *A = malloc(N * N * sizeof(double));
double *xtmp = malloc(N * N * sizeof(double));
printf(SEPARATOR);
printf("Matrix size: %dx%d\n", N, N);
printf("Maximum iterations: %d\n", MAX_ITERATIONS);
printf("Convergence threshold: %lf\n", CONVERGENCE_THRESHOLD);
printf(SEPARATOR);
for (int ii = 0; ii < N; ii++)
{
for (int jj = 0; jj < N; jj++)
{
double f;
fread(&f, sizeof(double), 1, data);
A[ii * N + jj] = f;
}
}
// Run Jacobi solver
start_timer();
int itr = run(A, xtmp);
stop_timer();
printf("Iterations = %d\n", itr);
printf("Solver runtime = %lf ms\n", elapsed_ns() / 1E6);
if (itr == MAX_ITERATIONS)
printf("WARNING: solution did not converge\n");
printf(SEPARATOR);
free(A);
free(xtmp);
fclose(data);
return 0;
}
int parse_int(const char *str)
{
char *next;
int value = strtoul(str, &next, 10);
return strlen(next) ? -1 : value;
}
double parse_double(const char *str)
{
char *next;
double value = strtod(str, &next);
return strlen(next) ? -1 : value;
}
void parse_arguments(int argc, char *argv[])
{
// Set default values
N = 500;
MAX_ITERATIONS = 2000;
CONVERGENCE_THRESHOLD = 0.001;
SEED = 0;
for (int i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "--convergence") || !strcmp(argv[i], "-c"))
{
if (++i >= argc || (CONVERGENCE_THRESHOLD = parse_double(argv[i])) < 0)
{
printf("Invalid convergence threshold\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--iterations") || !strcmp(argv[i], "-i"))
{
if (++i >= argc || (MAX_ITERATIONS = parse_int(argv[i])) < 0)
{
printf("Invalid number of iterations\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--norder") || !strcmp(argv[i], "-n"))
{
if (++i >= argc || (N = parse_int(argv[i])) < 0)
{
printf("Invalid matrix order\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--help") || !strcmp(argv[i], "-h"))
{
printf("\n");
printf("Usage: ./jacobi [OPTIONS]\n\n");
printf("Options:\n");
printf(" -h --help Print this message\n");
printf(" -c --convergence C Set convergence threshold\n");
printf(" -i --iterations I Set maximum number of iterations\n");
printf(" -n --norder N Set maxtrix order (500 or 1000)\n");
printf("\n");
exit(0);
}
else
{
printf("Unrecognized argument '%s' (try '--help')\n", argv[i]);
exit(1);
}
}
if (N == 1000)
data = fopen("data/jacobi-1000.bin", "rb");
else if (N == 500)
data = fopen("data/jacobi-500.bin", "rb");
else
{
printf("Invalid matrix order\n");
exit(1);
}
}

284
openmp/lab2/.solutions/jacobi-omp2.c Normal file
View file

@ -0,0 +1,284 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file jacobi.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief This code solves the steady state heat equation on a rectangular region.
* This code solves the steady state heat equation on a rectangular region.
* The sequential version of this program needs approximately
* 18/epsilon iterations to complete.
* The physical region, and the boundary conditions, are suggested
* by this diagram;
* W = 0
* +------------------+
* | |
* W = 100 | | W = 100
* | |
* +------------------+
* W = 100
* The region is covered with a grid of M by N nodes, and an N by N
* array W is used to record the temperature. The correspondence between
* array indices and locations in the region is suggested by giving the
* indices of the four corners:
* I = 0
* [0][0]-------------[0][N-1]
* | |
* J = 0 | | J = N-1
* | |
* [M-1][0]-----------[M-1][N-1]
* I = M-1
* The steady state solution to the discrete heat equation satisfies the
* following condition at an interior grid point:
* W[Central] = (1/4) * ( W[North] + W[South] + W[East] + W[West] )
* where "Central" is the index of the grid point, "North" is the index
* of its immediate neighbor to the "north", and so on.
*
* Given an approximate solution of the steady state heat equation, a
* "better" solution is given by replacing each interior point by the
* average of its 4 neighbors - in other words, by using the condition
* as an ASSIGNMENT statement:
* W[Central] <= (1/4) * ( W[North] + W[South] + W[East] + W[West] )
* If this process is repeated often enough, the difference between successive
* estimates of the solution will go to zero.
* This program carries out such an iteration, using a tolerance specified by
* the user, and writes the final estimate of the solution to a file that can
* be used for graphic processing.
* icensing:
* This code is distributed under the GNU LGPL license.
* odified:
* 18 October 2011
* uthor:
* Original C version by Michael Quinn.
* This C version by John Burkardt.
* eference:
* Michael Quinn,
* Parallel Programming in C with MPI and OpenMP,
* McGraw-Hill, 2004,
* ISBN13: 978-0071232654,
* LC: QA76.73.C15.Q55.
* ocal parameters:
* Local, double DIFF, the norm of the change in the solution from one iteration
* to the next.
* Local, double MEAN, the average of the boundary values, used to initialize
* the values of the solution in the interior.
* Local, double U[M][N], the solution at the previous iteration.
* Local, double W[M][N], the solution computed at the latest iteration.
*
*
* @see https://en.wikipedia.org/wiki/Jacobi_method
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include "utils.h"
static int N;
static int MAX_ITERATIONS;
static int SEED;
static double CONVERGENCE_THRESHOLD;
static FILE *data;
#define SEPARATOR "------------------------------------\n"
// Return the current time in seconds since the Epoch
double get_timestamp();
// Parse command line arguments to set solver parameters
void parse_arguments(int argc, char *argv[]);
// Run the Jacobi solver
// Returns the number of iterations performed
int run(double *A, double *xtmp)
{
int iter = 0, iterations_print = 1;
double err = 0.0;
#pragma omp parallel shared(err) firstprivate(iter, iterations_print) num_threads(NTHREADS)
{
do
{
#pragma omp barrier
#pragma omp single
err = 0.0;
#pragma omp for reduction(max \
: err) nowait
for (int i = 1; i < N - 1; i++)
{
for (int j = 1; j < N - 1; j++)
{
xtmp[i * N + j] = 0.25 * (A[(i - 1) * N + j] + A[(i + 1) * N + j] + A[i * N + j - 1] + A[i * N + j + 1]);
err = fmax(err, fabs(xtmp[i * N + j] - A[i * N + j]));
}
}
#pragma omp for
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
A[i * N + j] = xtmp[i * N + j];
}
}
iter++;
#ifdef DEBUG
if (iter == iterations_print)
{
printf(" %8d %f\n", iter, err);
iterations_print = 2 * iterations_print;
}
#endif
} while (err > CONVERGENCE_THRESHOLD && iter < MAX_ITERATIONS);
}
return iter;
}
int main(int argc, char *argv[])
{
parse_arguments(argc, argv);
double *A = malloc(N * N * sizeof(double));
double *xtmp = malloc(N * N * sizeof(double));
printf(SEPARATOR);
printf("Matrix size: %dx%d\n", N, N);
printf("Maximum iterations: %d\n", MAX_ITERATIONS);
printf("Convergence threshold: %lf\n", CONVERGENCE_THRESHOLD);
printf(SEPARATOR);
for (int ii = 0; ii < N; ii++)
{
for (int jj = 0; jj < N; jj++)
{
double f;
fread(&f, sizeof(double), 1, data);
A[ii * N + jj] = f;
}
}
// Run Jacobi solver
start_timer();
int itr = run(A, xtmp);
stop_timer();
printf("Iterations = %d\n", itr);
printf("Solver runtime = %lf ms\n", elapsed_ns() / 1E6);
if (itr == MAX_ITERATIONS)
printf("WARNING: solution did not converge\n");
printf(SEPARATOR);
free(A);
free(xtmp);
fclose(data);
return 0;
}
int parse_int(const char *str)
{
char *next;
int value = strtoul(str, &next, 10);
return strlen(next) ? -1 : value;
}
double parse_double(const char *str)
{
char *next;
double value = strtod(str, &next);
return strlen(next) ? -1 : value;
}
void parse_arguments(int argc, char *argv[])
{
// Set default values
N = 500;
MAX_ITERATIONS = 2000;
CONVERGENCE_THRESHOLD = 0.001;
SEED = 0;
for (int i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "--convergence") || !strcmp(argv[i], "-c"))
{
if (++i >= argc || (CONVERGENCE_THRESHOLD = parse_double(argv[i])) < 0)
{
printf("Invalid convergence threshold\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--iterations") || !strcmp(argv[i], "-i"))
{
if (++i >= argc || (MAX_ITERATIONS = parse_int(argv[i])) < 0)
{
printf("Invalid number of iterations\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--norder") || !strcmp(argv[i], "-n"))
{
if (++i >= argc || (N = parse_int(argv[i])) < 0)
{
printf("Invalid matrix order\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--help") || !strcmp(argv[i], "-h"))
{
printf("\n");
printf("Usage: ./jacobi [OPTIONS]\n\n");
printf("Options:\n");
printf(" -h --help Print this message\n");
printf(" -c --convergence C Set convergence threshold\n");
printf(" -i --iterations I Set maximum number of iterations\n");
printf(" -n --norder N Set maxtrix order (500 or 1000)\n");
printf("\n");
exit(0);
}
else
{
printf("Unrecognized argument '%s' (try '--help')\n", argv[i]);
exit(1);
}
}
if (N == 1000)
data = fopen("data/jacobi-1000.bin", "rb");
else if (N == 500)
data = fopen("data/jacobi-500.bin", "rb");
else
{
printf("Invalid matrix order\n");
exit(1);
}
}

286
openmp/lab2/.solutions/jacobi-omp3.c Normal file
View file

@ -0,0 +1,286 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file jacobi.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief This code solves the steady state heat equation on a rectangular region.
* This code solves the steady state heat equation on a rectangular region.
* The sequential version of this program needs approximately
* 18/epsilon iterations to complete.
* The physical region, and the boundary conditions, are suggested
* by this diagram;
* W = 0
* +------------------+
* | |
* W = 100 | | W = 100
* | |
* +------------------+
* W = 100
* The region is covered with a grid of M by N nodes, and an N by N
* array W is used to record the temperature. The correspondence between
* array indices and locations in the region is suggested by giving the
* indices of the four corners:
* I = 0
* [0][0]-------------[0][N-1]
* | |
* J = 0 | | J = N-1
* | |
* [M-1][0]-----------[M-1][N-1]
* I = M-1
* The steady state solution to the discrete heat equation satisfies the
* following condition at an interior grid point:
* W[Central] = (1/4) * ( W[North] + W[South] + W[East] + W[West] )
* where "Central" is the index of the grid point, "North" is the index
* of its immediate neighbor to the "north", and so on.
*
* Given an approximate solution of the steady state heat equation, a
* "better" solution is given by replacing each interior point by the
* average of its 4 neighbors - in other words, by using the condition
* as an ASSIGNMENT statement:
* W[Central] <= (1/4) * ( W[North] + W[South] + W[East] + W[West] )
* If this process is repeated often enough, the difference between successive
* estimates of the solution will go to zero.
* This program carries out such an iteration, using a tolerance specified by
* the user, and writes the final estimate of the solution to a file that can
* be used for graphic processing.
* icensing:
* This code is distributed under the GNU LGPL license.
* odified:
* 18 October 2011
* uthor:
* Original C version by Michael Quinn.
* This C version by John Burkardt.
* eference:
* Michael Quinn,
* Parallel Programming in C with MPI and OpenMP,
* McGraw-Hill, 2004,
* ISBN13: 978-0071232654,
* LC: QA76.73.C15.Q55.
* ocal parameters:
* Local, double DIFF, the norm of the change in the solution from one iteration
* to the next.
* Local, double MEAN, the average of the boundary values, used to initialize
* the values of the solution in the interior.
* Local, double U[M][N], the solution at the previous iteration.
* Local, double W[M][N], the solution computed at the latest iteration.
*
*
* @see https://en.wikipedia.org/wiki/Jacobi_method
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include "utils.h"
static int N;
static int MAX_ITERATIONS;
static int SEED;
static double CONVERGENCE_THRESHOLD;
static FILE *data;
#define SEPARATOR "------------------------------------\n"
// Return the current time in seconds since the Epoch
double get_timestamp();
// Parse command line arguments to set solver parameters
void parse_arguments(int argc, char *argv[]);
// Run the Jacobi solver
// Returns the number of iterations performed
int run(double *A, double *xtmp)
{
int iter = 0, iterations_print = 1;
double err = 0.0;
#pragma omp parallel shared(err) firstprivate(iter, iterations_print) num_threads(NTHREADS)
{
do
{
#pragma omp barrier
#pragma omp single
err = 0.0;
#pragma omp for reduction(max \
: err) nowait
for (int i = 1; i < N - 1; i++)
{
#pragma omp simd
for (int j = 1; j < N - 1; j++)
{
xtmp[i * N + j] = 0.25 * (A[(i - 1) * N + j] + A[(i + 1) * N + j] + A[i * N + j - 1] + A[i * N + j + 1]);
err = fmax(err, fabs(xtmp[i * N + j] - A[i * N + j]));
}
}
#pragma omp for
for (int i = 0; i < N; i++)
{
#pragma omp simd
for (int j = 0; j < N; j++)
{
A[i * N + j] = xtmp[i * N + j];
}
}
iter++;
#ifdef DEBUG
if (iter == iterations_print)
{
printf(" %8d %f\n", iter, err);
iterations_print = 2 * iterations_print;
}
#endif
} while (err > CONVERGENCE_THRESHOLD && iter < MAX_ITERATIONS);
}
return iter;
}
int main(int argc, char *argv[])
{
parse_arguments(argc, argv);
double *A = malloc(N * N * sizeof(double));
double *xtmp = malloc(N * N * sizeof(double));
printf(SEPARATOR);
printf("Matrix size: %dx%d\n", N, N);
printf("Maximum iterations: %d\n", MAX_ITERATIONS);
printf("Convergence threshold: %lf\n", CONVERGENCE_THRESHOLD);
printf(SEPARATOR);
for (int ii = 0; ii < N; ii++)
{
for (int jj = 0; jj < N; jj++)
{
double f;
fread(&f, sizeof(double), 1, data);
A[ii * N + jj] = f;
}
}
// Run Jacobi solver
start_timer();
int itr = run(A, xtmp);
stop_timer();
printf("Iterations = %d\n", itr);
printf("Solver runtime = %lf ms\n", elapsed_ns() / 1E6);
if (itr == MAX_ITERATIONS)
printf("WARNING: solution did not converge\n");
printf(SEPARATOR);
free(A);
free(xtmp);
fclose(data);
return 0;
}
int parse_int(const char *str)
{
char *next;
int value = strtoul(str, &next, 10);
return strlen(next) ? -1 : value;
}
double parse_double(const char *str)
{
char *next;
double value = strtod(str, &next);
return strlen(next) ? -1 : value;
}
void parse_arguments(int argc, char *argv[])
{
// Set default values
N = 500;
MAX_ITERATIONS = 1000;
CONVERGENCE_THRESHOLD = 0.001;
SEED = 0;
for (int i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "--convergence") || !strcmp(argv[i], "-c"))
{
if (++i >= argc || (CONVERGENCE_THRESHOLD = parse_double(argv[i])) < 0)
{
printf("Invalid convergence threshold\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--iterations") || !strcmp(argv[i], "-i"))
{
if (++i >= argc || (MAX_ITERATIONS = parse_int(argv[i])) < 0)
{
printf("Invalid number of iterations\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--norder") || !strcmp(argv[i], "-n"))
{
if (++i >= argc || (N = parse_int(argv[i])) < 0)
{
printf("Invalid matrix order\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--help") || !strcmp(argv[i], "-h"))
{
printf("\n");
printf("Usage: ./jacobi [OPTIONS]\n\n");
printf("Options:\n");
printf(" -h --help Print this message\n");
printf(" -c --convergence C Set convergence threshold\n");
printf(" -i --iterations I Set maximum number of iterations\n");
printf(" -n --norder N Set maxtrix order (500 or 1000)\n");
printf("\n");
exit(0);
}
else
{
printf("Unrecognized argument '%s' (try '--help')\n", argv[i]);
exit(1);
}
}
if (N == 1000)
data = fopen("data/jacobi-1000.bin", "rb");
else if (N == 500)
data = fopen("data/jacobi-500.bin", "rb");
else
{
printf("Invalid matrix order\n");
exit(1);
}
}

106
openmp/lab2/.solutions/pi-omp1.c Normal file
View file

@ -0,0 +1,106 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file exercise7.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Exercise 8
*
* Pi calculation
* @see https://dolly.fim.unimore.it/2019/course/view.php?id=152
*/
#include <stdio.h>
#include <omp.h>
#include "utils.h"
/**
* @brief EX 8- Pi Calculation
*
* This program computes pi as
* \pi = 4 arctan(1)
* = 4 \int _0 ^1 \frac{1} {1 + x^2} dx
*
* @return void
*/
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include "utils.h"
#if !defined(ITERS)
#define ITERS (4)
#endif
#define NSTEPS 134217728
void exercise(){
long i;
double dx = 1.0 / NSTEPS;
double pi = 0.0;
double start_time = omp_get_wtime();
#pragma omp parallel
for (i = 0; i < NSTEPS; i++)
{
double x = (i + 0.5) * dx;
#pragma omp critical
pi += 1.0 / (1.0 + x * x);
}
pi *= 4.0 * dx;
double run_time = omp_get_wtime() - start_time;
double ref_pi = 4.0 * atan(1.0);
printf("pi with %d steps is %.10f in %.6f seconds (error=%e)\n",
NSTEPS, pi, run_time, fabs(ref_pi - pi));
}
int
main(int argc, char** argv)
{
for(int i=0; i<ITERS; i++){
printf("\n\n");
printf("============================\n");
printf("Test - Iteration %d...\n", i);
printf("============================\n");
start_stats();
exercise();
collect_stats();
}
printf("\n\n");
printf("============================\n");
printf("Statistics\n");
printf("============================\n");
print_stats();
return 0;
}

107
openmp/lab2/.solutions/pi-omp2.c Normal file
View file

@ -0,0 +1,107 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file exercise7.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Exercise 8
*
* Pi calculation
* @see https://dolly.fim.unimore.it/2019/course/view.php?id=152
*/
#include <stdio.h>
#include <omp.h>
#include "utils.h"
/**
* @brief EX 8- Pi Calculation
*
* This program computes pi as
* \pi = 4 arctan(1)
* = 4 \int _0 ^1 \frac{1} {1 + x^2} dx
*
* @return void
*/
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include "utils.h"
#if !defined(ITERS)
#define ITERS (4)
#endif
#define NSTEPS 134217728
void exercise(){
long i;
double dx = 1.0 / NSTEPS;
double pi = 0.0;
double start_time = omp_get_wtime();
#pragma omp parallel for
for (i = 0; i < NSTEPS; i++)
{
double x = (i + 0.5) * dx;
double tmp = 1.0 / (1.0 + x * x);
#pragma omp atomic
pi = pi + tmp;
}
pi *= 4.0 * dx;
double run_time = omp_get_wtime() - start_time;
double ref_pi = 4.0 * atan(1.0);
printf("pi with %d steps is %.10f in %.6f seconds (error=%e)\n",
NSTEPS, pi, run_time, fabs(ref_pi - pi));
}
int
main(int argc, char** argv)
{
for(int i=0; i<ITERS; i++){
printf("\n\n");
printf("============================\n");
printf("Test - Iteration %d...\n", i);
printf("============================\n");
start_stats();
exercise();
collect_stats();
}
printf("\n\n");
printf("============================\n");
printf("Statistics\n");
printf("============================\n");
print_stats();
return 0;
}

105
openmp/lab2/.solutions/pi-omp3.c Normal file
View file

@ -0,0 +1,105 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file exercise7.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Exercise 8
*
* Pi calculation
* @see https://dolly.fim.unimore.it/2019/course/view.php?id=152
*/
#include <stdio.h>
#include <omp.h>
#include "utils.h"
/**
* @brief EX 8- Pi Calculation
*
* This program computes pi as
* \pi = 4 arctan(1)
* = 4 \int _0 ^1 \frac{1} {1 + x^2} dx
*
* @return void
*/
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include "utils.h"
#if !defined(ITERS)
#define ITERS (4)
#endif
#define NSTEPS 134217728
void exercise(){
long i;
double dx = 1.0 / NSTEPS;
double pi = 0.0;
double start_time = omp_get_wtime();
#pragma omp parallel for reduction(+:pi)
for (i = 0; i < NSTEPS; i++)
{
double x = (i + 0.5) * dx;
pi += 1.0 / (1.0 + x * x);
}
pi *= 4.0 * dx;
double run_time = omp_get_wtime() - start_time;
double ref_pi = 4.0 * atan(1.0);
printf("pi with %d steps is %.10f in %.6f seconds (error=%e)\n",
NSTEPS, pi, run_time, fabs(ref_pi - pi));
}
int
main(int argc, char** argv)
{
for(int i=0; i<ITERS; i++){
printf("\n\n");
printf("============================\n");
printf("Test - Iteration %d...\n", i);
printf("============================\n");
start_stats();
exercise();
collect_stats();
}
printf("\n\n");
printf("============================\n");
printf("Statistics\n");
printf("============================\n");
print_stats();
return 0;
}

80
openmp/lab2/.solutions/pi.c Normal file
View file

@ -0,0 +1,80 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file exercise7.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Exercise 8
*
* Pi calculation
* @see https://dolly.fim.unimore.it/2019/course/view.php?id=152
*/
#include <stdio.h>
#include <omp.h>
#include "utils.h"
/**
* @brief EX 8- Pi Calculation
*
* This program computes pi as
* \pi = 4 arctan(1)
* = 4 \int _0 ^1 \frac{1} {1 + x^2} dx
*
* @return void
*/
#include <stdio.h>
#include <math.h>
#include <omp.h>
#define NSTEPS 134217728
void exercise()
{
long i;
double dx = 1.0 / NSTEPS;
double pi = 0.0;
double start_time = omp_get_wtime();
#pragma omp parallel for reduction(+ \
: pi)
for (i = 0; i < NSTEPS; i++)
{
double x = (i + 0.5) * dx;
pi += 1.0 / (1.0 + x * x);
}
pi *= 4.0 * dx;
double run_time = omp_get_wtime() - start_time;
double ref_pi = 4.0 * atan(1.0);
printf("pi with %ld steps is %.10f in %.6f seconds (error=%e)\n",
NSTEPS, pi, run_time, fabs(ref_pi - pi));
}

29
openmp/lab2/Makefile Normal file
View file

@ -0,0 +1,29 @@
ifndef EXERCISE
EXERCISE=jacobi.c
endif
CC=gcc
LD=ld
OBJDUMP=objdump
OPT=-O3 -g -fopenmp
CFLAGS=$(OPT) -I. $(EXT_CFLAGS)
LDFLAGS=-lm $(EXT_LDFLAGS)
SRCS=utils.c
OBJS=$(SRCS:.c=.o) $(EXERCISE:.c=.o)
EXE=$(EXERCISE:.c=.exe)
$(EXE): $(OBJS)
$(CC) $(CFLAGS) $(OBJS) -o $@ $(LDFLAGS)
all: $(EXE)
.PHONY: run clean
run: $(EXE)
./$(EXE) $(EXT_ARGS)
clean:
rm -f $(OBJS) *.o *.exe *.out *~

BIN
openmp/lab2/data/jacobi-1000.bin Normal file
View file

Binary file not shown.

BIN
openmp/lab2/data/jacobi-500.bin Normal file
View file

Binary file not shown.

167
openmp/lab2/fibonacci.c Normal file
View file

@ -0,0 +1,167 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file fibonacci.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Recursive computation of Fibonacci
*
* @see https://en.wikipedia.org/wiki/Fibonacci_number
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "utils.h"
#define F_30 832040LL
#define F_40 102334155LL
#define F_50 12586269025LL
#define F_60 1548008755920LL
static int N;
static int CUTOFF;
#define SEPARATOR "------------------------------------\n"
// Parse command line arguments to set solver parameters
void parse_arguments(int argc, char *argv[]);
// Fibonacci Golden Model - DO NOT CHANGE!
unsigned long long fibonacci_g(unsigned long long n)
{
if (n < 2)
return n;
return fibonacci_g(n - 2) + fibonacci_g(n - 1);
}
// Run the Fibonacci
unsigned long long fib(unsigned long long n)
{
if (n < 2)
return n;
return fib(n - 2) + fib(n - 1);
}
int main(int argc, char *argv[])
{
parse_arguments(argc, argv);
printf(SEPARATOR);
printf("Number: %d\n", N);
printf("Cutoff: %d\n", CUTOFF);
printf(SEPARATOR);
// Run Jacobi solver
start_timer();
unsigned long long f_n = fib(N);
stop_timer();
// Check error of final solution
unsigned long long g_n;
if (N == 30)
g_n = F_30;
else if (N == 40)
g_n = F_40;
else if (N == 50)
g_n = F_50;
else if (N == 60)
g_n = F_60;
else
g_n = fibonacci_g(N);
unsigned long long err = f_n - g_n;
printf(SEPARATOR);
printf("F(%d) = %llu\n", N, f_n);
printf("Error = %llu\n", err);
printf("Runtime = %lf ms\n", elapsed_ns() / 1E6);
printf(SEPARATOR);
return 0;
}
int parse_int(const char *str)
{
char *next;
int value = strtoul(str, &next, 10);
return strlen(next) ? -1 : value;
}
double parse_double(const char *str)
{
char *next;
double value = strtod(str, &next);
return strlen(next) ? -1 : value;
}
void parse_arguments(int argc, char *argv[])
{
// Set default values
N = 30;
CUTOFF = 20;
for (int i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "--number") || !strcmp(argv[i], "-n"))
{
if (++i >= argc || (N = parse_int(argv[i])) < 0)
{
printf("Invalid matrix order\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--cutoff") || !strcmp(argv[i], "-c"))
{
if (++i >= argc || (CUTOFF = parse_int(argv[i])) < 0)
{
printf("Invalid seed\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--help") || !strcmp(argv[i], "-h"))
{
printf("\n");
printf("Usage: ./jacobi [OPTIONS]\n\n");
printf("Options:\n");
printf(" -h --help Print this message\n");
printf(" -c --cutoff C Set task cutoff\n");
printf(" -n --number N Set the Fibonacci number\n");
printf("\n");
exit(0);
}
else
{
printf("Unrecognized argument '%s' (try '--help')\n", argv[i]);
exit(1);
}
}
}

276
openmp/lab2/jacobi.c Normal file
View file

@ -0,0 +1,276 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file jacobi.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief This code solves the steady state heat equation on a rectangular region.
* This code solves the steady state heat equation on a rectangular region.
* The sequential version of this program needs approximately
* 18/epsilon iterations to complete.
* The physical region, and the boundary conditions, are suggested
* by this diagram;
* W = 0
* +------------------+
* | |
* W = 100 | | W = 100
* | |
* +------------------+
* W = 100
* The region is covered with a grid of M by N nodes, and an N by N
* array W is used to record the temperature. The correspondence between
* array indices and locations in the region is suggested by giving the
* indices of the four corners:
* I = 0
* [0][0]-------------[0][N-1]
* | |
* J = 0 | | J = N-1
* | |
* [M-1][0]-----------[M-1][N-1]
* I = M-1
* The steady state solution to the discrete heat equation satisfies the
* following condition at an interior grid point:
* W[Central] = (1/4) * ( W[North] + W[South] + W[East] + W[West] )
* where "Central" is the index of the grid point, "North" is the index
* of its immediate neighbor to the "north", and so on.
*
* Given an approximate solution of the steady state heat equation, a
* "better" solution is given by replacing each interior point by the
* average of its 4 neighbors - in other words, by using the condition
* as an ASSIGNMENT statement:
* W[Central] <= (1/4) * ( W[North] + W[South] + W[East] + W[West] )
* If this process is repeated often enough, the difference between successive
* estimates of the solution will go to zero.
* This program carries out such an iteration, using a tolerance specified by
* the user, and writes the final estimate of the solution to a file that can
* be used for graphic processing.
* icensing:
* This code is distributed under the GNU LGPL license.
* odified:
* 18 October 2011
* uthor:
* Original C version by Michael Quinn.
* This C version by John Burkardt.
* eference:
* Michael Quinn,
* Parallel Programming in C with MPI and OpenMP,
* McGraw-Hill, 2004,
* ISBN13: 978-0071232654,
* LC: QA76.73.C15.Q55.
* ocal parameters:
* Local, double DIFF, the norm of the change in the solution from one iteration
* to the next.
* Local, double MEAN, the average of the boundary values, used to initialize
* the values of the solution in the interior.
* Local, double U[M][N], the solution at the previous iteration.
* Local, double W[M][N], the solution computed at the latest iteration.
*
*
* @see https://en.wikipedia.org/wiki/Jacobi_method
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include "utils.h"
static int N;
static int MAX_ITERATIONS;
static int SEED;
static double CONVERGENCE_THRESHOLD;
static FILE *data;
#define SEPARATOR "------------------------------------\n"
// Return the current time in seconds since the Epoch
double get_timestamp();
// Parse command line arguments to set solver parameters
void parse_arguments(int argc, char *argv[]);
// Run the Jacobi solver
// Returns the number of iterations performed
int run(double *A, double *xtmp)
{
int iter = 0, iterations_print = 1;
double err = 0.0;
do
{
err = 0.0;
for (int i = 1; i < N - 1; i++)
{
for (int j = 1; j < N - 1; j++)
{
xtmp[i * N + j] = 0.25*(A[(i - 1) * N + j] + A[(i + 1) * N + j] + A[i * N + j - 1] + A[i * N + j + 1]);
err = fmax(err, fabs(xtmp[i * N + j] - A[i * N + j]));
}
}
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
A[i * N + j] = xtmp[i * N + j];
}
}
iter++;
#ifdef DEBUG
if (iter == iterations_print)
{
printf(" %8d %f\n", iter, err);
iterations_print = 2 * iterations_print;
}
#endif
} while (err > CONVERGENCE_THRESHOLD && iter < MAX_ITERATIONS);
return iter;
}
int main(int argc, char *argv[])
{
parse_arguments(argc, argv);
double *A = malloc(N * N * sizeof(double));
double *xtmp = malloc(N * N * sizeof(double));
printf(SEPARATOR);
printf("Matrix size: %dx%d\n", N, N);
printf("Maximum iterations: %d\n", MAX_ITERATIONS);
printf("Convergence threshold: %lf\n", CONVERGENCE_THRESHOLD);
printf(SEPARATOR);
for (int ii = 0; ii < N; ii++)
{
for (int jj = 0; jj < N; jj++)
{
double f;
fread(&f, sizeof(double), 1, data);
A[ii * N + jj] = f;
}
}
// Run Jacobi solver
start_timer();
int itr = run(A, xtmp);
stop_timer();
printf("Iterations = %d\n", itr);
printf("Solver runtime = %lf ms\n", elapsed_ns() / 1E6);
if (itr == MAX_ITERATIONS)
printf("WARNING: solution did not converge\n");
printf(SEPARATOR);
free(A);
free(xtmp);
fclose(data);
return 0;
}
int parse_int(const char *str)
{
char *next;
int value = strtoul(str, &next, 10);
return strlen(next) ? -1 : value;
}
double parse_double(const char *str)
{
char *next;
double value = strtod(str, &next);
return strlen(next) ? -1 : value;
}
void parse_arguments(int argc, char *argv[])
{
// Set default values
N = 500;
MAX_ITERATIONS = 2000;
CONVERGENCE_THRESHOLD = 0.001;
SEED = 0;
for (int i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "--convergence") || !strcmp(argv[i], "-c"))
{
if (++i >= argc || (CONVERGENCE_THRESHOLD = parse_double(argv[i])) < 0)
{
printf("Invalid convergence threshold\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--iterations") || !strcmp(argv[i], "-i"))
{
if (++i >= argc || (MAX_ITERATIONS = parse_int(argv[i])) < 0)
{
printf("Invalid number of iterations\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--norder") || !strcmp(argv[i], "-n"))
{
if (++i >= argc || (N = parse_int(argv[i])) < 0)
{
printf("Invalid matrix order\n");
exit(1);
}
}
else if (!strcmp(argv[i], "--help") || !strcmp(argv[i], "-h"))
{
printf("\n");
printf("Usage: ./jacobi [OPTIONS]\n\n");
printf("Options:\n");
printf(" -h --help Print this message\n");
printf(" -c --convergence C Set convergence threshold\n");
printf(" -i --iterations I Set maximum number of iterations\n");
printf(" -n --norder N Set maxtrix order (500 or 1000)\n");
printf("\n");
exit(0);
}
else
{
printf("Unrecognized argument '%s' (try '--help')\n", argv[i]);
exit(1);
}
}
if (N == 1000)
data = fopen("data/jacobi-1000.bin", "rb");
else if (N == 500)
data = fopen("data/jacobi-500.bin", "rb");
else
{
printf("Invalid matrix order\n");
exit(1);
}
}

104
openmp/lab2/pi.c Normal file
View file

@ -0,0 +1,104 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file exercise7.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief Exercise 8
*
* Pi calculation
* @see https://dolly.fim.unimore.it/2019/course/view.php?id=152
*/
#include <stdio.h>
#include <omp.h>
#include "utils.h"
/**
* @brief EX 8- Pi Calculation
*
* This program computes pi as
* \pi = 4 arctan(1)
* = 4 \int _0 ^1 \frac{1} {1 + x^2} dx
*
* @return void
*/
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include "utils.h"
#if !defined(ITERS)
#define ITERS (4)
#endif
#define NSTEPS 134217728
void exercise(){
long i;
double dx = 1.0 / NSTEPS;
double pi = 0.0;
double start_time = omp_get_wtime();
for (i = 0; i < NSTEPS; i++)
{
double x = (i + 0.5) * dx;
pi += 1.0 / (1.0 + x * x);
}
pi *= 4.0 * dx;
double run_time = omp_get_wtime() - start_time;
double ref_pi = 4.0 * atan(1.0);
printf("pi with %d steps is %.10f in %.6f seconds (error=%e)\n",
NSTEPS, pi, run_time, fabs(ref_pi - pi));
}
int
main(int argc, char** argv)
{
for(int i=0; i<ITERS; i++){
printf("\n\n");
printf("============================\n");
printf("Test - Iteration %d...\n", i);
printf("============================\n");
start_stats();
exercise();
collect_stats();
}
printf("\n\n");
printf("============================\n");
printf("Statistics\n");
printf("============================\n");
print_stats();
return 0;
}

150
openmp/lab2/utils.c Normal file
View file

@ -0,0 +1,150 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file utils.c
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief File containing utilities functions for HPC Unimore Class
*
* Utilities for OpenMP lab.
*
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#define _POSIX_C_SOURCE 199309L
#include <time.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include "utils.h"
#define MAX_ITERATIONS 100
static struct timespec timestampA, timestampB;
static unsigned long long statistics[MAX_ITERATIONS];
static int iterations = 0;
static unsigned long long __diff_ns(struct timespec start, struct timespec end)
{
struct timespec temp;
if ((end.tv_nsec - start.tv_nsec) < 0)
{
temp.tv_sec = end.tv_sec - start.tv_sec - 1;
temp.tv_nsec = 1000000000ULL + end.tv_nsec - start.tv_nsec;
}
else
{
temp.tv_sec = end.tv_sec - start.tv_sec;
temp.tv_nsec = end.tv_nsec - start.tv_nsec;
}
return temp.tv_nsec + temp.tv_sec * 1000000000ULL;
}
void start_timer()
{
asm volatile("" ::
: "memory");
clock_gettime(CLOCK_MONOTONIC_RAW, &timestampA);
asm volatile("" ::
: "memory");
}
void stop_timer()
{
unsigned long long elapsed = 0ULL;
asm volatile("" ::
: "memory");
clock_gettime(CLOCK_MONOTONIC_RAW, &timestampB);
asm volatile("" ::
: "memory");
}
unsigned long long elapsed_ns()
{
return __diff_ns(timestampA, timestampB);
}
void start_stats()
{
start_timer();
}
void collect_stats()
{
assert(iterations < MAX_ITERATIONS);
stop_timer();
statistics[iterations++] = elapsed_ns(timestampA, timestampB);
}
void print_stats()
{
unsigned long long min = ULLONG_MAX;
unsigned long long max = 0LL;
double average = 0.0;
double std_deviation = 0.0;
double sum = 0.0;
/* Compute the sum of all elements */
for (int i = 0; i < iterations; i++)
{
if (statistics[i] > max)
max = statistics[i];
if (statistics[i] < min)
min = statistics[i];
sum = sum + statistics[i] / 1E6;
}
average = sum / (double)iterations;
/* Compute variance and standard deviation */
for (int i = 0; i < iterations; i++)
{
sum = sum + pow((statistics[i] / 1E6 - average), 2);
}
std_deviation = sqrt(sum / (double)iterations);
printf("AvgTime\tMinTime\tMaxTime\tStdDev\n");
printf("%.4f ms\t%.4f ms\t%.4f ms\t%.4f\n", (double)average, (double)min / 1E6, (double)max / 1E6, (double)std_deviation);
}
#if defined(__GNUC__)
#pragma GCC push_options
#pragma GCC optimize("O0")
void work(unsigned long num)
#else
void work __attribute__((optnone)) (unsigned long num)
#endif
{
volatile int cnt = 0;
for (int i = 0; i < num; i++)
cnt += i;
}
#if defined(__GNUC__)
#pragma GCC pop_options
#endif

158
openmp/lab2/utils.h Normal file
View file

@ -0,0 +1,158 @@
/*
* BSD 2-Clause License
*
* Copyright (c) 2020, Alessandro Capotondi
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/**
* @file utils.h
* @author Alessandro Capotondi
* @date 27 Mar 2020
* @brief File containing utilities functions for HPC Unimore Class
*
* The header define time functions and dummy workload used on the example tests.
*
* @see http://algo.ing.unimo.it/people/andrea/Didattica/HPC/index.html
*/
#ifndef __UTILS_H__
#define __UTILS_H__
#include <stdarg.h>
#if defined(VERBOSE)
#define DEBUG_PRINT(x, ...) printf((x), ##__VA_ARGS__)
#else
#define DEBUG_PRINT(x, ...)
#endif
#if !defined(NTHREADS)
#define NTHREADS (4)
#endif
/**
* @brief The function set the timestampA
*
* The function is used to measure elapsed time between two execution points.
* The function start_timer() sets the starting point timestamp, while the function
* stop_timer() sets the termination timestamp. The elapsed time, expressed in nanoseconds,
* between the two points can be retrieved using the function elapsed_ns().
*
* Example usage:
* @code
* start_timer(); // Point A
* //SOME CODE HERE
* stop_timer(); // Point B
* printf("Elapsed time = %llu ns\n", elapsed_ns())); //Elapsed time between A and B
* //SOME OTHER CODE HERE
* stop_timer(); // Point C
* printf("Elapsed time = %llu ns\n", elapsed_ns())); //Elapsed time between A and C
* @endcode
*
* @return void
* @see start_timer()
* @see stop_timer()
* @see elapsed_ns()
*/
void start_timer();
/**
* @brief The function set the second timestamps
*
* The function is used to measure elapsed time between two execution points.
* The function start_timer() sets the starting point timestamp, while the function
* stop_timer() returns the elapsed time, expressed in nanoseconds between the last call
* of start_timer() and the current execution point.
*
* Example usage:
* @code
* start_timer(); // Point A
* //SOME CODE HERE
* stop_timer(); // Point B
* printf("Elapsed time = %llu ns\n", elapsed_ns())); //Elapsed time between A and B
* //SOME OTHER CODE HERE
* stop_timer(); // Point C
* printf("Elapsed time = %llu ns\n", elapsed_ns())); //Elapsed time between A and C
* @endcode
*
* @return void
* @see start_timer()
* @see stop_timer()
* @see elapsed_ns()
*/
void stop_timer();
/**
* @brief Elapsed nano seconds between start_timer() and stop_timer().
*
* @return Elapsed nano seconds
* @see start_timer()
* @see stop_timer()
*/
unsigned long long elapsed_ns();
/**
* @brief The function init the starting point of stat measurement.
*
* The function is similar to start_timer().
*
* @return void
* @see start_timer
*/
void start_stats();
/**
* @brief The function collects the elapsed time between the current exeuction point and the
* last call of start_stats().
*
* @return void
*/
void collect_stats();
/**
* @brief The function display the collected statistics.
* @return void
*/
void print_stats();
/**
* @brief The dummy work function
*
* The function is used to emulate some usefull workload.
*
* @param @num work duration in terms of loop iterations.
* @return void
*/
#if defined(__GNUC__)
#pragma GCC push_options
#pragma GCC optimize("O0")
void work(unsigned long num);
#else
void work __attribute__((optnone)) (unsigned long num);
#endif
#if defined(__GNUC__)
#pragma GCC pop_options
#endif
#endif /*__UTILS_H__*/