piputil.c

   1 #include <polylib/polylibgmp.h>
   2 #include <stdio.h>  /* needed for piplib ! */
   3 #include <piplib/piplibMP.h>
   4 #include <assert.h>
   5
   6 static PipMatrix *poly2pip(Polyhedron *P, int exist, int nparam)
   7 {
   8     PipMatrix *     matrix;
   9     int i, j;
  10     unsigned int nvar = P->Dimension - exist - nparam;
  11
  12     matrix = pip_matrix_alloc(P->NbConstraints, P->Dimension+2);
  13     for (i = 0; i < P->NbConstraints; ++i) {
  14         value_assign(matrix->p[i][0], P->Constraint[i][0]);
  15         for (j = 0; j < exist; ++j)
  16             value_assign(matrix->p[i][1+j], P->Constraint[i][1+nvar+j]);
  17         for (j = 0; j < nvar; ++j)
  18             value_assign(matrix->p[i][1+exist+j], P->Constraint[i][1+j]);
  19         for (j = 1+exist+nvar; j < P->Dimension+2; ++j)
  20             value_assign(matrix->p[i][j], P->Constraint[i][j]);
  21     }
  22
  23     return matrix;
  24 }
  25
  26 static int max_new(PipQuast *q, int max, int d, int *maxd)
  27 {
  28     PipNewparm *p;
  29
  30     for (p = q->newparm; p; p = p->next)
  31         if (p->rank > max)
  32             max = p->rank;
  33     if (q->condition) {
  34         if (++d > *maxd)
  35             *maxd = d;
  36         max = max_new(q->next_else, max, d, maxd);
  37         max = max_new(q->next_then, max, d, maxd);
  38     }
  39     return max;
  40 }
  41
  42 static int vectorpos2cpos(int v, int n, int e, int d, int j)
  43 {
  44     int np = d-v-n-e;
  45     int l;
  46
  47     if (j < v)
  48         l = 1+j;
  49     else if (j < np+v)
  50         l = 1+v+n+e+j-v;
  51     else
  52         l = 1+v+n+j-v-np;
  53
  54     return l;
  55 }
  56
  57 /*
  58  * i: next row
  59  * n: dimension of space in which minimum was taken
  60  *    i.e., number of existential vars in original problem
  61  * v: number of variables in original problem
  62  * e: number of extra existential variables
  63  * d: total dimension
  64  *
  65  * Dimensions in the resulting domain are:
  66  *      v n e (d-v-n-e)
  67  */
  68 static void add_quast(Polyhedron**D, Matrix* M, PipQuast *q,
  69                       int i, int v, int n, int e, int d)
  70 {
  71     PipNewparm *p;
  72
  73     for (p = q->newparm; p; p = p->next) {
  74         int j, l;
  75         PipVector *vec = p->vector;
  76
  77         value_set_si(M->p[i][0], 1);
  78         Vector_Set(M->p[i]+1, 0, d);
  79
  80         for (j = 0; j < vec->nb_elements-1; ++j) {
  81             l = vectorpos2cpos(v, n, e, d, j);
  82             value_assign(M->p[i][l], vec->the_vector[j]);
  83         }
  84         l = vectorpos2cpos(v, n, e, d, p->rank);
  85         value_oppose(M->p[i][l], p->deno);
  86         value_assign(M->p[i][1+d], vec->the_vector[j]);
  87         ++i;
  88
  89         value_set_si(M->p[i][0], 1);
  90         Vector_Set(M->p[i]+1, 0, d);
  91
  92         for (j = 0; j < vec->nb_elements-1; ++j) {
  93             l = vectorpos2cpos(v, n, e, d, j);
  94             value_oppose(M->p[i][l], vec->the_vector[j]);
  95         }
  96         l = vectorpos2cpos(v, n, e, d, p->rank);
  97         value_assign(M->p[i][l], p->deno);
  98         value_oppose(M->p[i][1+d], vec->the_vector[j]);
  99         value_addto(M->p[i][1+d], M->p[i][1+d], p->deno);
 100         value_decrement(M->p[i][1+d], M->p[i][1+d]);
 101         ++i;
 102     }
 103
 104     if (q->condition) {
 105         int j;
 106         value_set_si(M->p[i][0], 1);
 107         Vector_Set(M->p[i]+1, 0, d);
 108
 109         for (j = 0; j < q->condition->nb_elements-1; ++j) {
 110             int l = vectorpos2cpos(v, n, e, d, j);
 111             value_assign(M->p[i][l], q->condition->the_vector[j]);
 112         }
 113         value_assign(M->p[i][1+d], q->condition->the_vector[j]);
 114         add_quast(D, M, q->next_then, i+1, v, n, e, d);
 115
 116         for (j = 0; j < q->condition->nb_elements-1; ++j) {
 117             int l = vectorpos2cpos(v, n, e, d, j);
 118             value_oppose(M->p[i][l], q->condition->the_vector[j]);
 119         }
 120         value_oppose(M->p[i][1+d], q->condition->the_vector[j]);
 121         value_decrement(M->p[i][1+d], M->p[i][1+d]);
 122         add_quast(D, M, q->next_else, i+1, v, n, e, d);
 123     } else if (q->list) {
 124         PipList *l;
 125         Matrix *C;
 126         Polyhedron *P;
 127         int j, k;
 128         for (j = 0, l = q->list; l; ++j, l = l->next) {
 129             Vector_Set(M->p[i], 0, d+1);
 130             value_set_si(M->p[i][1+v+j], -1);
 131
 132             for (k = 0; k < l->vector->nb_elements-1; ++k) {
 133                 int ll = vectorpos2cpos(v, n, e, d, k);
 134                 value_assign(M->p[i][ll], l->vector->the_vector[k]);
 135             }
 136             value_assign(M->p[i][1+d], l->vector->the_vector[k]);
 137
 138             ++i;
 139         }
 140         C = Matrix_Alloc(i, 1+d+1);
 141         Vector_Copy(M->p[0], C->p[0], i * (1+d+1));
 142         P = Constraints2Polyhedron(C, 0);
 143         Matrix_Free(C);
 144         P->next = *D;
 145         *D = P;
 146     }
 147 }
 148
 149 static Polyhedron *quast2poly(PipQuast *q, int nvar, int nset)
 150 {
 151     int                 nparam, nexist;
 152     PipList*            l;
 153     int                 i, d, rows;
 154     Matrix*             M;
 155     Polyhedron*         D;
 156
 157     assert(nset != 0);  /* required ? */
 158     nparam = q->newparm ? q->newparm->rank :
 159              q->condition ? q->condition->nb_elements-1 :
 160                             q->list->vector->nb_elements-1;
 161     d = 0;
 162     nexist = max_new(q, nparam-1, 0, &d) - nparam+1;
 163     rows = nparam + 2 * nexist + d + nset;
 164     M = Matrix_Alloc(rows, 1+nset+nexist+nparam+1);
 165     /* All parameters are/should be positive */
 166     for (i = 0; i < nvar; ++i) {
 167         value_set_si(M->p[i][0], 1);
 168         value_set_si(M->p[i][1+i], 1);
 169     }
 170     for ( ; i < nparam; ++i) {
 171         value_set_si(M->p[i][0], 1);
 172         value_set_si(M->p[i][1+nvar+nset+nexist+i-nvar], 1);
 173     }
 174     D = 0;
 175     add_quast(&D, M, q, nparam, nvar, nset, nexist, nset+nexist+nparam);
 176     Matrix_Free(M);
 177     return D;
 178 }
 179
 180 Polyhedron *pip_lexmin(Polyhedron *P, int exist, int nparam)
 181 {
 182     PipOptions  *options;
 183     PipMatrix * domain;
 184     unsigned int nvar = P->Dimension - exist - nparam;
 185     PipMatrix   *context = pip_matrix_alloc(0, nvar+nparam+2);
 186     PipQuast    *sol;
 187     Polyhedron  *min;
 188
 189     domain = poly2pip(P, exist, nparam);
 190     /*
 191     pip_matrix_print(stderr, domain);
 192     pip_matrix_print(stderr, context);
 193     */
 194
 195     options = pip_options_init();
 196     sol = pip_solve(domain, context, -1, options);
 197
 198     /*
 199     pip_quast_print(stderr, sol, 0);
 200     */
 201
 202     min = quast2poly(sol, nvar, exist);
 203
 204     pip_quast_free(sol);
 205     pip_matrix_free(context);
 206     pip_matrix_free(domain);
 207     pip_options_free(options);
 208
 209     return min;
 210 }