kmazouzi / FormationMeso

Download as
Browse Code ยป

Commit 9cad13733896913c2918ed58ee8f7815593f7448

Authored by kmazouzi
1 parent 6c0f7fa151
Exists in master

MAJ

Showing 5 changed files with 7 additions and 7 deletions Inline Diff

lab1/pi_mpi.c View file @ 9cad137
#include <stdio.h> 1 1 #include <stdio.h>
#include <string.h> 2 2 #include <string.h>
#include "mpi/mpi.h" 3 3 #include "mpi.h"
#define N 100000000 4 4 #define N 100000000
5 5
int main (int argc, char** argv) 6 6 int main (int argc, char** argv)
{ 7 7 {
int rank; 8 8 int rank;
int size; 9 9 int size;
long long int n; 10 10 long long int n;
long long int i; 11 11 long long int i;
12 12
double l_sum,total_sum, x, h; 13 13 double l_sum,total_sum, x, h;
14 14
MPI_Init(NULL, NULL); 15 15 MPI_Init(NULL, NULL);
16 16
MPI_Comm_size(MPI_COMM_WORLD, &size); 17 17 MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank); 18 18 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
19 19
n=N; 20 20 n=N;
21 21
if(rank==0) 22 22 if(rank==0)
{ 23 23 {
if(argc==2) 24 24 if(argc==2)
{ 25 25 {
n = atoll(argv[1]); 26 26 n = atoll(argv[1]);
} 27 27 }
28 28
printf("MPI version with process = %d\n", size); 29 29 printf("MPI version with process = %d\n", size);
printf("Number of intervals: %lld\n", n); 30 30 printf("Number of intervals: %lld\n", n);
} 31 31 }
32 32
33 33
MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD); 34 34 MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
35 35
36 36
h = 1.0/n; 37 37 h = 1.0/n;
38 38
l_sum = 0.0; 39 39 l_sum = 0.0;
40 40
for (i = rank; i < n; i += size) 41 41 for (i = rank; i < n; i += size)
{ 42 42 {
lab1/pi_ser.c View file @ 9cad137
#include <stdio.h> 1 1 #include <stdio.h>
#include <stdlib.h> 2 2 #include <stdlib.h>
#include <string.h> 3 3 #include <string.h>
#define N 100000000 4 4 #define N 100000000
int main (int argc, char** argv) 5 5 int main (int argc, char** argv)
{ 6 6 {
long long int n; 7 7 long long int n;
long long int i; 8 8 long long int i;
9 9
double l_sum, x, h; 10 10 double l_sum, x, h;
11 11
n=N; 12 12 n=N;
13 13
if(argc==2) 14 14 if(argc==2)
{ 15 15 {
n=atol(argv[1]); 16 16 n=atol(argv[1]);
} 17 17 }
18 18
19 19
h = 1.0/n; 20 20 h = 1.0/n;
21 21
l_sum = 0.0; 22 22 l_sum = 0.0;
23 23
for (i = 0; i < n; i ++) 24 24 for (i = 0; i < n; i ++)
{ 25 25 {
x = (i + 0.5)*h; 26 26 x = (i+0.5)*h;
lab3/mpi_heat2D.c View file @ 9cad137
/**************************************************************************** 1 1 /****************************************************************************
* DESCRIPTION: 2 2 * DESCRIPTION:
* Serial HEAT2D Example - C Version 3 3 * Serial HEAT2D Example - C Version
* This example is based on a simplified 4 4 * This example is based on a simplified
* two-dimensional heat equation domain decomposition. The initial 5 5 * two-dimensional heat equation domain decomposition. The initial
* temperature is computed to be high in the middle of the domain and 6 6 * temperature is computed to be high in the middle of the domain and
* zero at the boundaries. The boundaries are held at zero throughout 7 7 * zero at the boundaries. The boundaries are held at zero throughout
* the simulation. During the time-stepping, an array containing two 8 8 * the simulation. During the time-stepping, an array containing two
* domains is used; these domains alternate between old data and new data. 9 9 * domains is used; these domains alternate between old data and new data.
* 10 10 *
* The physical region, and the boundary conditions, are suggested 11 11 * The physical region, and the boundary conditions, are suggested
by this diagram; 12 12 by this diagram;
13 13
u = 0 14 14 u = 0
+------------------+ 15 15 +------------------+
| | 16 16 | |
u = 100 | | u = 100 17 17 u = 100 | | u = 100
| | 18 18 | |
| | 19 19 | |
| | 20 20 | |
| | 21 21 | |
+------------------+ 22 22 +------------------+
u = 100 23 23 u = 100
24 24
Interrior point : 25 25 Interrior point :
u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] ) 26 26 u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
27 27
28 28
PARALLEL MPI VERSION : 29 29 PARALLEL MPI VERSION :
30 30
+-------------------+ 31 31 +-------------------+
| | P0 m=(n-2)/P +2 32 32 | | P0 m=(n-2)/P +2
+-------------------+ 33 33 +-------------------+
| | P1 34 34 | | P1
+-------------------+ 35 35 +-------------------+
n | | .. 36 36 n | | ..
+-------------------+ 37 37 +-------------------+
| | Pq 38 38 | | Pq
+-------------------+ 39 39 +-------------------+
40 40
<-------- n --------> 41 41 <-------- n -------->
<-------n-2 ------> 42 42 <-------n-2 ------>
43 43
44 44
****************************************************************************/ 45 45 ****************************************************************************/
#include <stdio.h> 46 46 #include <stdio.h>
#include <stdlib.h> 47 47 #include <stdlib.h>
#include <math.h> 48 48 #include <math.h>
#include <mpi/mpi.h> 49 49 #include <mpi.h>
#define NN 50 50 50 #define NN 50
#define MM 50 51 51 #define MM 50
52 52
#define RING 100 53 53 #define RING 100
#define ITER_PRINT 100 54 54 #define ITER_PRINT 100
#define PRINT_DATA 1 55 55 #define PRINT_DATA 1
#define MAX_ITER 1000 56 56 #define MAX_ITER 1000
#define _EPSILON 0.01 57 57 #define _EPSILON 0.01
58 58
59 59
float update(int rank,int size, int nx,int ny, float *u, float *unew); 60 60 float update(int rank,int size, int nx,int ny, float *u, float *unew);
void inidat(int rank, int size, int nx, int ny, float *u, float *unew); 61 61 void inidat(int rank, int size, int nx, int ny, float *u, float *unew);
void prtdat(int rank, int size, int nx, int ny, float *u,const char *fnam); 62 62 void prtdat(int rank, int size, int nx, int ny, float *u,const char *fnam);
63 63
64 64
65 65
66 66
int main(int argc, char *argv[]) 67 67 int main(int argc, char *argv[])
{ 68 68 {
int N=NN,M=MM; 69 69 int N=NN,M=MM;
70 70
int rank,size; 71 71 int rank,size;
72 72
float EPSILON=_EPSILON; 73 73 float EPSILON=_EPSILON;
74 74
75 75
/* INITIALIZE MPI */ 76 76 /* INITIALIZE MPI */
MPI_Init(&argc, &argv); 77 77 MPI_Init(&argc, &argv);
78 78
/* GET THE PROCESSOR ID AND NUMBER OF PROCESSORS */ 79 79 /* GET THE PROCESSOR ID AND NUMBER OF PROCESSORS */
MPI_Comm_rank(MPI_COMM_WORLD, &rank); 80 80 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size); 81 81 MPI_Comm_size(MPI_COMM_WORLD, &size);
82 82
//Only Rank 0 read application parameters 83 83 //Only Rank 0 read application parameters
if(rank==0) { 84 84 if(rank==0) {
85 85
if(argc !=3) 86 86 if(argc !=3)
{ 87 87 {
fprintf(stderr,"usage %s N EPSILON\n ", argv[0]); 88 88 fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);
fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n"); 89 89 fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");
fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n"); 90 90 fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");
return -1; 91 91 return -1;
} 92 92 }
93 93
N = M = atoi(argv[1]); 94 94 N = M = atoi(argv[1]);
EPSILON = atof(argv[2]); 95 95 EPSILON = atof(argv[2]);
96 96
if(N % size!=0) 97 97 if(N % size!=0)
{ 98 98 {
fprintf(stderr,"Grid Size MUST be divisible by the number of processors !"); 99 99 fprintf(stderr,"Grid Size MUST be divisible by the number of processors !");
return -1; 100 100 return -1;
} 101 101 }
102 102
} 103 103 }
104 104
//Wait for rank 0 , all process start here 105 105 //Wait for rank 0 , all process start here
MPI_Barrier(MPI_COMM_WORLD); 106 106 MPI_Barrier(MPI_COMM_WORLD);
107 107
//Exchange N 108 108 //Exchange N
MPI_Bcast(&N , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD); 109 109 MPI_Bcast(&N , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD);
//Exchange EPSILON 110 110 //Exchange EPSILON
MPI_Bcast(&EPSILON , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD); 111 111 MPI_Bcast(&EPSILON , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD);
112 112
//local size 113 113 //local size
M = (N-2) / size + 2 ; 114 114 M = (N-2) / size + 2 ;
115 115
float *u = (float *)malloc(N * M * sizeof(float)); 116 116 float *u = (float *)malloc(N * M * sizeof(float));
float *unew = (float *)malloc(N * M * sizeof(float)); 117 117 float *unew = (float *)malloc(N * M * sizeof(float));
118 118
if(u==0 || unew ==0) 119 119 if(u==0 || unew ==0)
{ 120 120 {
perror("Can't allocated data\n"); 121 121 perror("Can't allocated data\n");
return -1; 122 122 return -1;
} 123 123 }
124 124
if(rank==0) { 125 125 if(rank==0) {
126 126
printf ( "\n" ); 127 127 printf ( "\n" );
printf ( "HEATED_PLATE\n" ); 128 128 printf ( "HEATED_PLATE\n" );
printf ( " Parallel MPI version using %d processors \n",size ); 129 129 printf ( " Parallel MPI version using %d processors \n",size );
printf ( " A program to solve for the steady state temperature distribution\n" ); 130 130 printf ( " A program to solve for the steady state temperature distribution\n" );
printf ( " over a rectangular plate.\n" ); 131 131 printf ( " over a rectangular plate.\n" );
printf ( "\n" ); 132 132 printf ( "\n" );
printf ( " Spatial grid of %d by %d points.\n", N, N ); 133 133 printf ( " Spatial grid of %d by %d points.\n", N, N );
printf ( " Each processor will use grid of %d by %d points.\n", M, N ); 134 134 printf ( " Each processor will use grid of %d by %d points.\n", M, N );
printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON); 135 135 printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON);
136 136
} 137 137 }
138 138
/* Initialize grid and create input file 139 139 /* Initialize grid and create input file
* each process initialize its part 140 140 * each process initialize its part
* */ 141 141 * */
142 142
inidat(rank,size,M,N,u,unew); 143 143 inidat(rank,size,M,N,u,unew);
144 144
prtdat(rank,size,M,N,u, "initial.dat"); 145 145 prtdat(rank,size,M,N,u, "initial.dat");
146 146
147 147
float diff, globaldiff=1.0; 148 148 float diff, globaldiff=1.0;
int iter=0; 149 149 int iter=0;
150 150
/* 151 151 /*
* iterate (JACOBI ITERATION) until the new solution unew differs from the old solution u 152 152 * iterate (JACOBI ITERATION) until the new solution unew differs from the old solution u
* by no more than EPSILON. 153 153 * by no more than EPSILON.
**/ 154 154 **/
155 155
while(globaldiff> EPSILON) { 156 156 while(globaldiff> EPSILON) {
157 157
diff= update(rank,size,M,N, u, unew); 158 158 diff= update(rank,size,M,N, u, unew);
159 159
/** 160 160 /**
* COMPUTE GLOBAL CONVERGENCE 161 161 * COMPUTE GLOBAL CONVERGENCE
* */ 162 162 * */
163 163
MPI_Allreduce(&diff, &globaldiff , 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD); 164 164 MPI_Allreduce(&diff, &globaldiff , 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
165 165
166 166
if(rank==0) 167 167 if(rank==0)
if(iter%ITER_PRINT==0) 168 168 if(iter%ITER_PRINT==0)
printf("Processor #%d, iteration %d, epsilon = %f\n ", rank,iter,globaldiff); 169 169 printf("Processor #%d, iteration %d, epsilon = %f\n ", rank,iter,globaldiff);
iter++; 170 170 iter++;
} 171 171 }
172 172
prtdat(rank,size,M,N, u, "final.dat"); 173 173 prtdat(rank,size,M,N, u, "final.dat");
free(u); 174 174 free(u);
free(unew); 175 175 free(unew);
MPI_Finalize(); 176 176 MPI_Finalize();
} 177 177 }
178 178
179 179
180 180
/**************************************************************************** 181 181 /****************************************************************************
* subroutine update 182 182 * subroutine update
****************************************************************************/ 183 183 ****************************************************************************/
float update(int rank, int size, int nx,int ny, float *u, float *unew){ 184 184 float update(int rank, int size, int nx,int ny, float *u, float *unew){
int ix, iy; 185 185 int ix, iy;
float diff=0.0; 186 186 float diff=0.0;
MPI_Status status; 187 187 MPI_Status status;
188 188
189 189
/* 190 190 /*
* EXCHANGE GHOST CELL 191 191 * EXCHANGE GHOST CELL
*/ 192 192 */
if (rank > 0 && rank< size-1) 193 193 if (rank > 0 && rank< size-1)
{ 194 194 {
MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0, 195 195 MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0,
&u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status); 196 196 &u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status);
MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1, 197 197 MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1,
&u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status); 198 198 &u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status);
} 199 199 }
200 200
else if (rank == 0) 201 201 else if (rank == 0)
MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0, 202 202 MPI_Sendrecv(&u[ny*(nx-2)], ny, MPI_FLOAT, rank+1, 0,
&u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status); 203 203 &u[ny*(nx-1)], ny, MPI_FLOAT, rank+1, 1, MPI_COMM_WORLD, &status);
else if (rank == size-1) 204 204 else if (rank == size-1)
MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1, 205 205 MPI_Sendrecv(&u[ny*1], ny, MPI_FLOAT, rank-1, 1,
&u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status); 206 206 &u[ny*0], ny, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &status);
207 207
208 208
209 209
/** 210 210 /**
* PERFORM LOCAL COMPUTATION 211 211 * PERFORM LOCAL COMPUTATION
* */ 212 212 * */
213 213
for (ix = 1; ix < nx-1; ix++) { 214 214 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 215 215 for (iy = 1; iy < ny-1; iy++) {
unew[ix*ny+iy] = (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0 216 216 unew[ix*ny+iy] = (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0
; 217 217 ;
if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) 218 218 if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
{ 219 219 {
diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] ); 220 220 diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
} 221 221 }
} 222 222 }
223 223
} 224 224 }
225 225
226 226
for (ix = 1; ix < nx-1; ix++) { 227 227 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 228 228 for (iy = 1; iy < ny-1; iy++) {
u[ix*ny+iy] = unew[ix*ny+iy]; 229 229 u[ix*ny+iy] = unew[ix*ny+iy];
} 230 230 }
} 231 231 }
232 232
233 233
return diff; 234 234 return diff;
} 235 235 }
236 236
/***************************************************************************** 237 237 /*****************************************************************************
* Initialize Data 238 238 * Initialize Data
*****************************************************************************/ 239 239 *****************************************************************************/
240 240
void inidat(int rank, int size,int nx, int ny, float *u, float *unew) 241 241 void inidat(int rank, int size,int nx, int ny, float *u, float *unew)
{ 242 242 {
int ix, iy; 243 243 int ix, iy;
lab3/plot.py View file @ 9cad137
#!/usr/bin/python 1 1 #!/usr/bin/python
#-*-coding:utf8-*- 2 2 #-*-coding:utf8-*-
from sys import argv 3 3 from sys import argv
from pylab import * 4 4 from pylab import *
5 5
if __name__ == "__main__": 6 6 if __name__ == "__main__":
if len(argv) != 2: 7 7 if len(argv) != 2:
print "Usage: python", argv[0], " matrix.dat" 8 8 print "Usage: python", argv[0], " matrix.dat"
exit (0) 9 9 exit (0)
10 10
11 11
f = open ( argv[1] , 'r') 12 12 f = open ( argv[1] , 'r')
A = [ map(float,line.split()) for line in f ] 13 13 A = [ map(float,line.split()) for line in f ]
14 14
figure(1) 15 15 figure(1)
lab3/ser_heat2D.c View file @ 9cad137
/**************************************************************************** 1 1 /****************************************************************************
* DESCRIPTION: 2 2 * DESCRIPTION:
* Serial HEAT2D Example - C Version 3 3 * Serial HEAT2D Example - C Version
* This example is based on a simplified 4 4 * This example is based on a simplified
* two-dimensional heat equation domain decomposition. The initial 5 5 * two-dimensional heat equation domain decomposition. The initial
* temperature is computed to be high in the middle of the domain and 6 6 * temperature is computed to be high in the middle of the domain and
* zero at the boundaries. The boundaries are held at zero throughout 7 7 * zero at the boundaries. The boundaries are held at zero throughout
* the simulation. During the time-stepping, an array containing two 8 8 * the simulation. During the time-stepping, an array containing two
* domains is used; these domains alternate between old data and new data. 9 9 * domains is used; these domains alternate between old data and new data.
* 10 10 *
* The physical region, and the boundary conditions, are suggested 11 11 * The physical region, and the boundary conditions, are suggested
by this diagram; 12 12 by this diagram;
13 13
u = 0 14 14 u = 0
+------------------+ 15 15 +------------------+
| | 16 16 | |
u = 100 | | u = 100 17 17 u = 100 | | u = 100
| | 18 18 | |
+------------------+ 19 19 +------------------+
u = 100 20 20 u = 100
21 21
Interrior point : 22 22 Interrior point :
u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] ) 23 23 u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
24 24
****************************************************************************/ 25 25 ****************************************************************************/
#include <stdio.h> 26 26 #include <stdio.h>
#include <stdlib.h> 27 27 #include <stdlib.h>
#include <math.h> 28 28 #include <math.h>
29 29
#define NN 50 30 30 #define NN 50
#define MM 50 31 31 #define MM 50
32 32
#define ITER_PRINT 100 33 33 #define ITER_PRINT 100
#define PRINT_DATA 1 34 34 #define PRINT_DATA 1
35 35
#define _EPSILON 0.01 36 36 #define _EPSILON 0.01
37 37
38 38
float update(int nx,int ny, float *u, float *unew); 39 39 float update(int nx,int ny, float *u, float *unew);
void inidat(int nx, int ny, float *u, float *unew); 40 40 void inidat(int nx, int ny, float *u, float *unew);
void prtdat(int nx, int ny, float *u,const char *fnam); 41 41 void prtdat(int nx, int ny, float *u,const char *fnam);
42 42
43 43
44 44
45 45
int main(int argc, char *argv[]) 46 46 int main(int argc, char *argv[])
{ 47 47 {
int N=NN,M=MM; 48 48 int N=NN,M=MM;
49 49
float EPSILON=_EPSILON; 50 50 float EPSILON=_EPSILON;
51 51
if(argc !=3) 52 52 if(argc !=3)
{ 53 53 {
fprintf(stderr,"usage %s N EPSILON\n ", argv[0]); 54 54 fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);
fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n"); 55 55 fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");
fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n"); 56 56 fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");
return -1; 57 57 return -1;
} 58 58 }
59 59
N = M = atoi(argv[1]); 60 60 N = M = atoi(argv[1]);
EPSILON = atof(argv[2]); 61 61 EPSILON = atof(argv[2]);
62 62
float diff=1.0; 63 63 float diff=1.0;
64 64
float *u = (float *)malloc(N * M * sizeof(float)); 65 65 float *u = (float *)malloc(N * M * sizeof(float));
float *unew = (float *)malloc(N * M * sizeof(float)); 66 66 float *unew = (float *)malloc(N * M * sizeof(float));
67 67
if(u==0 || unew ==0) 68 68 if(u==0 || unew ==0)
{ 69 69 {
perror("Can't allocated data\n"); 70 70 perror("Can't allocated data\n");
return -1; 71 71 return -1;
} 72 72 }
73 73
printf ( "\n" ); 74 74 printf ( "\n" );
printf ( "HEATED_PLATE\n" ); 75 75 printf ( "HEATED_PLATE\n" );
printf ( " Serial version\n" ); 76 76 printf ( " Serial version\n" );
printf ( " A program to solve for the steady state temperature distribution\n" ); 77 77 printf ( " A program to solve for the steady state temperature distribution\n" );
printf ( " over a rectangular plate.\n" ); 78 78 printf ( " over a rectangular plate.\n" );
printf ( "\n" ); 79 79 printf ( "\n" );
printf ( " Spatial grid of %d by %d points.\n", M, N ); 80 80 printf ( " Spatial grid of %d by %d points.\n", M, N );
printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON); 81 81 printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON);
82 82
/* Initialize grid and create input file */ 83 83 /* Initialize grid and create input file */
printf("Initializing grid\n"); 84 84 printf("Initializing grid\n");
85 85
inidat(N, M,u,unew); 86 86 inidat(N, M,u,unew);
87 87
prtdat(N, M,u, "initial.dat"); 88 88 prtdat(N, M,u, "initial.dat");
89 89
printf("Start computing\n\n"); 90 90 printf("Start computing\n\n");
91 91
int iter=0; 92 92 int iter=0;
93 93
/* 94 94 /*
* iterate until the new solution unew differs from the old solution u 95 95 * iterate until the new solution unew differs from the old solution u
* by no more than EPSILON. 96 96 * by no more than EPSILON.
* */ 97 97 * */
98 98
while(diff> EPSILON) { 99 99 while(diff> EPSILON) {
100 100
diff= update(N, M, u, unew); 101 101 diff= update(N, M, u, unew);
102 102
if(iter%ITER_PRINT==0) 103 103 if(iter%ITER_PRINT==0)
104 104
printf("Iteration %d, diff = %f\n ", iter,diff); 105 105 printf("Iteration %d, diff = %e\n ", iter,diff);
106 106
iter++; 107 107 iter++;
} 108 108 }
109 109
prtdat(N, M, u, "final.dat"); 110 110 prtdat(N, M, u, "final.dat");
111 111
free(u); 112 112 free(u);
free(unew); 113 113 free(unew);
} 114 114 }
115 115
116 116
117 117
/**************************************************************************** 118 118 /****************************************************************************
* subroutine update 119 119 * subroutine update
****************************************************************************/ 120 120 ****************************************************************************/
float update(int nx,int ny, float *u, float *unew) 121 121 float update(int nx,int ny, float *u, float *unew)
{ 122 122 {
int ix, iy; 123 123 int ix, iy;
float diff=0.0; 124 124 float diff=0.0;
125 125
for (ix = 1; ix < nx-1; ix++) { 126 126 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 127 127 for (iy = 1; iy < ny-1; iy++) {
unew[ix*ny+iy] = (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0 128 128 unew[ix*ny+iy] = (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0
; 129 129 ;
if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) 130 130 if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
{ 131 131 {
diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] ); 132 132 diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
} 133 133 }
} 134 134 }
135 135
} 136 136 }
137 137
/** 138 138 /**
* COPY OLD DATA 139 139 * COPY OLD DATA
*/ 140 140 */
for (ix = 1; ix < nx-1; ix++) { 141 141 for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) { 142 142 for (iy = 1; iy < ny-1; iy++) {
u[ix*ny+iy] = unew[ix*ny+iy]; 143 143 u[ix*ny+iy] = unew[ix*ny+iy];
} 144 144 }
} 145 145 }
146 146
return diff; 147 147 return diff;
} 148 148 }
149 149
/***************************************************************************** 150 150 /*****************************************************************************
* Initialize Data 151 151 * Initialize Data
*****************************************************************************/ 152 152 *****************************************************************************/
void inidat(int nx, int ny, float *u, float *unew) 153 153 void inidat(int nx, int ny, float *u, float *unew)
{ 154 154 {
int ix, iy; 155 155 int ix, iy;
156 156
/* 157 157 /*
*Set boundary data and interrior values 158 158 *Set boundary data and interrior values
* */ 159 159 * */
160 160
161 161
// interior points 162 162 // interior points
for (ix = 1; ix < nx-1; ix++) 163 163 for (ix = 1; ix < nx-1; ix++)