ranking.c

Go to the documentation of this file.

/*
    This file is part of PolyLib.

    PolyLib 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
    (at your option) any later version.

    PolyLib 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 PolyLib.  If not, see <http://www.gnu.org/licenses/>.
*/

/**
 * Tools to compute the ranking function of an iteration J: the number of
 * integer points in P that are lexicographically inferior to J 
 * B. Meister
 * 6/2005
 * LSIIT-ICPS, UMR 7005 CNRS Universit� Louis Pasteur
 * HiPEAC Network
 */

#include <polylib/polylib.h>
#include <polylib/ranking.h>

/**
 * Returns a list of polytopes needed to compute
 * the number of points in P that are lexicographically
 * smaller than a given point in D.
 * Only the first dim dimensions are taken into account
 * for computing the lexsmaller relation.
 * The remaining variables are assumed to be extra
 * existential/control variables.
 * When P == D, this is the conventional ranking function.
 * P and D are assumed to have the same parameter domain C.
 *
 * The first polyhedron in the list returned is the
 * updated context: a combination of D and C or an extended C.
 *
 * The order of the variables in the remaining polyhedra is
 * - first dim variables of P
 * - existential variables of P
 * - existential variables of D
 * - first dim variables of D
 * - the parameters
 */
Polyhedron *LexSmaller(Polyhedron *P, Polyhedron *D, unsigned dim,
                        Polyhedron *C, unsigned MAXRAYS)
{
  unsigned i, j, k, r;
  unsigned nb_parms = C->Dimension;
  unsigned nb_vars = dim;
  unsigned P_extra = P->Dimension - nb_vars - nb_parms;
  unsigned D_extra = D->Dimension - nb_vars - nb_parms;
  unsigned nb_new_parms;
  unsigned ncons;
  Matrix * cur_element, * C_times_J, * Klon;
  Polyhedron * P1, *C1;
  Polyhedron * lexico_lesser_union = NULL;

  POL_ENSURE_INEQUALITIES(C);
  POL_ENSURE_INEQUALITIES(D);
  POL_ENSURE_INEQUALITIES(P);

  assert(P->Dimension >= C->Dimension + dim);
  assert(D->Dimension >= C->Dimension + dim);
  nb_new_parms = nb_vars;

  /* the number of variables must be positive */
  if (nb_vars<=0) {
    printf("\nRanking > No variables, returning NULL.\n"); 
    return NULL;
  }
  /*
   * if D has extra variables, then we can't squeeze the contraints
   * of D in the new context, so we simply add them to each element.
   */
  if (D_extra)
    cur_element = Matrix_Alloc(P->NbConstraints+D->NbConstraints+nb_new_parms, 
                               P->Dimension+D_extra+nb_new_parms+2);
  else
    cur_element = Matrix_Alloc(P->NbConstraints+nb_new_parms, 
                               P->Dimension+D_extra+nb_new_parms+2);


  /* 0- Put P in the first rows of cur_element */
  for (i=0; i < P->NbConstraints; i++) {
    Vector_Copy(P->Constraint[i], cur_element->p[i], nb_vars+P_extra+1);
    Vector_Copy(P->Constraint[i]+1+nb_vars+P_extra, 
                cur_element->p[i]+1+nb_vars+P_extra+D_extra+nb_new_parms, 
                nb_parms+1);
  }
  ncons = P->NbConstraints;
  if (D_extra) {
    for (i=0; i < D->NbConstraints; i++) {
      r = P->NbConstraints + i;
      Vector_Copy(D->Constraint[i], cur_element->p[r], 1);
      Vector_Copy(D->Constraint[i]+1, 
                  cur_element->p[r]+1+nb_vars+P_extra+D_extra, nb_new_parms);
      Vector_Copy(D->Constraint[i]+1+nb_new_parms, 
                  cur_element->p[r]+1+nb_vars+P_extra, D_extra);
      Vector_Copy(D->Constraint[i]+1+nb_new_parms+D_extra, 
                  cur_element->p[r]+1+nb_vars+P_extra+D_extra+nb_new_parms, 
                  nb_parms+1);
    }
    ncons += D->NbConstraints;
  }

  /* 1- compute the Ehrhart polynomial of each disjoint polyhedron defining the
     lexicographic order */
  for (k=0, r = ncons; k < nb_vars; k++, r++) {

    /* a- build the corresponding matrix
     *  the nb of rows of cur_element is fake, so that we do not have to
     *  re-allocate it. */
    cur_element->NbRows = r+1;

    /* convert the previous (strict) inequality into an equality */
    if (k>=1) {
      value_set_si(cur_element->p[r-1][0], 0);
      value_set_si(cur_element->p[r-1][cur_element->NbColumns-1], 0);
    }
    /* build the k-th inequality from P */
    value_set_si(cur_element->p[r][0], 1);
    value_set_si(cur_element->p[r][k+1], -1);
    value_set_si(cur_element->p[r][nb_vars+P_extra+D_extra+k+1], 1);
    /* we want a strict inequality */
    value_set_si(cur_element->p[r][cur_element->NbColumns-1], -1);
#ifdef ERDEBUG
    show_matrix(cur_element);
#endif

    /* b- add it to the current union
       as Constraints2Polyhedron modifies its input, we must clone cur_element */
    Klon = Matrix_Copy(cur_element);
    P1 = Constraints2Polyhedron(Klon, MAXRAYS);
    Matrix_Free(Klon);
    P1->next = lexico_lesser_union;
    lexico_lesser_union = P1;
  }
  
  /* 2- as we introduce n parameters, we must introduce them into the context
   * as well. 
   * The added constraints are P.M.(J N 1 )^T >=0 */
  if (D_extra)
    C_times_J = Matrix_Alloc(C->NbConstraints, nb_new_parms+nb_parms+2);
  else
    C_times_J = Matrix_Alloc(C->NbConstraints + D->NbConstraints, D->Dimension+2);
  /* copy the initial context while adding the new parameters */
  for (i = 0; i < C->NbConstraints; i++) {
    value_assign(C_times_J->p[i][0], C->Constraint[i][0]);
    Vector_Copy(C->Constraint[i]+1, C_times_J->p[i]+1+nb_new_parms, nb_parms+1);
  }

  /* copy constraints from evaluation domain */
  if (!D_extra)
    for (i = 0; i < D->NbConstraints; i++)
      Vector_Copy(D->Constraint[i], C_times_J->p[C->NbConstraints+i], 
                  D->Dimension+2);

#ifdef ERDEBUG
  show_matrix(C_times_J);
#endif
  C1 = Constraints2Polyhedron(C_times_J, POL_NO_DUAL);

  /* 4- clean up */
  Matrix_Free(cur_element);
  Matrix_Free(C_times_J);

  C1->next = P1;

  return C1;
} /* LexSmaller */


