diff --git a/lab3/Makefile b/lab3/Makefile index 7297cfc..e77d8a3 100644 --- a/lab3/Makefile +++ b/lab3/Makefile @@ -1,5 +1,6 @@ GCC = gcc -CFLAGS = -O3 -fopenmp +CFLAGS = -O3 -fopenmp +LDFLAG = -lm OMP_FLAG = -fopenmp RM = rm -rf MPI = mpicc @@ -12,13 +13,13 @@ all : $(EXE) ser_heat2D: ser_heat2D.o - $(GCC) $(CFLAGS) -o $@ $^ + $(GCC) $(CFLAGS) -o $@ $^ $(LDFLAG) omp_heat2D: omp_heat2D.o - $(GCC) $(CFLAGS) -o $@ $^ + $(GCC) $(CFLAGS) -o $@ $^ $(LDFLAG) mpi_heat2D: - $(MPI) $(MPI_FLAG) mpi_heat2D.c -o $@ + $(MPI) $(MPI_FLAG) mpi_heat2D.c -o $@ $(LDFLAG) diff --git a/lab3/mpi_heat2D.c b/lab3/mpi_heat2D.c index a08c7a0..5bb40a9 100644 --- a/lab3/mpi_heat2D.c +++ b/lab3/mpi_heat2D.c @@ -1,6 +1,6 @@ /**************************************************************************** * DESCRIPTION: - * Serial HEAT2D Example - C Version + * MPI HEAT2D Example - C Version * This example is based on a simplified * two-dimensional heat equation domain decomposition. The initial * temperature is computed to be high in the middle of the domain and @@ -53,10 +53,12 @@ PARALLEL MPI VERSION : #define RING 100 #define ITER_PRINT 100 #define PRINT_DATA 1 -#define MAX_ITER 1000 +#define _MAX_ITER 1000 #define _EPSILON 0.01 + + 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); void prtdat(int rank, int size, int nx, int ny, float *u,const char *fnam); @@ -65,14 +67,13 @@ void prtdat(int rank, int size, int nx, int ny, float *u,const char *fnam); int main(int argc, char *argv[]) -{ +{ int N=NN,M=MM; - - int rank,size; - float EPSILON=_EPSILON; - - + int MAX_ITER= _MAX_ITER; + + int size,rank; + /* INITIALIZE MPI */ MPI_Init(&argc, &argv); @@ -83,17 +84,16 @@ int main(int argc, char *argv[]) //Only Rank 0 read application parameters if(rank==0) { - if(argc !=3) - { - fprintf(stderr,"usage %s N EPSILON\n ", argv[0]); - fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n"); - fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n"); - return -1; - } - + if(argc !=4) + { + fprintf(stderr,"usage %s N EPSILON MAX_ITER\n ", argv[0]); + fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance, MAX_ITER is max iteration\n"); + fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1, MAX_ITER=1000\n"); + return -1; + } N = M = atoi(argv[1]); EPSILON = atof(argv[2]); - + MAX_ITER = atoi (argv[3]); if(N % size!=0) { fprintf(stderr,"Grid Size MUST be divisible by the number of processors !"); @@ -106,9 +106,12 @@ int main(int argc, char *argv[]) MPI_Barrier(MPI_COMM_WORLD); //Exchange N - MPI_Bcast(&N , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD); + MPI_Bcast(&N , 1, MPI_INT, 0 , MPI_COMM_WORLD); //Exchange EPSILON MPI_Bcast(&EPSILON , 1, MPI_FLOAT, 0 , MPI_COMM_WORLD); + + //Exchange MAX_ITER + MPI_Bcast(&MAX_ITER , 1, MPI_INT, 0 , MPI_COMM_WORLD); //local size M = (N-2) / size + 2 ; @@ -153,7 +156,7 @@ int main(int argc, char *argv[]) * by no more than EPSILON. **/ - while(globaldiff> EPSILON) { + while(globaldiff> EPSILON && iter < MAX_ITER) { diff= update(rank,size,M,N, u, unew); @@ -169,6 +172,7 @@ int main(int argc, char *argv[]) printf("Processor #%d, iteration %d, epsilon = %f\n ", rank,iter,globaldiff); iter++; } + prtdat(rank,size,M,N, u, "final.dat"); free(u); @@ -215,10 +219,7 @@ float update(int rank, int size, int nx,int ny, float *u, float *unew){ 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 ; - if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) - { - diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] ); - } + diff = fmax( diff, fabs(unew[ix*ny+iy] - u[ix*ny+iy])); } } diff --git a/lab3/omp_heat2D.c b/lab3/omp_heat2D.c index ae85c38..99e9a63 100644 --- a/lab3/omp_heat2D.c +++ b/lab3/omp_heat2D.c @@ -1,6 +1,6 @@ /**************************************************************************** * DESCRIPTION: - * Serial HEAT2D Example - C Version + * OpenMP HEAT2D Example - C Version * This example is based on a simplified * two-dimensional heat equation domain decomposition. The initial * temperature is computed to be high in the middle of the domain and @@ -23,18 +23,18 @@ Interrior point : u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] ) ****************************************************************************/ -#include + +#include #include #include -#include #define NN 50 #define MM 50 #define ITER_PRINT 100 +#define _MAX_ITER 400 #define PRINT_DATA 1 - -#define _EPSILON 0.001 +#define _EPSILON 0.01 float update(int nx,int ny, float *u, float *unew); @@ -46,21 +46,23 @@ void prtdat(int nx, int ny, float *u,const char *fnam); int main(int argc, char *argv[]) { - - float diff=1.0; - float EPSILON=_EPSILON; int N=NN,M=MM; + float EPSILON=_EPSILON; + int MAX_ITER= _MAX_ITER; - if(argc !=3) + if(argc !=4) { - fprintf(stderr,"usage %s N EPSILON\n ", argv[0]); - fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n"); - fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n"); + fprintf(stderr,"usage %s N EPSILON MAX_ITER\n ", argv[0]); + fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance, MAX_ITER is max iteration\n"); + fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1, MAX_ITER=1000\n"); return -1; } N = M = atoi(argv[1]); EPSILON = atof(argv[2]); + MAX_ITER = atoi (argv[3]); + + float diff=1.0; float *u = (float *)malloc(N * M * sizeof(float)); float *unew = (float *)malloc(N * M * sizeof(float)); @@ -73,11 +75,12 @@ int main(int argc, char *argv[]) printf ( "\n" ); printf ( "HEATED_PLATE\n" ); - printf ( " Parallel OpenMP version, using %d Threads\n",omp_get_max_threads() ); + printf ( " OpenMP version\n" ); printf ( " A program to solve for the steady state temperature distribution\n" ); printf ( " over a rectangular plate.\n" ); - printf ( " Spatial grid of %d by %d points.\n\n", M, N ); - + printf ( "\n" ); + printf ( " Spatial grid of %d by %d points.\n", M, N ); + printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON); /* Initialize grid and create input file */ printf("Initializing grid\n"); @@ -86,7 +89,7 @@ int main(int argc, char *argv[]) prtdat(N, M,u, "initial.dat"); - printf("Start computing\n"); + printf("Start computing\n\n"); int iter=0; @@ -95,16 +98,20 @@ int main(int argc, char *argv[]) * by no more than EPSILON. * */ - while(diff> EPSILON) { + while(diff> EPSILON && iter EPSILON) + printf("***Not converged MAX_ITERATION %d reached\n***", MAX_ITER); + prtdat(N, M, u, "final.dat"); free(u); @@ -121,57 +128,27 @@ float update(int nx,int ny, float *u, float *unew) int ix, iy; float diff=0.0; -#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy) +#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy) reduction(max:diff) for (ix = 1; ix < nx-1; ix++) { 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; - + u[ix*ny+iy+1] + u[ix*ny+iy-1] )*0.25; + diff = fmax( diff, fabs(unew[ix*ny+iy] - u[ix*ny+iy])); } } - - //compute reduction - - - float mydiff; - -/** - * IMPLEMENT OMP REDUCE MAX - */ - -#pragma omp parallel shared(nx,ny,u,unew, diff) private (ix,iy,mydiff) - { - mydiff=0.0; -#pragma omp for - for (ix = 1; ix < nx-1; ix++) { - for (iy = 1; iy < ny-1; iy++) { - if (mydiff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) - { - mydiff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] ); - } - } - } - - -#pragma omp critical - { - if (diff < mydiff ) - { - diff = mydiff; - } - } + /* * COPY OLD DATA */ -#pragma omp for +#pragma omp parallel for for (ix = 1; ix < nx-1; ix++) { for (iy = 1; iy < ny-1; iy++) { u[ix*ny+iy] = unew[ix*ny+iy]; } } - } + return diff; } @@ -194,7 +171,7 @@ void inidat(int nx, int ny, float *u, float *unew) #pragma omp for for (ix = 1; ix < nx-1; ix++) for (iy = 1; iy < ny-1; iy++) { - u[ix*ny+iy]=5.0; + u[ix*ny+iy]=0.0; } //boundary left @@ -221,7 +198,7 @@ void inidat(int nx, int ny, float *u, float *unew) //boundary top #pragma omp for for (iy = 0; iy < ny; iy++){ - u[iy]=0.0; + u[iy]=100.0; } diff --git a/lab3/ser_heat2D.c b/lab3/ser_heat2D.c index 3cf927c..0cfd729 100644 --- a/lab3/ser_heat2D.c +++ b/lab3/ser_heat2D.c @@ -31,8 +31,8 @@ Interrior point : #define MM 50 #define ITER_PRINT 100 +#define _MAX_ITER 400 #define PRINT_DATA 1 - #define _EPSILON 0.01 @@ -46,47 +46,48 @@ void prtdat(int nx, int ny, float *u,const char *fnam); int main(int argc, char *argv[]) { int N=NN,M=MM; - float EPSILON=_EPSILON; + int MAX_ITER= _MAX_ITER; - if(argc !=3) - { - fprintf(stderr,"usage %s N EPSILON\n ", argv[0]); - fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n"); - fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n"); - return -1; - } + if(argc !=4) + { + fprintf(stderr,"usage %s N EPSILON MAX_ITER\n ", argv[0]); + fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance, MAX_ITER is max iteration\n"); + fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1, MAX_ITER=1000\n"); + return -1; + } - N = M = atoi(argv[1]); - EPSILON = atof(argv[2]); + N = M = atoi(argv[1]); + EPSILON = atof(argv[2]); + MAX_ITER = atoi (argv[3]); float diff=1.0; float *u = (float *)malloc(N * M * sizeof(float)); float *unew = (float *)malloc(N * M * sizeof(float)); - + if(u==0 || unew ==0) { perror("Can't allocated data\n"); return -1; } - printf ( "\n" ); - printf ( "HEATED_PLATE\n" ); - printf ( " Serial version\n" ); - printf ( " A program to solve for the steady state temperature distribution\n" ); - printf ( " over a rectangular plate.\n" ); - printf ( "\n" ); - printf ( " Spatial grid of %d by %d points.\n", M, N ); - printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON); + printf ( "\n" ); + printf ( "HEATED_PLATE\n" ); + printf ( " Serial version\n" ); + printf ( " A program to solve for the steady state temperature distribution\n" ); + printf ( " over a rectangular plate.\n" ); + printf ( "\n" ); + printf ( " Spatial grid of %d by %d points.\n", M, N ); + printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON); /* Initialize grid and create input file */ printf("Initializing grid\n"); - + inidat(N, M,u,unew); prtdat(N, M,u, "initial.dat"); - + printf("Start computing\n\n"); int iter=0; @@ -95,20 +96,23 @@ int main(int argc, char *argv[]) * iterate until the new solution unew differs from the old solution u * by no more than EPSILON. * */ - - while(diff> EPSILON) { - diff= update(N, M, u, unew); - + while(diff> EPSILON && iter EPSILON) + printf("***Not converged MAX_ITERATION %d reached\n***", MAX_ITER); + prtdat(N, M, u, "final.dat"); - + free(u); free(unew); } @@ -125,12 +129,9 @@ float update(int nx,int ny, float *u, float *unew) for (ix = 1; ix < nx-1; ix++) { 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 - ; - if (diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] )) - { - diff = fabs ( unew[ix*ny+iy] - u[ix*ny+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] )*0.25; + diff = fmax( diff, fabs(unew[ix*ny+iy] - u[ix*ny+iy])); } } @@ -162,10 +163,10 @@ void inidat(int nx, int ny, float *u, float *unew) // interior points for (ix = 1; ix < nx-1; ix++) for (iy = 1; iy < ny-1; iy++) { - u[ix*ny+iy]=5.0; + u[ix*ny+iy]=0.0; } - //boundary left + //boundary left for (ix = 1; ix < nx-1; ix++){ u[ix*ny]=100.0; @@ -178,13 +179,13 @@ void inidat(int nx, int ny, float *u, float *unew) } //boundary down - for (iy = 0; iy < ny; iy++){ + for (iy = 0; iy < ny; iy++){ u[(nx-1)*(ny)+iy]=100.0; } //boundary top - for (iy = 0; iy < ny; iy++){ + for (iy = 0; iy < ny; iy++){ u[iy]=100.0; } @@ -199,26 +200,26 @@ void inidat(int nx, int ny, float *u, float *unew) **************************************************************************/ void prtdat(int nx, int ny, float *u,const char *fname) { - + int ix, iy; FILE *fp; if(ITER_PRINT==0)return; - + fp = fopen(fname, "w"); - // fprintf ( fp, "%d\n", M ); - // fprintf ( fp, "%d\n", N ); + // fprintf ( fp, "%d\n", M ); + // fprintf ( fp, "%d\n", N ); for (ix = 0 ; ix < nx; ix++) { for (iy =0; iy < ny; iy++) { fprintf(fp, "%6.2f ", u[ix*ny+iy]); } - fputc ( '\n', fp); + fputc ( '\n', fp); } - printf (" Data written to the output file %s\n", fname); + printf (" Data written to the output file %s\n", fname); fclose(fp); }