/*

  $Id$

  Copyright 1989-2016 MINES ParisTech

  This file is part of PIPS.

  PIPS is free software: you can redistribute it and/or modify it
  under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  any later version.

  PIPS is distributed in the hope that it will be useful, but WITHOUT ANY
  WARRANTY; without even the implied warranty of MERCHANTABILITY or
  FITNESS FOR A PARTICULAR PURPOSE.

  See the GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with PIPS.  If not, see <http://www.gnu.org/licenses/>.

*/
#ifdef HAVE_CONFIG_H
    #include "pips_config.h"
#endif
/* package tiling
 *
 * 1. Why?
 *    - memory hierarchy (registers, caches L1/L2/L3, memory, virtual memory, out-of-core,...)
 *    - granularity (amortization of fix costs: synchronization, communication, control,...)
 * 2. Legality
 *    - TO BE IMPLEMENTED
 * 3. Selection
 *    - directions (e.g. communication minimization): Darte/Robert, Hoegsted
 *    - ratios (e.g. critical path minimization)
 *    - volume (fix cost amortization) under memory constraints
 *    - mapping on finite hardware (Robert/...)
 *    restriction to orthogonal tilings: ratios, volume, mapping (Andonov/Rajopadhye)
 * 4. Code Generation (Xue?)
 *    - control and memory addressing overheads
 * 5. Hierarchical Tiling (Ferrante/Carter,...)
 * 6. Data vs Control Tiling
 * 7. Extensions
 *    - perfectly nested loops (IMPLEMENTED)
 *    - non perfectly nested loops (e.g. matrix multiply)
 *    - general nested loops (e.g. ADI)
 *    - sequence of loop nests (signal processing, Thomson-CSF)
 *    - ...
 * See: http://www.irisa.fr/api/Rajopadhye/tiling
 */

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

#include "genC.h"
#include "linear.h"
#include "ri.h"
#include "effects.h"
#include "misc.h"
#include "text.h"
#include "pipsdbm.h"
#include "boolean.h"
#include "vecteur.h"
#include "contrainte.h"
#include "sc.h"
#include "matrice.h"
#include "matrix.h"
#include "sparse_sc.h"
#include "ri-util.h"
#include "effects-util.h"
#include "conversion.h"
#include "dg.h"
typedef dg_arc_label arc_label;
typedef dg_vertex_label vertex_label;
#include "graph.h"
#include "ricedg.h"

#include "graph.h"
#include "transformations.h"
#include "properties.h"
#include "tiling.h"
#include "movements.h"

#include "hyperplane.h"

/* Create a new entity for tile index. Because of the module name, it is easier to postfix
 * by "_t" than to prefix by "t_".
 */
static entity 
make_tile_index_entity(entity old_index)
{
  return make_new_index_entity(old_index, "_t");
}
static entity 
make_local_tile_index_entity(entity old_index)
{
  return make_new_index_entity(old_index, "_l");
}

bool
static_partitioning_matrix(matrice P, int n, const char* serialized_matrix)
{
    pips_assert("interactive_partitioning_matrix", n>=1);
    bool status = false;
    if( serialized_matrix && !string_undefined_p(serialized_matrix) && !empty_string_p(serialized_matrix))
    {
        int row,col;
	string ctxt0 = string_undefined, ctxt1 = string_undefined;
        string saved_ptr,buffer,elem = NULL;
        saved_ptr= buffer = strdup(serialized_matrix);
        string line = strtok_r(buffer,",",&ctxt0);
        DENOMINATOR(P) = VALUE_ONE;
        for(row=1;row<=n;row++){
            elem = strtok_r(line," ",&ctxt1);

            for(col=1;col<=n;col++)
            {
                if(!elem) elem = col==row ? "1" : "0";
                ACCESS(P, n, row, col)=atoi(elem);
                elem = strtok_r(NULL," ",&ctxt1);
            }
            line = strtok_r(NULL,",",&ctxt0);
        }
        status= ( line == NULL ) && (elem == NULL );
        free(saved_ptr);
    }
    return status;
}

/* Query the user for a partitioning matrix P
 */

bool
interactive_partitioning_matrix(matrice P, int n)
{
    int n_read;
    string resp = string_undefined;
    string cn = string_undefined;
    bool return_status = false;
    int row;
    int col;

    DENOMINATOR(P) = VALUE_ONE;
    /* Query the user for P's components */
    pips_assert("interactive_partitioning_matrix", n>=1);
    debug(8, "interactive_partitioning_matrix", "Reading P\n");

    for(row=1; row<=n; row++) {
	resp = user_request("Partitioning matrix (%dx%d)?\n"
			    "(give all its integer coordinates on one line of %d per row): ",
			    n, n, n);
	if (resp[0] == '\0') {
	    user_log("Tiling loop transformation has been cancelled.\n");
	    return_status = false;
	}
	else {    
	    cn = strtok(resp, " \t");

	    return_status = true;
	    for( col = 1; col<=n; col++) {
		if(cn==NULL) {
		    user_log("Too few coordinates. "
			     "Tiling loop transformation has been cancelled.\n");
		    return_status = false;
		    break;
		}
		n_read = sscanf(cn," " VALUE_FMT, &ACCESS(P, n, row, col));
		if(n_read!=1) {
		    user_log("Too few coordinates. "
			     "Hyperplane loop transformation has been cancelled.\n");
		    return_status = false;
		    break;
		}
		cn = strtok(NULL, " \t");
	    }
	}

	if(cn!=NULL) {
	    user_log("Too many coordinates. "
		     "Tiling loop transformation has been cancelled.\n");
	    return_status = false;
	}
    }

    ifdebug(8) {
	if(return_status) {
	    pips_debug(8, "Partitioning matrix:\n");
	    matrice_fprint(stderr, P, n, n);
	    (void) fprintf(stderr,"\n");
	    pips_debug(8, "End\n");
	}
	else {
	    pips_debug(8, "Ends with failure\n");
	}
    }

    return return_status;
}


