diff --git a/README.md b/README.md index 5bc3275..9b28b59 100644 --- a/README.md +++ b/README.md @@ -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*. diff --git a/openmp/lab2/.solutions/fib-omp1.c b/openmp/lab2/.solutions/fib-omp1.c new file mode 100644 index 0000000..642e17c --- /dev/null +++ b/openmp/lab2/.solutions/fib-omp1.c @@ -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 +#include +#include +#include + +#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); + } + } +} diff --git a/openmp/lab2/.solutions/fib-omp2.c b/openmp/lab2/.solutions/fib-omp2.c new file mode 100644 index 0000000..8c7f392 --- /dev/null +++ b/openmp/lab2/.solutions/fib-omp2.c @@ -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 +#include +#include +#include + +#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); + } + } +} diff --git a/openmp/lab2/.solutions/fib-omp3.c b/openmp/lab2/.solutions/fib-omp3.c new file mode 100644 index 0000000..60e3122 --- /dev/null +++ b/openmp/lab2/.solutions/fib-omp3.c @@ -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 +#include +#include +#include + +#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); + } + } +} diff --git a/openmp/lab2/.solutions/jacobi-omp1.c b/openmp/lab2/.solutions/jacobi-omp1.c new file mode 100644 index 0000000..e78312c --- /dev/null +++ b/openmp/lab2/.solutions/jacobi-omp1.c @@ -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 +#include +#include +#include +#include + +#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); + } +} diff --git a/openmp/lab2/.solutions/jacobi-omp2.c b/openmp/lab2/.solutions/jacobi-omp2.c new file mode 100644 index 0000000..fec9efd --- /dev/null +++ b/openmp/lab2/.solutions/jacobi-omp2.c @@ -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 +#include +#include +#include +#include + +#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); + } +} diff --git a/openmp/lab2/.solutions/jacobi-omp3.c b/openmp/lab2/.solutions/jacobi-omp3.c new file mode 100644 index 0000000..eff31f9 --- /dev/null +++ b/openmp/lab2/.solutions/jacobi-omp3.c @@ -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 +#include +#include +#include +#include + +#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); + } +} diff --git a/openmp/lab2/.solutions/pi-omp1.c b/openmp/lab2/.solutions/pi-omp1.c new file mode 100644 index 0000000..09a9025 --- /dev/null +++ b/openmp/lab2/.solutions/pi-omp1.c @@ -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 +#include + +#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 +#include +#include +#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 +#include + +#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 +#include +#include +#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 +#include + +#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 +#include +#include +#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 +#include + +#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 +#include +#include + +#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)); +} diff --git a/openmp/lab2/Makefile b/openmp/lab2/Makefile new file mode 100644 index 0000000..52896c0 --- /dev/null +++ b/openmp/lab2/Makefile @@ -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 *~ + diff --git a/openmp/lab2/data/jacobi-1000.bin b/openmp/lab2/data/jacobi-1000.bin new file mode 100644 index 0000000..3e46b9c Binary files /dev/null and b/openmp/lab2/data/jacobi-1000.bin differ diff --git a/openmp/lab2/data/jacobi-500.bin b/openmp/lab2/data/jacobi-500.bin new file mode 100644 index 0000000..e62657f Binary files /dev/null and b/openmp/lab2/data/jacobi-500.bin differ diff --git a/openmp/lab2/fibonacci.c b/openmp/lab2/fibonacci.c new file mode 100644 index 0000000..06985bf --- /dev/null +++ b/openmp/lab2/fibonacci.c @@ -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 +#include +#include +#include + +#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); + } + } +} diff --git a/openmp/lab2/jacobi.c b/openmp/lab2/jacobi.c new file mode 100644 index 0000000..feca2e7 --- /dev/null +++ b/openmp/lab2/jacobi.c @@ -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 +#include +#include +#include +#include + +#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); + } +} diff --git a/openmp/lab2/pi.c b/openmp/lab2/pi.c new file mode 100644 index 0000000..d23253d --- /dev/null +++ b/openmp/lab2/pi.c @@ -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 +#include + +#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 +#include +#include +#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 +#include +#include +#include +#include + +#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, ×tampA); + asm volatile("" :: + : "memory"); +} + +void stop_timer() +{ + unsigned long long elapsed = 0ULL; + asm volatile("" :: + : "memory"); + clock_gettime(CLOCK_MONOTONIC_RAW, ×tampB); + 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 diff --git a/openmp/lab2/utils.h b/openmp/lab2/utils.h new file mode 100644 index 0000000..b552540 --- /dev/null +++ b/openmp/lab2/utils.h @@ -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 + +#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__*/