source/tab.c

   1 /******************************************************************************
   2  *                     PIP : Parametric Integer Programming                   *
   3  ******************************************************************************
   4  *                                   tab.h                                    *
   5  ******************************************************************************
   6  *                                                                            *
   7  * Copyright Paul Feautrier, 1988, 1993, 1994, 1996, 2002                     *
   8  *                                                                            *
   9  * This is free software; you can redistribute it and/or modify it under the  *
  10  * terms of the GNU General Public License as published by the Free Software  *
  11  * Foundation; either version 2 of the License, or (at your option) any later *
  12  * version.                                                                   *
  13  *                                                                            *
  14  * This software is distributed in the hope that it will be useful, but       *
  15  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
  16  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
  17  * for more details.                                                          *
  18  *                                                                            *
  19  * You should have received a copy of the GNU General Public License along    *
  20  * with software; if not, write to the Free Software Foundation, Inc.,        *
  21  * 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA                     *
  22  *                                                                            *
  23  * Written by Paul Feautrier and Cedric Bastoul                               *
  24  *                                                                            *
  25  *****************************************************************************/
  26
  27
  28 #include <stdio.h>
  29 #include <stdlib.h>
  30
  31 #include "pip.h"
  32
  33 #define TAB_CHUNK 4096*sizeof(Entier)
  34
  35 static char *tab_free, *tab_top;
  36 static struct A *tab_base;
  37
  38 extern int allocation;
  39 extern long int cross_product, limit;
  40 static int chunk_count;
  41
  42 int dgetc(FILE *);
  43 #if defined(LINEAR_VALUE_IS_MP)
  44 int dscanf(FILE *, Entier);
  45 #else
  46 int dscanf(FILE *, Entier *);
  47 #endif
  48
  49 extern FILE * dump;
  50
  51 #define sizeof_struct_A ((sizeof(struct A) % sizeof(Entier)) ?              \
  52                          (sizeof(struct A) + sizeof(Entier)                 \
  53                                 - (sizeof(struct A) % sizeof(Entier))) :    \
  54                           sizeof(struct A))
  55
  56 void tab_init(void)
  57 {
  58  tab_free = malloc(sizeof_struct_A);
  59  if(tab_free == NULL)
  60      {fprintf(stderr, "Your computer doesn't have enough memory\n");
  61       exit(1);
  62      }
  63  allocation = 1;
  64  tab_top = tab_free + sizeof_struct_A;
  65  tab_base = (struct A *)tab_free;
  66  tab_free += sizeof_struct_A;
  67  tab_base->precedent = NULL;
  68  tab_base->bout = tab_top;
  69  tab_base->free = tab_free;
  70  chunk_count = 1;
  71 }
  72
  73
  74 void tab_close(void)
  75 {
  76   if (tab_base) free(tab_base);
  77 }
  78
  79
  80 struct high_water_mark tab_hwm(void)
  81 {struct high_water_mark p;
  82  p.chunk = chunk_count;
  83  p.top = tab_free;
  84  return p;
  85 }
  86
  87
  88 #if defined(LINEAR_VALUE_IS_MP)
  89 /* the clear_tab routine clears the GMP objects which may be referenced
  90    in the given Tableau.
  91 */
  92 void tab_clear(Tableau *tp)
  93 {
  94   int i, j;
  95   /* clear the determinant */
  96   mpz_clear(tp->determinant);
  97
  98   for(i=0; i<tp->height; i++){
  99     /* clear the denominator */
 100     mpz_clear(Denom(tp, i));
 101     if((Flag(tp, i) & Unit) == 0)
 102       for(j=0; j<tp->width; j++)
 103         mpz_clear(Index(tp,i,j));
 104   }
 105 }
 106 #endif
 107
 108 void tab_reset(struct high_water_mark by_the_mark)
 109
 110 {struct A *g;
 111  char *p;
 112  while(chunk_count > by_the_mark.chunk)
 113      {
 114       g = tab_base->precedent;
 115
 116       #if defined(LINEAR_VALUE_IS_MP)
 117       /* Before actually freeing the memory, one has to clear the
 118        * included Tableaux. If this is not done, the GMP objects
 119        * referenced in the Tableaux will be orphaned.
 120        */
 121
 122       /* Enumerate the included tableaux. */
 123       p = (char *)tab_base + sizeof_struct_A;
 124       while(p < tab_base->free){
 125         Tableau *pt;
 126         pt = (Tableau *) p;
 127         tab_clear(pt);
 128         p += pt->taille;
 129       }
 130       #endif
 131
 132       free(tab_base);
 133       tab_base = g;
 134       tab_top = tab_base->bout;
 135       chunk_count--;
 136      }
 137  if(chunk_count > 0) {
 138      #if defined(LINEAR_VALUE_IS_MP)
 139      /* Do not forget to clear the tables in the current chunk above the
 140         high water mark */
 141      p = (char *)by_the_mark.top;
 142      while(p < tab_base->free) {
 143         Tableau *pt;
 144         pt = (Tableau *) p;
 145         tab_clear(pt);
 146         p += pt->taille;
 147         }
 148      #endif
 149      tab_free = by_the_mark.top;
 150      tab_base->free = tab_free;
 151      }
 152  else {
 153      fprintf(stderr, "Syserr: tab_reset : error in memory allocation\n");
 154      exit(1);
 155      }
 156 }
 157
 158 Tableau * tab_alloc(int h, int w, int n)
 159
 160 /* h : le nombre de ligne reelles;
 161    n : le nombre de lignes virtuelles
 162 */
 163 {
 164  char *p; Tableau *tp;
 165  Entier *q;
 166  unsigned long taille;
 167  int i, j;
 168  taille = sizeof(struct T) + (h+n-1) * sizeof (struct L)
 169           + h * w * sizeof (Entier);
 170  if(tab_free + taille >= tab_top)
 171      {struct A * g;
 172       unsigned long d;
 173       d = taille + sizeof_struct_A;
 174       if(d < TAB_CHUNK) d = TAB_CHUNK;
 175       tab_free = malloc(d);
 176       if(tab_free == NULL)
 177           {printf("Memory overflow\n");
 178            exit(23);
 179           }
 180       chunk_count++;
 181       g = (struct A *)tab_free;
 182       g->precedent = tab_base;
 183       tab_top = tab_free + d;
 184       tab_free += sizeof_struct_A;
 185       tab_base = g;
 186       g->bout = tab_top;
 187      }
 188  p = tab_free;
 189  tab_free += taille;
 190  tab_base->free = tab_free;
 191  tp = (Tableau *)p;
 192  q = (Entier *)(p +  sizeof(struct T) + (h+n-1) * sizeof (struct L));
 193  #if defined(LINEAR_VALUE_IS_MP)
 194  mpz_init_set_ui(tp->determinant,1);
 195  #else
 196  tp->determinant[0] = (Entier) 1;
 197  tp->l_determinant = 1;
 198  #endif
 199  for(i = 0; i<n ; i++){
 200    tp->row[i].flags = Unit;
 201    tp->row[i].objet.unit = i;
 202    #if defined(LINEAR_VALUE_IS_MP)
 203    mpz_init_set_ui(Denom(tp, i), 1);
 204    #else
 205    Denom(tp, i) = UN ;
 206    #endif
 207  }
 208  for(i = n; i < (h+n); i++){
 209    tp->row[i].flags = 0;
 210    tp->row[i].objet.val = q;
 211    for(j = 0; j < w; j++)
 212    #if defined(LINEAR_VALUE_IS_MP)
 213    mpz_init_set_ui(*q++, 0); /* loop body. */
 214    mpz_init_set_ui(Denom(tp, i), 0);
 215    #else
 216    *q++ = 0;                 /* loop body. */
 217    Denom(tp, i) = ZERO ;
 218    #endif
 219  }
 220  tp->height = h + n; tp->width = w;
 221  #if defined(LINEAR_VALUE_IS_MP)
 222  tp->taille = taille ;
 223  #endif
 224
 225  return(tp);
 226 }
 227
 228 Tableau * tab_get(foo, h, w, n)
 229 FILE * foo;
 230 int h, w, n;
 231 {
 232  Tableau *p;
 233  int i, j, c;
 234  Entier x;
 235  #if defined(LINEAR_VALUE_IS_MP)
 236  mpz_init(x);
 237  #endif
 238
 239  p = tab_alloc(h, w, n);
 240  while((c = dgetc(foo)) != EOF)
 241       if(c == '(')break;
 242  for(i = n; i<h+n; i++)
 243      {p->row[i].flags = Unknown;
 244       #if defined(LINEAR_VALUE_IS_MP)
 245       mpz_set_ui(Denom(p, i), 1);
 246       #else
 247       Denom(p, i) = UN;
 248       #endif
 249       while((c = dgetc(foo)) != EOF)if(c == '[')break;
 250       for(j = 0; j<w; j++){
 251         #if defined(LINEAR_VALUE_IS_MP)
 252         if(dscanf(foo, x) < 0)
 253           return NULL;
 254         else
 255           mpz_set(p->row[i].objet.val[j], x);
 256         #else
 257         if(dscanf(foo, &x) < 0)
 258           return NULL;
 259         else
 260           p->row[i].objet.val[j] = x;
 261         #endif
 262         }
 263       }
 264       while((c = dgetc(foo)) != EOF)if(c == ']')break;
 265
 266  #if defined(LINEAR_VALUE_IS_MP)
 267  mpz_clear(x);
 268  #endif
 269
 270  return(p);
 271 }
 272
 273
 274 /* Fonction tab_Matrix2Tableau :
 275  * Cette fonction effectue la conversion du format de matrice de la polylib
 276  * vers le format de traitement de Pip. matrix est la matrice a convertir.
 277  * Nineq est le nombre d'inequations necessaires (dans le format de la
 278  * polylib, le premier element d'une ligne indique si l'equation decrite
 279  * est une inequation ou une egalite. Pip ne gere que les inequations. On
 280  * compte donc le nombre d'inequations total pour reserver la place
 281  * necessaire, et on scinde toute egalite p(x)=0 en p(x)>=0 et -p(x)>=0).
 282  * Nv est le nombre de variables dans la premiere serie de variables (c'est
 283  * a dire que si les premiers coefficients dans les lignes de la matrice
 284  * sont ceux des inconnues, Nv est le nombre d'inconnues, resp. parametres).
 285  * n est le nombre de lignes 'virtuelles' contenues dans la matrice (c'est
 286  * a dire en fait le nombre d'inconnues). Si Shift vaut 0, on va rechercher
 287  * le minimum lexicographique non-negatif, sinon on recherche le maximum
 288  * (Shift = 1) ou bien le minimum tout court (Shift = -1). La fonction
 289  * met alors en place le bignum s'il n'y est pas deja et prepare les
 290  * contraintes au calcul du maximum lexicographique.
 291  * 27 juillet 2001 : Premiere version, Ced.
 292  * 30 juillet 2001 : Nombreuses modifications. Le calcul du nombre total
 293  *                   d'inequations (Nineq) se fait a present a l'exterieur.
 294  *  3 octobre 2001 : Pas mal d'ameliorations.
 295  * 18 octobre 2003 : Mise en place de la possibilite de calculer le
 296  *                   maximum lexicographique (parties 'if (Max)').
 297  */
 298 Tableau * tab_Matrix2Tableau(matrix, Nineq, Nv, n, Shift, Bg, Urs_parms)
 299 PipMatrix * matrix ;
 300 int Nineq, Nv, n, Shift, Bg, Urs_parms;
 301 { Tableau * p ;
 302   unsigned i, j, k, current, new, nb_columns, decal=0, bignum_is_new ;
 303   int inequality;
 304   Entier bignum;
 305
 306   entier_init(bignum);
 307   nb_columns = matrix->NbColumns - 1 ;
 308   /* S'il faut un BigNum et qu'il n'existe pas, on lui reserve sa place. */
 309   bignum_is_new = Shift && (Bg > (matrix->NbColumns - 2));
 310   if (bignum_is_new)
 311     nb_columns++;
 312   /* Ce sont juste des parametres. */
 313   if (Bg <= Nv)
 314     Shift = 0;
 315
 316   p = tab_alloc(Nineq,nb_columns+Urs_parms,n) ;
 317
 318   /* La variable decal sert a prendre en compte les lignes supplementaires
 319    * issues des egalites.
 320    */
 321   for (i = 0; i < matrix->NbRows; i++) {
 322     current = i + n + decal;
 323     Flag(p,current) = Unknown ;
 324     entier_set_si(Denom(p,current), 1);
 325     if (Shift)
 326       entier_set_si(bignum, 0);
 327     /* Pour passer l'indicateur d'egalite/inegalite. */
 328     inequality = entier_notzero_p(matrix->p[i][0]);
 329
 330     /* Dans le format de la polylib, l'element constant est place en
 331      * dernier. Dans le format de Pip, il se trouve apres la premiere
 332      * serie de variables (inconnues ou parametres). On remet donc les
 333      * choses dans l'ordre de Pip. Ici pour p(x) >= 0.
 334      */
 335     for (j=0;j<Nv;j++) {
 336       if (bignum_is_new && 1+j == Bg)
 337         continue;
 338       if (Shift)
 339         entier_addto(bignum, bignum, matrix->p[i][1+j]);
 340       if (Shift > 0)
 341         entier_oppose(p->row[current].objet.val[j], matrix->p[i][1+j]);
 342       else
 343         entier_assign(p->row[current].objet.val[j], matrix->p[i][1+j]);
 344     }
 345     for (k=j=Nv+1;j<nb_columns;j++) {
 346         if (bignum_is_new && j == Bg)
 347           continue;
 348         entier_assign(p->row[current].objet.val[j], matrix->p[i][k++]);
 349     }
 350     for (j=0; j < Urs_parms; ++j) {
 351         int pos = nb_columns - Urs_parms + j;
 352         if (pos <= Nv)
 353             --pos;
 354         if (pos <= Bg)
 355             --pos;
 356         entier_oppose(p->row[current].objet.val[nb_columns+j],
 357                      p->row[current].objet.val[pos]);
 358     }
 359     entier_assign(p->row[current].objet.val[Nv],
 360                  matrix->p[i][matrix->NbColumns-1]);
 361     if (Shift) {
 362       if (Shift < 0)
 363         entier_oppose(bignum, bignum);
 364
 365       if (bignum_is_new)
 366         entier_assign(p->row[current].objet.val[Bg], bignum);
 367       else
 368         entier_addto(p->row[current].objet.val[Bg],
 369                     p->row[current].objet.val[Bg], bignum);
 370     }
 371
 372     /* Et ici lors de l'ajout de -p(x) >= 0 quand on traite une egalite. */
 373     if (!inequality) {
 374       decal ++ ;
 375       new = current + 1 ;
 376       Flag(p,new)= Unknown ;
 377       entier_set_si(Denom(p,new), 1);
 378
 379       for (j=0;j<nb_columns+Urs_parms;j++)
 380         entier_oppose(p->row[new].objet.val[j], p->row[current].objet.val[j]);
 381     }
 382   }
 383   entier_clear(bignum);
 384
 385   return(p);
 386 }
 387
 388
 389 char *Attr[] = {"Unit", "+", "-", "0", "*", "?"};
 390
 391 void tab_display(p, foo)
 392 FILE *foo;
 393 Tableau *p;
 394 {
 395
 396  int i, j, ff, fff, n;
 397  Entier x, d;
 398  #if defined(LINEAR_VALUE_IS_MP)
 399  mpz_init(d);
 400  #endif
 401
 402  fprintf(foo, "%ld/[%d * %d]\n", cross_product, p->height, p->width);
 403  for(i = 0; i<p->height; i++){
 404    fff = ff = p->row[i].flags;
 405    /* if(fff ==0) continue; */
 406    #if defined(LINEAR_VALUE_IS_MP)
 407    mpz_set(d, Denom(p, i));
 408    #else
 409    d = Denom(p, i);
 410    #endif
 411    n = 0;
 412    while(fff){
 413      if(fff & 1) fprintf(foo, "%s ",Attr[n]);
 414      n++; fff >>= 1;
 415    }
 416    fprintf(foo, "%f #[", p->row[i].size);
 417    if(ff & Unit)
 418      for(j = 0; j<p->width; j++)
 419        fprintf(foo, " /%d/",(j == p->row[i].objet.unit)? 1: 0);
 420    else
 421      for(j = 0; j<p->width; j++){
 422        #if defined(LINEAR_VALUE_IS_MP)
 423        mpz_out_str(foo, 10, Index(p, i, j));
 424        putc(' ', foo);
 425        #else
 426        x = Index(p,i,j);
 427        fprintf(foo, FORMAT, x);
 428        fprintf(foo, " ");
 429        #endif
 430      }
 431    fprintf(foo, "]/");
 432    #if defined(LINEAR_VALUE_IS_MP)
 433    mpz_out_str(foo, 10, d);
 434    #else
 435    fprintf(foo, "%d", (int)d);
 436    #endif
 437    putc('\n', foo);
 438  }
 439  #if defined(LINEAR_VALUE_IS_MP)
 440  mpz_clear(d);
 441  #endif
 442 }