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

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

#define DEBUG 0 // Display all error messages
//#define NX 256 // number of cells in the x-direction
//#define NY 256 // number of cells in the y-direction
#define L 1.0 // domain length
#define W 1.0 // domain width
#define C 1.0 // c, material conductivity. Uniform assumption.
#define TEND 1.0 // tEnd, output time
#define DX (L/256) // dx, cell size
#define DY (W/256) // 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 BLOCK_SIZE_X 16
#define BLOCK_SIZE_Y 16

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


    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) un[i+j*NX] = o + KX*(e-2*o+w) + KY*(n-2*o+s);
}

/***********************************/
/* JACOBI ITERATION FUNCTION - CPU */
/***********************************/
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];
	}
    }
  } 
}


/******************************/
/* TEMPERATURE INITIALIZATION */
/******************************/
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() {

  // --- 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
  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));
