kmazouzi / FormationMeso

Commit 3604dfae01896fa69a8327b313cf2f5da11bc4c1

Authored by kmazouzi 2015-04-08 20:10:43 +0200

Exists in master

MPI version

Showing 4 changed files with 451 additions and 63 deletions Inline Diff

lab3/Makefile
lab3/mpi_heat2D.c
lab3/omp_heat2D.c
lab3/ser_heat2D.c

GCC = gcc	1	1	GCC = gcc
CFLAGS = -O3 -fopenmp	2	2	CFLAGS = -O3 -fopenmp
OMP_FLAG = -fopenmp	3	3	OMP_FLAG = -fopenmp
RM = rm -rf	4	4	RM = rm -rf
		5	MPI = mpicc
		6	MPI_FLAG = -O1 -g
		7	EXE = omp_heat2D ser_heat2D mpi_heat2D
	5	8
	6
EXE = omp_heat2D ser_heat2D	7
	8
all : $(EXE)	9	9	all : $(EXE)
	10	10
#.PHONY: all clean purge	11	11	#.PHONY: all clean purge
	12	12
	13	13
pi_ser: ser_heat2D.o	14	14	ser_heat2D: ser_heat2D.o
$(GCC) $(CFLAGS) -o $@ $^	15	15	$(GCC) $(CFLAGS) -o $@ $^
	16	16
pi_task: omp_heat2D.o	17	17	omp_heat2D: omp_heat2D.o
$(GCC) $(CFLAGS) -o $@ $^	18	18	$(GCC) $(CFLAGS) -o $@ $^
		19
		20	mpi_heat2D:
		21	$(MPI) $(MPI_FLAG) mpi_heat2D.c -o $@
		22
	19	23

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

/****************************************************************************	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>
#include <omp.h>	29	29	#include <omp.h>
	30	30
#define NN 50	31	31	#define NN 50
#define MM 50	32	32	#define MM 50
	33	33
#define ITER_PRINT 100	34	34	#define ITER_PRINT 100
#define PRINT_DATA 1	35	35	#define PRINT_DATA 1
	36	36
#define _EPSILON 0.001	37	37	#define _EPSILON 0.001
	38	38
	39	39
void update(int nx,int ny, float u, float unew, float * diff);	40	40	float update(int nx,int ny, float u, float unew);
void inidat(int nx, int ny, float u, float unew);	41	41	void inidat(int nx, int ny, float u, float unew);
void prtdat(int nx, int ny, float u,const char fnam);	42	42	void prtdat(int nx, int ny, float u,const char fnam);
	43	43
	44	44
	45	45
	46	46
int main(int argc, char *argv[])	47	47	int main(int argc, char *argv[])
{	48	48	{
	49	49
float diff=1.0;	50	50	float diff=1.0;
float EPSILON=_EPSILON;	51	51	float EPSILON=_EPSILON;
int N=NN,M=MM;	52	52	int N=NN,M=MM;
	53	53
if(argc !=3)	54	54	if(argc !=3)
{	55	55	{
fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);	56	56	fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);
fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");	57	57	fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");
fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");	58	58	fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");
return -1;	59	59	return -1;
}	60	60	}
	61	61
N = M = atoi(argv[1]);	62	62	N = M = atoi(argv[1]);
EPSILON = atof(argv[2]);	63	63	EPSILON = atof(argv[2]);
	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 ( " Parallel OpenMP version, using %d Threads\n",omp_get_max_threads() );	76	76	printf ( " Parallel OpenMP version, using %d Threads\n",omp_get_max_threads() );
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 ( " Spatial grid of %d by %d points.\n\n", M, N );	79	79	printf ( " Spatial grid of %d by %d points.\n\n", M, N );
	80	80
	81	81
/* Initialize grid and create input file */	82	82	/* Initialize grid and create input file */
printf("Initializing grid\n");	83	83	printf("Initializing grid\n");
	84	84
inidat(N, M,u,unew);	85	85	inidat(N, M,u,unew);
	86	86
prtdat(N, M,u, "initial.dat");	87	87	prtdat(N, M,u, "initial.dat");
	88	88
printf("Start computing\n");	89	89	printf("Start computing\n");
	90	90
int iter=0;	91	91	int iter=0;
	92	92
/*	93	93	/*
* iterate until the new solution unew differs from the old solution u	94	94	* iterate until the new solution unew differs from the old solution u
* by no more than EPSILON.	95	95	* by no more than EPSILON.
* */	96	96	* */
	97	97
