#include <stdio.h>
#include <time.h>
#if !(defined(_WIN32) || defined(__WIN32__) || defined(WIN32))
#include <sys/time.h>
#endif
#include <vector>
using namespace std;


#if (defined(_WIN32) || defined(__WIN32__) || defined(WIN32))
clock_t global_time;
#else
struct timeval time_start,time_stop;
#endif

//*** IMPORTANT KERNEL STATIC OPTION.
//*** UnComment/Comment the following line
// In the first case (parameter defined) it simply assigns the block variables.
// In the second case it performs fast unit propagation after any assignment
#define SIMPLEBLOCK 0

//**** STATIC SIZE FOR block_vars and delta_vars
//**** ARRAYS USED IN PARALLELIZING SEARCH
#define  BV_SIZE 101 
#define  DV_SIZE 100
#define  FN_SIZE 10000 

//*** Comment/Uncomment the following line for information
//*** on memory allocation during CUDA execution
//#define MYTRACE 0

// The following PARAMETERS decide the parallelism degree
// and the complexity of the formula to be delagetd to CUDA
// We allow two versions, with static (#define) and dynamic values

// *** PARAMETER #1 (argv 1)
// *** Modes for using GPU. Currently 0-3
int USEGPU = 0;   // The default value 

// *** PARAMETER #2 (argv 2)
// NUMBER OF BLOCKS = CUDABLOCKS
// LOG_BLKS is the number of vars fixed by the block address
int LOG_BLKS=2;
int CUDABLOCKS=(1<<LOG_BLKS); 

// *** PARAMETER #3 (argv 3)
// NUMBER OF THREADS per BLOCK  = CUDATHREADS 
// LOG_THRDS is the number of vars fixed by the thread
// In the current GPUs 10 is the maximum value (1024 threads per block)
int LOG_THRDS=5; 
int CUDATHREADS=(1<<LOG_THRDS); 

// *** PARAMETER #4 (argv 4)
// MaxV is the max number of non ground variables delegated to the GPU
int MaxV=50;  

// Delta is the number of variables treated ND by a CUDA thread. 
// The "+1" is needed since delta_vars[0] is not used!
int Delta= MaxV - LOG_BLKS + 1;   

// *** STATIC LIMITS TO FORMULA SIZE (before and after learning)

int MAX_NL   = 8000000;  
int MAX_NC   = 550000;
int nclausel = 0;

// *** PARAMETERS FOR FOR PARALLELIZING UNIT PROPAGATION
//   THREADS  (typically, 256, 512, 1024)
#define THREADS 512
//  Number of blocks used (assigned later using instance and HW data)
int PARblocks=1;          



//*********************************************
// GLOBAL ARRAYS USED
//*********************************************

int * formula;            // Formula and vars in the host    
int * clause_pointer;     // 
int * filtered_formula;   // Reduced formula in the tail 
int * filtered_clause_pointer; // stage
int * dev_formula;        // Formula to be passed in the device
int * dev_cp;             // in the UP 
int * dev_mapped_formula; // Formula in the device
int * dev_mapped_cp;      // in the unit propagation
int * host_vars;          // variables in the host
int * filtered_vars;      // Device vars (tail)
int * dev_vars;           // Device vars (u.p.)
int * dev_parma_vars;     // Device vars (u.p.)

//*********************************************

int * host_strat_vars; // Can be used for strategy 
int * mask_data; //*** unificato per unico trasferimento
int * h_mask_data; 
int * maps_vars;


// *** stack/trail information

int * level;
int * learnt_clause;
int * seen;
int * refs;
int * trail;

int learning      =1;
int level_to_jump = -2; // in caso di backtracking/jumping

// *** GLOBAL VARIABLES

int NV,NC,NL;        // Formula parameters
int CUDA_count=0;    // It counts the calls to a device
int CUDA_count_lower=0;    // It counts the calls to a device (option 1 and 3)
int backtracking=0;  // It counts the backtrackings in the host
int backjump=0;  // It counts the backtrackings in the host
int sat_val,clauind,selected_var; // Again, for backtracking purposes
float deltatime;
cudaEvent_t start, stop;

//*** Error Handling
__global__ void parmask_prop( int* mask_data, int* vars,int* dev_mapped_formula,int* dev_mapped_cp, int g_tn, int dbg);
static void HandleError(cudaError_t err,const char *file,int line){
    if (err != cudaSuccess){
      printf("%s in %s at line %d\n",cudaGetErrorString(err),file,line);
      printf("ncl %d at lev %d\n",nclausel,level[0]);
      for(int i=1;i<trail[0];i++)
           printf("%d@%dr%d,",trail[i],level[abs(trail[i])],refs[abs(trail[i])]);
      printf("\n");
      exit( EXIT_FAILURE );}}
#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))
#define HANDLE_NULL( a ){if (a == NULL) { \
     printf( "Host memory failed in %s at line %d\n", __FILE__, __LINE__ ); \
     exit( EXIT_FAILURE );}}
#define IMIN(a,b) (a<b?a:b)
#define IMAX(a,b) (a>b?a:b)

// *** ALLOCATE / DEALLOCATE

__host__ void allocate_first(){
//*** SIZE INDEPENDENT ON THE FORMULA  
  if (USEGPU>0){
    HANDLE_ERROR( cudaHostAlloc( (void**)&formula,         MAX_NL * sizeof(int), cudaHostAllocDefault ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_mapped_formula, MAX_NL * sizeof(int)));
    HANDLE_ERROR( cudaMalloc((void**)&dev_formula,         MAX_NL * sizeof(int)));
    HANDLE_ERROR( cudaHostAlloc( (void**)&clause_pointer, (MAX_NC+1) * sizeof(int), cudaHostAllocDefault ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_mapped_cp,     (MAX_NC+1) * sizeof(int) ));
    HANDLE_ERROR( cudaMalloc((void**)&dev_cp,             (MAX_NC+1) * sizeof(int) ));
  }
  else{
    formula=(int*)malloc(MAX_NL * sizeof(int));
    clause_pointer=(int*)malloc((MAX_NC+1) * sizeof(int));
  }
  filtered_clause_pointer=(int*) malloc((MAX_NC+1) * sizeof(int));
  filtered_formula       =(int*) malloc(MAX_NL * sizeof(int));
}    
    
__host__ void allocate_second(){
#ifdef MYTRACE
    printf("Memory allocation (in main):\n");
    printf("    (cuda)                  phi:\t%d bytes for %d literals\n", NL * sizeof(int), NL);
    printf("    (cuda)                   cp:\t%d bytes for %d clauses\n", (NC + 1) * sizeof(int), NC);
    printf("    (cuda)            host_vars:\t%d bytes for %d vars\n", NV * sizeof(int), NV);
    printf("    (cuda)       dev_parma_vars:\t%d bytes for %d vars\n", NV * sizeof(int), NV);
    printf(" (host)              h_mask_num:\t%d bytes for %d blocks\n", PARblocks * sizeof(int), PARblocks);
    printf(" (host)               h_mask_id:\t%d bytes for %d blocks\n", PARblocks * sizeof(int), PARblocks);
    fflush(stdout);
#endif
//*** SIZE DEPENDING ON THE FORMULA
    host_strat_vars=(int*)malloc(2*NV * sizeof(int));
    host_vars=(int*)malloc(NV * sizeof(int)); 

    if (USEGPU>0){
      HANDLE_ERROR( cudaHostAlloc( (void**)&h_mask_data, (3* PARblocks+NV) * sizeof(int),  cudaHostAllocDefault ) ); //questo non serve su solo cpu 
      HANDLE_ERROR( cudaMalloc((void**)&dev_vars,              NV * sizeof( int) )) ;
      HANDLE_ERROR( cudaMalloc( (void**)&dev_parma_vars,       NV * sizeof(int) ) );
      HANDLE_ERROR( cudaMalloc( (void**)&mask_data,      (3* PARblocks+NV) * sizeof(int) ) );
    }
    filtered_vars = (int*) malloc(NV * sizeof( int));
    maps_vars     = (int*) malloc(NV * sizeof( int));
    level         = (int*) malloc(NV * sizeof(int));
    refs          = (int*) malloc(NV * sizeof(int));
    trail         = (int*) malloc(NV * sizeof(int));
    seen          = (int*) malloc(NV * sizeof(int));
    learnt_clause = (int*) malloc(NV * sizeof(int)); // AGO" giusto NV?
    // *** Initialization:
    if (USEGPU>0){
    HANDLE_ERROR( cudaMemcpy( dev_mapped_formula, formula,    NL * sizeof(int), cudaMemcpyHostToDevice ) );
    HANDLE_ERROR( cudaMemcpy( dev_mapped_cp, clause_pointer, (NC+1) * sizeof(int), cudaMemcpyHostToDevice ) );
    }
    for(int i=0;i<NV;i++){host_vars[i]=-1; level[i] = -1; refs[i] = -1; seen[i]=0;}
    trail[0] = 1; 
    level[0] = 0;
}