/* Generate the tile membership constraints between a tile coordinates and
 an iteration coordinate
 */

static Psysteme
tile_membership_constraints(Pbase initial_basis,
			    Pbase tile_basis,
			    matrice HT,
			    Pvecteur tiling_offset)
{
    Psysteme mc = sc_new();
    int dim = base_dimension(initial_basis);
    int row;
    int col;
    Value k = DENOMINATOR(HT);
    Value up = k - VALUE_ONE;
    Pbase civ = BASE_UNDEFINED;
    Pbase ctv = BASE_UNDEFINED;

    ifdebug(8) {
	debug(8, "tile_membership_constraints", "Begin with Matrix HT:\n");
	matrice_fprint(stderr, HT, dim, dim);
    }

    pips_assert("The two bases have the same dimension", dim == base_dimension(tile_basis));

    for(row = 1; row <= dim; row++) {
	Pvecteur upper = VECTEUR_NUL;
	Pvecteur lower = VECTEUR_NUL;
	Pcontrainte cupper = CONTRAINTE_UNDEFINED;
	Pcontrainte clower = CONTRAINTE_UNDEFINED;

	for(col = 1, civ = initial_basis, ctv = tile_basis;
	    col <= dim;
	    col++, civ = vecteur_succ(civ), ctv = vecteur_succ(ctv)) {
	    if(ACCESS(HT, dim, row, col)!=VALUE_ZERO) {
		Value coeff = ACCESS(HT, dim, row, col);
		Value offset = vect_coeff(vecteur_var(civ), tiling_offset);

		vect_add_elem(&upper, vecteur_var(civ), coeff);
		vect_add_elem(&upper, TCST, value_uminus(offset*coeff));
	    }
	    if(col==row) {
		vect_add_elem(&upper, vecteur_var(ctv), value_uminus(k));
	    }
	}
	lower = vect_dup(upper);
	vect_chg_sgn(lower);
	vect_add_elem(&upper, TCST, value_uminus(up));
	cupper = contrainte_make(upper);
	clower = contrainte_make(lower);
	sc_add_inegalite(mc, cupper);
	sc_add_inegalite(mc, clower);
    }

    sc_creer_base(mc);

    ifdebug(8) {
	pips_debug(8, "End with constraint system mc:\n");
	sc_fprint(stderr, mc, (get_variable_name_t)entity_local_name);
    }

    return mc;
}

/* Generate the tile membership constraints between a tile coordinates and
 an iteration coordinate
*/

Psysteme
tile_hyperplane_constraints(Pbase initial_basis,
			     Pbase tile_basis,
			     matrice HT,
			     Pvecteur tiling_offset)
{
  Psysteme mc = sc_new();
  int dim = base_dimension(initial_basis);
  int row;
  int col,col2;
  Value k = DENOMINATOR(HT);
  Value up = k - VALUE_ONE;
  Pbase civ = BASE_UNDEFINED;
  Pbase ctv = BASE_UNDEFINED;
  matrice  G_tile;
  Value *h_tile; // hyperplane direction
  ifdebug(8) {
    debug(8, "tile_membership_constraints", "Begin with Matrix HT:\n");
    matrice_fprint(stderr,HT, dim, dim);
  }
  int n = dim; 
  pips_assert("The two bases have the same dimension", dim == base_dimension(tile_basis));

  ifdebug(8) {
    pips_debug(8,"hyperplane matrix HT:");
    matrice_fprint(stderr, HT, n, n);
    (void) fprintf(stderr,"\n");
  }

  /* computation of the hyperplane tile direction in the tile basis: = sum(hi)= 1H */
  h_tile = (Value*)(malloc(n*sizeof(Value)));
  G_tile = matrice_new(n,n);
 
  string tile_direction = (string) get_string_property("TILE_DIRECTION");
  if(!empty_string_p(tile_direction) && strcmp(tile_direction,"TP") == 0) {
    for(col=0; col<n; col++) {
      h_tile[col]=VALUE_ONE;
    }
    
    ifdebug(8) {
      (void) fprintf(stdout,"The hyperplane direction h_tile is :");
      for(col=0; col<n; col++) 
	(void) printf("%d ;",(int)h_tile[col]);
      (void) printf("\n");
    }
    /* computation of the tile scanning base G.  
     */
   scanning_base_hyperplane(h_tile, n, G_tile);
    ifdebug(8) {
      (void) fprintf(stderr,"The tile scanning base G_tile is :");
      matrice_fprint(stderr, G_tile, n, n);
    }
  }
  else matrice_identite(G_tile,n,0);
  
    /* Build the constraints
       0 <= det(HT).[HT. (i -i0)]<det(HT).1 with 
       i0=P.G_tile.i_tile
     <==> 0 <= det(HT).[HT. (i -P.G_tile.i_tile)]<det(HT).1
     <==> 0 <= det(HT).[HT.i -G_tile.i_tile]<= det(HT)-1
  */
  for(row = 1; row <= dim; row++) {
    Pvecteur upper = VECTEUR_NUL;
    Pvecteur lower = VECTEUR_NUL;
    Pcontrainte cupper = CONTRAINTE_UNDEFINED;
    Pcontrainte clower = CONTRAINTE_UNDEFINED;

    for(col = 1, civ = initial_basis; col <= dim; col++, civ = vecteur_succ(civ)) {
      if(ACCESS(HT, dim, row, col)!=VALUE_ZERO) {
	Value coeff = ACCESS(HT, dim, row, col);
	Value offset = vect_coeff(vecteur_var(civ), tiling_offset);

	vect_add_elem(&upper, vecteur_var(civ), coeff);
	vect_add_elem(&upper, TCST, value_uminus(offset*coeff));
     }
    }
    for(col2 = 1, ctv = tile_basis; col2 <= dim; col2++,ctv = vecteur_succ(ctv)) {

      /* Find the origin of the iteration domain. Use 0 as default coordinate */
      if(ACCESS(G_tile, dim, row, col2)!=VALUE_ZERO) {
	Value coeff2 = ACCESS(G_tile, dim,row, col2);
	vect_add_elem(&upper, vecteur_var(ctv), value_uminus(k) *coeff2);	
      }
    }
  
    lower = vect_dup(upper);
    vect_chg_sgn(lower);
    vect_add_elem(&upper, TCST, value_uminus(up));
    cupper = contrainte_make(upper);
    clower = contrainte_make(lower);
    sc_add_inegalite(mc, cupper);
    sc_add_inegalite(mc, clower);
  }
  sc_creer_base(mc);

  ifdebug(8) {
    pips_debug(8, "End with constraint system mc:\n");
    sc_fprint(stderr, mc, (get_variable_name_t)entity_local_name);
  }

  return mc;
}

