#include #include #ifdef _OPENMP #include #endif /* ================================================================================= The distribution of heat over time is described by the so called heat equation: df/dt = alpha*(d2f/dx2 + d2f/dy2) for a function f(x,y,t). This formula may be discretized in a regular grid G(:,:) by computing the new value G1(x,y) in a point (x,y) at each time step as: G1(x,y) = G(x,y) + Ax * ( G(x+1,y) + G(x-1,y) - 2.0 * G(x,y)) + Ay * ( G(x,y+1) + G(x,y-1) - 2.0 * G(x,y)) For each point in the grid the next value depends on the values of the four up and down, left and right adjacent points. This solution is based on a source code written by Blaise Barney, D. Turner, George L. Gusciora **************************************************************************** Explanation of constants and variables Nrows = x dimension of problem grid Ncols = y dimension of problem grid Tsteps = number of time Tsteps MAXWORKER = maximum number of workers tasks MINWORKER = minimum number of workers tasks IniDataTag, LTAG, RTAG, DONE = message tags NoNeighbor = indicates no neighbor Ax, Ay = used in heat equation Grid = array for grids my_rank,Root = Process ranks numprocs = number of processes --------------------------------------------------------------------------- */ #define Nrows 200 #define Ncols 200 int Nrow0, Nrow1, Ncol0, Ncol1, l_Ncol0, l_Ncol1, g_Ncol0, g_Ncol1, sp; int my_Ncols; float *Grid; #define Ax 0.1 #define Ay 0.1 // to address C arrays in Fortran-like column order #define index2D(a,b,LD1) a + (b*LD1) #define index3D(a,b,c,LD1,LD2) a + (b*LD1) + (c*LD1*LD2) void compute(int y0, int y1, float *g1, float *g2) { /* Compute next step using matrices g1, g2 of dimension (nr,nc) */ int i, j,j0,j1,nr=Nrows; if (y0<=0) j0=1; else j0=y0; if (y1>=Ncols-1) j1=Ncols-2; else j1=y1; for ( j = j0; j < j1; j++ ) { for ( i = 1; i < Nrows-1; i++ ) { g2[index2D(i,j,nr)] = g1[index2D(i,j,nr)] + Ax * ( g1[index2D((i+1),j,nr)] + g1[index2D((i-1),j,nr)] - 2.0 * g1[index2D(i,j,nr)]) + Ay * ( g1[index2D(i,(j+1),nr)] + g1[index2D(i,(j-1),nr)] - 2.0 * g1[index2D(i,j,nr)]); } } } void inidat() { /* Initialize Grid data */ int i,j; for ( i = 0; i < Nrows*Ncols; i++ ) Grid[i] = 0.0; for ( i = 0; i < Nrows; i++ ) { for ( j = 1; j <= 2; j++ ) { Grid[index3D(i,j,0,Nrows,Ncols)] = (float)i*(float)(Nrows-i-1) + (float)j*(float)(Ncols-j-1); } } } int MINVAL(int s, float* a) { int e; float v; v = a[0]; for ( e = 0; e < s; e++ ) { if ( v > a[e] ) v = a[e]; } return(v); } int MAXVAL(int s, float* a) { int e; float v; v = a[0]; for ( e = 0; e < s; e++ ) { if ( v < a[e] ) v = a[e]; } return(v); } void RealData2pgm(int s1, int s2, float *rdata, char* name) { /* Simple subroutine to dump integer data in a PGM format */ int i, j; FILE* ouni; int vp, vs; float rmin, rmax; char fname[80]; /* Write on unit 700 with PGM format */ strcpy(fname,name); strcat(fname,".pgm\0"); ouni = fopen(fname,"w"); if ( ! ouni ) { fprintf(stderr,"!!!! Error write access to file %s\n",fname); } /* Magic code */ fprintf(ouni,"P2\n"); /* Dimensions */ fprintf(ouni,"%d %d\n",s1, s2); /* Maximum value */ fprintf(ouni,"255\n"); /* Values from 0 to 255 */ rmin = MINVAL(s1*s2,rdata); rmax = MAXVAL(s1*s2,rdata); vs = 0; for ( i = 0; i < s1; i++ ) { for ( j = 0; j < s2; j++ ) { vp = (int) ( (rdata[i+(j*s1)] - rmin) * 255.0 / (rmax - rmin) ); vs++; fprintf(ouni,"%4d ", vp); if ( vs >= 10 ) { fprintf(ouni," \n"); vs = 0; } } fprintf(ouni," "); vs = 0; } fclose(ouni); return; } int ProcessBand(int length, int* start, int* end, int intervals, int n) { int delta; delta = ( length + (intervals - 1) ) / intervals; *start = delta * n; *end = delta * ( n + 1) - 1; if ( *end > length-1 ) *end = length-1; return(0); } int main( int argc, char *argv[] ) /* Program Heat2D */ { int Tsteps=1000, MAXWORKER=Ncols, MINWORKER=1, IniDataTag=100, LTAG=200, RTAG=300, DONE=400, NoNeighbor=-1, Root=0, lborder, rborder; int my_rank=0,numprocs=1,avecol,cols,offset,extra, p, dest,source,left,right,msgtype, rc,start,end,i,ix,iy,it,ierr, st; char fname[80]; int l1, l2, ls; double time1, time2; // Initialize grid - root allocates entire Grid Grid = (float*)malloc(sizeof((float)1.0)*Nrows*Ncols*2); if ( Grid == NULL ) { fprintf(stderr," Error allocating Grid(%d,%d,2)",Nrows,Ncols); ierr = -1; return(-1); } // Insert parallel region, pay attention to shared or private variables // This two statements are executed only if openmp is defined //numprocs=omp_get_num_threads(); //my_rank=omp_get_thread_num(); printf("my_rank,numprocs= %d %d \n",my_rank,numprocs); // Decide which process or processes execute the check and the call to inidat and RealData2pgm // Check if numprocs are too many if ( numprocs > Ncols ) { fprintf(stderr," NUMPROCS must not exceed Ncols = %d\n",Ncols); } inidat(); // Dump initial configuration RealData2pgm(Nrows,Ncols,Grid,"tmap0000"); // This statement is executed only if openmp is defined //time1 = omp_get_wtime(); // Calculate portion of columns ProcessBand(Ncols,&Ncol0,&Ncol1,numprocs,my_rank); // // Compute iterations: before calculating communicate borders to neighbors. // It is enough to exchange only vertical borders sp=0; for ( it = 1; it <= Tsteps; it++ ) { // Check for possible syncronization issues // Calculate grid points compute(Ncol0,Ncol1, &Grid[index3D(0,0,sp,Nrows,Ncols)], &Grid[index3D(0,0,(1-sp),Nrows,Ncols)]); sp=1-sp; } // This statement is executed only if openmp is defined //time2 = omp_get_wtime(); // Decide which process or processes execute the following statements fprintf(stdout," Elapsed time with %d processes = %lf\n",numprocs,(time2-time1)); // Print and show results it = Tsteps; sprintf(fname,"tmap%4.4d",it); RealData2pgm(Nrows,Ncols,&Grid[index3D(0,0,sp,Nrows,my_Ncols)],fname); return(0); }