__host__ void deallocate() {      
  free( host_vars );  //AGO: solo queste deallochiamo?
  if (USEGPU>0){
    HANDLE_ERROR( cudaFree( mask_data ));      //
    HANDLE_ERROR( cudaFree(dev_vars));
    HANDLE_ERROR( cudaFree(dev_formula));
    HANDLE_ERROR( cudaFree(dev_cp));
    cudaEventDestroy( start );
    cudaEventDestroy( stop );
  }
}
    

//**************************************************
// INPUT: 
// Data Structure: the formula is stored in two vectors.
// The first stores the literals in sequence. 
// Its length is NL (0..NL-1)
// The second stores the beginning of the various clauses
// The number of clauses is NC. The vector is 0..NC
// The last (redundant) cell stores the value NL
// It also returns NV = 1+number of vars in the formula
// Remark: the formula is read twice 
//**************************************************

__host__ int simple_parsing(char c, FILE *miofile){
  int sum=c-'0';
  char t = fgetc(miofile);
  while (47 <t  && t < 58){
       sum=sum*10+t-48;
       t = fgetc(miofile);
  }
  return sum;   
}

__host__ void leggilinea(FILE *miofile){
  char t = fgetc(miofile);
  while (t != 13 && t != 10) // Look for CR or NL
       t = fgetc(miofile);
}

__host__ void load_formula(char *filename, int *NV, int *NC,int *NL){
  int i=0,j=1,max = 0,flag=0;
  char t;
  FILE *miofile;

  clause_pointer[0]=0;  
  if((miofile = fopen(filename,"r"))==NULL)
     printf("File not found\n");
  else { 
    while(!feof(miofile)){
      t = fgetc(miofile);
      if (t=='c' || t=='p' || t=='%')
         leggilinea(miofile); // Skip comments
      else if (t=='-'){       
         flag=1;
         formula[i] = -simple_parsing(48,miofile);
         max = IMAX(max,-formula[i]); 
         i++;
      }      
      else if ('0' <t  && t <= '9'){
         flag=1;
         formula[i] = simple_parsing(t, miofile);
         max = IMAX(max,-formula[i]); 
         i++;
      } 
      else if (t == '0' && flag==1){
         clause_pointer[j]=i; j++; flag=0;
      }   
    }            
    printf("Num clauses: %d, Num vars: %d, Num literals: %d \n",j-1,max,i);       
    *NV = max+1;   *NC = j-1;    *NL = i;
    fclose(miofile);     
  } // else  
}

//**************************************************
// OUTPUT functions: 
//**************************************************

__host__ void  print_info(cudaDeviceProp prop){
    printf("INFO: multiproc %d\n",prop.multiProcessorCount);
    printf("INFO: sh mem per block %d\n",(int)prop.sharedMemPerBlock);
    printf("INFO: threads per block %d\n",prop.maxThreadsPerBlock);
    printf("INFO: kernelExecTimeoutEnabled %d\n",prop.kernelExecTimeoutEnabled);
    printf("INFO: max tex 1d %d\n",prop.maxTexture1D);
    printf("INFO: max tex 2d %d %d\n",prop.maxTexture2D[0],prop.maxTexture2D[1]);
    printf("INFO: tex align %d\n",(int)prop.textureAlignment);
    printf("INFO: Block Allocation ");
    #ifdef SIMPLEBLOCK
        printf("simple\n");
    #else
        printf("with unit propagation\n");
    #endif
}

__host__ void print_instance( int * phi, int * cp, int NC){
      int i=0,j=1; 
      printf("SAT INSTANCE: \n");
      while(j <= NC){
         printf("%d ",phi[i]);
         i++;
         if(i == cp[j]){printf("\n");j++;}
      }   
      printf("*****************\n"); 
}       

__host__ void print_result(int *vars, int NV){
  printf("\n %d GPU mask prop %d GPU lower tree search, %d backtracks, %d backjumps, %d clauses learnt\n",CUDA_count,CUDA_count_lower,
	 backtracking,backjump,nclausel);
    if (vars[0]){
         for(int i=1;i<NV;i++){
           printf( "V[%d]=%d", i, vars[i]);
           if (i < NV-1) printf(", ");
           if (i % 7 == 0) printf("\n");
        }
        printf("\n");
    }    
    else if  (vars[0]==-1)
         printf( "Computation ended by timeout\n");
    else printf( "No solution\n");
}    

// print_time prints a time expressed in ms 

__host__ void print_time(float time){
    int MIN;
    if (time>=60000){
          MIN = (int)time/60000;
          time = (time - MIN*60000)/1000;
          printf("Execution Time: %dm %.5f sec\n",MIN,time);
    }      
    else printf("Execution Time: %.5f sec\n",time/1000);
}


//*****************************************************
//*** PROPAGATION PHASE
//*****************************************************

// unit propagation is implemented by mask_prop (in one
// of its variants).
// Given a (partial) assignment "vars", for each clause i, 
// mask[i]= 0 => clause i is true 
// mask[i]=-1 => clause i is false 
// mask[i]> 0 => clause i can be satisfied and there are mask[i] 
//               unground variables
// *****************************************************************
// return  0 if all clauses are satisfied
// return -1 if there is a false clause
// return  the smallest index of a clause with fewer unground literals otherwise
// **** Moreover,  *point returns a pointer to ONE of these unassigned literal
// ******************************************************************
// *** There are versions for the device and for the host
// ******************************************************************

__host__ int mask_prop(int* sel_var, int* minmask){
  int i,j,min=0,t,mask=0,pt,bestpt=0,bestcl=0;
  bool initial=true;
  int maybe_sat; // 2 surely sat, 1 maybe, 0 unsat
  
  //***********************************************************************
  //****** SIMPLE CASE (coherent with mask_prop_delta and PARMASKprop)
  //***********************************************************************
  // *** AGO: some little speed-up can be obtained storing pt=t instead of
  // *** pt = j (and changing all tests accordingly - basically removing formula[...]). 
  // *** Left this way to be fully consistent with mask_prop_delta
  //***********************************************************************
  i=0;    
  while((mask >= 0) && (i <NC)){
    mask = 0;   
    j=clause_pointer[i];
    maybe_sat = 0;
    pt = 0; // Greater than number of vars
    while((maybe_sat < 2) && (j<clause_pointer[i+1])) {
      t=formula[j]; // t = literal in position j 
      if ((t<0 && host_vars[-t]==0)||(t>0 && host_vars[t]==1)){
	mask      = 0;
	maybe_sat = 2; 
      }  
      else if (host_vars[abs(t)] < 0) { // Literal Unknown 
	mask++; 
	if (pt == 0 || (pt > 0 && abs(t) < abs(formula[pt]))) 
	  pt=j;  // ==> Here the remark above
	j++;
	maybe_sat = 1;
      } 
      else j++;
    }  // while
    if (j==clause_pointer[i+1] && maybe_sat==0) mask = -1; // EOC: UNSAT!
    if (initial) {
      if (mask != 0)
	{min=mask; bestpt=pt; initial=false; bestcl=i;} 
    }       
    else if (mask!=0){ // not initial
      if  (mask<min || (mask==min && abs(formula[pt]) < abs(formula[bestpt])) ||
	   (mask==min && abs(formula[pt]) == abs(formula[bestpt]) && i < bestcl))
	{min=mask; bestpt=pt; bestcl=i;}
    }
    i++;     
  } // while principale
  //*** OUTPUT 
  *sel_var  =  formula[bestpt];
  *minmask =  min;
  return    bestcl;        
} // mask_prop


// ******************************************************************
// Device version without the array and with careful treatment of vars
// vars[i]<0 means still unlabeled,
// in this case, a=vars[i] is -index of deltavars for variable i