static void compute_local_change_of_basis( Value *h_tile, matrice P, matrice G, int option, int dim, int selected_Pdim)
{
  int row;
  int col;
  if (option ==1) {
    /* computation of the hyperplane tile direction in the tile basis: = sum(hi)= 1H */
    for(row=0; row<dim; row++) {
      h_tile[row]=VALUE_ZERO;
      for(col=0; col<dim; col++) {
	h_tile[row]+=ACCESS(P, dim, row+1, col+1);
      }
    }
  }
  else if (option==2)
    { /* The local tile direction is colinear to the (selected_Pdim) partitioning vecteur  */
      for(col=0; col<dim; col++) 
	h_tile[col]=ACCESS(P, dim, col+1, selected_Pdim);
    }
  ifdebug(8) {
    (void) fprintf(stdout,"The first scanning direction h_tile in %s mode is :",(option==1)?"LP":"LS");
    for(col=0; col<dim; col++) 
      (void) printf("%d ;",(int)h_tile[col]);
    (void) printf("\n");
  }
  
  /* computation of the local tile scanning base G according to h_tile 
   */
  scanning_base_hyperplane(h_tile, dim, G);
}

Psysteme
local_tile_constraints(Pbase initial_basis,Pbase local_tile_basis, Psysteme sc, matrice P, Pvecteur *pvch, statement st)
{
  Psysteme mc = sc_dup(sc),mc2;
  int dim = base_dimension(initial_basis);
  int row;
  int col;
  Pbase civ = BASE_UNDEFINED;
  Pbase ctv = BASE_UNDEFINED;
  Value *h_tile; // hyperplane direction
  int n = dim; 
  matrice V = matrice_new(n,n);
  matrice G_tile = matrice_new(n,n);
  bool legal=false;
  
  pips_assert("The two bases have the same dimension", dim == base_dimension(local_tile_basis));

  string local_tile_direction = (string) get_string_property("LOCAL_TILE_DIRECTION");
  if (!empty_string_p(local_tile_direction) && strcmp(local_tile_direction,"LI") != 0) {
    h_tile = (Value*)(malloc(n*sizeof(Value)));
    if(strcmp(local_tile_direction,"LP") == 0) {
      /* computation of the hyperplane tile direction in the tile basis: = sum(hi)= 1H */
      compute_local_change_of_basis(h_tile, P, G_tile, 1 , n, 1);
      matrice_general_inversion(G_tile,V,n);
      legal = check_tiling_legality(st,V,n);
      if (!legal)
	matrice_identite(G_tile,n,0); 
    }
    else // search a valid direction colinear to one of the partitioning vectors
      if(strcmp(local_tile_direction,"LS") == 0) {
	int pdim;
	legal=false;
	for (pdim=1; !legal && pdim<=n; pdim++) {
	  /* The local tile direction is legal and colinear to the partitioning vecteur (pdim) */
	  compute_local_change_of_basis(h_tile, P, G_tile, 2, n, pdim);
	  ifdebug(8) {
	    fprintf(stderr,"Check the legality of the chosen LOCAL basis");
	    (void) fprintf(stderr,"Temporary local tile scanning base G_tile obtained with %d-th Pvecteur :",pdim);
	    matrice_fprint(stderr, G_tile, n, n);
	  }
	  matrice_general_inversion(G_tile,V,n);
	  legal = check_tiling_legality(st,V,n);
	}
	if (pdim>n && !legal)
	  matrice_identite(G_tile,n,0);  
      }
  }
  else // the local tile direction is the original basis
    matrice_identite(G_tile,n,0);
  
  ifdebug(8) {
    (void) fprintf(stderr,"The local tile scanning base G_tile - i = G_tile.i_prime is :");
    matrice_fprint(stderr, G_tile, n, n);
  }
  /* Build the constraints
     0 <= det(HT).[HT. (i -i0)]<det(HT).1 with 
     i0=P.G_tile.i_tile
     <==> 0 <= det(HT).[HT. (i -P.G_tile.i_tile)]<det(HT).1
     <==> 0 <= det(HT).[HT.i -G_tile.i_tile]<= det(HT)-1
  */
  for(ctv = local_tile_basis; !VECTEUR_NUL_P(ctv) ; ctv = vecteur_succ(ctv)) {
    base_add_variable(mc->base,vecteur_var(ctv));
    sc_dimension(mc)++;
  }
  for(row = 1, civ = initial_basis;  row <= dim; row++, civ = vecteur_succ(civ)) {
    Pvecteur veq = vect_new(vecteur_var(civ),VALUE_ONE);
    Pcontrainte ceq = CONTRAINTE_UNDEFINED;

    for(col = 1, ctv = local_tile_basis; col <= dim ; col++, ctv = vecteur_succ(ctv)) {
      if(ACCESS(G_tile, dim, row, col)!=VALUE_ZERO) {
	Value coeff = ACCESS(G_tile, dim, row, col);
	vect_add_elem(&veq, vecteur_var(ctv), -1*coeff);
      }
    }  
    ceq = contrainte_make(veq);
    sc_add_egalite(mc, ceq);
  }

  ifdebug(8) {
    pips_debug(8, "End with constraint system mc:\n");
    sc_fprint(stderr, mc, (get_variable_name_t)entity_local_name);
    fprintf(stderr, "sc dimension = %d \n",mc->dimension);
  }
  mc2=sc_dup(mc);
  sc_projection_along_variables_ofl_ctrl(&mc2,initial_basis , OFL_CTRL);
  ifdebug(8) {
    pips_debug(8, "End with constraint system mc after change of basis:\n");
    sc_fprint(stderr, mc, (get_variable_name_t)entity_local_name);
  }
  scanning_base_to_vect(G_tile, n, local_tile_basis, pvch);
  return mc2;
}


