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

rSCINTROPARALLEL/poisson.c
7176566415eaimplementation/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>
#include <omp.h>
#define epsilon 0.005
static void allocate_grids(int n, float *** u, float *** uo, float *** f) {
int i;
// Allocation of the data storage
float * u_data = (float *)malloc(n * n * sizeof(float));
float * uo_data = (float *)malloc(n * n * sizeof(float));
float * f_data = (float *)malloc(n * n * sizeof(float));
// Allocation of arrays of pointers for rows
*u = (float **) malloc(n * sizeof(float *));
*uo = (float **) malloc(n * sizeof(float *));
*f = (float **) malloc(n * sizeof(float *));
// set the row pointers in the memory
for (i = 0; i < n; i++) {
(*u)[i] = u_data + i * n;
(*uo)[i] = uo_data + i * n;
(*f)[i] = f_data + i * n;
}
}
static void deallocate_grids(float ***uo, float ***u, float ***f) {
// de-allocate the data
free((*u)[0]);
free((*uo)[0]);
free((*f)[0]);
// de-allocate the rows pointers
free(*u);
free(*uo);
free(*f);
*u = NULL;
*uo = NULL;
*f = NULL;
}
static void swap_grids(float *** uo, float *** u) {
float ** tmp = *uo;
*uo = *u;
*u = tmp;
}
static inline float compute_row(int i, int n, float ** uo, float ** u,
float ** f, float h) {
int j;
float l2 = 0.;
for (j = 1; j < n - 1; 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]);
}
return l2;
}
static float compute_step(int n, float ** uo, float ** u, float ** f, float h) {
float l2 = 0.;
int i;
#pragma omp parallel for private(i), reduction(+:l2)
for (i = 1; i < n - 1; i++) {
l2 += compute_row(i, n, uo, u, f, h);
}
return l2;
}
int main(int argc, char * argv[]) {
int i, j, k = 0;
float l2;
float **u, **uo, **f;
float h;
int n = 512;
if (argc == 2) {
n = atoi(argv[1]);
}
h = 1. / n;
allocate_grids(n, &uo, &u, &f);
// initialization of u0 and f
for (i = 0; i < n; i++) {
for (j = 0; j < n; 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 = compute_step(n, uo, u, f, h);
// copy new grid in old grid
swap_grids(&uo, &u);
// write the new old grid to disk
// write_to_file(n, uo, k, -1., 1.);
++k;
} while (l2 > epsilon);
double ttime = second() - start;
printf("%d %d %f %f\n", n, k, l2, ttime);
deallocate_grids(&uo, &u, &f);
return 0;
}
c4science · Help