__device__ inline int mask_prop_delta(
         int * formula, int * cp,
         int * block_vars,    
         int NC, int * lit, int *cpidx){
	
    int i=0,t,min=0,mask=0;
    int j;
    int pt, bestpt=0 ;
    int bestcpidx=0;
    bool initial=true;
    int maybe_sat=0; // 2 surely sat, 1 maybe, 0 unsat
    
    while((mask >= 0) && (i <NC)){
  
      mask = 0;   
      j=cp[i];
      maybe_sat = 0;
      pt = 0;
    // if((flag==1)&&(NC>0)){printf("cp[%d]: %d,cp[%d]: %d\n", i, cp[i], i+1,cp[i+1]);}
      while((maybe_sat < 2) && (j<cp[i+1])) {
            t=formula[j];         // Temp variables are used to reduce 
              
              if((t<0 && block_vars[-t]==0)||(t>0 && block_vars[t]==1)){
				mask = 0;
				maybe_sat = 2; 
               
			  }  
			   // Literal Unknown  
			  else if(block_vars[abs(t)]<0){
               mask++; 
               if (pt == 0 || (pt > 0 && abs(t) < abs(formula[pt]))) 
                    pt=j; // CAMBIATO QUI - 22/5
               j++;
               maybe_sat = 1;
             } 
            else j++;
      }  // while
      if (j==cp[i+1] && maybe_sat==0) mask = -1; // EOC: UNSAT!
      if (initial) {
         if (mask != 0)
             {min=mask; bestpt=pt; initial=false; bestcpidx=i; } 
      }       
      else if (mask!=0){ // not initial
           if  ( mask< min || (mask==min && abs(formula[pt]) < abs(formula[bestpt])) ||
				(mask==min && abs(formula[pt])== abs(formula[bestpt]) && i<bestcpidx))
               {min=mask; bestpt=pt; bestcpidx=i;}
           }
      i++;     
   } // while principale 
   *lit = bestpt;
   *cpidx=bestcpidx;
   return min;      
} 

//*********************************************************************
//*********************************************************************
//*********************************************************************
//*********************************************************************


__global__ void parmask_prop( int* mask_data, int* vars,int* dev_mapped_formula,int* dev_mapped_cp, int g_tn){
  int stride = blockDim.x * gridDim.x;
  int g_tid = threadIdx.x + blockIdx.x * blockDim.x;
  int tid = threadIdx.x;    
  __shared__ int s_mask_num[THREADS];
  __shared__ int s_mask_id[THREADS];   
  __shared__ int s_mask_cl[THREADS];   // memorizza la clausola che fa unit prop  

  s_mask_num[tid] = 0;  

  while(g_tid < g_tn){
    int l,v;
    short maybe_sat; // maybe_sat == 2 surely sat, 1 maybe, 0 unsat
    int iter,end;
         
    iter = dev_mapped_cp[g_tid];
    end =  dev_mapped_cp[g_tid + 1];      
	
    maybe_sat = 0;
    int num=0; // portato a registri, cosi' posso tenere s_mask per tutte le clausole
    int id=2147483647; // MAXINT <==========
    while((maybe_sat < 2) && (iter < end)) { // Per ogni lett non sat nella cl
        l = dev_mapped_formula[iter];
        v = vars[abs(l)];      
        iter++;                            // passo al prossimo lett
        if ((( l < 0 ) && ( v == 0 )) ||   // se ha segno concorde
	        (( l > 0 ) && ( v == 1 ))) {   // con la variabile ground
	         num = 0;                      // la clausola è vera
	         maybe_sat = 2;                // sicuramente sat
        }  
        else if (v == -1){ // Literal non-ground
	    // Nota. s_mask_num e s_mask_id sono shared, dunque ad accesso veloce
	         num++;          // +1 al num di lett non-ground in clausola
	         if (abs(l)<abs(id))
	             id = l;      // salvo l'indirizzo del lett 
	         maybe_sat = 1;              // maybe sat
        } 
    }  // while
    if (iter == end && maybe_sat == 0) 
      num = -1; // EndOfClause, clausola falsa.

    // aggiorna analisi clausola nello store shared
    int bestnum=s_mask_num[tid];
    int bestid=s_mask_id[tid];
    int bestcl=s_mask_cl[tid];

    // se nuovo calcolato e' piu' interessante -> memorizza!
    if( num == -1 || bestnum == -1) {
        if (num==-1 && bestnum!=-1){ // la nuova cl e' conflitto
	         s_mask_num[tid] = -1; // copia solo se il precedente non era fallito -> tengo g_tid minimo per cella tid
	         s_mask_cl[tid] = g_tid; // segno la clausola che fa fallire
        }
    }
    else if (bestnum == 0) {
        s_mask_num[tid] = num;
        s_mask_id[tid] = id;
        s_mask_cl[tid] = g_tid;
    }
    else if (num > 0 && 
	     (bestnum > num || 
         (bestnum == num && abs(bestid) > abs(id)) || 
         (bestnum == num && abs(bestid) == abs(id) && bestcl > g_tid)
	      )) {
        s_mask_num[tid] = num;
        s_mask_id[tid] = id;
        s_mask_cl[tid] = g_tid;
    }

    g_tid+=stride;
  } /// while ogni pezzo

  //*************************************************************************************
  __syncthreads();

  // DISTRIBUTED algorithm for implementing first-fail choice (and minimum)      
  short tn = blockDim.x;          
  while(tn > 1){ 
     short half = (tn >> 1); // stesso che tn = tn/2
     short tid_2 = tid + half;
     if(tid < half){
        if( s_mask_num[tid_2] == -1 || s_mask_num[tid] == -1) {
	        if ((s_mask_num[tid_2] == -1 && s_mask_num[tid]>-1 ) || 
                (s_mask_num[tid_2] == -1 && s_mask_cl[tid_2]<s_mask_cl[tid] )){ // copio se il piu' alto e' -1 e scrivo il minimo (a parita' di -1)
	           s_mask_num[tid] = -1;
	           s_mask_cl[tid] = s_mask_cl[tid_2];
	        //printf("riduzione cl %d (%d) num %d (%d)-> da %d a %d\n",s_mask_cl[tid_2],s_mask_cl[tid],s_mask_num[tid_2],s_mask_num[tid],tid_2,tid);
	        }
        }
//*** AGO GIU12: i seguenti due bodies sono uguali. Bug?
    else if (s_mask_num[tid] == 0) {
	    s_mask_num[tid] = s_mask_num[tid_2];
	    s_mask_id[tid] = s_mask_id[tid_2];
	    s_mask_cl[tid] = s_mask_cl[tid_2];
    }
    else if  (s_mask_num[tid_2] > 0 && 
		(s_mask_num[tid] > s_mask_num[tid_2] || 
          (s_mask_num[tid] == s_mask_num[tid_2] && abs(s_mask_id[tid]) > abs(s_mask_id[tid_2]))  || 
          (s_mask_num[tid] == s_mask_num[tid_2] && abs(s_mask_id[tid]) == abs(s_mask_id[tid_2]) && 
           s_mask_cl[tid] > s_mask_cl[tid_2])
		 )) {
	    s_mask_num[tid] = s_mask_num[tid_2];
	    s_mask_id[tid] = s_mask_id[tid_2];
	    s_mask_cl[tid] = s_mask_cl[tid_2];
      }
    }
    
    //*************************************************************************************
    __syncthreads();
    tn = half;
  } // while (tn>1)              

  if(tid == 0){      
    mask_data[3*blockIdx.x] = s_mask_num[0];
    mask_data[3*blockIdx.x+1] = s_mask_id[0];
    mask_data[3*blockIdx.x+2] = s_mask_cl[0];
  }
      
} // parmask_prop


//****************************************************************
//*** FORMULA SIMPLIFICATION (ATL-VEN July 12, 2011)
// simplifies a (still satisfiable) formula by 
// removing true clauses and false literals
//*****************************************************************

__host__ int filter_formula(){ 

 int i,lp,cp,flp,fv,fnc,vera;
 
 filtered_vars[0]=-1;
 maps_vars[0]=0;
 
// Here only the unassigned (free) vars are considered. 
// maps_vars stores the correspondence new var / old var
 
 fv = 0;
 for(i=1;i<NV;i++){
    if (host_vars[i] == -1){
      fv++;
      maps_vars[i] = fv; // OLD NONGROUND VARS ARE ASSIGNED TO NEW
      filtered_vars[fv] = -1;
    }  
    // else maps_vars[i] = i; //(this value would be never used)
 }                         
 
// The formula is simplified and projected
// on the new variables.
     
 lp=0; fnc=0; 
 filtered_clause_pointer[0] = 0;
  
 for (cp=0; cp < NC; cp++){
    vera=1;
    flp = filtered_clause_pointer[fnc];
    while (lp < clause_pointer[cp+1]){
       // FIRST CASE: the literal, hence the clause C_cp is satisfied. Remove C_cp
       if ((formula[lp] > 0) && (host_vars[formula[lp]] == 1) ||
           (formula[lp] < 0) && (host_vars[-formula[lp]] == 0)){
          vera = 1;
          lp = clause_pointer[cp+1]; break; // Exit the loop
       }  // SECOND CASE: the literal is not ground. Copy it, renaming the var
       else if (host_vars[abs(formula[lp])] == -1){
          filtered_formula[flp] =  maps_vars[abs(formula[lp])];  
          if (formula[lp] <= 0)  filtered_formula[flp]=-filtered_formula[flp];
          flp++;
          lp++;
          vera=0;
       }  //THIRD CASE: the literal is false. Just skip it            
       else lp++;
   } // while
   if (!vera) { // NOTA: se e' falsa, la tolgo come se fosse vera
       fnc++;   // Non dovrebbe mai essere chiamata in tal caso.
       filtered_clause_pointer[fnc] = flp;
   }                       
  } // for  
  
  return fnc;
  
} // function filter_formula


