     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);
+  if (i>0 && i<NX-1 && j>0 && j<NY-1) un[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);
+void Laplace2d(float * __restrict u,float * __restrict un){
+  // Using (i,j) = [i+N*j] indexes
+  int i, j, o, n, s, e, w;
+  for (j = 0; j < NY; j++) {
+    for (i = 0; i < NX; i++) {
+
+        o =  i + NX*j ; // node( j,i )     n
+	n = i+NX*(j+1); // node(j+1,i)     |
+	s = i+NX*(j-1); // node(j-1,i)  w--o--e
+	e = (i+1)+NX*j; // node(j,i+1)     |
+	w = (i-1)+NX*j; // node(j,i-1)     s
+
+	// only update "interior" nodes
+	if(i>0 && i<NX-1 && j>0 && j<NY-1) {
+	  un[o] = u[o] + KX*(u[e]-2*u[o]+u[w]) + KY*(u[n]-2*u[o]+u[s]);
+	} else {
+	  un[o] = u[o];
 	}
-      
-      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;
-}
-
+void Set_IC(float * __restrict u0){
+  int i, j, o, IC; 
+
+  // select IC
+  IC=2;
+
+  switch (IC) {
+  case 1: {
+    for (j = 0; j < NY; j++) {
+      for (i = 0; i < NX; i++) {
+	// set all domain's cells equal to zero
+	o = i+NX*j;  u0[o] = 0.0;
+	// set BCs in the domain 
+	if (j==0)    u0[o] = 0.0; // bottom
+	if (i==0)    u0[o] = 0.0; // left
+	if (j==NY-1) u0[o] = 1.0; // top
+	if (i==NX-1) u0[o] = 1.0; // right
+      }
+    }
+    break;
+  }
+  case 2: {
+    float u_bl = 0.7f;
+    float u_br = 1.0f;
+    float u_tl = 0.7f;
+    float u_tr = 1.0f;
+
+    for (j = 0; j < NY; j++) {
+      for (i = 0; i < NX; i++) {
+	// set all domain's cells equal to zero
+	o = i+NX*j;  u0[o] = 0.0;
+	// set BCs in the domain 
+	if (j==0)    u0[o] = u_bl + (u_br-u_bl)*i/(NX+1); // bottom
+	if (j==NY-1) u0[o] = u_tl + (u_tr-u_tl)*i/(NX+1); // top
+	if (i==0)    u0[o] = u_bl + (u_tl-u_bl)*j/(NY+1); // left
+	if (i==NX-1) u0[o] = u_br + (u_tr-u_br)*j/(NY+1); // right
+      }
+    }
+    break;
+  }
+  case 3: {
+    for (j = 0; j < NY; j++) {
+      for (i = 0; i < NX; i++) {
+	// set all domain's cells equal to zero
+	o = i+NX*j;  u0[o] = 0.0;
+	// set left wall to 1
+	if (i==NX-1) u0[o] = 1.0;
+      }
+    }
+    break;
+  }
+    // here to add another IC
+  }
+}
 
 /********/
 /* 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
+int main() {
+
+  // --- host 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
+  // --- Set initial condition
+  Set_IC(h_T);
+  Set_IC(h_T_old);
+
+  // --- device 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));

/*********************************************/
/* 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

