   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++)

/*********************************************/
/* JACOBI ITERATION FUNCTION - GPU - TEXTURE */
/*********************************************/
__global__ void Laplace2d_texture(float * __restrict__ T_new, const bool flag, const int NX, const int NY) {
    
  float o, n, s, e, w;
  //Threads id

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define DEBUG 0 // Display all error messages
#define NX 1024 // number of cells in the x-direction
#define NY 1024 // number of cells in the y-direction
#define L 10.0 // domain length
#define W 10.0 // domain width
#define C 1.0 // c, material conductivity. Uniform assumption.
#define TEND 1.0 // tEnd, output time
#define DX (L/NX) // dx, cell size
#define DY (W/NY) // dy, cell size
#define DT (1/(2*C*(1/DX/DX+1/DY/DY))) // dt, fix time step size
#define KX (C*DT/(DX*DX)) // numerical conductivity
#define KY (C*DT/(DY*DY)) // numerical conductivity
#define NO_STEPS (TEND/DT) // No. of time steps
#define PI 3.1415926535897932f
#define COMPARE 0 // compare to CPU solution

#define BLOCK_SIZE_X 32
#define BLOCK_SIZE_Y 32

#define DIVIDE_INTO(x,y) (((x)+(y)-1)/(y)) // define No. of blocks/warps