//***************************************************************
//***** KERNEL FUNCTION.
//***** DEVICE VERSION OF DPLL
//***** Without recursion
//***** Use block/thread address to guide the search
//***** Use shared block memories
//***************************************************************

// __shared__ short int* block_varsE;

//**************************learning su GPU****************************
   
__device__ int clause_learning_gpu(
		       int * dev_level, int * dev_refs, int * dev_trail, int* dev_seen,
		       int* formula, int* new_formula, int* cp, int* new_cp, int old_NC, int* NC, int* new_NC, int NV) {
		       	
	//imposto il numero di clausole
	int NCP=*new_NC;

	
	int learnt_clause[50];
	
	// sets jump and level_to_jump
	int level_to_jump;

	// inizializzazioni
	int i=0, c=0, p=0, q, var_p=0, var_q, clause, btlevel=0, tptr=dev_trail[0];
	int lp=0;

	for (int i=0;i<NV;i++){
	   dev_seen[i]=0;  
	}

	// *** Trovo first-UIP e creo la clausola
	do {
		if (p==0) {				// inizio con la conflict clause 
			clause = dev_refs[0];	// conflict clause
			
		}else{
			clause = dev_refs[var_p]; // reason clause
				
		}	
		//printf("refs[%d]: %d\n", var_p, dev_refs[var_p]);	
			
		if(dev_refs[var_p]>=old_NC){
			clause=clause-old_NC;
			
			//printf("cp[%d]:%d\n", clause,cp[clause]);	
			//printf("cp[%d+1]:%d\n", clause+1,cp[clause+1]);
			for( i=new_cp[clause];i<new_cp[clause+1];i++){
		
				//controllo il valore di refs per vedere se devo lavorare sulla formula nuova 
				//printf("i:%d\n",i);
				
				var_q=abs(new_formula[i]);
				q=-new_formula[i];
				
				//printf("q:%d\n",q);
				
				if((p==0)||(var_q != var_p)){ // se sono al primo giro (p==0) devo processare tutta la clausola, altrimenti skippo la unit propagation
				
					//printf("dev_level[0] %d --- dev_level[%d]: %d \n",dev_level[0], var_q, dev_level[var_q]);
				
					if(!dev_seen[var_q]){
						dev_seen[var_q]=1;
							
						if(dev_level[var_q]==dev_level[0]){ // (q == decision level corrente)
							c++;
							//printf("c:%d\n",c);
						}
						else 
							if(dev_level[var_q]>0){   // se q è di un decision level inferiore, ma non il root level (al momento il root level è inutilizzato)
												
								learnt_clause[lp++]= -q;
									
								// calcolo level to backtrack, il massimo della learnt_clause, escluso il first-uip che verrà aggiunto alla fine
								btlevel=max(btlevel,dev_level[var_q]);
							}
					}
				}
			}//chiudo il for
		}//ref[var_p]>old_NC
		else{
			for( i=cp[clause];i<cp[clause+1];i++){
		
				//controllo il valore di refs per vedere se devo lavorare sulla formula nuova 
				//printf("i:%d\n",i);
				
				var_q=abs(formula[i]);
				q=-formula[i];
			
				//printf("q:%d\n",q);
				
				if((p==0)||(var_q != var_p)){ // se sono al primo giro (p==0) devo processare tutta la clausola, altrimenti skippo la unit propagation
				
					//printf("dev_level[0] %d --- dev_level[%d]: %d \n",dev_level[0], var_q, dev_level[var_q]);
				
					if(!dev_seen[var_q]){
						dev_seen[var_q]=1;
							
						if(dev_level[var_q]==dev_level[0]){ // (q == decision level corrente)
							c++;
							//printf("c:%d\n",c);
						}
						else 
							if(dev_level[var_q]>0){   // se q è di un decision level inferiore, ma non il root level (al momento il root level è inutilizzato)				
								learnt_clause[lp++]= -q;
						
								// calcolo level to backtrack, il massimo della learnt_clause, escluso il first-uip che verrà aggiunto alla fine
								btlevel=max(btlevel,dev_level[var_q]);
							}
					}
				}
			}//chiudo il for
		}//chiudo else ref[var_p]<old_NC
	
		do{
			p=dev_trail[--tptr];
			var_p=abs(p);
		}
		while(!dev_seen[var_p]);
		c--;
	}while(c>0);
	
	
	  // Blocca i thread del blocco
    // finchè non tutti hanno terminato l'operazione di scrittura
    
    __syncthreads();
	
	//ALE: 8 giu commentata per avere backjumping
	//btlevel = max(btlevel,dev_level[abs(p)]-1); // in caso di solo UIP senza reasons, comunque devo annullare la decisione del livello corrente (dato che poi sommo 1, qui lo sottraggo)
	
	//in p ora ho il first UIP, ci va il suo negato nella learnt clause
		learnt_clause[lp++]=-p;
	
	//condizioni per memorizzare solo una determinata categoria di clausole
	//if (lp <16 ){
	
		for(i=0;i<lp;i++){
			new_formula[new_cp[NCP]+i]=learnt_clause[i];
			//printf("new_formula[new_cp[%d]+%i]: %d\n", NCP, i, new_formula[new_cp[NCP]+i]); 
		}
		new_cp[NCP+1]=new_cp[NCP]+lp;
		//printf("new_cp[%d]: %d,new_cp[%d]: %d\n", NCP+1, new_cp[NCP+1], NCP,new_cp[NCP]); 
		NCP++;
		*new_NC=NCP;
		*NC++;
	
		// se ho abbastanza memoria... (non è detto che debba fare forgetting!)
		// mi arriva:	learnt_clause[] di interi con i numeri delle host_vars della clausola
		//	 lp = intero con la lunghezza della clausola
		// ..e la aggiunge alla formula
		// - aggiorno il clause pointer finale	  
		// nclausel++;
		
		//aggiorno il livello di salto per il backjumping
		level_to_jump=btlevel;	// non chronological backtracking
		
	
/*	}else{ //aggiorno il livello di salto nel caso non abbia imparato la clausola
	
	    level_to_jump=dev_level[0];	// chronological backtracking
	    
	}*/
	
	//printf("level_to_jump %d\n",level_to_jump);
	return level_to_jump;


}




//***************************************************************
//***** KERNEL FUNCTION.
//***** DEVICE VERSION OF DPLL
//***** Without recursion
//***** Use block/thread address to guide the search
//***** Use shared block memories
//***************************************************************

// __shared__ short int* block_varsE;

