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