-
Notifications
You must be signed in to change notification settings - Fork 0
/
orig3d.c
99 lines (90 loc) · 1.92 KB
/
orig3d.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#define Max(a,b) ((a)>(b)?(a):(b))
#define N 256
double maxeps = 0.1e-7;
int itmax = 100;
int i,j,k;
double eps;
double A [N][N][N];
void relax();
void init();
void verify();
void wtime();
int main(int an, char **as)
{
int it;
double time0, time1;
init();
/* time0=omp_get_wtime (); */
wtime(&time0);
for(it=1; it<=itmax; it++)
{
eps = 0.;
relax();
printf( "it=%4i eps=%f\n", it,eps);
if (eps < maxeps) break;
}
wtime(&time1);
/* time1=omp_get_wtime (); */
printf("Time in seconds=%gs\t",time1-time0);
verify();
return 0;
}
void init()
{
for(i=0; i<=N-1; i++)
for(j=0; j<=N-1; j++)
for(k=0; k<=N-1; k++)
{
if(i==0 || i==N-1 || j==0 || j==N-1 || k==0 || k==N-1)
A[i][j][k]= 0.;
else A[i][j][k]= ( 4. + i + j + k) ;
}
}
void relax()
{
for(i=1; i<=N-2; i++)
for(j=1; j<=N-2; j++)
for(k=1; k<=N-2; k++)
{
A[i][j][k] = (A[i-1][j][k]+A[i+1][j][k])/2.;
}
for(i=1; i<=N-2; i++)
for(j=1; j<=N-2; j++)
for(k=1; k<=N-2; k++)
{
A[i][j][k] =(A[i][j-1][k]+A[i][j+1][k])/2.;
}
for(i=1; i<=N-2; i++)
for(j=1; j<=N-2; j++)
for(k=1; k<=N-2; k++)
{
double e;
e=A[i][j][k];
A[i][j][k] = (A[i][j][k-1]+A[i][j][k+1])/2.;
eps=Max(eps,fabs(e-A[i][j][k]));
}
}
void verify()
{
double s;
s=0.;
for(i=0; i<=N-1; i++)
for(j=0; j<=N-1; j++)
for(k=0; k<=N-1; k++)
{
s=s+A[i][j][k]*(i+1)*(j+1)*(k+1)/(N*N*N);
}
printf(" S = %f\n",s);
}
void wtime(double *t)
{
static int sec = -1;
struct timeval tv;
gettimeofday(&tv, (void *)0);
if (sec < 0) sec = tv.tv_sec;
*t = (tv.tv_sec - sec) + 1.0e-6*tv.tv_usec;
}