__global__ void CUDADPLL(int* formula, int* cp,  int *vars, int NC, int NV, 
              int CUDATHREADS, int LOG_BLKS, int Delta, int LOG_THRDS){

    
  //*** If some other thread already found solution, skip the thread

  //if (threadIdx.x==0)  printf("cuda %d-%d ",blockIdx.x,threadIdx.x);
if (vars[0]<1){
	// *** SHARED ARRAY OF THE BLOCK
	__shared__ int block_vars[BV_SIZE];
	// extern
	//__shared__ 
	int new_formula[FN_SIZE]; 
	//__shared__ 
	int new_NC;
   
	//__shared__ 
	int new_cp[1000];
   
   
	new_NC=0;
	
	int old_NC=NC;
    
	int dev_trail[DV_SIZE];
	int dev_seen[DV_SIZE];
   
	int dev_refs[DV_SIZE];
	int dev_level[DV_SIZE];
   
	for(int i=0;i<NV;i++){dev_level[i] = -1; dev_refs[i] = -1; dev_seen[i]=0;}
	dev_trail[0]=1;
	dev_level[0]=0;
          
	int addr=blockIdx.x;
	int addrTh=threadIdx.x;
	int pos,i;
	
	int lit, new_lit, old_lit;
	int cpidx, new_cpidx, old_cpidx;
	bool end=false;
	int min_mask=0;
	int new_min_mask=0;
	int old_min_mask=0;
	int level_to_jump=0;
	int isbackjumping=0;
   
	//****************************************************************
	// The first LOG_BLKS vars are assigned using block coordinates
	// The others are delegated to delta_vars
	//****************************************************************


	// *************************************************************************
	// **** SIMPLE CODE. GOOD FOR DATA PARALLELISM. SOME DETERMINISM IS
	// **** NOT WELL-EXPLOITED
	// *************************************************************************
 	
	block_vars[0]=0;
		for(i=1;i<NV;i++){
			if (i <= LOG_BLKS){ 
				block_vars[i] =  addr % 2; 
				addr=addr/2;
			}         
			else {
				if (i <= LOG_THRDS+LOG_BLKS){
					int j=i-LOG_BLKS-1;
					//printf("j:%d\n",j);
					if(addrTh%2==0){
						new_formula[j] = i;
						
					}else{ 
						new_formula[j] = -i;
						
					}
					new_cp[new_NC]=j+1;
					//printf("new_cp[%d]:%d\n",new_NC,new_cp[new_NC]);
					
					//printf("new_formula[%d]:%d\n",j,new_formula[j]); 
					addrTh=addrTh/2; 
					block_vars[i] = -1;
				} 
				else{   
					block_vars[i] = -1;
				} 
			}      
		} // for i ...
	
	

	// Use the (block,thread) coordinate to guide the successive 
	// LOG_THRDS ND choices
	// top is now a pointer to the top of the stack 
	// lit points to the selected unknown literal
	// pos is its sign


   
       
   do{
	isbackjumping=0;
       
		//*****************************************************
		//*** Partial substitution evaluation here:
		//*****************************************************       
		int flag=1;
		new_min_mask = mask_prop_delta(new_formula, new_cp, block_vars,new_NC,&new_lit, &new_cpidx);
		//printf("new_min_mask: %d, new_formula[new_lit]: %d, new_cpidx:%d \n", new_min_mask, new_formula[new_lit], new_cpidx);
		flag=0;
		old_min_mask = mask_prop_delta(formula, cp, block_vars,old_NC,&old_lit, &old_cpidx);
    
		//printf("min_mask: %d, formula[lit]: %d, cpidx:%d \n", old_min_mask,formula[old_lit], old_cpidx);
		
		
		
		
		
	if((old_min_mask==-1)||(new_min_mask==-1)){
		if(old_min_mask==-1){
			min_mask=-1;
			flag=0;
			lit=old_lit;
			cpidx=old_cpidx;
		}else{
			min_mask=new_min_mask;
			flag=1;
			lit=new_lit;
			cpidx=new_cpidx+NC;
		}
	}else{ 
		if((old_min_mask==0)&&(new_min_mask==0)){
			min_mask=0;
		}else{
			if((old_min_mask==0)&&(new_min_mask > 0)){
				min_mask=new_min_mask;
				flag=1;
				lit=new_lit;
				cpidx=new_cpidx+NC;
			}else{ 
				if((old_min_mask>0)&&(new_min_mask == 0)){
					min_mask=old_min_mask;
					flag=0;
					lit=old_lit;
					cpidx=old_cpidx;
				}else{
					if(old_min_mask < new_min_mask){
						min_mask=old_min_mask;
						flag=0;
						lit=old_lit;
						cpidx=old_cpidx;
					}else{
						min_mask=new_min_mask;
						flag=1;
						lit=new_lit;
						cpidx=new_cpidx+NC;
					}
				}
			}
		}  	
	}//chiudo il controllo per scegliere min_mask
		
		
		
		
		
       
		if (min_mask == 0) {
			end=true; break;  // logically superfluous "break", but it speeds it up...
		}else 
			if (min_mask > 0) { // min_mask > 0
				if(flag==0){
					pos = (formula[lit] > 0); // sign of the unknown literal
				}else{
					pos = (new_formula[lit] > 0); // sign of the unknown literal
				}
				if(min_mask == 1){ // There is ONE non-ground literal: determinism 
					
					int VAR;
					if(flag==0){
						VAR=abs(formula[lit]);
						dev_trail[dev_trail[0]] = formula[lit];
					}else{
						VAR=abs(new_formula[lit]);
						dev_trail[dev_trail[0]] = new_formula[lit];
					}
					
					block_vars[VAR]=pos;
					dev_level[VAR] = dev_level[0];
					
					dev_refs[VAR] = cpidx; // devo sapere la clausola	
		 			dev_trail[0]++;	
		 			 
				}//chiudo if (min_mask==1)
				else {
					//if (min_mask>=2){
					int selected;
					dev_level[0]++; // apro nuovo ramo
					
					if(flag==0){
						selected=abs(formula[lit]);
						dev_trail[dev_trail[0]]=formula[lit];
						
					}else{
						selected=abs(new_formula[lit]);
						dev_trail[dev_trail[0]]=new_formula[lit];
					}
	
					block_vars[selected]=pos;
					dev_level[selected] = dev_level[0];		
					
					dev_trail[0]++;
					
				} //else(min_mask>1)
		
			} // if  min_mask > 0  
      
			else {   // if (min_mask < 0) *** Failure: Backtracking
				
				//********************************************************************
				dev_refs[0]=cpidx;
				
				//faccio il learning
				level_to_jump= clause_learning_gpu( dev_level, dev_refs, dev_trail, dev_seen, formula, new_formula, cp, new_cp, old_NC, &NC, &new_NC, NV);
				int loop=dev_trail[0]>0;
				
				while (loop){ // questo mette a posto trail e variabili. level sistemato sulla ricorsione
					int VARP=abs(dev_trail[dev_trail[0]-1]);
					//printf("VARP: %d\n",VARP);
					loop=dev_trail[0]>0 && ((dev_level[VARP]>level_to_jump)   ||  (dev_refs[VARP]>=0)); //<------------------ da casini
				
					if ((dev_level[VARP]>level_to_jump) || (dev_refs[VARP]>=0)){
						
						block_vars[VARP]=-1;
						dev_level[VARP] = -1;
						//printf("level[%d]:%d\n",VARP,dev_level[VARP]);
						dev_refs[VARP] = -1;
						dev_trail[0]--;
					}
					else{ // refs<0 && level=leveltojump
						if (dev_level[VARP]<=level_to_jump && dev_refs[VARP]==-1){ // variabile gia' testata -> metto caso opposto
							
							dev_level[0]=dev_level[VARP];
						}
					}
				} //chiudo while
		  		// caso speciale se backjump annulla anche primo livello
				if (level_to_jump==0){
					dev_level[0]=0;
				}
				
			isbackjumping=1;
			
				if((dev_level[0]<0)&& (min_mask<0)){
					vars[0]=0;
				}
			
       //********************************************************************

      }//chiudo l'else(min_mask<0)
    } while( isbackjumping || (dev_level[0]>0 && !end));

//*************************************************************************************
	__syncthreads();

	//printf("exit %d-%d-%d ",blockIdx.x,threadIdx.x,min_mask);
	
	if (vars[0]<1 && (min_mask==0)){ // *** A solution: export it on vars
		//printf("yes\n");
		vars[0]=1;
		for(pos=1; pos<NV; pos++) 
				vars[pos] = block_vars[pos];
	} 
 } // if (vars[0]...
} 


//**************************************************************************
//*********    CUDA_caller                         *************************
//**************************************************************************

