1 /* automatically generated by sed scripts from the c source named below: */
2 /* cddproj.c:  Polyhedral Projections in cddlib
3    written by Komei Fukuda, fukuda@cs.mcgill.ca
4    Version 0.94, Aug. 4, 2005
5 */
6 
7 /* cddlib : C-library of the double description method for
8    computing all vertices and extreme rays of the polyhedron
9    P= {x :  b - A x >= 0}.
10    Please read COPYING (GNU General Public Licence) and
11    the manual cddlibman.tex for detail.
12 */
13 
14 #include "setoper.h"  /* set operation library header (Ver. June 1, 2000 or later) */
15 #include "cdd_f.h"
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <time.h>
19 #include <math.h>
20 #include <string.h>
21 
ddf_BlockElimination(ddf_MatrixPtr M,ddf_colset delset,ddf_ErrorType * error)22 ddf_MatrixPtr ddf_BlockElimination(ddf_MatrixPtr M, ddf_colset delset, ddf_ErrorType *error)
23 /* Eliminate the variables (columns) delset by
24    the Block Elimination with ddf_DoubleDescription algorithm.
25 
26    Given (where y is to be eliminated):
27    c1 + A1 x + B1 y >= 0
28    c2 + A2 x + B2 y =  0
29 
30    1. First construct the dual system:  z1^T B1 + z2^T B2 = 0, z1 >= 0.
31    2. Compute the generators of the dual.
32    3. Then take the linear combination of the original system with each generator.
33    4. Remove redundant inequalies.
34 
35 */
36 {
37   ddf_MatrixPtr Mdual=NULL, Mproj=NULL, Gdual=NULL;
38   ddf_rowrange i,h,m,mproj,mdual,linsize;
39   ddf_colrange j,k,d,dproj,ddual,delsize;
40   ddf_colindex delindex;
41   myfloat temp,prod;
42   ddf_PolyhedraPtr dualpoly;
43   ddf_ErrorType err=ddf_NoError;
44   ddf_boolean localdebug=ddf_FALSE;
45 
46   *error=ddf_NoError;
47   m= M->rowsize;
48   d= M->colsize;
49   delindex=(long*)calloc(d+1,sizeof(long));
50   ddf_init(temp);
51   ddf_init(prod);
52 
53   k=0; delsize=0;
54   for (j=1; j<=d; j++){
55     if (set_member(j, delset)){
56       k++;  delsize++;
57       delindex[k]=j;  /* stores the kth deletion column index */
58     }
59   }
60   if (localdebug) ddf_WriteMatrix(stdout, M);
61 
62   linsize=set_card(M->linset);
63   ddual=m+1;
64   mdual=delsize + m - linsize;  /* #equalitions + dimension of z1 */
65 
66   /* setup the dual matrix */
67   Mdual=ddf_CreateMatrix(mdual, ddual);
68   Mdual->representation=ddf_Inequality;
69   for (i = 1; i <= delsize; i++){
70     set_addelem(Mdual->linset,i);  /* equality */
71     for (j = 1; j <= m; j++) {
72       ddf_set(Mdual->matrix[i-1][j], M->matrix[j-1][delindex[i]-1]);
73     }
74   }
75 
76   k=0;
77   for (i = 1; i <= m; i++){
78     if (!set_member(i, M->linset)){
79       /* set nonnegativity for the dual variable associated with
80          each non-linearity inequality. */
81       k++;
82       ddf_set(Mdual->matrix[delsize+k-1][i], ddf_one);
83     }
84   }
85 
86   /* 2. Compute the generators of the dual system. */
87   dualpoly=ddf_DDMatrix2Poly(Mdual, &err);
88   Gdual=ddf_CopyGenerators(dualpoly);
89 
90   /* 3. Take the linear combination of the original system with each generator.  */
91   dproj=d-delsize;
92   mproj=Gdual->rowsize;
93   Mproj=ddf_CreateMatrix(mproj, dproj);
94   Mproj->representation=ddf_Inequality;
95   set_copy(Mproj->linset, Gdual->linset);
96 
97   for (i=1; i<=mproj; i++){
98     k=0;
99     for (j=1; j<=d; j++){
100       if (!set_member(j, delset)){
101         k++;  /* new index of the variable x_j  */
102         ddf_set(prod, ddf_purezero);
103         for (h = 1; h <= m; h++){
104           ddf_mul(temp,M->matrix[h-1][j-1],Gdual->matrix[i-1][h]);
105           ddf_add(prod,prod,temp);
106         }
107         ddf_set(Mproj->matrix[i-1][k-1],prod);
108       }
109     }
110   }
111   if (localdebug) printf("Size of the projection system: %ld x %ld\n", mproj, dproj);
112 
113   ddf_FreePolyhedra(dualpoly);
114   free(delindex);
115   ddf_clear(temp);
116   ddf_clear(prod);
117   ddf_FreeMatrix(Mdual);
118   ddf_FreeMatrix(Gdual);
119   return Mproj;
120 }
121 
122 
ddf_FourierElimination(ddf_MatrixPtr M,ddf_ErrorType * error)123 ddf_MatrixPtr ddf_FourierElimination(ddf_MatrixPtr M,ddf_ErrorType *error)
124 /* Eliminate the last variable (column) from the given H-matrix using
125    the standard Fourier Elimination.
126  */
127 {
128   ddf_MatrixPtr Mnew=NULL;
129   ddf_rowrange i,inew,ip,in,iz,m,mpos=0,mneg=0,mzero=0,mnew;
130   ddf_colrange j,d,dnew;
131   ddf_rowindex posrowindex, negrowindex,zerorowindex;
132   myfloat temp1,temp2;
133   ddf_boolean localdebug=ddf_FALSE;
134 
135   *error=ddf_NoError;
136   m= M->rowsize;
137   d= M->colsize;
138   if (d<=1){
139     *error=ddf_ColIndexOutOfRange;
140     if (localdebug) {
141       printf("The number of column is too small: %ld for Fourier's Elimination.\n",d);
142     }
143     goto _L99;
144   }
145 
146   if (M->representation==ddf_Generator){
147     *error=ddf_NotAvailForV;
148     if (localdebug) {
149       printf("Fourier's Elimination cannot be applied to a V-polyhedron.\n");
150     }
151     goto _L99;
152   }
153 
154   if (set_card(M->linset)>0){
155     *error=ddf_CannotHandleLinearity;
156     if (localdebug) {
157       printf("The Fourier Elimination function does not handle equality in this version.\n");
158     }
159     goto _L99;
160   }
161 
162   /* Create temporary spaces to be removed at the end of this function */
163   posrowindex=(long*)calloc(m+1,sizeof(long));
164   negrowindex=(long*)calloc(m+1,sizeof(long));
165   zerorowindex=(long*)calloc(m+1,sizeof(long));
166   ddf_init(temp1);
167   ddf_init(temp2);
168 
169   for (i = 1; i <= m; i++) {
170     if (ddf_Positive(M->matrix[i-1][d-1])){
171       mpos++;
172       posrowindex[mpos]=i;
173     } else if (ddf_Negative(M->matrix[i-1][d-1])) {
174       mneg++;
175       negrowindex[mneg]=i;
176     } else {
177       mzero++;
178       zerorowindex[mzero]=i;
179     }
180   }  /*of i*/
181 
182   if (localdebug) {
183     ddf_WriteMatrix(stdout, M);
184     printf("No of  (+  -  0) rows = (%ld, %ld, %ld)\n", mpos,mneg, mzero);
185   }
186 
187   /* The present code generates so many redundant inequalities and thus
188      is quite useless, except for very small examples
189   */
190   mnew=mzero+mpos*mneg;  /* the total number of rows after elimination */
191   dnew=d-1;
192 
193   Mnew=ddf_CreateMatrix(mnew, dnew);
194   ddf_CopyArow(Mnew->rowvec, M->rowvec, dnew);
195 /*  set_copy(Mnew->linset,M->linset);  */
196   Mnew->numbtype=M->numbtype;
197   Mnew->representation=M->representation;
198   Mnew->objective=M->objective;
199 
200 
201   /* Copy the inequalities independent of x_d to the top of the new matrix. */
202   for (iz = 1; iz <= mzero; iz++){
203     for (j = 1; j <= dnew; j++) {
204       ddf_set(Mnew->matrix[iz-1][j-1], M->matrix[zerorowindex[iz]-1][j-1]);
205     }
206   }
207 
208   /* Create the new inequalities by combining x_d positive and negative ones. */
209   inew=mzero;  /* the index of the last x_d zero inequality */
210   for (ip = 1; ip <= mpos; ip++){
211     for (in = 1; in <= mneg; in++){
212       inew++;
213       ddf_neg(temp1, M->matrix[negrowindex[in]-1][d-1]);
214       for (j = 1; j <= dnew; j++) {
215         ddf_LinearComb(temp2,M->matrix[posrowindex[ip]-1][j-1],temp1,\
216           M->matrix[negrowindex[in]-1][j-1],\
217           M->matrix[posrowindex[ip]-1][d-1]);
218         ddf_set(Mnew->matrix[inew-1][j-1],temp2);
219       }
220       ddf_Normalize(dnew,Mnew->matrix[inew-1]);
221     }
222   }
223 
224 
225   free(posrowindex);
226   free(negrowindex);
227   free(zerorowindex);
228   ddf_clear(temp1);
229   ddf_clear(temp2);
230 
231  _L99:
232   return Mnew;
233 }
234 
235 
236 /* end of cddproj.c */
237