while(diff> EPSILON) {	98	98	while(diff> EPSILON) {
	99	99
update(N, M, u, unew,&diff);	100	100	diff = update(N, M, u, unew);
	101	101
if(iter%ITER_PRINT==0)	102	102	if(iter%ITER_PRINT==0)
printf("Iteration %d, diff = %f\n ", iter,diff);	103	103	printf("Iteration %d, diff = %f\n ", iter,diff);
	104	104
iter++;	105	105	iter++;
}	106	106	}
	107	107
prtdat(N, M, u, "final.dat");	108	108	prtdat(N, M, u, "final.dat");
	109	109
free(u);	110	110	free(u);
free(unew);	111	111	free(unew);
}	112	112	}
	113	113
	114	114
	115	115
/****************************************************************************	116	116	/****************************************************************************
* subroutine update	117	117	* subroutine update
****************************************************************************/	118	118	****************************************************************************/
void update(int nx,int ny, float u, float unew, float * diff)	119	119	float update(int nx,int ny, float u, float unew)
{	120	120	{
int ix, iy;	121	121	int ix, iy;
*diff=0.0;	122	122	float diff=0.0;
	123	123
#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy)	124	124	#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy)
for (ix = 1; ix < nx-1; ix++) {	125	125	for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) {	126	126	for (iy = 1; iy < ny-1; iy++) {
unew[ix*ny+iy] =	127	127	unew[ix*ny+iy] =
(u[(ix+1)ny+iy] + u[(ix-1)ny+iy] +	128	128	(u[(ix+1)ny+iy] + u[(ix-1)ny+iy] +
u[ixny+iy+1] + u[ixny+iy-1] )/4.0;	129	129	u[ixny+iy+1] + u[ixny+iy-1] )/4.0;
	130	130
}	131	131	}
}	132	132	}
	133	133
//compute reduction	134	134	//compute reduction
	135	135
	136	136
float mydiff;	137	137	float mydiff;
	138	138
		139	/**
		140	* IMPLEMENT OMP REDUCE MAX
		141	*/
		142
#pragma omp parallel shared(nx,ny,u,unew, diff) private (ix,iy,mydiff)	139	143	#pragma omp parallel shared(nx,ny,u,unew, diff) private (ix,iy,mydiff)
{	140	144	{
mydiff=0.0;	141	145	mydiff=0.0;
#pragma omp for	142	146	#pragma omp for
for (ix = 1; ix < nx-1; ix++) {	143	147	for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) {	144	148	for (iy = 1; iy < ny-1; iy++) {
if (mydiff < fabs (unew[ixny+iy] - u[ixny+iy] ))	145	149	if (mydiff < fabs (unew[ixny+iy] - u[ixny+iy] ))
{	146	150	{
mydiff = fabs ( unew[ixny+iy] - u[ixny+iy] );	147	151	mydiff = fabs ( unew[ixny+iy] - u[ixny+iy] );
}	148	152	}
}	149	153	}
}	150	154	}
	151	155
	152	156
# pragma omp critical	153	157	#pragma omp critical
{	154	158	{
if (*diff < mydiff )	155	159	if (diff < mydiff )
{	156	160	{
*diff = mydiff;	157	161	diff = mydiff;
}	158	162	}
}	159	163	}
	160	164
	161	165	/*
		166	* COPY OLD DATA
		167	*/
#pragma omp for	162	168	#pragma omp for
for (ix = 1; ix < nx-1; ix++) {	163	169	for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) {	164	170	for (iy = 1; iy < ny-1; iy++) {
u[ixny+iy] = unew[ixny+iy];	165	171	u[ixny+iy] = unew[ixny+iy];
}	166	172	}
}	167	173	}
}	168	174	}
		175
		176	return diff;
}	169	177	}
	170	178
/*****************************************************************************	171	179	/*****************************************************************************
* Initialize Data	172	180	* Initialize Data
*****************************************************************************/	173	181	*****************************************************************************/
		182