__host__ int CUDA_caller(int FV){

   int unk_clauses;
   
   CUDA_count_lower++;                       
    
   // *** filter_formula call - global variables are used
   
    unk_clauses= filter_formula();
    
   //    printf("Called GPU with %d free vars\n",FV);
   //    printf("Before: V=%d, C=%d, L=%d, After: V=%d, C=%d, L=%d\n",
   //           NV-1,NC,clause_pointer[NC],FV,unk_clauses,filtered_clause_pointer[unk_clauses]);

//   DEBUG PRINT
/*
    printf("Vars Mapping:\n");
    for(int i=1; i < NV; i++) if(maps_vars[i] < i)
             printf("var[%d] -> %d \t",maps_vars[i],i);   
    printf("\n");
    print_instance(filtered_formula,filtered_clause_pointer,unk_clauses);  
    
    for(int i=1; i < NV; i++) 
             printf("%d -> var[%d]\t",host_vars[i],i);   
    printf("\n");
*/

    //printf("unk_clau %d nl %d, fv %d\n",unk_clauses,filtered_clause_pointer[unk_clauses],FV);
    //print_instance(filtered_formula,filtered_clause_pointer,unk_clauses);  
        
    // copy the value of filtered vars, formula, and clause pointers 
    // to device global variables

    HANDLE_ERROR(cudaMemcpy(dev_vars, filtered_vars, (1+FV)*sizeof( int), cudaMemcpyHostToDevice ));
    HANDLE_ERROR(cudaMemcpy(dev_formula, filtered_formula, filtered_clause_pointer[unk_clauses]*sizeof( int), cudaMemcpyHostToDevice ));
    HANDLE_ERROR(cudaMemcpy(dev_cp, filtered_clause_pointer, (unk_clauses+1)*sizeof(int), cudaMemcpyHostToDevice ));
    
#ifdef MYTRACE
  printf("\nMore memory allocation (in CUDA_caller):\n");
  printf("    (cuda)             dev_vars:\t%d bytes for %d vars\n", (1+FV)* sizeof(int), FV);
  printf("    (cuda)          dev_formula:\t%d bytes for %d vars\n", filtered_clause_pointer[unk_clauses]* sizeof(int), filtered_clause_pointer[unk_clauses]);
  printf("    (cuda)               dev_cp:\t%d bytes for %d clauses\n", (unk_clauses+1) * sizeof(int), (unk_clauses+1));
  printf(" (host)           filtered_vars:\t%d bytes for %d vars\n", (1+FV) * sizeof(int), FV);
  printf(" (host)               maps_vars:\t%d bytes for %d vars\n", NV * sizeof(int), NV);
  printf(" (host)        filtered_formula:\t%d bytes for %d vars\n", clause_pointer[NC] * sizeof(int), clause_pointer[NC]);
  printf(" (host) filtered_clause_pointer:\t%d bytes for %d clauses\n", (NC+1) * sizeof(int), NC);
  printf("CUDA_caller calling CUDADPPL3: blocks = %d  threads per block  = %d\n",CUDABLOCKS,CUDATHREADS);
  printf("          (shared) block_vars:\t%d bytes for %d vars. Total of %d bytes per block\n", MaxV * sizeof(int), MaxV, MaxV * sizeof(int));
  printf("          (thread) delta_vars:\t%d bytes for %d vars. Total of %d bytes per block\n", Delta * sizeof(int), Delta, CUDATHREADS*Delta * sizeof(int));
  printf("            (thread) varstack:\t%d bytes for %d vals. Total of %d bytes per block\n", 2*Delta * sizeof(int), 2*Delta, CUDATHREADS*2*Delta * sizeof(int));
  fflush(stdout);
#endif

    // ***************   DEVICE COMPUTATION: *************************
    //printf("Delta %d (Check if it is less than  %d)\n",Delta,DV_SIZE);
    CUDADPLL<<<CUDABLOCKS,CUDATHREADS>>>(dev_formula,dev_cp,dev_vars,unk_clauses,1+FV,CUDATHREADS,LOG_BLKS,Delta,LOG_THRDS); 

    // ***************   Results' analysis   *************************
    // Only the first value (a flag) is copied back
    HANDLE_ERROR(cudaMemcpy(host_vars,dev_vars, sizeof( int), cudaMemcpyDeviceToHost));
    //  If there is no solution, skip, otherwise copy the whole
    //  assignment back from DEVICE to HOST
    //  Then assembly the solution from vars and filtered_vars
    
    if (host_vars[0]==1){ 
      HANDLE_ERROR(cudaMemcpy(filtered_vars,dev_vars, (1+FV)*sizeof( int), cudaMemcpyDeviceToHost));
       for(int i = 1; i < NV; i++)  
          if (host_vars[i] < 0)          
             host_vars[i] = filtered_vars[maps_vars[i]];              
    }  
    
   
   // In realta' le prox quattro potrebbero anche essere globali o statiche,
   // comunque grandi (stanno su host) e si  evita  di riallocare
   
    return (host_vars[0]==1)*2; // 2 if good, 0 if bad
}


// **************************************************************************************************
// **************************************************************************************************
// *******************     clause learning **********************************************************
// **************************************************************************************************
// **************************************************************************************************

__host__ void clause_learning() {
  // sets jump and level_to_jump

	// inizializzazioni
	int i, c=0, p=0, q, var_p=0, var_q, clause, lp=0, btlevel=0, tptr=trail[0],
		verbose=0;

	if (verbose>1) {
		printf("Current level:\t\t%d\n",level[0]);
		printf("Conflict clause:\t%d\n",refs[0]);
		int i;
		printf("trail[%d] = [",trail[0]);
		for(i=1;i<trail[0]-1;i++)
			printf("%d@%d,",trail[i],level[abs(trail[i])]);
		printf("%d@%d]\n",trail[i],level[abs(trail[i])]);
	}
	//for(int u=0;u<NV;u++)
	//	if(host_vars[u]!=-1)
	//		printf("refs[%d] = %d\n",u,refs[u]);

	// *** Trovo first-UIP e creo la clausola

	//printf("\ntrail:\n");
	//for(int b=1;b<trail[0];b++)
	//	printf("%-5d @ %-5d r %-5d\n",trail[b],level[abs(trail[b])],refs[abs(trail[b])]);
	//printf("\n\n");

	memset(seen,0,NV*sizeof(int));
	do {
	  if (p==0) {				// inizio con la conflict clause 
	    clause = refs[0];	// conflict clause
	  } else
	    clause = refs[var_p];	// reason clause
	  if (verbose>1) {
	    printf("p:\t%d\t",p);
	    printf("var_p:\t%d\t",var_p);
	    printf("clause:\t%d\t",clause);
	    printf("c_start:\t%d\t",clause_pointer[clause]);
	    printf("end: %d\n",clause_pointer[clause+1]);
	    printf("reason (clause %d) is:\t",clause);
	    for(i=clause_pointer[clause];i<clause_pointer[clause+1];i++) {
	      printf("%d\t",formula[i]);
	    }
	    printf("\n");
	  }
	  //printf("clause %-6d cp[clause] %-10 clause_pointer[clause+1] %-10d var_p %-6d refs[var_p] %-10d\n",clause,clause_pointer[clause],clause_pointer[clause+1],var_p,refs[var_p]);
	  for(i=clause_pointer[clause];i<clause_pointer[clause+1];i++) {
	    var_q = abs(formula[i]);
	    q = -formula[i];
	    if ((p==0) || (var_q!=var_p)) { // se sono al primo giro (p==0) devo processare tutta la clausola, altrimenti skippo la unit propagation
	      if (verbose>1) {
		printf("q:%d\n",q);
		printf("lp at %d, q is %d at level %d, ref %d\n",lp,q,level[var_q],refs[var_q]);
	      }
	      if (!seen[var_q]) {
		if (verbose>1) printf("q never seen! ");
		seen[var_q] = 1;
		if (level[var_q] == level[0]) { // (q == decision level corrente)
		  c++;
		  if (verbose>1) printf("and of same current level (%d==%d) so c++.\n",level[var_q],level[0]);
		}
		else if (level[var_q] > 0) { // se q è di un decision level inferiore, ma non il root level (al momento il root level è inutilizzato)

		  learnt_clause[lp++] = -q;		    

		  // calcolo level to backtrack, il massimo della learnt_clause, escluso il first-uip che verrà aggiunto alla fine
		  btlevel = max(btlevel,level[var_q]);
		  if (verbose>1) {
		    printf("and of minor level (%d<%d) so add %d.\n",level[var_q],level[0],-trail[i]);
		    if (lp>0) {
		      printf("current learnt clause is: ");
		      int j;
		      for(j=0;j<lp-1;j++)
			printf("%d,",learnt_clause[j]);
		      printf("%d.\n",learnt_clause[j]);
		    }
		  }
		}
	      } else if (verbose>1) printf("q is seen\n");
	    }
	  }
	  do {
	    p = trail[--tptr];
	    var_p = abs(p);
	    if (verbose>1) {
	      printf("c:%d p%d\n",c,p);
	    }
	    if (verbose>1) {
	      printf("scan p=%d\ttptr=%d\tseen[var_p]=%d\n",p,tptr,seen[var_p]);
	    }	
	  }
	  while (!seen[var_p]);
	  c--;
	  //seen[var_p] = 0;
	  if (verbose>1) printf("c--, process var %d, c:%d, ref %d val %d\n",var_p,c,refs[var_p],host_vars[var_p]);	
	}
	while (c>0);

	//ALE: 8 giu commentata per avere backjumping
	//btlevel = max(btlevel,level[abs(p)]-1); // in caso di solo UIP senza reasons, comunque devo annullare la decisione del livello corrente (dato che poi sommo 1, qui lo sottraggo)

	// in p ora ho il first-UIP, ci va il suo negato nella learnt clause
	learnt_clause[lp++] = -p;

	// controllo inserimento clausole doppie: todo

	// debug
	if (verbose>0) {
	  printf("Add First-UIP: %d lev %d\nLearn clause %d: ",-p,level[abs(p)],NC);
	  //		printf("* LEARNT CLAUSE: ");
		for(i=0;i<lp-1;i++) {
			//printf("%d ",learnt_clause[i]); // stampa per .cnf
		  //			printf("%d@%d, ",learnt_clause[i],level[abs(learnt_clause[i])]);
		  printf("%d ",learnt_clause[i]);
		}
		//printf("%d 0\n",learnt_clause[i]);// stampa per .cnf
		//		printf("%d@%d.\n",learnt_clause[i],level[abs(learnt_clause[i])]);
		printf("%d 0\n",learnt_clause[i]);
	}
	
	// *** aggiungo la clausola imparata
	  if (//lp <10 &&
	    ((NL+lp)<=MAX_NL) && (NC+1<=MAX_NC)) { // se ho abbastanza memoria... (non è detto che debba fare forgetting!)
	  // mi arriva:	learnt_clause[] di interi con i numeri delle host_vars della clausola
	  //	 lp = intero con la lunghezza della clausola
	  // ..e la aggiunge alla formula
	  // - aggiorno il clause pointer finale	  
	  nclausel++;

	  clause_pointer[NC+1] = clause_pointer[NC]+lp;
	  // - aggiungo la clausola alla formula
	  for(i=0;i<lp;i++)
	    formula[clause_pointer[NC]+i] = learnt_clause[i];
	  // - aggiorno il numero dei letterali e il numero di clausole

	  if (USEGPU%2==1){ // solo se uso gpu
	    HANDLE_ERROR( cudaMemcpy( dev_mapped_formula+NL, formula+NL, lp * sizeof(int), cudaMemcpyHostToDevice ) );
	    HANDLE_ERROR( cudaMemcpy( dev_mapped_cp+NC+1, clause_pointer+NC+1, 1 * sizeof(int), cudaMemcpyHostToDevice ) );
	  }

	  NL += lp;
	  NC += 1;

	  if (verbose>1) printf("backjump from level	%d to level %d\n\n\n",level[0],btlevel);
	  level_to_jump=btlevel;	// non chronological backtracking
	  }
	  else {
	    if (verbose>1) printf("backtrack to level %d\n",level[0]);
	    level_to_jump=level[0];	// chronological backtracking
	  }
}


