Homec4science
Diffusion Poisson code for introduction to parallelism (implementation/c++) Browse

rSCINTROPARALLEL/poisson.c
652c5ad5225dimplementation/c++

Last Change
Raw File

poisson.c
View Options

#include "io_binary.h"
#include "timing.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define epsilon 0.005
static inline void compute_row(int i, int n, float ** uo, float ** u,
float ** f, float * l2, float h) {
int j;
for (j = 1; j < n; j++) {
// computation of the new step
u[i][j] = 0.25 * (uo[i - 1][j] + uo[i + 1][j] + uo[i][j - 1] +
uo[i][j + 1] - f[i][j] * h * h);
// L2 norm
*l2 += (uo[i][j] - u[i][j]) * (uo[i][j] - u[i][j]);
}
}
static inline void swap_grids(int n, float ** uo, float ** u) {
int i, j;
for (i = 0; i < n + 1; i++) {
for (j = 0; j < n + 1; j++) {
uo[i][j] = u[i][j];
}
}
}
int main(int argc, char * argv[]) {
int i, j, k = 0;
float l2;
float **u, **uo, **f;
float *u_data, *uo_data, *f_data;
float h;
int N = 256;
if(argc == 2) {
N = atoi(argv[1]);
}
h = 1. / (N + 1);
// Allocation of the data storage
u_data = (float *) malloc((N + 1) * (N + 1) * sizeof(float));
uo_data = (float *) malloc((N + 1) * (N + 1) * sizeof(float));
f_data = (float *) malloc((N + 1) * (N + 1) * sizeof(float));
// Allocation of arrays of pointers for rows
u = (float **) malloc((N + 1) * sizeof(float *));
uo = (float **) malloc((N + 1) * sizeof(float *));
f = (float **) malloc((N + 1) * sizeof(float *));
// set the row pointers in the memory
for (i = 0; i < N + 1; i++) {
u[i] = u_data + i * (N + 1);
uo[i] = uo_data + i * (N + 1);
f[i] = f_data + i * (N + 1);
}
// initialization of u0 and f
for (i = 0; i < N + 1; i++) {
for (j = 0; j < N + 1; j++) {
u[i][j] = 0;
uo[i][j] = 0;
f[i][j] = -2. * 100. * M_PI * M_PI * sin(10. * M_PI * i * h) *
sin(10. * M_PI * j * h);
}
}
double start = second();
do {
l2 = 0.;
for (i = 1; i < N; i++) {
compute_row(i, N, uo, u, f, &l2, h);
}
// copy new grid in old grid
swap_grids(N, uo, u);
write_to_file(N + 1, u, k, -1., 1.);
++k;
} while (l2 > epsilon);
double ttime = second() - start;
printf("Execution time = %f [s] (%f [ms/step]) - nb steps %d [l2 = %g]\n",
ttime, (ttime / k) * 1e3, k, l2);
// de-allocate the data
free(u_data);
free(uo_data);
free(f_data);
// de-allocate the rows pointers
free(u);
free(uo);
free(f);
return 0;
}
c4science · Help