JacobiHPC/src/impl/sequential.c

51 lines
1.3 KiB
C

/*
* Sequential version.
* No paralleliization used.
*/
#include <stdio.h>
#include <math.h>
#include <mpi.h>
#include "../config.h"
#include "../utils.h"
double *compute_jacobi(int rank, int numprocs, int n, double init_value, double threshold, borders b, int *iterations) {
double *x;
double max_diff, new_x;
int i, j;
int nb = n + 2; // n plus the border
if (numprocs != 1) {
MPI_Abort(MPI_COMM_WORLD, 1);
}
/* Initialize boundary regions */
x = create_sa_matrix(n + 2, n + 2);
for (i = 1; i <= n; i++) {
x[IDX(nb, 0, i)] = b.north;
x[IDX(nb, n + 1, i)] = b.south;
x[IDX(nb, i, 0)] = b.west;
x[IDX(nb, i, n + 1)] = b.east;
}
/* Initialize the rest of the matrix */
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
x[IDX(nb, i, j)] = init_value;
}
}
/* Iterative refinement of x until values converge */
*iterations = 0;
do {
max_diff = 0;
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
new_x = 0.25 * (x[IDX(nb, i - 1, j)] + x[IDX(nb, i, j + 1)] + x[IDX(nb, i + 1, j)] + x[IDX(nb, i, j - 1)]);
max_diff = (double) fmax(max_diff, fabs(new_x - x[IDX(nb, i, j)]));
x[IDX(nb, i, j)] = new_x;
}
}
(*iterations)++;
} while (max_diff > threshold);
return x;
}