//**************************************************************************
//**************************************************************************
//****** AUXILIARY
//**************************************************************************
//**************************************************************************

__host__ inline void backtrackingf(int level_to_jump){

  // annulla trail stack e riferimenti a variabili
  //  printf("backtrack to %d\n",level_to_jump);

  backtracking++;
  int loop=trail[0]>0;
  while (loop){ // questo mette a posto trail e variabili. level sistemato sulla ricorsione
    int VARP=abs(trail[trail[0]-1]);
    loop=trail[0]>0 &&       
      ((level[VARP]>level_to_jump)   ||  (refs[VARP]>=0)  ||
       (level[VARP]<=level_to_jump && refs[VARP]==-2));
     
    if (level[VARP]>level_to_jump || refs[VARP]>=0 || refs[VARP]==-2){
      //printf("back var %d lev %d\n",VARP,level[VARP]);
      host_vars[VARP] = -1;		// restore "unknown" status
      level[VARP] = -1;
      refs[VARP] = -1;
      trail[0]--;
    }
    else{
      if (level[VARP]<=level_to_jump && refs[VARP]==-1){ // variabile gia' testata -> metto caso opposto
	//printf("switch var %d at lev %d -> %d\n",VARP,level[VARP],-trail[trail[0]-1]);
	refs[VARP] = -2;
	host_vars[VARP] = 1-host_vars[VARP]; // inverse status      
	trail[trail[0]-1]=-trail[trail[0]-1]; // inverto valore su trail
	level[0]=level[VARP];
      }    
    }
  }
} 

__host__ inline void backjumpingf(int level_to_jump){

  // riparte con ultima scelta fatta a level_to_jump
  // annulla trail stack e riferimenti a variabili
  //  printf("backtrack to %d\n",level_to_jump);

  backjump++;
  int loop=trail[0]>0;
  while (loop){ // questo mette a posto trail e variabili. level sistemato sulla ricorsione
    int VARP=abs(trail[trail[0]-1]);
    loop=trail[0]>0 &&       
      ((level[VARP]>level_to_jump)   ||  (refs[VARP]>=0)
       //|| (level[VARP]<=level_to_jump && refs[VARP]==-2)
       );
     
    if (level[VARP]>level_to_jump || refs[VARP]>=0){
      //printf("back var %d lev %d\n",VARP,level[VARP]);
      host_vars[VARP] = -1;		// restore "unknown" status
      level[VARP] = -1;
      refs[VARP] = -1;
      trail[0]--;
    }
    else{ // refs<0 && level=leveltojump
      if (level[VARP]<=level_to_jump && refs[VARP]==-1){ // variabile gia' testata -> metto caso opposto
	//printf("switch var %d at lev %d -> %d\n",VARP,level[VARP],-trail[trail[0]-1]);
	//refs[VARP] = -2;
	//host_vars[VARP] = 1-host_vars[VARP]; // inverse status      
	//trail[trail[0]-1]=-trail[trail[0]-1]; // inverto valore su trail
	level[0]=level[VARP];
      }    
    }
  }
  if (level_to_jump==0)// caso speciale se backjump annulla anche primo livello
    level[0]=0;
} 


//**************************************************************************
//**************************************************************************
//********** MAIN PROCEDURE: twolevel_DPLL *********************************
//********** a part of DPLL (recursive) is handled by the host.  ***********
//********** The final part is made in the device(s) ***********************
//**************************************************************************
//**************************************************************************
//**************************************************************************


__host__  void mask_propagation_cpu(){
  int dbg=0;
    do{ // CPU mask_prop
      clauind = mask_prop(&selected_var,&sat_val);
      if (dbg) printf("prop %d %d %d: ",sat_val,selected_var,clauind);
      if (dbg)
        for (int i=clause_pointer[clauind];i<clause_pointer[clauind+1];i++)
	  printf("%d ",formula[i]);
      if (dbg) printf("\n");

      if (sat_val==1){
	//	unit++;
	int VAR = abs(selected_var);		// Look for its variable and sign
	if (dbg) printf("var prop %d, lev %d\n",selected_var,level[0]);
	host_vars[VAR] = selected_var > 0;
	level[VAR] = level[0];
	refs[VAR] = clauind; // devo sapere la clausola
	trail[trail[0]] = selected_var;
	trail[0]++;
      }    
    }
    while (sat_val==1);
}

__host__ void mask_propagation_gpu(){
  int dbg=0;
    if (1==0 && dbg){
      printf("at lev %d\n",level[0]);
      for(int i=1;i<trail[0];i++)
	printf("%d@%dr%d,",trail[i],level[abs(trail[i])],refs[abs(trail[i])]);
      printf("\n");
    }
    if (dbg) printf("bl %d, ncl %d NC %d NL %d\n",PARblocks,nclausel,NC,NL);
    do{ // GPU mask_prop
      if (dbg) printf(".\n");
      HANDLE_ERROR( cudaMemcpy( dev_parma_vars, host_vars, NV * sizeof(int), cudaMemcpyHostToDevice ) );
      parmask_prop<<<PARblocks, THREADS>>>( mask_data, dev_parma_vars, dev_mapped_formula, dev_mapped_cp,NC);
      HANDLE_ERROR( cudaMemcpy( h_mask_data, mask_data, (3*PARblocks)*sizeof(int), cudaMemcpyDeviceToHost)); // puo' copiare anche 3*blocks+NV per risparmiare chiamate (se servono dati su vars)
      if (dbg) printf("..\n");
      CUDA_count++;
      // colleziona dati da blocchi
      int bestnum=0;
      int bestid=0;                         
      int bestcl=0;                         
      for (int i=0;i<PARblocks;i++){                                
	if (dbg) printf("%d: %d,%d,%d (%d)",i,h_mask_data[3*i],h_mask_data[3*i+1],h_mask_data[3*i+2],bestid);

	if ((h_mask_data[3*i]==-1 && bestnum!=-1) ||
	    (h_mask_data[3*i]==-1 && bestnum==-1 && bestcl>h_mask_data[3*i+2]) ||
	    bestnum==0 ||
	    (bestnum>0 && h_mask_data[3*i]>0 && 
	     (h_mask_data[3*i]<bestnum ||
	      (h_mask_data[3*i]==bestnum && abs(h_mask_data[3*i+1])<abs(bestid)) ||
	      (h_mask_data[3*i]==bestnum && abs(h_mask_data[3*i+1])==abs(bestid) && h_mask_data[3*i+2]<bestcl)
	      ))){
	  bestnum=h_mask_data[3*i];
	  bestid=h_mask_data[3*i+1];
	  bestcl=h_mask_data[3*i+2];
	}    

	/*
	 // TODO: aumenta propagazione in parallelo (una per ogni blocco possibile)
	if (h_mask_data[3*i]==1){
	  int VAR = abs(h_mask_data[3*i+1]);		// Look for its variable and sign
	  if (level[VAR]==-1){ // se non e' gia'stata dedotta da un blocco parallelo (e' sempre -1 grazie a backtracking)
	    if (dbg) printf("var prop %d, lev %d\n",h_mask_data[3*i+1],level[0]);
	    host_vars[VAR] = h_mask_data[3*i+1] > 0;
	    level[VAR] = level[0];
	    refs[VAR] = h_mask_data[3*i+2]; // devo sapere la clausola
	    trail[trail[0]] = h_mask_data[3*i+1];
	    trail[0]++;	
	  }	  
	}
	*/

      }
	if (dbg) printf("\nbest %d %d %d\n",bestnum,bestid,bestcl);
      //printf("\n");
      sat_val  =bestnum;
      selected_var=bestid;
      clauind=bestcl;
      if (dbg) printf("prop %d %d %d: ",sat_val,selected_var,clauind);
      if (dbg)
        for (int i=clause_pointer[clauind];i<clause_pointer[clauind+1];i++)
	  printf("%d ",formula[i]);
      if (dbg) printf("\n");

      if (sat_val==1){
	int VAR = abs(selected_var);		// Look for its variable and sign
	if (dbg) printf("var prop %d, lev %d\n",selected_var,level[0]);
	host_vars[VAR] = selected_var > 0;
	level[VAR] = level[0];
	refs[VAR] = clauind; // devo sapere la clausola
	trail[trail[0]] = selected_var;
	trail[0]++;		  
      }	
    }
    while (sat_val==1);
}


