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