Pvecteur
loop_nest_to_offset(list lls)
{
    Pvecteur origin = VECTEUR_NUL;
    list cs = list_undefined;

    for (cs = lls; cs != NIL; cs = CDR(cs)){
	loop l = instruction_loop(statement_instruction(STATEMENT(CAR(cs))));
	entity ind = loop_index(l);
	range r = loop_range(l);
	expression lower = range_lower(r);
	intptr_t val;

	if(expression_integer_value(lower, &val)) {
	    vect_chg_coeff(&origin, (Variable) ind, (Value) val);
	}
    }

    return origin;
}

/* Generate tiled code for a loop nest, PPoPP'91, p. 46, Figure 15.
 *
 * The row-echelon algorithm is called from new_loop_bound().
 */

statement tiling_transformation(list lls, _UNUSED_ bool (*u)(statement))
{
    Psysteme sci;			/* iteration domain */
    Psysteme sc_tile_scan;
    Psysteme sc_tile;
    Psysteme mc = SC_UNDEFINED; /* Tile membership constraints */
    Psysteme sc_B_prime = SC_UNDEFINED;
    Psysteme sc_B_second = SC_UNDEFINED;
    Pbase initial_basis = NULL;
    Pbase tile_basis = NULL;
    Pbase reverse_tile_basis = NULL;
    /* Pbase local_basis = NULL; */
    Pbase new_basis = NULL;
    matrice P; /* Partitioning matrix */
    matrice HT; /* Transposed matrix of the inverse of P */
    matrice G; /* Change of basis in the tile space to use vector 1 as hyperplane direction */
    int n;				/* number of indices, i.e. loop nest depth */
    Value *h;
    statement s_lhyp;
    Pvecteur *pvg;
    Pbase pb;
    expression lower, upper;
    int col;
    Pvecteur to = VECTEUR_NUL; /* Tiling offset: 0 by default */

    debug_on("TILING_DEBUG_LEVEL");

    pips_debug(8, "Begin with iteration domain:\n");

    /* make the constraint system for the iteration space and find a good
       origin for the tiling */

    const char* smatrix = get_string_property("LOOP_TILING_MATRIX");

    if (get_bool_property("PARTIAL_LOOP_TILING")) {
      int ndim=0;
      if (smatrix && !string_undefined_p(smatrix) && !empty_string_p(smatrix))
	for (int k=1; k <=(int)strlen(smatrix); k++)
	  if (smatrix[k]==',') ndim++;
      if (ndim!=0 && ndim<(int)gen_length(lls)-1) {
	for (int k =1; k< (int)gen_length(lls)-ndim && lls != NIL; k++)
	  lls=CDR(lls);
      }
    }
    
    sci = loop_iteration_domaine_to_sc(lls, &initial_basis);
    n = base_dimension(initial_basis);
    to = loop_nest_to_offset(lls);
    ifdebug(8) {
      sc_fprint(stderr, sci, (get_variable_name_t)entity_local_name);
      pips_debug(8,"And with origin:\n");
      vect_fprint(stderr, to, (get_variable_name_t)entity_local_name);
    }

    /* computation of the partitioning matrix P and its inverse HT */

    P = matrice_new(n, n);
    HT = matrice_new(n, n);
    if(
       !static_partitioning_matrix(P,n,smatrix) &&
       !interactive_partitioning_matrix(P, n)
       ) {
      pips_user_error("A proper partitioning matrix was not provided\n");
    }

    ifdebug(8) {
      pips_debug(8,"Partitioning matrix P:");
      matrice_fprint(stderr, P, n, n);
      (void) fprintf(stderr,"\n");
    }

    matrice_general_inversion(P, HT, n);
    
    ifdebug(8) {
      pips_debug(8,"Inverse partitioning matrix HT:");
      matrice_fprint(stderr, HT, n, n);
      (void) fprintf(stderr,"\n");
    }

    /* Compute B': each iteration i in the iteration space is linked to its tile s */

    derive_new_basis(initial_basis, &tile_basis, make_tile_index_entity);
     mc = tile_membership_constraints(initial_basis, tile_basis, HT, to);
    mc = sc_normalize(mc);
    ifdebug(8) {
	(void) fprintf(stderr,"Tile membership constraints:\n");
	sc_fprint(stderr, mc, (get_variable_name_t)entity_local_name);
    }
    /* mc and SC_B_prime are aliased after this call */
    sc_B_prime = sc_append(mc, sci);
    ifdebug(8) {
	(void) fprintf(stderr,"sc_B_prime after call to sc_append (is the basis ok?):\n");
	sc_fprint(stderr, sc_B_prime, (get_variable_name_t)entity_local_name);
    }

  
    mc = SC_UNDEFINED;
    /* Save a copy to compute B" later */
    sc_B_second = sc_dup(sc_B_prime);

    /* Get constraints on tile coordinates */

    sc_projection_along_variables_ofl_ctrl(&sc_B_prime, initial_basis, OFL_CTRL);
    ifdebug(8) {
	(void) fprintf(stderr,"Tile domain:\n");
	sc_fprint(stderr, sc_B_prime, (get_variable_name_t)entity_local_name);
    }

    /* Build the constraint system to scan the set of tiles */
    sc_tile_scan = new_loop_bound(sc_B_prime, tile_basis);
    ifdebug(8) {
	(void) fprintf(stderr,"Tile domain in echelon format:\n");
	sc_fprint(stderr, sc_tile_scan, (get_variable_name_t)entity_local_name);
    }

    /* CA: Build the new basis (tile_basis+initial_basis)*/
    /* base It, Jt, I, J  pour notre exemple */
    new_basis = vect_add(vect_dup(initial_basis),vect_dup(tile_basis));
    ifdebug(8) {
	(void) fprintf(stderr,"new_basis\n");
	vect_fprint(stderr, new_basis, (get_variable_name_t)entity_local_name);
    }

    /* Build the constraint system sc_tile to scan one tile (BS IN
       PPoPP'91 paper) */
    ifdebug(8) {
	(void) fprintf(stderr,"sc_B_second:\n");
	sc_fprint(stderr, sc_B_second, (get_variable_name_t)entity_local_name);
    }
    sc_tile = new_loop_bound(sc_B_second, new_basis);
    ifdebug(8) {
	(void) fprintf(stderr,"Iteration domain for one tile:\n");
	sc_fprint(stderr, sc_tile, (get_variable_name_t)entity_local_name);
    }

    /* computation of the hyperplane tile direction: let's use the
       default 1 vector */
    h = (Value*)(malloc(n*sizeof(Value)));
    for(col=0; col<n; col++) {
	h[col] = VALUE_ONE;
    }
    /* computation of the tile scanning base G: right now, let's
     * assume it's Id.  This is OK to tile parallel loops... or to
     * scan tiles sequentially on a monoprocessor.
     */
    G = matrice_new(n,n);
    scanning_base_hyperplane(h, n, G);
    ifdebug(8) {
	(void) fprintf(stderr,"The tile scanning base G is before ID:");
	matrice_fprint(stderr, G, n, n);
    }
    matrice_identite(G, n, 0);
    ifdebug(8) {
	(void) fprintf(stderr,"The tile scanning base G is:");
	matrice_fprint(stderr, G, n, n);
    }

    /* generation of code for scanning one tile */

    /* Compute the local coordinate changes: there should be none for the time being
     * because we keep the initial basis to scan iterations within one tile, i.e
     * G must be the identity matrix
     */
    pvg = (Pvecteur *)malloc((unsigned)n*sizeof(Svecteur));
    scanning_base_to_vect(G, n, initial_basis, pvg);
     ifdebug(8) {
      (void) fprintf(stderr,"The local coordinate change vector is:");
      for (int i=1; i<=n; i++)
	vect_fprint(stderr, pvg[i], (get_variable_name_t)entity_local_name);
    }
    /* generation of code to scan one tile, and update of loop body using pvg */

    s_lhyp = code_generation(lls, pvg, initial_basis,
			     new_basis, sc_tile, false);

    /* generation of code for scanning all tiles */

    reverse_tile_basis = base_reversal(tile_basis);
    for (pb = reverse_tile_basis; pb!=NULL; pb=pb->succ) {
	loop tl = loop_undefined;

	make_bound_expression(pb->var, tile_basis, sc_tile_scan, &lower, &upper);
	tl = make_loop((entity) vecteur_var(pb),
		       make_range(copy_expression(lower), copy_expression(upper),
				  int_to_expression(1)),
		       s_lhyp,
		       entity_empty_label(),
		       make_execution(is_execution_sequential, UU),
		       NIL);
	s_lhyp = instruction_to_statement(make_instruction(is_instruction_loop, tl));
    }

    pips_debug(8,"End\n");

    debug_off();

    return(s_lhyp);
}

