   dim3 dimBlock(BLOCK_SIZE_X, BLOCK_SIZE_Y);
   dim3 dimGrid (DIVIDE_INTO(NX, BLOCK_SIZE_X), DIVIDE_INTO(NY, BLOCK_SIZE_Y));
 
+  // Request computer current time
+  time_t t = clock();
+
   // --- Jacobi iterations on the device - texture case
-  for (int k=0; k<MAX_ITER; k=k+2) {
-    // Update d_T_tex starting from data stored in d_T_old_tex
-    Laplace2d_texture<<<dimGrid, dimBlock>>>(d_T_old_tex,0,NX,NY); cudaDeviceSynchronize();
-    // Update d_T_old_tex starting from data stored in d_T_tex
-    Laplace2d_texture<<<dimGrid, dimBlock>>>(d_T_tex,1,NX,NY); cudaDeviceSynchronize();
+  printf("Using GPU-Texture solver\n");
+  for (int step=0; step < NO_STEPS; step+=2) {
+    if (step%100000==0) printf("Step %d of %d\n",step,(int)NO_STEPS);
+    Laplace2d_texture<<<dimGrid, dimBlock>>>(d_T_old_tex,0); cudaDeviceSynchronize();
+    Laplace2d_texture<<<dimGrid, dimBlock>>>( d_T_tex , 1 ); cudaDeviceSynchronize();
   }	
+    // --- Copy results from device to host
+  cudaMemcpy(h_T_GPU_tex_result,d_T_tex,NX*NY*sizeof(float),cudaMemcpyDeviceToHost);
+
+  // Measure and Report computation time
+  t = clock()-t; printf("Computing time (%f seconds).\n",((float)t)/CLOCKS_PER_SEC);
 
   // --- Unbind textures
   cudaUnbindTexture(tex_T);
   cudaUnbindTexture(tex_T_old);
 
+  /***********************************/
+  /* JACOBI ITERATION FUNCTION - CPU */
+  /***********************************/
+  if (COMPARE) {
+  // Request computer current time
+  t = clock();
+
   // --- Jacobi iterations on the host
-  Laplace2d(h_T,h_T_old,NX,NY,MAX_ITER);
-
-  // --- Copy results from device to host
-  cudaMemcpy(h_T_GPU_tex_result,d_T_tex,NX*NY*sizeof(float),cudaMemcpyDeviceToHost);
-	
+  printf("Using CPU solver\n");
+  for (int step=0; step < NO_STEPS; step+=2) {
+    if (step%100000==0) printf("Step %d of %d\n",step,(int)NO_STEPS);
+    Laplace2d(h_T,h_T_old); 
+    Laplace2d(h_T_old,h_T);
+  }
+
+  // Measure and Report computation time
+  t = clock()-t; printf("Computing time (%f seconds).\n",((float)t)/CLOCKS_PER_SEC);
+  }
+
+  /*******************/
+  /* POST-PROCESSING */
+  /*******************/
+  if (COMPARE) {
   // --- Calculate percentage root mean square error between host and device results
   float sum_tex = 0.f, sum_ref = 0.f;
   for (int j=0; j<NY; j++)

    w = tex2D(tex_T,i-1,j); // node(i-1,j)     s
  }
  // --- Only update "interior" (not boundary) node points
  if (i>0 && i<NX-1 && j>0 && j<NY-1) T_new[i+j*NX] = o + KX*(e-2*o+w) + KY*(n-2*o+s);
}

/***********************************/
/* JACOBI ITERATION FUNCTION - CPU */
/***********************************/
void Laplace2d(float * __restrict T, float * __restrict T_new, const int NX, const int NY, const int MAX_ITER)
{
  for(int iter=0; iter<MAX_ITER; iter=iter+2)
    {
      // --- Only update "interior" (not boundary) node points
      for(int j=1; j<NY-1; j++) 
	for(int i=1; i<NX-1; i++) {
	  float T_E = T[(i+1) + NX*j];
	  float T_W = T[(i-1) + NX*j];
	  float T_N = T[i + NX*(j+1)];
	  float T_S = T[i + NX*(j-1)];
	  T_new[i+NX*j] = 0.25*(T_E + T_W + T_N + T_S);
	}
      
      for(int j=1; j<NY-1; j++) 
	for(int i=1; i<NX-1; i++) {
	  float T_E = T_new[(i+1) + NX*j];
	  float T_W = T_new[(i-1) + NX*j];
	  float T_N = T_new[i + NX*(j+1)];
	  float T_S = T_new[i + NX*(j-1)];
	  T[i+NX*j] = 0.25*(T_E + T_W + T_N + T_S);
	}
    }
}

/******************************/
/* TEMPERATURE INITIALIZATION */
/******************************/
void Initialize(float * __restrict h_T, const int NX, const int NY)
{
    // --- Set left wall to 1
    for(int j=0; j<NY; j++) h_T[j * NX] = 1.0;
}


/********/
/* MAIN */
/********/
int main()
{
  const int NX = 256;		// --- Number of discretization points along the x axis
  const int NY = 256;		// --- Number of discretization points along the y axis
  const int MAX_ITER = 4000;	// --- Number of Jacobi iterations

  // --- CPU temperature distributions
  float *h_T	 = (float *)calloc(NX*NY,sizeof(float));
  float *h_T_old = (float *)calloc(NX*NY,sizeof(float));
  float *h_T_GPU_tex_result= (float *)malloc(NX*NY*sizeof(float));

  // --- set initial condition
  Initialize(h_T,     NX, NY);
  Initialize(h_T_old, NX, NY);

  // --- GPU temperature distribution
  float *d_T;		cudaMalloc((void**)&d_T,	NX*NY*sizeof(float));
  float *d_T_old;	cudaMalloc((void**)&d_T_old,	NX*NY*sizeof(float));
  float *d_T_tex;	cudaMalloc((void**)&d_T_tex,	NX*NY*sizeof(float));

      sum_ref = sum_ref + h_T[j*NX+i]*h_T[j*NX+i];
    }
  printf("Percentage root mean square error texture   = %f\n", 100.*sqrt(sum_tex / sum_ref));
  }
	
  // print result to txt file
  FILE *pFile = fopen("result.txt", "w");
