-
Notifications
You must be signed in to change notification settings - Fork 2
/
cdd.h
338 lines (307 loc) · 12.8 KB
/
cdd.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
/* cdd.h: Header file for cdd.C
written by Komei Fukuda, fukuda@ifor.math.ethz.ch
Version 0.77, August 19, 2003
*/
/* cdd.C : C++-Implementation of the double description method for
computing all vertices and extreme rays of the polyhedron
P= {x : b - A x >= 0}.
Please read COPYING (GNU General Public Licence) and
the manual cddman.tex for detail.
*/
#define COPYRIGHT "Copyright (C) 1999, Komei Fukuda, fukuda@ifor.math.ethz.ch"
#define DDVERSION "Version 0.77(August 19, 2003)"
#ifdef RATIONAL
#ifdef GMP
#define ARITHMETIC "Compiled for Rational Exact Arithmetic with GMP"
#else
#define ARITHMETIC "Compiled for Rational Exact Arithmetic with G++"
#endif
#else
#define ARITHMETIC "Compiled for Floating-Point Arithmetic"
#endif
#include <ctime>
#include <iostream>
#include <fstream>
using std::string;
using std::ifstream;
using std::ofstream;
#ifdef RATIONAL
typedef Rational myTYPE;
#else
typedef double myTYPE;
#endif // RATIONAL
typedef Rational myRational;
typedef int boolean;
typedef long rowrange;
typedef long colrange;
typedef set_type rowset; /* set_type defined in setoper.h */
typedef set_type colset;
typedef long *rowindex;
/* rowindex should be intialized to be an array of [mm+1] components */
typedef long colindex[NMAX+1];
typedef myTYPE *Amatrix[MMAX];
typedef myTYPE *Arow;
typedef myTYPE *Bmatrix[NMAX];
typedef set_type Aincidence[MMAX];
typedef int *SignAmatrix[MMAX]; /* Sign (+1, 0 ,-1)-matrix of Amatrix type */
typedef char DataFileType[filenamelen];
typedef char LineType[linelenmax];
typedef char WordType[wordlenmax];
typedef struct RayRecord {
myTYPE *Ray;
rowset ZeroSet;
rowrange FirstInfeasIndex; /* the first inequality the ray violates */
boolean feasible; /* flag to store the feasibility */
myTYPE *ARay; /* temporary area to store some row of A*Ray */
struct RayRecord *Next;
} RayRecord;
typedef struct AdjacencyRecord {
RayRecord *Ray1, *Ray2;
struct AdjacencyRecord *Next;
} AdjacencyRecord;
typedef struct node {long key; struct node *next;} node;
typedef enum {
Combinatorial, Algebraic
} AdjacencyTestType;
typedef enum {
MaxIndex, MinIndex, MinCutoff, MaxCutoff, MixCutoff,
LexMin, LexMax, RandomRow, LineShelling
} HyperplaneOrderType;
typedef enum {
Real, Rational, Integer, Unknown
} NumberType;
typedef enum {
ZeroRHS, NonzeroRHS
} InequalityType;
typedef enum {
IneToExt, ExtToIne, Projection,
LPmax, LPmin, FacetListing, VertexListing, TopeListing, InteriorFind,
FacetListingExternal,VertexListingExternal
} ConversionType;
typedef enum {
RowDecomposition, RowSubproblemSolve, Nothing
} SpecialConversionType;
typedef enum {
CrissCross,DualSimplex,CombMaxImprove
} LPsolverType;
typedef enum {
IncOff=0, IncCardinality, InputIncidence, OutputIncidence, IOIncidence
} IncidenceOutputType;
typedef enum {
AdjOff=0, OutputAdjacency, InputAdjacency, IOAdjacency
} AdjacencyOutputType;
typedef enum {
Auto, SemiAuto, Manual
} FileInputModeType;
/* Auto if a input filename is specified by command arguments */
typedef enum {
DimensionTooLarge, LowColumnRank, ImproperInputFormat, DependentMarkedSet,
ImproperExecutable, FileNotFound, None
} ErrorType;
typedef enum {
InProgress, AllFound, RegionEmpty
} CompStatusType;
typedef enum {
LPSundecided, Optimal, Inconsistent, DualInconsistent, Unbounded, DualUnbounded
} LPStatusType;
extern long minput, ninput; /*size of input data [b -A] */
extern long mm, nn; /*size of the homogenous system to be solved by dd*/
extern long projdim; /*dimension of orthogonal preprojection */
extern colset projvars; /*set of variables spanning the space of preprojection,
i.e. the remaining variables are to be removed*/
extern rowset EqualitySet, NonequalitySet, GroundSet, Face, Face1,SubproblemRowSet, RedundantRowSet;
extern rowrange Iteration, hh;
extern rowindex OrderVector;
extern rowindex EqualityIndex;
extern rowset AddedHyperplanes, WeaklyAddedHyperplanes, InitialHyperplanes;
extern long RayCount, FeasibleRayCount, WeaklyFeasibleRayCount,
TotalRayCount, VertexCount,ZeroRayCount;
extern long EdgeCount,TotalEdgeCount;
extern long count_int,count_int_good,count_int_bad;
extern boolean DynamicWriteOn, DynamicRayWriteOn, LogWriteOn,
ShowSignTableauOn, PostAnalysisOn,
CondensedListOn, SignPivotOn, ManualPivotOn, debug;
extern Amatrix AA;
extern Bmatrix InitialRays;
extern colindex InitialRayIndex;
extern Arow LPcost;
extern LPStatusType LPStatus;
extern colrange RHScol; /* LP RHS column */
extern rowrange OBJrow; /* LP OBJ row */
extern RayRecord *ArtificialRay, *FirstRay, *LastRay;
extern RayRecord *PosHead, *ZeroHead, *NegHead, *PosLast, *ZeroLast, *NegLast;
extern AdjacencyRecord *Edges[MMAX]; /* adjacency relation storage for iteration k */
extern boolean RecomputeRowOrder, inputsuccessful;
extern HyperplaneOrderType HyperplaneOrder;
extern AdjacencyTestType AdjacencyTest;
extern NumberType Number;
extern string InputNumberString, OutputNumberString;
extern InequalityType Inequality;
extern boolean NondegAssumed; /* Nondegeneacy preknowledge flag */
extern boolean InitBasisAtBottom; /* if it is on, the initial Basis will be selected at bottom */
extern boolean RestrictedEnumeration; /* Restricted enumeration Switch (True if it is restricted on the intersection of EqualitySet hyperplanes) */
extern boolean RelaxedEnumeration; /* Relaxed enumeration Switch (True if NonequalitySet inequalities must be satisfied with strict inequality) */
extern boolean VerifyInput; /* Verification switch for the input data */
extern boolean PreOrderedRun;
extern CompStatusType CompStatus; /* Computation Status */
extern ConversionType Conversion;
extern SpecialConversionType SpecialConversion; /* rowdecomposition, rowsubproblem */
extern LPsolverType LPsolver;
extern IncidenceOutputType IncidenceOutput;
extern AdjacencyOutputType AdjacencyOutput;
extern ErrorType Error;
extern FileInputModeType FileInputMode;
extern DataFileType inputfile,ifilehead,ifiletail,
outputfile,projfile,incfile,adjfile,logfile,dexfile,verfile,xtnfile;
extern time_t starttime, endtime;
extern unsigned int rseed;
extern myTYPE zero; /*Rational or floating zero*/
extern boolean Round_Output; /* rounding option for floating-point output. */
extern int output_digits; /* Float digits for output. Does not affect the computation. */
void SetInputFile(boolean *);
void SetWriteFileName(DataFileType, char, char *);
void SetReadFileName(DataFileType, char, char *);
myTYPE FABS(myTYPE);
void SetNumberType(string);
void ProcessCommandLine(ifstream &, string);
void AmatrixInput(boolean *);
void SetInequalitySets(rowindex);
void RandomPermutation(rowindex, long, unsigned int);
void QuickSort(rowindex, long, long, Amatrix, long);
myTYPE AValue(myTYPE *, rowrange );
void WriteIncidence(ostream &, RayRecord *);
void StoreRay1(myTYPE *, RayRecord *, boolean *);
void StoreRay2(myTYPE *, RayRecord *, boolean *, boolean *);
void AddRay(myTYPE *);
void AddArtificialRay(void);
void ConditionalAddEdge(RayRecord *Ray1, RayRecord *Ray2, RayRecord *ValidFirstRay);
void CreateInitialEdges(void);
void UpdateEdges(RayRecord *RRbegin, RayRecord *RRend);
void FreeDDMemory(void);
void Normalize(myTYPE *);
void ZeroIndexSet(myTYPE *, rowset);
void CopyBmatrix(Bmatrix T, Bmatrix TCOPY);
void SelectPivot1(Amatrix, HyperplaneOrderType,
rowrange, rowset, colset, rowrange *, colrange *,boolean *);
myTYPE TableauEntry(Amatrix, Bmatrix T, rowrange, colrange);
void WriteTableau(ostream &,Amatrix, Bmatrix T, InequalityType);
char Sign(myTYPE);
int myTYPE2sign(myTYPE);
int long2sign(long);
void WriteSignAmatrix(ostream &f, SignAmatrix X,
rowindex OV, long bflag[], rowrange objrow, colrange rhscol);
void WriteSignTableau(ostream &f, Amatrix X, Bmatrix T,
rowindex OV, long bflag[], rowrange objrow, colrange rhscol);
void WriteDictionary(ostream &f, Amatrix X, Bmatrix T,
rowindex OV, long bflag[], colindex NBIndex, rowrange objrow, colrange rhscol);
void OutputTableau(Amatrix, Bmatrix T,InequalityType);
void SelectPivot2(Amatrix, Bmatrix T,
HyperplaneOrderType,rowrange, rowset, colset,
rowrange *, colrange *,boolean *);
void GausianColumnPivot1(Amatrix, Bmatrix, SignAmatrix, rowrange, colrange);
void GausianColumnPivot2(Amatrix, Bmatrix, rowrange, colrange);
void InitializeBmatrix(Bmatrix T);
void SetToIdentity(Bmatrix T);
void free_Bmatrix(Bmatrix T);
void WriteBmatrix(ostream &, Bmatrix T);
void ReduceAA(rowset, colset);
void DualizeAA(Bmatrix T);
void EnlargeAAforInteriorFinding(void);
void EnlargeAAforZeroRHSLP(void);
void RecoverAAafterInteriorFinding(void);
void ShiftPointsAroundOrigin(ostream &, ostream &, Arow);
boolean RowEquivalent_Q(Arow a1, Arow a2, colrange n);
void FindRowEquivalenceClasses(long *classno, rowindex);
void WriteSubMatrixOfAA(ostream &, rowset, colset, InequalityType);
void WriteAmatrix(ostream &, Amatrix, long, long, InequalityType);
void WriteErrorMessages(ostream &);
void ComputeRank(Amatrix, unsigned long *, long *);
void ComputeBInverse(Amatrix, long, Bmatrix InvA1, long *);
void FindBasis(Amatrix, HyperplaneOrderType, rowset, long *,
Bmatrix BasisInverse, long *);
void SelectCrissCrossPivot(Amatrix, Bmatrix, rowindex,
long, rowrange,colrange,rowrange *,colrange *,
boolean *, LPStatusType *);
void SelectDualSimplexPivot(boolean, Amatrix, Bmatrix, rowindex,
colindex, long, rowrange, colrange,
rowrange *, colrange *, boolean *, LPStatusType *);
void CrissCrossMinimize(ostream &, ostream &, Amatrix,Bmatrix,
rowrange, colrange, boolean,
LPStatusType *, myTYPE *, Arow, Arow, colindex,
rowrange *, colrange *, long *);
void CrissCrossMaximize(ostream &, ostream &, Amatrix,Bmatrix,
rowrange, colrange, boolean,
LPStatusType *, myTYPE *, Arow, Arow, colindex,
rowrange *, colrange *, long *);
void DualSimplexMinimize(ostream &, ostream &, Amatrix,Bmatrix,
rowrange, colrange, boolean,
LPStatusType *, myTYPE *, Arow, Arow, colindex,
rowrange *, colrange *, long *);
void DualSimplexMaximize(ostream &, ostream &, Amatrix,Bmatrix,
rowrange, colrange, boolean,
LPStatusType *, myTYPE *, Arow, Arow, colindex,
rowrange *, colrange *, long *);
void ManualPivot(ostream &, ostream &, Amatrix,Bmatrix BasisInverse,
rowrange, colrange,
LPStatusType *, myTYPE *optvalue, Arow, Arow, colindex,
rowrange *, colrange *, long *);
void WriteLPResult(ostream &, LPStatusType, myTYPE,
Arow, Arow, colindex, rowrange, colrange, long);
void FindInitialRays(rowset, Bmatrix, colindex, boolean *);
void CheckAdjacency1(RayRecord **, RayRecord **,boolean *);
void CheckAdjacency2(RayRecord **, RayRecord **,boolean *);
void CheckEquality(RayRecord **, RayRecord **, boolean *);
void Eliminate(RayRecord **);
void CreateNewRay(RayRecord *, RayRecord *, rowrange);
void EvaluateARay1(rowrange);
void EvaluateARay2(rowrange);
void FeasibilityIndices(long *, long *, rowrange);
boolean LexSmaller(myTYPE *, myTYPE *, long);
boolean LexLarger(myTYPE *, myTYPE *, long);
void LineShellingOrder(rowindex, myTYPE *, myTYPE *);
void CompileDecompResult(ofstream &);
void ReadExtFile(ifstream &);
void CopyArow(myTYPE *, myTYPE *, long);
void ComputeRowOrderVector(rowindex OV, HyperplaneOrderType ho);
void UpdateRowOrderVector(rowset PriorityRows);
void SelectNextHyperplane(HyperplaneOrderType, unsigned long *, rowrange *, boolean *);
void SelectPreorderedNext(long *excluded, rowindex, rowrange *hnext);
void AddNewHyperplane1(rowrange);
void AddNewHyperplane2(rowrange);
void WriteRunningMode(ostream &);
void WriteRunningMode0(ostream &);
void WriteRunningMode2(ostream &);
void WriteCompletionStatus(ostream &);
void WriteTimes(ostream &);
void WriteProgramDescription(ostream &);
void WriteSolvedProblem(ostream &);
void WriteNumber(ostream &, myTYPE);
void WriteRayRecord(ostream &, RayRecord *);
void WriteRayRecord2(ostream &, RayRecord *);
void WriteExtFile(ostream &, ostream &);
void WriteIncidenceFile(ostream &);
void WriteInputIncidenceFile(ostream &);
void WriteInputAdjacencyFile(ostream &);
void WriteAdjacencyFile(ostream &);
void WriteProjResult(ostream &, ostream &, long *);
void WriteProjRayRecord(ostream &, RayRecord *, long *);
void InitialWriting(ostream &, ostream &);
void WriteDecompResult(ostream &f, ostream &f_log);
void WriteRowEquivalence(ostream &f, long classno, rowindex rowequiv);
void OutputHeading(void);
void InitialDataSetup(void);
void LPInit(void);
void DDInit(void);
void LPMain(ostream &, ostream &);
void PostAnalysisMain(ifstream &, ostream &);
void InteriorFindMain(ostream &, ostream &, boolean *);
void WriteCurrentSolution(ostream &, Amatrix, Bmatrix, rowrange, colrange, colindex);
void CheckConversionConsistency(ostream &, ostream &);
// procedures in cddrevs.C
boolean Facet_Q(topeOBJECT, rowrange);
boolean Facet_Q2(topeOBJECT, rowrange,colindex,Arow);
void ReverseSearch(ostream &, topeOBJECT, long delta);
void FacetandVertexListMain(ostream &, ostream &);
void FacetandVertexExternalListMain(ostream &, ostream &);
void TopeListMain(ostream &, ostream &);
/* end of cdd.h */