bool loop_tiling(const char* module_name)
{
    bool return_status = false;

    return_status = interactive_loop_transformation(module_name,  (statement (*)(list, bool (*)(statement)))tiling_transformation);
    
    return return_status;
}



static bool statement_in_loopnest_p = false;
static statement test_statement_of_reference;
static void statement_in_loopnest(statement s)
{
  if (statement_number(test_statement_of_reference) == statement_number(s))
    statement_in_loopnest_p= true;
}

/*
   Check that H.v >0
   v is a dependence vecteur  or a dependence vertex expressed in base b
   Matrice H  is the tiling scanning base, expressed in base b.
 */
static void legal_point_p(Pvecteur v, Ptsg gs, matrice H,int dim, bool *legal)
{
  Pbase b=gs->base;
  Pvecteur coord;
  int row, col, accu;
  if(vect_in_basis_p(v, b)) {
    for (row=1; row<=dim; row++) {
      accu = 0;
      for(coord = b, col=1; !VECTEUR_NUL_P(coord) && col<=dim; coord = coord->succ,col++) {
	Variable var = vecteur_var(coord);
	if(VARIABLE_DEFINED_P(var)) {
	  accu += vect_coeff(var, v)*ACCESS(H, dim, row, col);
	}
      }
      if (accu <0) {
	ifdebug(8) {
	  pips_debug(8,"Tile scanning matrix H:");
	  fprintf(stdout,"NON LEGAL TILING: H(%d,*)*SG =%d <0\n", row,accu);
	  sg_fprint_as_dense(stdout, gs, gs->base );
	}
	*legal=(*legal) && false;
      }
      else  {
	ifdebug(8) {
	  pips_debug(8,"Tile scanning matrix H:");
	  fprintf(stdout,"LEGAL TILING: H(%d,*)*SG =%d >=0\n", row,accu);
	  sg_fprint_as_dense(stdout, gs, gs->base );
	}
      }			
    }		      
  }
}