//**********************************************************************************
//**********************************************************************************
//**********************************************************************************
//******************   TWO LEVEL DPLL (no recursive)   *****************************
//**********************************************************************************
//**********************************************************************************
//**********************************************************************************


__host__ int twolevel_DPLL(){
  int selected;
  int pos;
  int FV;
  short good=0;
  int dbg=0;  
  int isbackjumping=0; // a 1 se fatto backj e quindi devo continuare (anche se il bj ha portato il livello a 0 = ho imparato clausola unitaria, e con una UP posso ripartire da livello 1 con l'altro valore)
  do {    

    isbackjumping=0;

    if (dbg) printf("twolevel_DPLL lev %d\n ",level[0]);
    
    // mask_propagation
    if (USEGPU %2 == 0)
      mask_propagation_cpu();
    else 
      mask_propagation_gpu();

  if (dbg){ 
    printf("dopo maskprop lev %d host_vars: ",level[0]);
      for(int i=1;i<trail[0];i++)
	printf("%d@%dr%d,",trail[i],level[abs(trail[i])],refs[abs(trail[i])]);
      printf("\n");
      /*    for (int i=0;i<NV;i++)
      printf("%d>%d(%d,%d) ",i,host_vars[i],level[i],refs[i]);
    printf("\n");
      */
    printf("sat_val %d, selvar %d\n",sat_val,selected_var);     
  }   

  // new case
  if (level[0]<=0 && sat_val<0){ // failed, UP was enough, no need to learn, return fail!
    return 0;
  }
  
  //************* UNSATISFIABLE ASSIGNMENT
  if (sat_val < 0){ // At least one clause is false
    good = 0;
    refs[0] = clauind;	// conflict clause
    if (learning) {// CLAUSE LEARNING
      clause_learning(); 
      if (NC%1000 ==0) printf("NC: %d\n",NC);
      if (dbg) printf("cl: level to backjump to %d\n",level_to_jump);      
      if (dbg) printf("backjump %d\n",level_to_jump);
      backjumpingf(level_to_jump);
      isbackjumping=1;
    }
    else{
      if (dbg) printf("backtrack %d\n",level[0]);
      backtrackingf(level[0]); // se non imparo, apro fratello (tengo livello corrente)
    }

    // END CLAUSE LEARNING
  }  
  else if (sat_val == 0){  //************* FOUND A SOLUTION
      good = 1;  
  }      
  else    
  { // There is a non-ground literal - sat_val > 0
    pos = (selected_var > 0);  
    selected = abs(selected_var);  // Look for its variable and sign 
    FV = NV-trail[0]; // free vars
    if (USEGPU >= 2 && FV <= MaxV)   { //*** GPU CALL: 
      //printf("CUDA caller with %d learnt clauses\n",nclausel);
      good=CUDA_caller(FV);
      if (!good)
	backtrackingf(level[0]); //backtrack
    }   
    else { ///////// vado con CPU
      level[0]++; // apro nuovo ramo
      if (dbg)
	printf("lev %d, 2 scelte: var %d val %d\n",level[0],selected,pos);
      level[selected] = level[0];
      host_vars[selected] = pos;
      trail[trail[0]++] = selected_var;
    } // CPU
  } //sat_val >2


  }
  while(isbackjumping || (level[0]>0 && !good));
  return good;
}


//************************************************************************
//************************************************************************
//*******************       main        **********************************
//************************************************************************
//************************************************************************

__host__ int  main(int argc, char** argv) {
    int retv = 0;
    cudaDeviceProp prop;
    int whichDevice;
    
    if (argc!=7){
      printf("usage: %s mode N_vars_blocks(2LOG_BLKS) N_vars_thread(LOG_THRDS) Maxv learning(0/1) filename\n",argv[0]);
      return -1;
    }    
    
    // *** PARAMETER SETTINGS
    USEGPU      = atoi(argv[1]);
    LOG_BLKS    = atoi(argv[2]);
    LOG_THRDS   = atoi(argv[3]);
    CUDABLOCKS  = (1<<LOG_BLKS);   //*** pow(2,LOG_BLKS); 
    CUDATHREADS = (1<<LOG_THRDS);  //*** pow(2,LOG_THRDS);  
    MaxV        = atoi(argv[4]);  
    Delta       = max(1,MaxV - LOG_BLKS + 1); // *** In case MaxV < LOG_BLKS use default "1"  
    learning    = atoi(argv[5]);

    if (USEGPU>0){
      // CUDA timings initialization
      cudaEventCreate( &start );
      cudaEventCreate( &stop );
      // Checking of CUDA HW capabilities
      HANDLE_ERROR( cudaGetDevice( &whichDevice ) );
      HANDLE_ERROR( cudaGetDeviceProperties( &prop, whichDevice ) );
      if (prop.canMapHostMemory != 1) {
	            printf( "Device cannot map memory.\n" );
	  return 0;
      }
    }

    //*** Some "warning" and exit
    if(USEGPU>0 && CUDATHREADS > prop.maxThreadsPerBlock){
      printf("Cant handle so many vars per block (max %d < %d threads per block)\n",CUDATHREADS,prop.maxThreadsPerBlock);
      return -1;
    }
    printf("using LOG_BLKS=%d, LOG_THRDS=%d MaxV=%d \n",LOG_BLKS,LOG_THRDS,MaxV);
    if (Delta > DV_SIZE){
        printf("Delta too large, recompile kernel with %d in the size of arrays\n",Delta);
        return -1;
    }    
    
    //**** CORE COMPUTATION
     allocate_first();
    // *** Load a SAT formula in DIMACS format
     load_formula(argv[6],&NV,&NC,&NL);
    //*** Print and use some HW and instance info     
     if (USEGPU>0){
       print_info(prop);
       PARblocks = IMIN( ( NC + THREADS - 1) / THREADS, 2*prop.multiProcessorCount );
     } 
    
     allocate_second();
    //*** Start running time (formula loading is not counted)    
     if (USEGPU %2 == 1)
         printf("blocks %d, threads %d\n",PARblocks, THREADS);
     srand ( time(NULL) );
     printf("start\n");
     
     if (USEGPU>0){
       cudaEventRecord(start,0);
       cudaEventSynchronize(start);    
     }else{
#if (defined(_WIN32) || defined(__WIN32__) || defined(WIN32))
       global_time=clock();
#else
       gettimeofday(&time_start,NULL);
#endif
     }

//*********************************************************************
//******* Call the main procedure (store the result in vars[0]): 
//*********************************************************************
      
     retv = twolevel_DPLL();  // <======== HERE! 
     host_vars[0] = retv;  
    
    //*** CUDA timings:
    if (USEGPU>0){
      cudaEventRecord( stop, 0 );    
      cudaEventSynchronize( stop );    
      cudaEventElapsedTime( &deltatime, start, stop );
    }
    else{
#if (defined(_WIN32) || defined(__WIN32__) || defined(WIN32))
     clock_t final=clock()-global_time;
     deltatime=((float)final/CLK_TCK)*1000;
#else
     gettimeofday(&time_stop,NULL);
     deltatime=(0.0+time_stop.tv_sec+time_stop.tv_usec/1000000.0) - (0.0+time_start.tv_sec+time_start.tv_usec/1000000.0);
     deltatime*=1000;
#endif
    }

    /* DISPLAY RESULT AND RUNNING TIME: */
    print_result( host_vars, NV );    
    print_time( deltatime );
             
    deallocate();
    return 0;
}            

//********* END *******************************************************