/**
 * Returns the number of points in P that are lexicographically
 * smaller than a given point in D.
 * Only the first dim dimensions are taken into account
 * for computing the lexsmaller relation.
 * The remaining variables are assumed to be extra
 * existential/control variables.
 * When P == D, this is the conventional ranking function.
 * P and D are assumed to have the same parameter domain C.
 * The variables in the Enumeration correspond to the first dim variables
 * in D followed by the parameters of D (the variables of C).
 */
Enumeration *Polyhedron_LexSmallerEnumerate(Polyhedron *P, Polyhedron *D, 
                                            unsigned dim,
                                            Polyhedron *C, unsigned MAXRAYS)
{
  Enumeration * ranking;
  Polyhedron *RC, *RD;

  RC = LexSmaller(P, D, dim, C, MAXRAYS);
  RD = RC->next;
  RC->next = NULL;

  /* Compute the ranking, which is the sum of the Ehrhart polynomials of the n
     disjoint polyhedra we just put in P1. */
  /* OPT : our polyhdera are (already) disjoint, so Domain_Enumerate does
     probably too much work uselessly */
  ranking = Domain_Enumerate(RD, RC, MAXRAYS, NULL);

  Domain_Free(RD);
  Polyhedron_Free(RC);

  return ranking;
}


/*
 * Returns a function that assigns a unique number to each point in the
 * polytope P ranging from zero to (number of points in P)-1.
 * The order of the numbers corresponds to the lexicographical order.
 *
 * C is the parameter context of the polytope
 */
Enumeration *Polyhedron_Ranking(Polyhedron *P, Polyhedron *C, unsigned MAXRAYS)
{
    return Polyhedron_LexSmallerEnumerate(P, P, P->Dimension-C->Dimension, 
                                          C, MAXRAYS);
}

---