mirror of
https://github.com/Steffo99/unimore-hpc-assignments.git
synced 2024-11-22 16:14:24 +00:00
286 lines
9.2 KiB
C
286 lines
9.2 KiB
C
/*
|
|
* 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);
|
|
}
|
|
}
|