static void check_positive_dependence(Ptsg gs, matrice H,int dim, bool *legal)
{
  if( sg_nbre_rayons(gs) > 0) {
    for (Pray_dte e = sg_rayons(gs); e != NULL; e = e->succ) {
      legal_point_p(e->vecteur, gs, H,dim, legal);
    }
  }
  if( legal && sg_nbre_droites(gs) > 0) {
    for (Pray_dte e = sg_droites(gs); e != NULL; e = e->succ) {
      Pvecteur v=vect_copy(e->vecteur);
      legal_point_p(e->vecteur, gs, H,dim, legal);
      vect_chg_sgn(v);
      legal_point_p(e->vecteur, gs, H,dim, legal);
      vect_rm(v);
    }
  }
  if( legal && sg_nbre_sommets(gs) > 0) {
    for (Psommet e = sg_sommets(gs); e != NULL; e = e->succ) {
      legal_point_p(e->vecteur, gs, H,dim, legal);
    }
  }
}
/*
   Test whether H.D >0. 
   Returns TRUE if ok, FALSE otherwise
   D is dependence cone (vertex, rays, lines) in base b 
     among all the dependence conflicts belonging to loopnest
   Matrice H is a squared matrix of dimension dim. 
   It is the tiling scanning base, expressed in base b.
 */
bool check_tiling_legality(statement loopnest_st, matrice H,int dim)
{
  const char * module_name = get_current_module_name();
  graph dg = (graph) db_get_memory_resource(DBR_DG, module_name, true);
   int bl;
   matrice HC;
  // statement mod_stat = get_current_module_statement();
  bool legal=true;
  FOREACH(VERTEX,v,graph_vertices(dg)) {
    statement s = vertex_to_statement(v);
    // Verify s is a statement in the loop nest
    statement_in_loopnest_p= false;
    test_statement_of_reference = s;
    gen_recurse(loopnest_st, statement_domain, gen_true,statement_in_loopnest);
	
    if(statement_in_loopnest_p) {
      FOREACH(SUCCESSOR,su,vertex_successors(v)) {
	      
	vertex v2 = successor_vertex(su);
	statement s2 = vertex_to_statement( v2 );
	// Verify s2 is a statement in the loop nest      
	statement_in_loopnest_p= false;
	test_statement_of_reference = s2;
	gen_recurse(loopnest_st, statement_domain, gen_true,statement_in_loopnest);
	
	if(statement_in_loopnest_p) { // s and s2 are in the loop nest
	  FOREACH(CONFLICT,c,dg_arc_label_conflicts(successor_arc_label(su))) {
                       
	    if ( conflict_cone(c) != cone_undefined ) {
	      //test H.D >0 for sommets, rays, and lines	    
	      Ptsg gs = (Ptsg) cone_generating_system(conflict_cone(c));
	      if( !SG_UNDEFINED_P(gs)) {
		Pbase b=gs->base;
		bl=vect_size(b);
		
		if (bl==dim)
		  check_positive_dependence(gs, H,dim, &legal);
		else {
		  if (bl > dim) {
		    HC=matrice_new(bl,bl);
		    int row,col;
		    ifdebug(8) {
		      fprintf(stdout," Different SIZES for base %d and matrice %d, Partial tiling case ?\n",bl,dim);
		    }
		    // Complete matrix Ht for the first non tiled dimension of the loop nest
		    // for 1<=i<=bl-dim, if di!=0  hji>=loop increment
		    // only loop with increment 1 are tiled ==> hji =1
		    matrice_identite(HC,bl,0);
		    DENOMINATOR(HC)=DENOMINATOR(H);
		    for (row=1; row<=bl; row++) {
		      if (row>=bl-dim+1) {
			for(col=1;  col<bl-dim+1;col++) {
			  ACCESS(HC, bl, row, col)=DENOMINATOR(H) ;
			}
			for(col=bl-dim+1;  col<=bl;col++) {
			  ACCESS(HC, bl, row, col)= ACCESS(H, dim, row-bl+dim, col-bl+dim);
			}
		      }
		      else {
			ACCESS(HC, bl, row, row)=DENOMINATOR(H);
		      }
		    }
		    check_positive_dependence(gs, HC,bl, &legal);
		    matrice_free(HC);
		  }
		  else
		    pips_user_warning(" Different SIZES for dependence cone basis %d and partitioning matrice %d\n",bl,dim);
		}
	      }
	    }
	  }
	}
      }
    }
  }
  return(legal);
}

