Interleaving of row processing

Signed-off-by: Gaurav Kukreja <mailme.gaurav@gmail.com>

heat-seq/heat.c

@@ -12,7 +12,7 @@
 #include "heat.h"
 #include "timing.h"
 
-#define BLOCKSIZE 100
+#define INTERLEAVING_COUNT 10 
 
 void usage( char *s )
 {
@@ -27,8 +27,7 @@ int main( int argc, char *argv[] )
     FILE *infile, *resfile;
     char *resfilename;
 
-	unsigned BlockSize = BLOCKSIZE;
-
+	int Interleaving_Count;
 	double *tmp;
 	
     // algorithmic parameters
@@ -119,10 +118,12 @@ int main( int argc, char *argv[] )
 
 		case 0: // JACOBI
       
-		    residual = relax_jacobi_return_residual(param.u, param.uhelp, np, np, BlockSize);
-			tmp = param.u;
-			param.u = param.uhelp;
-			param.uhelp = tmp;
+			Interleaving_Count = INTERLEAVING_COUNT;
+
+			residual = relax_jacobi_return_residual(param.u, param.uhelp, np, np, Interleaving_Count);
+			tmp=param.u;
+			param.u=param.uhelp;
+			param.uhelp=tmp;
 		    break;
 
 		case 1: // GAUSS
@@ -132,7 +133,7 @@ int main( int argc, char *argv[] )
 		    break;
 	    }
 	    
-	    iter++;
+	    iter = iter + Interleaving_Count;
 
 	    // solution good enough ?
 	    if (residual < 0.000005) break;

heat-seq/heat.h

@@ -51,18 +51,18 @@ int coarsen(double *uold, unsigned oldx, unsigned oldy ,
 	    double *unew, unsigned newx, unsigned newy );
 
 // Gauss-Seidel: relax_gauss.c
+#if 0
 double residual_gauss( double *u, double *utmp,
 		       unsigned sizex, unsigned sizey );
+#endif
 void relax_gauss( double *u, 
 		  unsigned sizex, unsigned sizey  );
 
 // Jacobi: relax_jacobi.c
-#if 0
 double residual_jacobi( double *u,
 			unsigned sizex, unsigned sizey );
-#endif
 double relax_jacobi_return_residual( double *u, double *utmp,
-		   unsigned sizex, unsigned sizey, unsigned BlockSize ); 
+		   unsigned sizex, unsigned sizey, int Interleaving_Count ); 
 
 
 #endif // JACOBI_H_INCLUDED

heat-seq/misc.c

@@ -92,7 +92,7 @@ int initialize( algoparam_t *param )
 		    (param->heatsrcs[i].range-dist) / 
 		    param->heatsrcs[i].range *
 		    param->heatsrcs[i].temp;
-		(param->uhelp)[ j*np ] = (param->u)[ j*np ];
+			(param->uhelp)[j*np] = (param->u)[j*np];
 	    }
 	}
       
@@ -109,7 +109,7 @@ int initialize( algoparam_t *param )
 		    (param->heatsrcs[i].range-dist) /
 		    param->heatsrcs[i].range *
 		    param->heatsrcs[i].temp;
-		(param->uhelp)[ j*np+(np-1) ] = (param->u)[ j*np+(np-1) ];
+			(param->uhelp)[j*np+(np-1)] = (param->u)[j*np+(np-1)];
 	    }
 	}
     }

heat-seq/relax_jacobi.c

@@ -7,13 +7,13 @@
 
 #include "heat.h"
 
-#if 0
+#if 0			//residual_jacobi
 /*
  * Residual (length of error vector)
  * between current solution and next after a Jacobi step
  */
 double residual_jacobi( double *u,
-			unsigned sizex, unsigned sizey )
+		unsigned sizex, unsigned sizey)
 {
     unsigned i, j;
     double unew, diff, sum=0.0;
@@ -22,46 +22,54 @@ double residual_jacobi( double *u,
     {
 	    for( j=1; j<sizex-1; j++ )
 	    {
-	    unew = 0.25 * (u[ i*sizex     + (j-1) ]+  // left
-			   u[ i*sizex     + (j+1) ]+  // right
-			   u[ (i-1)*sizex + j     ]+  // top
-			   u[ (i+1)*sizex + j     ]); // bottom
+		    unsigned offset = i*sizex + j;
+		    unew = 0.25 * (u[ offset -1 ]+  // left
+				    u[ offset +1 ]+  // right
+				    u[ offset - sizex ]+  // top
+				    u[ offset + sizex ]); // bottom
 
-	    diff = unew - u[i*sizex + j];
+		    diff = unew - u[offset];
 		    sum += diff * diff;
 	    }
     }
 
     return sum;
+
 }
-#endif
+#endif			//residual_jacobi
 
 /*
  * One Jacobi iteration step
  */
 double relax_jacobi_return_residual( double *u, double *utmp,
-		unsigned sizex, unsigned sizey, unsigned BlockSize )
+		unsigned sizex, unsigned sizey, int Interleaving_Count )
 {
-	int i, j, k, l;
+	int i, j, k, l=0;
 	double unew, diff, sum=0.0;
 
-	unsigned BlockCountX = (sizex-2)/BlockSize;
-	unsigned BlockCountY = (sizey-2)/BlockSize;
-
-	for ( k = 0; k < BlockCountY; k++)
-		for ( l = 0; l < BlockCountX; l++)
+	//MODIFIED: exchanged the outer and the inner loop
+	for( i=1; i < sizey-2+Interleaving_Count; i++ )
 	{
-			for( i=1 + BlockSize*k; i <= BlockSize*(k+1); i++ )
-				for( j=1 + BlockSize*l; j <= BlockSize*(l+1); j++ )
+		for(k=i; k > i-Interleaving_Count; k--)
+		{
+			for( j=1; j<sizex-1; j++ )
 			{
-					utmp[i*sizex+j]= 0.25 * (u[ i*sizex     + (j-1) ]+  // left
-							u[ i*sizex     + (j+1) ]+  // right
-							u[ (i-1)*sizex + j     ]+  // top
-							u[ (i+1)*sizex + j     ]); // bottom
+				if(k == 0) break;
+				if(k > sizey-2) continue;
 
-					diff = utmp[i*sizex + j] - u[i*sizex +j];
+				//MODIFIED: storing the offset in a variable
+				u[k*sizex + j] = utmp[k*sizex + j];
+				utmp[k*sizex + j]= 0.25 * (u[ k*sizex + j -1 ]+  // left-old
+						utmp[ k*sizex + j +1 ]+  // right-new
+						u[ k*sizex + j - sizex ]+  // top-old
+						utmp[ k*sizex + j + sizex ]); // bottom-new
+				if(k == i-Interleaving_Count+1)
+				{
+					diff = utmp[k*sizex + j] - u[k*sizex + j];
 					sum += diff * diff;
 				}
 			}
+		}
+	}
 	return sum;
 }