void inidat(int nx, int ny, float u, float unew)	174	183	void inidat(int nx, int ny, float u, float unew)
{	175	184	{
int ix, iy;	176	185	int ix, iy;
	177	186
/*	178	187	/*
*Set boundary data and interrior values	179	188	*Set boundary data and interrior values
* */	180	189	* */
for (ix = 0; ix < nx; ix++)	181
for (iy = 0; iy < ny; iy++) {	182
	183	190
if(ix==0)	184	191	#pragma omp parallel private (ix,iy)
{	185	192	{
u[ix*ny+iy]=0.0;	186	193	// interior points
}	187	194	#pragma omp for
else	188	195	for (ix = 1; ix < nx-1; ix++)
if(iy==0 && ix!=0)	189	196	for (iy = 1; iy < ny-1; iy++) {
{	190	197	u[ix*ny+iy]=5.0;
u[ix*ny+iy]=100.0;	191	198	}
}else	192
	193	199
if(ix==nx-1)	194	200	//boundary left
{	195	201	#pragma omp for
u[ix*ny+iy]=100.0;	196	202	for (ix = 1; ix < nx-1; ix++){
}else	197	203	u[ix*ny]=100.0;
	198	204
if(iy==ny-1 && ix!=0)	199	205	}
{	200
u[ix*ny+iy]=100.0;	201
}else	202
	203	206
u[ix*ny+iy]=0.0;	204	207	//boundary right
}	205	208	#pragma omp for
		209	for (ix = 1; ix < nx-1; ix++){
		210	u[ix*ny+ (ny-1)]=100.0;
		211
		212	}
		213
		214	//boundary down
		215	#pragma omp for
		216	for (iy = 0; iy < ny; iy++){
		217	u[(nx-1)*(ny)+iy]=100.0;
		218
		219	}
		220
		221	//boundary top
		222	#pragma omp for
		223	for (iy = 0; iy < ny; iy++){
		224	u[iy]=0.0;
		225
		226	}
		227
		228	}
		229
}	206	230	}

/****************************************************************************	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
void update(int nx,int ny, float u, float unew, float * diff);	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
update(N, M, u, unew,&diff);	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 = %f\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	****************************************************************************/
void update(int nx,int ny, float u, float unew, float * diff)	121	121	float update(int nx,int ny, float u, float unew)
{	122	122	{
int ix, iy;	123	123	int ix, iy;
*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[ixny+iy] = (u[(ix+1)ny+iy] + u[(ix-1)ny+iy] + u[ixny+iy+1] + u[ix*ny+iy-1] )/4.0	128	128	unew[ixny+iy] = (u[(ix+1)ny+iy] + u[(ix-1)ny+iy] + u[ixny+iy+1] + u[ix*ny+iy-1] )/4.0
;	129	129	;
if (diff < fabs (unew[ixny+iy] - u[ix*ny+iy] ))	130	130	if (diff < fabs (unew[ixny+iy] - u[ixny+iy] ))
{	131	131	{
diff = fabs ( unew[ixny+iy] - u[ix*ny+iy] );	132	132	diff = fabs ( unew[ixny+iy] - u[ixny+iy] );
}	133	133	}
}	134	134	}
	135	135
}	136	136	}
	137	137
	138	138	/**
		139	* COPY OLD DATA
		140	*/
for (ix = 1; ix < nx-1; ix++) {	139	141	for (ix = 1; ix < nx-1; ix++) {
for (iy = 1; iy < ny-1; iy++) {	140	142	for (iy = 1; iy < ny-1; iy++) {
u[ixny+iy] = unew[ixny+iy];	141	143	u[ixny+iy] = unew[ixny+iy];
}	142	144	}
}	143	145	}
	144	146
		147	return diff;
}	145	148	}
	146	149
/*****************************************************************************	147	150	/*****************************************************************************
* Initialize Data	148	151	* Initialize Data
*****************************************************************************/	149	152	*****************************************************************************/
void inidat(int nx, int ny, float u, float unew)	150	153	void inidat(int nx, int ny, float u, float unew)
{	151	154	{
int ix, iy;	152	155	int ix, iy;
	153	156
/*	154	157	/*
*Set boundary data and interrior values	155	158	*Set boundary data and interrior values
* */	156	159	* */
for (ix = 0; ix < nx; ix++)	157
for (iy = 0; iy < ny; iy++) {	158
	159	160
if(ix==0)	160
{	161
u[ix*ny+iy]=0.0;	162
}	163
else	164
if(iy==0 && ix!=0)	165
{	166
u[ix*ny+iy]=100.0;	167
}else	168
	169
if(ix==nx-1)	170
{	171
u[ix*ny+iy]=100.0;	172
}else	173
	174	161
if(iy==ny-1 && ix!=0)	175	162	// interior points
{	176	163	for (ix = 1; ix < nx-1; ix++)
u[ix*ny+iy]=100.0;	177	164	for (iy = 1; iy < ny-1; iy++) {
}else	178	165	u[ix*ny+iy]=5.0;
	179
u[ix*ny+iy]=0.0;	180
}	181	166	}
		167
		168	//boundary left
		169	for (ix = 1; ix < nx-1; ix++){
		170	u[ix*ny]=100.0;
		171
		172	}
		173
		174	//boundary right
		175	for (ix = 1; ix < nx-1; ix++){
		176	u[ix*ny+ (ny-1)]=100.0;
		177
		178	}
		179
		180	//boundary down