Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F64656716
poisson.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, May 28, 11:58
Size
2 KB
Mime Type
text/x-c
Expires
Thu, May 30, 11:58 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17936220
Attached To
rSCINTROPARALLEL Poisson code for introduction to parallelism
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;
}
Event Timeline
Log In to Comment