Psysteme sc_projection_concat_proj_on_variables(
  Psysteme sc,
  Pbase index_base)
{
  Pvecteur lvar_proj;
  Psysteme sc1 = sc_dup(sc);
  int nb=vect_size(index_base);
  // Automatic variables read in a CATCH block need to be declared volatile as
  // specified by the documentation
  volatile Psysteme sc2,sc3;
   sc3=sc_init_with_sc(sc);
   
  if (!VECTEUR_NUL_P(index_base))
  {
    volatile Pvecteur pv1;
    lvar_proj = vect_copy(base_reversal(index_base));
    sc2 = sc_copy(sc);
    for (pv1 = lvar_proj;!VECTEUR_NUL_P(pv1) && nb>1; pv1=pv1->succ, nb--) {
	    
      CATCH(overflow_error) {
	sc_elim_var(sc2, vecteur_var(pv1));
      }
      TRY {// force the overflow forwarding in order to apply sc_elim_var in case of overflow
	sc_projection_along_variable_ofl_ctrl
	  (&sc2,vecteur_var(pv1), FWD_OFL_CTRL);
	UNCATCH(overflow_error);
      }
      sc2 = sc_normalize(sc2);
      if (SC_EMPTY_P(sc2)) {
	sc2 = sc_empty(base_copy(sc->base));
	break;
      }
      else {
        build_sc_nredund_1pass(&sc2);
      }
      // Concatenation of the intermediate system sc2 resulting
      // from the projection of pv1 and the previous concatenated system sc3
      sc3 = sc_append(sc3,sc2);
      if (SC_EMPTY_P(sc3)) {
	sc3 = sc_empty(base_copy(sc->base));
	break;
      }
    }
  }
  sc1 = sc_append(sc1,sc3);
  sc1 = sc_normalize(sc1);
  if (SC_EMPTY_P(sc1))
    sc1 = sc_empty(base_copy(sc->base));     
  return sc1;
}


