source/piplibMP.c

   1 /******************************************************************************
   2  *                     PIP : Parametric Integer Programming                   *
   3  ******************************************************************************
   4  *                                 piplib.c                                   *
   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 Cedric Bastoul                                                  *
  24  *                                                                            *
  25  ******************************************************************************/
  26
  27 /* Premiere version du 30 juillet 2001. */
  28
  29 # include <stdlib.h>
  30 # include <stdio.h>
  31 # include <ctype.h>
  32
  33 #include <piplib/piplib.h>
  34
  35 extern long int cross_product, limit ;
  36 extern int verbose ;
  37 extern FILE * dump ;
  38 extern int compa_count ;
  39 extern char dump_name[] ;
  40
  41
  42 /******************************************************************************
  43  *                    Fonctions d'affichage des structures                    *
  44  ******************************************************************************/
  45
  46
  47 /* Fonction pip_matrix_print :
  48  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
  49  * que contient la structure de type PipMatrix qu'elle recoit en parametre.
  50  * Premiere version : Ced. 29 juillet 2001.
  51  * 24 octobre 2002 : premiere version MP.
  52  */
  53 void pip_matrix_print(FILE * foo, PipMatrix * Mat)
  54 { Entier * p;
  55   int i, j ;
  56   unsigned NbRows, NbColumns ;
  57
  58   fprintf(foo,"%d %d\n", NbRows=Mat->NbRows, NbColumns=Mat->NbColumns) ;
  59   for (i=0;i<NbRows;i++)
  60   { p=*(Mat->p+i) ;
  61     for (j=0;j<NbColumns;j++)
  62     { fprintf(foo," ") ;
  63       mpz_out_str(foo,10,*p++) ;
  64     }
  65     fprintf(foo, "\n") ;
  66   }
  67 }
  68
  69
  70 /* Fonction pip_vector_print :
  71  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
  72  * que contient la structure de type PipVector qu'elle recoit en parametre.
  73  * Premiere version : Ced. 20 juillet 2001.
  74  * 24 octobre 2002 : premiere version MP.
  75  */
  76 void pip_vector_print(FILE * foo, PipVector * vector)
  77 { int i ;
  78
  79   if (vector != NULL)
  80   { fprintf(foo,"#[") ;
  81     for (i=0;i<vector->nb_elements;i++)
  82     { fprintf(foo," ") ;
  83       mpz_out_str(foo,10,vector->the_vector[i]) ;
  84       if (mpz_cmp(vector->the_deno[i],UN) != 0)
  85       { fprintf(foo,"/") ;
  86         mpz_out_str(foo,10,vector->the_deno[i]) ;
  87       }
  88     }
  89     fprintf(foo,"]") ;
  90   }
  91 }
  92
  93
  94 /* Fonction pip_newparm_print :
  95  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
  96  * que contient la structure de type PipNewparm qu'elle recoit en parametre.
  97  * Le parametre indent est le nombre d'espaces blancs en debut de chaque
  98  * ligne avant indentation. Une valeur negative de indent signifie qu'on ne
  99  * desire pas d'indentation.
 100  * Premiere version : Ced. 18 octobre 2001.
 101  * 24 octobre 2002 : premiere version MP.
 102  */
 103 void pip_newparm_print(FILE * foo, PipNewparm * newparm, int indent)
 104 { int i ;
 105
 106   if (newparm != NULL)
 107   { do
 108     { for (i=0;i<indent;i++) fprintf(foo," ") ;             /* Indent. */
 109       fprintf(foo,"(newparm ") ;
 110       fprintf(foo,"%d",newparm->rank) ;
 111       fprintf(foo," (div ") ;
 112       pip_vector_print(foo,newparm->vector) ;
 113       fprintf(foo," ") ;
 114       mpz_out_str(foo,10,newparm->deno) ;
 115       fprintf(foo,"))\n") ;
 116     }
 117     while ((newparm = newparm->next) != NULL) ;
 118   }
 119 }
 120
 121
 122 /* Fonction pip_list_print :
 123  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
 124  * que contient la structure de type PipList qu'elle recoit en parametre.
 125  * Le parametre indent est le nombre d'espaces blancs en debut de chaque
 126  * ligne avant indentation. Une valeur negative de indent signifie qu'on ne
 127  * desire pas d'indentation.
 128  * Premiere version : Ced. 18 octobre 2001.
 129  */
 130 void pip_list_print(FILE * foo, PipList * list, int indent)
 131 { int i ;
 132
 133   if (list == NULL)
 134   { for (i=0;i<indent;i++) fprintf(foo," ") ;               /* Indent. */
 135     fprintf(foo,"()\n") ;
 136   }
 137   else
 138   { for (i=0;i<indent;i++) fprintf(foo," ") ;               /* Indent. */
 139     fprintf(foo,"(list\n") ;
 140     do
 141     { for (i=0;i<indent+1;i++) fprintf(foo," ") ;           /* Indent. */
 142       pip_vector_print(foo,list->vector) ;
 143       fprintf(foo,"\n") ;
 144     }
 145     while ((list = list->next) != NULL) ;
 146     for (i=0;i<indent;i++) fprintf(foo," ") ;               /* Indent. */
 147     fprintf(foo,")\n") ;
 148   }
 149 }
 150
 151
 152 /* Fonction pip_quast_print :
 153  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
 154  * que contient la structure de type PipQuast qu'elle recoit en parametre.
 155  * Le parametre indent est le nombre d'espaces blancs en debut de chaque
 156  * ligne avant indentation. Une valeur negative de indent signifie qu'on ne
 157  * desire pas d'indentation.
 158  * 20 juillet 2001 : Premiere version, Ced.
 159  * 18 octobre 2001 : eclatement.
 160  */
 161 void pip_quast_print(FILE * foo, PipQuast * solution, int indent)
 162 { int i ;
 163   PipVector * vector ;
 164
 165   if (solution != NULL)
 166   { pip_newparm_print(foo,solution->newparm,indent) ;
 167     if (solution->condition == NULL)
 168     pip_list_print(foo,solution->list,indent) ;
 169     else
 170     { for (i=0;i<indent;i++) fprintf(foo," ") ;             /* Indent. */
 171       fprintf(foo,"(if ") ;
 172       pip_vector_print(foo,solution->condition) ;
 173       fprintf(foo,"\n") ;
 174       if (indent>=0)                                        /* Indent. */
 175       { pip_quast_print(foo,solution->next_then,indent+1) ;
 176         pip_quast_print(foo,solution->next_else,indent+1) ;
 177       }
 178       else
 179       { pip_quast_print(foo,solution->next_then,indent) ;
 180         pip_quast_print(foo,solution->next_else,indent) ;
 181       }
 182       for (i=0;i<indent;i++) fprintf(foo," ") ;             /* Indent. */
 183       fprintf(foo,")\n") ;
 184     }
 185   }
 186   else
 187   { for (i=0;i<indent;i++) fprintf(foo," ") ;               /* Indent. */
 188     fprintf(foo,"void\n") ;
 189   }
 190 }
 191
 192
 193 /******************************************************************************
 194  *                       Fonctions de liberation memoire                      *
 195  ******************************************************************************/
 196
 197
 198 /* Fonction pip_matrix_free :
 199  * Cette fonction libere la memoire reservee a la structure de type PipMatrix
 200  * que pointe son parametre.
 201  * Premiere version : Ced. 29 juillet 2001.
 202  * 24 octobre 2002 : premiere version MP.
 203  */
 204 void pip_matrix_free(PipMatrix * matrix)
 205 { int i, j ;
 206   Entier * p ;
 207
 208   p = matrix->p_Init ;
 209   for (i=0;i<matrix->NbRows;i++)
 210   for (j=0;j<matrix->NbColumns;j++)
 211   mpz_clear(*p++) ;
 212
 213   if (matrix != NULL)
 214   { free(matrix->p_Init) ;
 215     free(matrix->p) ;
 216     free(matrix) ;
 217   }
 218 }
 219
 220
 221 /* Fonction pip_vector_free :
 222  * Cette fonction libere la memoire reservee a la structure de type PipVector
 223  * que pointe son parametre.
 224  * 20 juillet 2001 : Premiere version, Ced.
 225  * 18 octobre 2001 : simplification suite a l'eclatement de PipVector.
 226  * 24 octobre 2002 : premiere version MP.
 227  */
 228 void pip_vector_free(PipVector * vector)
 229 { int i ;
 230
 231   for (i=0;i<vector->nb_elements;i++)
 232   { mpz_clear(vector->the_vector[i]);
 233     mpz_clear(vector->the_deno[i]);
 234   }
 235   free(vector->the_vector) ;
 236   free(vector->the_deno) ;
 237   free(vector) ;
 238 }
 239
 240
 241 /* Fonction pip_newparm_free :
 242  * Cette fonction libere la memoire reservee a la structure de type PipNewparm
 243  * que pointe son parametre. Sont liberes aussi tous les elements de la
 244  * liste chainee dont il pouvait etre le depart.
 245  * Premiere version : Ced. 18 octobre 2001.
 246  * 24 octobre 2002 : premiere version MP.
 247  */
 248 void pip_newparm_free(PipNewparm * newparm)
 249 { PipNewparm * next ;
 250
 251   while (newparm != NULL)
 252   { next = newparm->next ;
 253     mpz_clear(newparm->deno);
 254     pip_vector_free(newparm->vector) ;
 255     free(newparm) ;
 256     newparm = next ;
 257   }
 258 }
 259
 260
 261 /* Fonction pip_list_free :
 262  * Cette fonction libere la memoire reservee a la structure de type PipList
 263  * que pointe son parametre. Sont liberes aussi tous les elements de la
 264  * liste chainee dont il pouvait etre le depart.
 265  * Premiere version : Ced. 18 octobre 2001.
 266  */
 267 void pip_list_free(PipList * list)
 268 { PipList * next ;
 269
 270   while (list != NULL)
 271   { next = list->next ;
 272     pip_vector_free(list->vector) ;
 273     free(list) ;
 274     list = next ;
 275   }
 276 }
 277
 278
 279 /* Fonction pip_quast_free :
 280  * Cette fonction libere la memoire reservee a la structure de type
 281  * PipSolution que pointe son parametre. Sont liberees aussi toutes les
 282  * differentes listes chainees qui pouvaient en partir.
 283  * 20 juillet 2001 : Premiere version, Ced.
 284  * 18 octobre 2001 : simplification suite a l'eclatement de PipVector.
 285  */
 286 void pip_quast_free(PipQuast * solution)
 287 { if (solution != NULL)
 288   { if (solution->newparm != NULL)
 289     pip_newparm_free(solution->newparm) ;
 290
 291     if (solution->list != NULL)
 292     pip_list_free(solution->list) ;
 293
 294     if (solution->condition != NULL)
 295     { pip_vector_free(solution->condition) ;
 296       pip_quast_free(solution->next_then) ;
 297       pip_quast_free(solution->next_else) ;
 298     }
 299     free(solution) ;
 300   }
 301 }
 302
 303
 304 /******************************************************************************
 305  *                     Fonctions d'acquisition de matrices                    *
 306  ******************************************************************************/
 307
 308
 309 /* Fonction pip_matrix_alloc :
 310  * Fonction (tres) modifiee de Matrix_Alloc de la polylib. Elle alloue l'espace
 311  * memoire necessaire pour recevoir le contenu d'une matrice de NbRows lignes
 312  * et de NbColumns colonnes, et initialise les valeurs a 0. Elle retourne un
 313  * pointeur sur l'espace memoire alloue.
 314  * Premiere version : Ced. 18 octobre 2001.
 315  * 24 octobre 2002 : premiere version MP.
 316  */
 317 PipMatrix * pip_matrix_alloc(unsigned NbRows, unsigned NbColumns)
 318 { PipMatrix * matrix ;
 319   Entier ** p, * q ;
 320   int i, j ;
 321
 322   matrix = (PipMatrix *)malloc(sizeof(PipMatrix)) ;
 323   if (matrix == NULL)
 324   { fprintf(stderr, "Memory Overflow.\n") ;
 325     exit(1) ;
 326   }
 327   matrix->NbRows = NbRows ;
 328   matrix->NbColumns = NbColumns ;
 329   if (NbRows == 0)
 330   { matrix->p = NULL ;
 331     matrix->p_Init = NULL ;
 332   }
 333   else
 334   { if (NbColumns == 0)
 335     { matrix->p = NULL ;
 336       matrix->p_Init = NULL ;
 337     }
 338     else
 339     { p = (Entier **)malloc(NbRows*sizeof(Entier *)) ;
 340       if (p == NULL)
 341       { fprintf(stderr, "Memory Overflow.\n") ;
 342         exit(1) ;
 343       }
 344       q = (Entier *)malloc(NbRows * NbColumns * sizeof(Entier)) ;
 345       if (q == NULL)
 346       { fprintf(stderr, "Memory Overflow.\n") ;
 347         exit(1) ;
 348       }
 349       matrix->p = p ;
 350       matrix->p_Init = q ;
 351       for (i=0;i<NbRows;i++)
 352       { *p++ = q ;
 353         for (j=0;j<NbColumns;j++)
 354         mpz_init_set_si(*(q+j),0) ;
 355         q += NbColumns ;
 356       }
 357     }
 358   }
 359   return matrix ;
 360 }
 361
 362
 363 /* Fonction pip_matrix_read :
 364  * Adaptation de Matrix_Read de la polylib. Cette fonction lit les donnees
 365  * d'une matrice dans un fichier 'foo' et retourne un pointeur vers une
 366  * structure ou elle a copie les informations de cette matrice. Les donnees
 367  * doivent se presenter tq :
 368  * - des lignes de commentaires commencant par # (optionnelles),
 369  * - le nombre de lignes suivit du nombre de colonnes de la matrice, puis
 370  *   eventuellement d'un commentaire sur une meme ligne,
 371  * - des lignes de la matrice, chaque ligne devant etre sur sa propre ligne de
 372  *   texte et eventuellement suivies d'un commentaire.
 373  * Premiere version : Ced. 18 octobre 2001.
 374  * 24 octobre 2002 : premiere version MP, attention, uniquement capable de
 375  *                   lire des long long pour l'instant. On utilise pas
 376  *                   mpz_inp_str car on lit depuis des char * et non des FILE.
 377  */
 378 PipMatrix * pip_matrix_read(FILE * foo)
 379 { unsigned NbRows, NbColumns ;
 380   int i, j, n ;
 381   long long val ;
 382   char *c, s[1024], str[1024] ;
 383   PipMatrix * matrix ;
 384   Entier * p ;
 385
 386   while (fgets(s,1024,foo) == 0) ;
 387   while ((*s=='#' || *s=='\n') || (sscanf(s," %u %u",&NbRows,&NbColumns)<2))
 388   fgets(s, 1024, foo) ;
 389
 390   matrix = pip_matrix_alloc(NbRows,NbColumns) ;
 391
 392   p = matrix->p_Init ;
 393   for (i=0;i<matrix->NbRows;i++)
 394   { do
 395     { c = fgets(s,1024,foo) ;
 396       while (isspace(*c) && (*c != '\n'))
 397       c++ ;
 398     }
 399     while (c != NULL && (*c == '#' || *c == '\n'));
 400
 401     if (c == NULL)
 402     { fprintf(stderr, "Not enough rows.\n") ;
 403       exit(1) ;
 404     }
 405     for (j=0;j<matrix->NbColumns;j++)
 406     { if (c == NULL || *c == '#' || *c == '\n')
 407       { fprintf(stderr, "Not enough columns.\n") ;
 408         exit(1) ;
 409       }
 410       /* NdCed : Dans le n ca met strlen(str). */
 411       if (sscanf(c,"%s%n",str,&n) == 0)
 412       { fprintf(stderr, "Not enough rows.\n") ;
 413         exit(1) ;
 414       }
 415       sscanf(str,FORMAT,&val) ;
 416       mpz_init_set_si(*p++,val) ;
 417       c += n ;
 418     }
 419   }
 420   return matrix ;
 421 }
 422
 423
 424 /******************************************************************************
 425  *                           Fonction de resolution                           *
 426  ******************************************************************************/
 427
 428
 429 /* Fonction pip_solve :
 430  * Cette fonction fait un appel a Pip pour la resolution d'un probleme. Le
 431  * probleme est fourni dans les arguments. Deux matrices de la forme de celles
 432  * utilisees dans la Polylib forment les systemes d'equations/inequations :
 433  * un pour les inconnues, l'autre pour les parametres. Bg est le 'bignum', Nq
 434  * est un booleen renseignant si on cherche une solution entiere (vrai=1) ou
 435  * non (faux=0). Verbose est un booleen permettant de rendre Pip bavard
 436  * (Verbose a vrai=1), il imprimera alors la plupart de ses traitements dans le
 437  * fichier dont le nom est dans la variable d'environnement DEBUG, ou si DEBUG
 438  * n'est pas placee, dans un nouveau fichier de nom genere par mkstemp, si
 439  * Verbose est a faux=0, Pip restera muet. Simplify est un booleen permettant
 440  * de demander a Pip de simplifier sa solution (en eliminant les formes de type
 441  * 'if #[...] () ()') quand il est a vrai=1, ou non quand il est a faux=0. Max
 442  * n'est pas encore utilise et doit etre mis a 0. Cette fonction retourne la
 443  * solution sous la forme d'un arbre de structures PipQuast.
 444  * 30 juillet 2001 : Premiere version, Ced.
 445  * 18 octobre 2001 : suppression de l'argument Np, le nombre de parametres. Il
 446  *                   est a present deduit de ineqpar. Si ineqpar vaut NULL,
 447  *                   c'est que Np vaut 0. S'il y a des parametres mais pas de
 448  *                   contraintes dessus, ineqpar sera une matrice de 0 lignes
 449  *                   mais du bon nombre de colonnes (Np + 2).
 450  * 24 octobre 2002 : premiere version MP.
 451  */
 452 PipQuast * pip_solve(inequnk, ineqpar, Bg, Nq, Verbose, Simplify, Max)
 453 PipMatrix * inequnk, * ineqpar ;
 454 int Bg, Nq, Verbose, Simplify, Max ;
 455 { Tableau * ineq, * context, * ctxt ;
 456   int i, Np, Nn, Nl, Nm, p, q, xq, non_vide ;
 457   char * g ;
 458   struct high_water_mark hq ;
 459   PipQuast * solution ;
 460
 461   mpz_init_set_si(UN, 1);
 462   mpz_init_set_si(ZERO, 0);
 463
 464   /* initialisations diverses :
 465    * - la valeur de Verbose est placee dans sa variable globale. Dans le cas
 466    *   ou on doit etre en mode verbose, on ouvre le fichier dans lequel
 467    *   ecrire les tracages. Si la variable d'environnement DEBUG est placee,
 468    *   on ecrira dans le nom de fichier correspondant, sinon, dans un nouveau
 469    *   fichier de nom genere par mkstemp,
 470    * - limit est mis au 0 long int (sa valeur par defaut dans Pip original),
 471    * - on lance les initialisations pour tab et sol (autres mises en place
 472    *   de variables globales).
 473    */
 474   verbose = Verbose ;
 475   if (verbose)
 476   { g = getenv("DEBUG") ;
 477     if(g && *g)
 478     { dump = fopen(g, "w") ;
 479       if(dump == NULL)
 480       { fprintf(stderr,"%s unaccessible\n",g) ;
 481         verbose = 0 ;
 482       }
 483     }
 484     else
 485     { mkstemp(dump_name) ;
 486       dump = fopen(dump_name, "w") ;
 487     }
 488   }
 489   limit = 0LL ;
 490   sol_init() ;
 491   tab_init() ;
 492
 493   /* Si inequnk est NULL, la solution est automatiquement void (NULL). */
 494   if (inequnk != NULL)
 495   { /* Np vaut 0 si ineqpar vaut NULL, ineqpar->NbColumns - 2 sinon (on a -1
 496      * pour la constante et -1 pour le marqueur d'egalite/inegalite = -2).
 497      */
 498     Np = (ineqpar == NULL) ? 0 : ineqpar->NbColumns - 2 ;
 499     /* Calcul du nombre d'inconnues du probleme. Comme les matrices d'entree
 500      * sont sous la forme de la polylib.
 501      */
 502     Nn = inequnk->NbColumns - Np - 2 ;
 503     /* Calcul du nombre d'inequations du probleme. Le format de matrice de la
 504      * polylib permet les egalites, on doit donc les compter double quand il
 505      * y en a.
 506      */
 507     Nl = inequnk->NbRows ;
 508     for (i=0;i<inequnk->NbRows;i++)
 509     if (**(inequnk->p + i) == 0)
 510     Nl ++ ;
 511
 512     /* On prend les 'marques' de debut de traitement. */
 513     cross_product = 0 ;
 514     hq = tab_hwm() ;
 515     xq = p = sol_hwm();
 516
 517     /* On s'assure d'abord que le systeme pour le contexte n'est pas vide
 518      * avant de commencer le traitement. Si c'est le cas, la solution est
 519      * void (NULL).
 520      */
 521     if (ineqpar != NULL)
 522     { Nm = ineqpar->NbRows ;
 523       for (i=0;i<ineqpar->NbRows;i++)
 524       if (**(ineqpar->p + i) == 0)
 525       Nm ++ ;
 526       context = tab_Matrix2Tableau(ineqpar,Nm,Np,0) ;
 527       if (Nm)
 528       { /* Traduction du format de matrice de la polylib vers celui de
 529          * traitement de Pip. Puis traitement proprement dit.
 530          */
 531         ctxt = expanser(context, Np, Nm, Np+1, Np, 0, 0) ;
 532         traiter(ctxt, NULL, True, Np, 0, Nm, 0, -1) ;
 533         non_vide = is_not_Nil(p) ;
 534         sol_reset(p) ;
 535       }
 536       else
 537       non_vide = True ;
 538     }
 539     else
 540     { Nm = 0 ;
 541       context = NULL ; /* ATTENTION, en toute rigueur on devrait faire un
 542                         * tab_Matrix2Tableau d'une matrice 0*0.
 543                         */
 544       non_vide = True ;
 545     }
 546
 547     if (verbose)
 548     fprintf(dump, "%d %d %d %d %d %d\n",Nn, Np, Nl, Nm, Bg, Nq) ;
 549
 550     /* S'il est possible de trouver une solution, on passe au traitement. */
 551     if (non_vide)
 552     { ineq = tab_Matrix2Tableau(inequnk,Nl,Nn,Nn) ;
 553       compa_count = 0 ;
 554       traiter(ineq, context, Nq, Nn, Np, Nl, Nm, Bg) ;
 555
 556       if (Simplify)
 557       sol_simplify(xq) ;
 558       q = sol_hwm() ;
 559       /* On traduit la solution du format de solution de Pip vers un arbre
 560        * de structures de type PipQuast.
 561        */
 562       solution = sol_quast_edit(&xq,NULL) ;
 563       sol_reset(p) ;
 564     }
 565     else
 566     return NULL ;
 567     tab_reset(hq) ;
 568   }
 569   else
 570   return NULL ;
 571
 572   return(solution) ;
 573 }