1 /* automatically generated by sed scripts from the c source named below: */
2 /* cddlib.c: cdd library (library version of cdd)
3 written by Komei Fukuda, fukuda@ifor.math.ethz.ch
4 Version 0.94f, February 7, 2008
5 Standard ftp site: ftp.ifor.math.ethz.ch, Directory: pub/fukuda/cdd
6 */
7
8 /* cdd : C-Implementation of the double description method for
9 computing all vertices and extreme rays of the polyhedron
10 P= {x : b - A x >= 0}.
11 Please read COPYING (GNU General Public Licence) and
12 the manual cddlibman.tex for detail.
13 */
14
15 /* This program is free software; you can redistribute it and/or modify
16 it under the terms of the GNU General Public License as published by
17 the Free Software Foundation; either version 2 of the License, or
18 (at your option) any later version.
19
20 This program is distributed in the hope that it will be useful,
21 but WITHOUT ANY WARRANTY; without even the implied warranty of
22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 GNU General Public License for more details.
24
25 You should have received a copy of the GNU General Public License
26 along with this program; if not, write to the Free Software
27 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
28 */
29
30 /* The first version C0.21 was created on November 10,1993
31 with Dave Gillespie's p2c translator
32 from the Pascal program pdd.p written by Komei Fukuda.
33 */
34
35 #include "setoper.h"
36 /* set operation library header (Ver. June 1, 2000 or later) */
37 #include "cdd_f.h"
38 #include <stdio.h>
39 #include <stdlib.h>
40 #include <time.h>
41 #include <math.h>
42 #include <string.h>
43
44 /* Global Variables */
45 ddf_boolean ddf_debug =ddf_FALSE;
46 ddf_boolean ddf_log =ddf_FALSE;
47 /* GLOBAL CONSTANTS and STATICS VARIABLES (to be set by ddf_set_global_constants() */
48 myfloat ddf_zero;
49 myfloat ddf_one;
50 myfloat ddf_purezero;
51 myfloat ddf_minuszero;
52 myfloat ddf_minusone;
53
54 time_t ddf_statStartTime; /* cddlib starting time */
55 long ddf_statBApivots; /* basis finding pivots */
56 long ddf_statCCpivots; /* criss-cross pivots */
57 long ddf_statDS1pivots; /* phase 1 pivots */
58 long ddf_statDS2pivots; /* phase 2 pivots */
59 long ddf_statACpivots; /* anticycling (cc) pivots */
60 #ifdef ddf_GMPRATIONAL
61 long ddf_statBSpivots; /* basis status checking pivots */
62 #endif
63 ddf_LPSolverType ddf_choiceLPSolverDefault; /* Default LP solver Algorithm */
64 ddf_LPSolverType ddf_choiceRedcheckAlgorithm; /* Redundancy Checking Algorithm */
65 ddf_boolean ddf_choiceLexicoPivotQ; /* whether to use the lexicographic pivot */
66
67 /* #include <profile.h> THINK C PROFILER */
68 /* #include <console.h> THINK C PROFILER */
69
ddf_DDInit(ddf_ConePtr cone)70 void ddf_DDInit(ddf_ConePtr cone)
71 {
72 cone->Error=ddf_NoError;
73 cone->CompStatus=ddf_InProgress;
74 cone->RayCount = 0;
75 cone->TotalRayCount = 0;
76 cone->FeasibleRayCount = 0;
77 cone->WeaklyFeasibleRayCount = 0;
78 cone->EdgeCount=0; /* active edge count */
79 cone->TotalEdgeCount=0; /* active edge count */
80 ddf_SetInequalitySets(cone);
81 ddf_ComputeRowOrderVector(cone);
82 cone->RecomputeRowOrder=ddf_FALSE;
83 }
84
ddf_DDMain(ddf_ConePtr cone)85 void ddf_DDMain(ddf_ConePtr cone)
86 {
87 ddf_rowrange hh, itemp, otemp;
88 ddf_boolean locallog=ddf_log; /* if ddf_log=ddf_FALSE, no log will be written. */
89
90 if (cone->d<=0){
91 cone->Iteration=cone->m;
92 cone->FeasibleRayCount=0;
93 cone->CompStatus=ddf_AllFound;
94 goto _L99;
95 }
96 if (locallog) {
97 fprintf(stderr,"(Initially added rows ) = ");
98 set_fwrite(stderr,cone->InitialHalfspaces);
99 }
100 while (cone->Iteration <= cone->m) {
101 ddf_SelectNextHalfspace(cone, cone->WeaklyAddedHalfspaces, &hh);
102 if (set_member(hh,cone->NonequalitySet)){ /* Skip the row hh */
103 if (ddf_debug) {
104 fprintf(stderr,"*The row # %3ld should be inactive and thus skipped.\n", hh);
105 }
106 set_addelem(cone->WeaklyAddedHalfspaces, hh);
107 } else {
108 if (cone->PreOrderedRun)
109 ddf_AddNewHalfspace2(cone, hh);
110 else{
111 ddf_AddNewHalfspace1(cone, hh);
112 }
113 set_addelem(cone->AddedHalfspaces, hh);
114 set_addelem(cone->WeaklyAddedHalfspaces, hh);
115 }
116 if (!cone->PreOrderedRun){
117 for (itemp=1; cone->OrderVector[itemp]!=hh; itemp++);
118 otemp=cone->OrderVector[cone->Iteration];
119 cone->OrderVector[cone->Iteration]=hh;
120 /* store the dynamic ordering in ordervec */
121 cone->OrderVector[itemp]=otemp;
122 /* store the dynamic ordering in ordervec */
123 }
124 if (locallog){
125 fprintf(stderr,"(Iter, Row, #Total, #Curr, #Feas)= %5ld %5ld %9ld %6ld %6ld\n",
126 cone->Iteration, hh, cone->TotalRayCount, cone->RayCount,
127 cone->FeasibleRayCount);
128 }
129 if (cone->CompStatus==ddf_AllFound||cone->CompStatus==ddf_RegionEmpty) {
130 set_addelem(cone->AddedHalfspaces, hh);
131 goto _L99;
132 }
133 (cone->Iteration)++;
134 }
135 _L99:;
136 if (cone->d<=0 || cone->newcol[1]==0){ /* fixing the number of output */
137 cone->parent->n=cone->LinearityDim + cone->FeasibleRayCount -1;
138 cone->parent->ldim=cone->LinearityDim - 1;
139 } else {
140 cone->parent->n=cone->LinearityDim + cone->FeasibleRayCount;
141 cone->parent->ldim=cone->LinearityDim;
142 }
143 }
144
145
ddf_InitialDataSetup(ddf_ConePtr cone)146 void ddf_InitialDataSetup(ddf_ConePtr cone)
147 {
148 long j, r;
149 ddf_rowset ZSet;
150 static ddf_Arow Vector1,Vector2;
151 static ddf_colrange last_d=0;
152
153 if (last_d < cone->d){
154 if (last_d>0) {
155 for (j=0; j<last_d; j++){
156 ddf_init(Vector1[j]);
157 ddf_init(Vector2[j]);
158 }
159 free(Vector1); free(Vector2);
160 }
161 Vector1=(myfloat*)calloc(cone->d,sizeof(myfloat));
162 Vector2=(myfloat*)calloc(cone->d,sizeof(myfloat));
163 for (j=0; j<cone->d; j++){
164 ddf_init(Vector1[j]);
165 ddf_init(Vector2[j]);
166 }
167 last_d=cone->d;
168 }
169
170 cone->RecomputeRowOrder=ddf_FALSE;
171 cone->ArtificialRay = NULL;
172 cone->FirstRay = NULL;
173 cone->LastRay = NULL;
174 set_initialize(&ZSet,cone->m);
175 ddf_AddArtificialRay(cone);
176 set_copy(cone->AddedHalfspaces, cone->InitialHalfspaces);
177 set_copy(cone->WeaklyAddedHalfspaces, cone->InitialHalfspaces);
178 ddf_UpdateRowOrderVector(cone, cone->InitialHalfspaces);
179 for (r = 1; r <= cone->d; r++) {
180 for (j = 0; j < cone->d; j++){
181 ddf_set(Vector1[j], cone->B[j][r-1]);
182 ddf_neg(Vector2[j], cone->B[j][r-1]);
183 }
184 ddf_Normalize(cone->d, Vector1);
185 ddf_Normalize(cone->d, Vector2);
186 ddf_ZeroIndexSet(cone->m, cone->d, cone->A, Vector1, ZSet);
187 if (set_subset(cone->EqualitySet, ZSet)){
188 if (ddf_debug) {
189 fprintf(stderr,"add an initial ray with zero set:");
190 set_fwrite(stderr,ZSet);
191 }
192 ddf_AddRay(cone, Vector1);
193 if (cone->InitialRayIndex[r]==0) {
194 ddf_AddRay(cone, Vector2);
195 if (ddf_debug) {
196 fprintf(stderr,"and add its negative also.\n");
197 }
198 }
199 }
200 }
201 ddf_CreateInitialEdges(cone);
202 cone->Iteration = cone->d + 1;
203 if (cone->Iteration > cone->m) cone->CompStatus=ddf_AllFound; /* 0.94b */
204 set_free(ZSet);
205 }
206
ddf_CheckEmptiness(ddf_PolyhedraPtr poly,ddf_ErrorType * err)207 ddf_boolean ddf_CheckEmptiness(ddf_PolyhedraPtr poly, ddf_ErrorType *err)
208 {
209 ddf_rowset R, S;
210 ddf_MatrixPtr M=NULL;
211 ddf_boolean answer=ddf_FALSE;
212
213 *err=ddf_NoError;
214
215 if (poly->representation==ddf_Inequality){
216 M=ddf_CopyInequalities(poly);
217 set_initialize(&R, M->rowsize);
218 set_initialize(&S, M->rowsize);
219 if (!ddf_ExistsRestrictedFace(M, R, S, err)){
220 poly->child->CompStatus=ddf_AllFound;
221 poly->IsEmpty=ddf_TRUE;
222 poly->n=0;
223 answer=ddf_TRUE;
224 }
225 set_free(R);
226 set_free(S);
227 ddf_FreeMatrix(M);
228 } else if (poly->representation==ddf_Generator && poly->m<=0){
229 *err=ddf_EmptyVrepresentation;
230 poly->IsEmpty=ddf_TRUE;
231 poly->child->CompStatus=ddf_AllFound;
232 answer=ddf_TRUE;
233 poly->child->Error=*err;
234 }
235
236 return answer;
237 }
238
239
ddf_DoubleDescription(ddf_PolyhedraPtr poly,ddf_ErrorType * err)240 ddf_boolean ddf_DoubleDescription(ddf_PolyhedraPtr poly, ddf_ErrorType *err)
241 {
242 ddf_ConePtr cone=NULL;
243 ddf_boolean found=ddf_FALSE;
244
245 *err=ddf_NoError;
246
247 if (poly!=NULL && (poly->child==NULL || poly->child->CompStatus!=ddf_AllFound)){
248 cone=ddf_ConeDataLoad(poly);
249 /* create a cone associated with poly by homogenization */
250 time(&cone->starttime);
251 ddf_DDInit(cone);
252 if (poly->representation==ddf_Generator && poly->m<=0){
253 *err=ddf_EmptyVrepresentation;
254 cone->Error=*err;
255 goto _L99;
256 }
257 /* Check emptiness of the polyhedron */
258 ddf_CheckEmptiness(poly,err);
259
260 if (cone->CompStatus!=ddf_AllFound){
261 ddf_FindInitialRays(cone, &found);
262 if (found) {
263 ddf_InitialDataSetup(cone);
264 if (cone->CompStatus==ddf_AllFound) goto _L99;
265 ddf_DDMain(cone);
266 if (cone->FeasibleRayCount!=cone->RayCount) *err=ddf_NumericallyInconsistent; /* cddlib-093d */
267 }
268 }
269 time(&cone->endtime);
270 }
271
272 _L99: ;
273
274 return found;
275 }
276
ddf_DoubleDescription2(ddf_PolyhedraPtr poly,ddf_RowOrderType horder,ddf_ErrorType * err)277 ddf_boolean ddf_DoubleDescription2(ddf_PolyhedraPtr poly, ddf_RowOrderType horder, ddf_ErrorType *err)
278 {
279 ddf_ConePtr cone=NULL;
280 ddf_boolean found=ddf_FALSE;
281
282 *err=ddf_NoError;
283
284 if (poly!=NULL && (poly->child==NULL || poly->child->CompStatus!=ddf_AllFound)){
285 cone=ddf_ConeDataLoad(poly);
286 /* create a cone associated with poly by homogenization */
287 cone->HalfspaceOrder=horder; /* set the row order */
288 time(&cone->starttime);
289 ddf_DDInit(cone);
290 if (poly->representation==ddf_Generator && poly->m<=0){
291 *err=ddf_EmptyVrepresentation;
292 cone->Error=*err;
293 goto _L99;
294 }
295 /* Check emptiness of the polyhedron */
296 ddf_CheckEmptiness(poly,err);
297
298 if (cone->CompStatus!=ddf_AllFound){
299 ddf_FindInitialRays(cone, &found);
300 if (found) {
301 ddf_InitialDataSetup(cone);
302 if (cone->CompStatus==ddf_AllFound) goto _L99;
303 ddf_DDMain(cone);
304 if (cone->FeasibleRayCount!=cone->RayCount) *err=ddf_NumericallyInconsistent; /* cddlib-093d */
305 }
306 }
307 time(&cone->endtime);
308 }
309
310 _L99: ;
311
312 return found;
313 }
314
ddf_DDInputAppend(ddf_PolyhedraPtr * poly,ddf_MatrixPtr M,ddf_ErrorType * err)315 ddf_boolean ddf_DDInputAppend(ddf_PolyhedraPtr *poly, ddf_MatrixPtr M,
316 ddf_ErrorType *err)
317 {
318 /* This is imcomplete. It simply solves the problem from scratch. */
319 ddf_boolean found;
320 ddf_ErrorType error;
321
322 if ((*poly)->child!=NULL) ddf_FreeDDMemory(*poly);
323 ddf_AppendMatrix2Poly(poly, M);
324 (*poly)->representation=ddf_Inequality;
325 found=ddf_DoubleDescription(*poly, &error);
326 *err=error;
327 return found;
328 }
329
ddf_DDFile2File(char * ifile,char * ofile,ddf_ErrorType * err)330 ddf_boolean ddf_DDFile2File(char *ifile, char *ofile, ddf_ErrorType *err)
331 {
332 /* The representation conversion from an input file to an outfile. */
333 /* modified by D. Avis to allow stdin/stdout */
334 ddf_boolean found=ddf_TRUE;
335 FILE *reading=NULL,*writing=NULL;
336 ddf_PolyhedraPtr poly;
337 ddf_MatrixPtr M, A, G;
338
339 if (strcmp(ifile,"**stdin") == 0 )
340 reading = stdin;
341 else if ( ( reading = fopen(ifile, "r") )!= NULL) {
342 fprintf(stderr,"input file %s is open\n", ifile);
343 }
344 else{
345 fprintf(stderr,"The input file %s not found\n",ifile);
346 found=ddf_FALSE;
347 *err=ddf_IFileNotFound;
348 goto _L99;
349 }
350
351 if (found){
352 if (strcmp(ofile,"**stdout") == 0 )
353 writing = stdout;
354 else if ( (writing = fopen(ofile, "w") ) != NULL){
355 fprintf(stderr,"output file %s is open\n",ofile);
356 found=ddf_TRUE;
357 } else {
358 fprintf(stderr,"The output file %s cannot be opened\n",ofile);
359 found=ddf_FALSE;
360 *err=ddf_OFileNotOpen;
361 goto _L99;
362 }
363 }
364
365 M=ddf_PolyFile2Matrix(reading, err);
366 if (*err!=ddf_NoError){
367 goto _L99;
368 } poly=ddf_DDMatrix2Poly(M, err); /* compute the second representation */ ddf_FreeMatrix(M);
369
370 if (*err==ddf_NoError) {
371 ddf_WriteRunningMode(writing, poly);
372 A=ddf_CopyInequalities(poly);
373 G=ddf_CopyGenerators(poly);
374
375 if (poly->representation==ddf_Inequality) {
376 ddf_WriteMatrix(writing,G);
377 } else {
378 ddf_WriteMatrix(writing,A);
379 }
380
381 ddf_FreePolyhedra(poly);
382 ddf_FreeMatrix(A);
383 ddf_FreeMatrix(G);
384 }
385
386 _L99: ;
387 if (*err!=ddf_NoError) ddf_WriteErrorMessages(stderr,*err);
388 if (reading!=NULL) fclose(reading);
389 if (writing!=NULL) fclose(writing);
390 return found;
391 }
392
393 /* end of cddlib.c */
394