statement parallel_tiling(list lls, _UNUSED_ bool (*u)(statement))
{
    Psysteme sci;			/* iteration domain */
    Psysteme sc_tile;
    Psysteme sc_local_tile;
    Psysteme mc = SC_UNDEFINED; /* Tile membership constraints */
    Psysteme sc_B_prime = SC_UNDEFINED;
    Psysteme sc_B_prime_2 = SC_UNDEFINED;
    Pbase initial_basis = NULL;
    Pbase tile_basis = NULL;
    Pbase local_tile_basis = NULL;
    Pbase reverse_tile_basis = NULL;
    Pbase new_basis = NULL;
    matrice P; /* Partitioning matrix */
    matrice HT; /* Transposed matrix of the inverse of P */
    int n,i;				/* number of indices, i.e. loop nest depth */
    statement s_lhyp;
    Pvecteur *pvch;
    Pbase pb;
    expression lower, upper;
    Pvecteur to = VECTEUR_NUL; /* Tiling offset: 0 by default */
    statement stmf = statement_undefined;
    list lls2; // first loop in the loop nest;
#define maxscinfosize 100
    int sc_info[maxscinfosize][4];
    Psysteme sc_proj2,sc_tmp;
    int space;
    Psysteme *list_of_systems;
    statement mod_stat = get_current_module_statement();
    set_ordering_to_statement(mod_stat);

    debug_on("TILING_DEBUG_LEVEL");
    pips_debug(8, "Begin with iteration domain:\n");
    const char* smatrix = get_string_property("LOOP_TILING_MATRIX");

    if (get_bool_property("PARTIAL_LOOP_TILING")) {
      int ndim=0;
      if (smatrix && !string_undefined_p(smatrix) && !empty_string_p(smatrix))
	for (int k=1; k <=(int)strlen(smatrix); k++)
	  if (smatrix[k]==',') ndim++;
      if (ndim!=0 && ndim<(int)gen_length(lls)-1) {
	for (int k =1; k< (int)gen_length(lls)-ndim && lls != NIL; k++)
	  lls=CDR(lls);
      }
    }
    // keep the original statement in case of NOT legal tiling
    for (lls2=lls;(int)gen_length(lls2)>1 && lls2 != NIL; lls2=CDR(lls2));
    stmf=STATEMENT(CAR(lls2));

    sci = loop_iteration_domaine_to_sc(lls, &initial_basis);
    n = base_dimension(initial_basis);
    to = loop_nest_to_offset(lls);
    ifdebug(8) {
      sc_fprint(stderr, sci, (get_variable_name_t)entity_user_name);
      pips_debug(8,"And with origin:\n");
      vect_fprint(stderr, to, (get_variable_name_t)entity_user_name);
    }

    /* computation of the partitioning matrix P and its inverse HT */
    P = matrice_new(n, n);
    HT = matrice_new(n, n);
    if(
       !static_partitioning_matrix(P,n,smatrix) &&
       !interactive_partitioning_matrix(P, n)
       ) {
      pips_user_error("A proper partitioning matrix was not provided\n");
    }

    ifdebug(8) {
      pips_debug(8,"Partitioning matrix P:");
      matrice_fprint(stderr, P, n, n);
      (void) fprintf(stderr,"\n");
    }
    matrice_general_inversion(P, HT, n);  
    ifdebug(8) {
      pips_debug(8,"Inverse partitioning matrix HT = P^-1:");
      matrice_fprint(stderr, HT, n, n);
      (void) fprintf(stderr,"\n");
    }  
    /* make the constraint system for the iteration space and find a good
       origin for the tiling if the transformation is legal*/
    if (!get_bool_property("CHECK_TRANSFORMATION_LEGALITY") || check_tiling_legality(STATEMENT(CAR(lls2)),HT,n)) {
      
      derive_new_basis(initial_basis, &tile_basis, make_tile_index_entity);
      derive_new_basis(initial_basis, &local_tile_basis, make_local_tile_index_entity);
      /* Compute B': each iteration i in the iteration space is linked to its tile s */
     
      mc = tile_hyperplane_constraints(initial_basis, tile_basis, HT, to);
      mc = sc_normalize(mc);
      ifdebug(8) {
	(void) fprintf(stderr,"Tile membership constraints:\n");
	sc_fprint(stderr, mc, (get_variable_name_t)entity_user_name);
      }
      /* mc and SC_B_prime are aliased after this call */
      sc_B_prime = sc_append(mc, sci);
      ifdebug(8) {
	(void) fprintf(stderr,"sc_B_prime after call to sc_append :\n");
	sc_fprint(stderr, sc_B_prime, (get_variable_name_t)entity_user_name);
      }

      /* computation of the base scanning local tile elements and its constraints*/
      pvch = (Pvecteur *)malloc((unsigned)n*sizeof(Svecteur));
      sc_B_prime_2= local_tile_constraints(initial_basis,local_tile_basis,sc_B_prime, P, pvch, STATEMENT(CAR(lls2)));

      ifdebug(8) {
	(void) fprintf(stderr,"sc_B_prime2 after change of local basis:\n");
	sc_fprint(stderr, sc_B_prime_2, (get_variable_name_t)entity_user_name);
      }

      /* Build the new basis (tile_basis+local_tile_basis)*/
      new_basis = vect_add(vect_dup(local_tile_basis),vect_dup(tile_basis));
      ifdebug(8) {
	(void) fprintf(stderr,"new_basis\n");
	vect_fprint(stderr, new_basis, (get_variable_name_t)entity_user_name);
      }
      {
	sc_tmp= sc_dup(sc_B_prime_2);
	space = (vect_size(new_basis)+1) * sizeof(Ssysteme);
	list_of_systems = (Psysteme *) malloc((unsigned) space);
	 
	sc_proj2 = sc_dup(sc_B_prime_2);
	// Add redundant constraints on intermediate loop indices to improve the generated code
	sc_proj2 = sc_projection_concat_proj_on_variables(sc_proj2,new_basis);
	sc_integer_projection_information(sc_proj2,new_basis, sc_info,sc_proj2->dimension,sc_proj2->dimension);
	sc_proj2=build_integer_sc_nredund(sc_proj2,new_basis,sc_info,1,sc_proj2->dimension,sc_proj2->dimension);

	// Distribute constraints according to loop nest indices of new_basis
	for (i=1;i<=sc_proj2->dimension;i++)
	  list_of_systems[i] = sc_init_with_sc(sc_proj2);
	constraint_distribution(sc_proj2,list_of_systems,new_basis,sc_info);

	// Remove redundant constraints according to the cumulated constraints from outer to inner loop bounds
	sc_tmp=sc_dup(list_of_systems[1]);
	for (i = 2; sc_tmp != NULL && i <=vect_size(new_basis);i++) {
	  list_of_systems[i] = elim_redund_sc_with_sc(list_of_systems[i],
						      sc_tmp,
						      new_basis,
						      vect_size(new_basis));
	  ifdebug(8) {
	    fprintf(stdout," Constraints on the %d-th var. :\n",i);
	    sc_fprint(stdout,list_of_systems[i],(get_variable_name_t)entity_user_name);
	  }
	  sc_tmp=sc_intersection(sc_tmp,sc_tmp,list_of_systems[i]);	
	}
	sc_rm(sc_tmp);
     
      // Build  system sc_tile to scan one tile
	sc_tile =sc_dup(list_of_systems[1]);
	for (i=2; i<=vect_size(tile_basis);i++)
	  sc_tile =sc_intersection(sc_tile,sc_tile,list_of_systems[i]);
	ifdebug(8) {
	  (void) fprintf(stderr,"Tile domain in echelon format:\n");
	  sc_fprint(stderr, sc_tile, (get_variable_name_t)entity_user_name);
	}
	// Constraint system sc_tile to scan one tile
	sc_local_tile =sc_dup(list_of_systems[i]);
	i++;
	for (; i<=vect_size(tile_basis)+vect_size(local_tile_basis);i++)
	  sc_local_tile =sc_intersection(sc_local_tile,sc_local_tile,list_of_systems[i]);
	ifdebug(8) {
	  (void) fprintf(stderr,"Iteration domain for one tile:\n");
	  sc_fprint(stderr, sc_local_tile, (get_variable_name_t)entity_user_name);
	}
      }
      
      /* generation of code for scanning one tile */
      ifdebug(8) {
	(void) fprintf(stderr,"The local coordonate change vector is:");
	for (int i=1; i<=n; i++)
	  vect_fprint(stderr, pvch[i], (get_variable_name_t)entity_user_name);
      }
      /* generation of code to scan one tile and update of loop body using pvch */
      s_lhyp = code_generation(lls, pvch, initial_basis,
			       local_tile_basis, sc_local_tile, false);

     /* generation of code for scanning all tiles */
     reverse_tile_basis = base_reversal(tile_basis);
     for (pb = reverse_tile_basis; pb!=NULL; pb=pb->succ) {
       loop tl = loop_undefined;

       make_bound_expression(pb->var, tile_basis, sc_tile, &lower, &upper);
       tl = make_loop((entity) vecteur_var(pb),
		      make_range(copy_expression(lower), copy_expression(upper),
				 int_to_expression(1)),
		      s_lhyp,
		      entity_empty_label(),
		      make_execution(is_execution_sequential, UU),
		      NIL);
       s_lhyp = instruction_to_statement(make_instruction(is_instruction_loop, tl));
     }
     stmf =s_lhyp; 
    }
    else pips_user_warning("PIPS legality test for Tiling transformation NOT VERIFIED, Verify data dependencies or simplify array access functions\n");
    
    pips_debug(8,"End\n");

    debug_off();
    reset_ordering_to_statement();
    return(stmf);
}

bool parallel_loop_tiling(const char* module_name)
{
  bool return_status = false;
    return_status = interactive_loop_transformation(module_name, parallel_tiling);
    return return_status;
}