58 lines
1.6 KiB
C
58 lines
1.6 KiB
C
/*
|
|
* OpenMP version.
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <omp.h>
|
|
#include "../config.h"
|
|
#include "../utils.h"
|
|
|
|
double *compute_jacobi(int n, double init_value, double threshold, borders b, int *iterations) {
|
|
double *x;
|
|
double *new_x;
|
|
double *tmp_x;
|
|
double max_diff, new_value;
|
|
int i, j;
|
|
int nb = n + 2; // n plus the border
|
|
|
|
/* Initialize boundary regions */
|
|
x = create_sa_matrix(nb, nb);
|
|
new_x = create_sa_matrix(nb, nb);
|
|
for (i = 0; i < nb; 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;
|
|
new_x[IDX(nb, 0, i)] = b.north;
|
|
new_x[IDX(nb, n + 1, i)] = b.south;
|
|
new_x[IDX(nb, i, 0)] = b.west;
|
|
new_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 */
|
|
omp_set_num_threads(2);
|
|
*iterations = 0;
|
|
do {
|
|
max_diff = 0;
|
|
#pragma omp parallel for schedule(static, 20) reduction (max:max_diff) private(new_value, j) firstprivate(n, nb) shared(x, new_x)
|
|
for (i = 1; i <= n; i++) {
|
|
for (j = 1; j <= n; j++) {
|
|
new_value = 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_value - x[IDX(nb, i, j)]));
|
|
new_x[IDX(nb, i, j)] = new_value;
|
|
}
|
|
}
|
|
tmp_x = new_x;
|
|
new_x = x;
|
|
x = tmp_x;
|
|
(*iterations)++;
|
|
} while (max_diff > threshold);
|
|
return x;
|
|
}
|