/SqlSpatialIndexing/Lib/Cdd/IO.cs
C# | 2168 lines | 1751 code | 272 blank | 145 comment | 377 complexity | 7cf76b086103d7ab21ee0643618c47d2 MD5 | raw file
- using System;
- using System.Collections.Generic;
- using System.Text;
- using System.IO;
-
- namespace CddSharp
- {
- public partial class Cdd
- {
- /* cddio.c: Basic Input and Output Procedures for cddlib
- written by Komei Fukuda, fukuda@ifor.math.ethz.ch
- Version 0.94, Aug. 4, 2005
- */
-
- /* cddlib : C-library 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 cddlibman.tex for detail.
- */
-
- #if nemkell
- void SetInputFile(FILE **f,DataFileType inputfile,ErrorType *Error)
- {
- int opened=0,stop,quit=0;
- int i,dotpos=0,trial=0;
- char ch;
- char *tempname;
-
-
- *Error=NoError;
- while (!opened && !quit) {
- fprintf(stderr,"\n>> Input file: ");
- scanf("%s",inputfile);
- ch=getchar();
- stop=false;
- for (i=0; i<filenamelen && !stop; i++){
- ch=inputfile[i];
- switch (ch) {
- case '.':
- dotpos=i+1;
- break;
- case ';': case ' ': case '\0': case '\n': case '\t':
- stop=true;
- tempname=(char*)calloc(filenamelen,sizeof(ch));
- strncpy(tempname,inputfile,i);
- strcpy(inputfile,tempname);
- free(tempname);
- break;
- }
- }
- if ( ( *f = fopen(inputfile,"r") )!= null) {
- fprintf(stderr,"input file %s is open\n",inputfile);
- opened=1;
- *Error=NoError;
- }
- else{
- fprintf(stderr,"The file %s not found\n",inputfile);
- trial++;
- if (trial>5) {
- *Error=IFileNotFound;
- quit=1;
- }
- }
- }
- }
-
- void SetWriteFileName(DataFileType inputfile, DataFileType outfile, char cflag, RepresentationType rep)
- {
- char *extension;
- DataFileType ifilehead="";
- int i,dotpos;
-
- switch (cflag) {
- case 'o':
- switch (rep) {
- case Generator:
- extension=".ine"; break; /* output file for ine data */
- case Inequality:
- extension=".ext"; break; /* output file for ext data */
- default:
- extension=".xxx";break;
- }
- break;
-
- case 'a': /* decide for output adjacence */
- if (rep==Inequality)
- extension=".ead"; /* adjacency file for ext data */
- else
- extension=".iad"; /* adjacency file for ine data */
- break;
- case 'i': /* decide for output incidence */
- if (rep==Inequality)
- extension=".ecd"; /* ext incidence file */
- else
- extension=".icd"; /* ine incidence file */
- break;
- case 'n': /* decide for input incidence */
- if (rep==Inequality)
- extension=".icd"; /* ine incidence file */
- else
- extension=".ecd"; /* ext incidence file */
- break;
- case 'j': /* decide for input adjacency */
- if (rep==Inequality)
- extension=".iad"; /* ine adjacency file */
- else
- extension=".ead"; /* ext adjacency file */
- break;
- case 'l':
- extension=".ddl";break; /* log file */
- case 'd':
- extension=".dex";break; /* decomposition output */
- case 'p':
- extension="sub.ine";break; /* preprojection sub inequality file */
- case 'v':
- extension=".solved";break; /* verify_input file */
- case 's':
- extension=".lps";break; /* LP solution file */
- default:
- extension=".xxx";break;
- }
- dotpos=-1;
- for (i=0; i< strlen(inputfile); i++){
- if (inputfile[i]=='.') dotpos=i;
- }
- if (dotpos>1) strncpy(ifilehead, inputfile, dotpos);
- else strcpy(ifilehead,inputfile);
- if (strlen(inputfile)<=0) strcpy(ifilehead,"tempcdd");
- strcpy(outfile,ifilehead);
- strcat(outfile,extension);
- if (strcmp(inputfile, outfile)==0) {
- strcpy(outfile,inputfile);
- strcat(outfile,extension);
- }
- /* fprintf(stderr,"outfile name = %s\n",outfile); */
- }
-
-
- NumberType GetNumberType(char *line)
- {
- NumberType nt;
-
- if (strncmp(line, "integer", 7)==0) {
- nt = Integer;
- }
- else if (strncmp(line, "rational", 8)==0) {
- nt = Rational;
- }
- else if (strncmp(line, "real", 4)==0) {
- nt = Real;
- }
- else {
- nt=Unknown;
- }
- return nt;
- }
-
- void ProcessCommandLine(FILE *f, Matrix M, char *line)
- {
- char newline[linelenmax];
- int j;
- double value;
-
- init(value);
- if (strncmp(line, "hull", 4)==0) {
- M.representation = Generator;
- }
- if (strncmp(line, "debug", 5)==0) {
- debug = true;
- //#ifdef GMPRATIONAL
- ddf_debug = ddf_TRUE;
- //#endif
- }
- if (strncmp(line, "partial_enum", 12)==0 ||
- strncmp(line, "equality", 8)==0 ||
- strncmp(line, "linearity", 9)==0 ) {
- fgets(newline,linelenmax,f);
- SetLinearity(M,newline);
- }
- if (strncmp(line, "maximize", 8)==0 ||
- strncmp(line, "minimize", 8)==0) {
- if (strncmp(line, "maximize", 8)==0) M.objective=LPmax;
- else M.objective=LPmin;
- for (j = 1; j <= M.colsize; j++) {
- if (M.numbtype==Real) {
- //#if !defined(GMPRATIONAL)
- double rvalue;
- fscanf(f, "%lf", &rvalue);
- set_d(value, rvalue);
- //#endif
- } else {
- fread_rational_value (f, value);
- }
- set(M.rowvec[j - 1],value);
- if (debug) {fprintf(stderr,"cost(%5ld) =",j); WriteNumber(stderr,value);}
- } /*of j*/
- }
- clear(value);
- }
- #endif
-
- public bool AppendMatrix2Poly(Polyhedra poly, Matrix M)
- {
- bool success = false;
- Matrix Mpoly, Mnew = null;
- ErrorType err;
-
- if (poly != null && poly.m >= 0 && poly.d >= 0 &&
- poly.d == M.colsize && M.rowsize > 0)
- {
- Mpoly = CopyInput(poly);
- Mnew = AppendMatrix(Mpoly, M);
-
- poly = DDMatrix2Poly(Mnew, out err);
-
- if (err == ErrorType.NoError) success = true;
- }
- return success;
- }
-
- public Matrix MatrixCopy(Matrix M)
- {
- Matrix Mcopy = null;
- int m;
- int d;
-
- m = M.rowsize;
- d = M.colsize;
- if (m >= 0 && d >= 0)
- {
- Mcopy = CreateMatrix(m, d);
- CopyAmatrix(Mcopy.matrix, M.matrix, m, d);
- CopyArow(Mcopy.rowvec, M.rowvec, d);
- M.linset.CopyTo(Mcopy.linset);
- Mcopy.numbtype = M.numbtype;
- Mcopy.representation = M.representation;
- Mcopy.objective = M.objective;
- }
- return Mcopy;
- }
-
- public Matrix CopyMatrix(Matrix M)
- {
- return MatrixCopy(M);
- }
-
- public Matrix MatrixNormalizedCopy(Matrix M)
- {
- Matrix Mcopy = null;
- int m;
- int d;
-
- m = M.rowsize;
- d = M.colsize;
- if (m >= 0 && d >= 0)
- {
- Mcopy = CreateMatrix(m, d);
- CopyNormalizedAmatrix(Mcopy.matrix, M.matrix, m, d);
- CopyArow(Mcopy.rowvec, M.rowvec, d);
- M.linset.CopyTo(Mcopy.linset);
- Mcopy.numbtype = M.numbtype;
- Mcopy.representation = M.representation;
- Mcopy.objective = M.objective;
- }
- return Mcopy;
- }
-
-
- public Matrix MatrixAppend(Matrix M1, Matrix M2)
- {
- Matrix M = null;
- int i, m, m1, m2;
- int j, d, d1, d2;
-
- m1 = M1.rowsize;
- d1 = M1.colsize;
- m2 = M2.rowsize;
- d2 = M2.colsize;
-
- m = m1 + m2;
- d = d1;
-
- if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)
- {
- M = CreateMatrix(m, d);
- CopyAmatrix(M.matrix, M1.matrix, m1, d);
- CopyArow(M.rowvec, M1.rowvec, d);
- for (i = 0; i < m1; i++)
- {
- if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);
- }
- for (i = 0; i < m2; i++)
- {
- for (j = 0; j < d; j++)
- M.matrix[m1 + i][j] = M2.matrix[i][j];
- /* append the second matrix */
- if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);
- }
- M.numbtype = M1.numbtype;
- }
- return M;
- }
-
- public Matrix MatrixNormalizedSortedCopy(Matrix M, out int[] newpos) /* 094 */
- {
- /* Sort the rows of double[][] lexicographically, and return a link to this sorted copy.
- The vector newpos is allocated, where newpos[i] returns the new row index
- of the original row i (i=1,...,M.rowsize). */
- Matrix Mcopy = null, Mnorm = null;
- int m, i;
- int d;
- int[] roworder;
-
- /* if (newpos!=null) free(newpos); */
- m = M.rowsize;
- d = M.colsize;
- roworder = new int[m + 1];
- newpos = new int[m + 1];
- if (m >= 0 && d >= 0)
- {
- Mnorm = MatrixNormalizedCopy(M);
- Mcopy = CreateMatrix(m, d);
- for (i = 1; i <= m; i++) roworder[i] = i;
-
- RandomPermutation(roworder, m, 123);
- QuickSort(roworder, 1, m, Mnorm.matrix, d);
-
- PermuteCopyAmatrix(Mcopy.matrix, Mnorm.matrix, m, d, roworder);
- CopyArow(Mcopy.rowvec, M.rowvec, d);
- for (i = 1; i <= m; i++)
- {
- if (M.linset.Get(roworder[i])) Mcopy.linset.Set(i, true);
- newpos[roworder[i]] = i;
- }
- Mcopy.numbtype = M.numbtype;
- Mcopy.representation = M.representation;
- Mcopy.objective = M.objective;
-
- }
-
- return Mcopy;
- }
-
- public Matrix MatrixUniqueCopy(Matrix M, out int[] newpos)
- {
- /* Remove row duplicates, and return a link to this sorted copy.
- Linearity rows have priority over the other rows.
- It is better to call this after sorting with MatrixNormalizedSortedCopy.
- The vector newpos is allocated, where *newpos[i] returns the new row index
- of the original row i (i=1,...,M.rowsize). *newpos[i] is negative if the original
- row is dominated by -*newpos[i] and eliminated in the new copy.
- */
- Matrix Mcopy = null;
- int m, i, uniqrows;
- BitArray preferredrows;
- int d;
- int[] roworder;
-
- /* if (newpos!=null) free(newpos); */
- m = M.rowsize;
- d = M.colsize;
- preferredrows = M.linset;
- roworder = new int[m + 1];
- if (m >= 0 && d >= 0)
- {
- for (i = 1; i <= m; i++) roworder[i] = i;
- UniqueRows(roworder, 1, m, M.matrix, d, preferredrows, out uniqrows);
-
- Mcopy = CreateMatrix(uniqrows, d);
- PermutePartialCopyAmatrix(Mcopy.matrix, M.matrix, m, d, roworder, 1, m);
- CopyArow(Mcopy.rowvec, M.rowvec, d);
- for (i = 1; i <= m; i++)
- {
- if (roworder[i] > 0 && M.linset.Get(i)) Mcopy.linset.Set(roworder[i], true);
- }
- Mcopy.numbtype = M.numbtype;
- Mcopy.representation = M.representation;
- Mcopy.objective = M.objective;
- }
- newpos = roworder;
- return Mcopy;
- }
-
-
- public Matrix MatrixNormalizedSortedUniqueCopy(Matrix M, out int[] newpos) /* 094 */
- {
- /* Sort and remove row duplicates, and return a link to this sorted copy.
- Linearity rows have priority over the other rows.
- It is better to call this after sorting with MatrixNormalizedSortedCopy.
- The vector newpos is allocated, where *newpos[i] returns the new row index
- of the original row i (i=1,...,M.rowsize). *newpos[i] is negative if the original
- row is dominated by -*newpos[i] and eliminated in the new copy.
- */
- Matrix M1 = null, M2 = null;
- int m, i;
- int d;
- int[] newpos1 = null, newpos1r = null, newpos2 = null;
-
- /* if (newpos!=null) free(newpos); */
- m = M.rowsize;
- d = M.colsize;
- newpos = new int[m + 1];
- newpos1r = new int[m + 1];
- if (m >= 0 && d >= 0)
- {
- M1 = MatrixNormalizedSortedCopy(M, out newpos1);
- for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i; /* reverse of newpos1 */
- M2 = MatrixUniqueCopy(M1, out newpos2);
- M2.linset.SetAll(false); //emptyset
- for (i = 1; i <= m; i++)
- {
- if (newpos2[newpos1[i]] > 0)
- {
- //printf("newpos1[%ld]=%ld, newpos2[newpos1[%ld]]=%ld\n",i,newpos1[i], i,newpos2[newpos1[i]]);
- if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);
- newpos[i] = newpos2[newpos1[i]];
- }
- else
- {
- newpos[i] = -newpos1r[-newpos2[newpos1[i]]];
- }
- }
-
- }
-
- return M2;
- }
-
- public Matrix MatrixSortedUniqueCopy(Matrix M, out int[] newpos) /* 094 */
- {
- /* Same as MatrixNormalizedSortedUniqueCopy except that it returns a unnormalized origial data
- with original ordering.
- */
- Matrix M1 = null, M2 = null;
- int m, i, k, ii;
- int d;
- int[] newpos1 = null, newpos1r = null, newpos2 = null;
-
- /* if (newpos!=null) free(newpos); */
- m = M.rowsize;
- d = M.colsize;
- newpos = new int[m + 1];
- newpos1r = new int[m + 1];
- if (m >= 0 && d >= 0)
- {
- M1 = MatrixNormalizedSortedCopy(M, out newpos1);
- for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i; /* reverse of newpos1 */
- M2 = MatrixUniqueCopy(M1, out newpos2);
-
- M2.linset.SetAll(false);
- for (i = 1; i <= m; i++)
- {
- if (newpos2[newpos1[i]] > 0)
- {
- if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);
- newpos[i] = newpos2[newpos1[i]];
- }
- else
- {
- newpos[i] = -newpos1r[-newpos2[newpos1[i]]];
- }
- }
-
- ii = 0;
- M2.linset.SetAll(false);
- for (i = 1; i <= m; i++)
- {
- k = newpos[i];
- if (k > 0)
- {
- ii += 1;
- newpos[i] = ii;
- CopyArow(M2.matrix[ii - 1], M.matrix[i - 1], d);
- if (M.linset.Get(i)) M2.linset.Set(ii, true);
- }
- }
-
- }
-
- return M2;
- }
-
- public Matrix AppendMatrix(Matrix M1, Matrix M2)
- {
- return MatrixAppend(M1, M2);
- }
-
- public bool MatrixAppendTo(ref Matrix M1, Matrix M2)
- {
- Matrix M = null;
- int i, m, m1, m2;
- int j, d, d1, d2;
- bool success = false;
-
- m1 = M1.rowsize;
- d1 = M1.colsize;
- m2 = M2.rowsize;
- d2 = M2.colsize;
-
- m = m1 + m2;
- d = d1;
-
- if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)
- {
- M = CreateMatrix(m, d);
- CopyAmatrix(M.matrix, M1.matrix, m1, d);
- CopyArow(M.rowvec, M1.rowvec, d);
- for (i = 0; i < m1; i++)
- {
- if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);
- }
- for (i = 0; i < m2; i++)
- {
- for (j = 0; j < d; j++)
- M.matrix[m1 + i][j] = M2.matrix[i][j];
- /* append the second matrix */
- if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);
- }
- M.numbtype = M1.numbtype;
-
- M1 = M;
- success = true;
- }
- return success;
- }
-
- public bool MatrixRowRemove(Matrix M, int r) /* 092 */
- {
- int i, m;
- int d;
- bool success = false;
-
- m = M.rowsize;
- d = M.colsize;
-
- if (r >= 1 && r <= m)
- {
- M.rowsize = m - 1;
-
- M.linset.Set(r, false);
- /* slide the row headers */
- for (i = r; i < m; i++)
- {
- M.matrix[i - 1] = M.matrix[i];
- if (M.linset.Get(i + 1))
- {
- M.linset.Set(i + 1, false);
- M.linset.Set(i, true);
- }
- }
- success = true;
- }
- return success;
- }
-
- public bool MatrixRowRemove2(Matrix M, int r, out int[] newpos) /* 094 */
- {
- int i, m;
- int d;
- bool success = false;
- int[] roworder;
-
- newpos = null;
-
- m = M.rowsize;
- d = M.colsize;
-
- if (r >= 1 && r <= m)
- {
- roworder = new int[m + 1];
- M.rowsize = m - 1;
-
- M.linset.Set(r, false);
- /* slide the row headers */
- for (i = 1; i < r; i++) roworder[i] = i;
- roworder[r] = 0; /* meaning it is removed */
- for (i = r; i < m; i++)
- {
- M.matrix[i - 1] = M.matrix[i];
- roworder[i + 1] = i;
- if (M.linset.Get(i + 1))
- {
- M.linset.Set(i + 1, false);
- M.linset.Set(i, true);
- }
- }
- success = true;
- }
- return success;
- }
-
- public Matrix MatrixSubmatrix(Matrix M, BitArray delset) /* 092 */
- {
- Matrix Msub = null;
- int i, isub = 1, m, msub;
- int d;
-
- m = M.rowsize;
- d = M.colsize;
- msub = m;
- if (m >= 0 && d >= 0)
- {
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i)) msub -= 1;
- }
- Msub = CreateMatrix(msub, d);
- for (i = 1; i <= m; i++)
- {
- if (!delset.Get(i))
- {
- CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);
- if (M.linset.Get(i))
- {
- Msub.linset.Set(isub, true);
- }
- isub++;
- }
- }
- CopyArow(Msub.rowvec, M.rowvec, d);
- Msub.numbtype = M.numbtype;
- Msub.representation = M.representation;
- Msub.objective = M.objective;
- }
- return Msub;
- }
-
- public Matrix MatrixSubmatrix2(Matrix M, BitArray delset, out int[] newpos) /* 092 */
- { /* returns a pointer to a new matrix which is a submatrix of M with rows in delset
- removed. *newpos[i] returns the position of the original row i in the new matrix.
- It is -1 if and only if it is deleted.
- */
-
- Matrix Msub = null;
- int i, isub = 1, m, msub;
- int d;
- int[] roworder;
-
- newpos = null;
-
- m = M.rowsize;
- d = M.colsize;
- msub = m;
- if (m >= 0 && d >= 0)
- {
- roworder = new int[m + 1];
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i)) msub -= 1;
- }
- Msub = CreateMatrix(msub, d);
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i))
- {
- roworder[i] = 0; /* zero means the row i is removed */
- }
- else
- {
- CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);
- if (M.linset.Get(i))
- {
- Msub.linset.Set(isub, true);
- }
- roworder[i] = isub;
- isub++;
- }
- }
- newpos = roworder;
- CopyArow(Msub.rowvec, M.rowvec, d);
- Msub.numbtype = M.numbtype;
- Msub.representation = M.representation;
- Msub.objective = M.objective;
- }
- return Msub;
- }
-
-
- public Matrix MatrixSubmatrix2L(Matrix M, BitArray delset, out int[] newpos) /* 094 */
- { /* This is same as MatrixSubmatrix2 except that the linearity rows will be shifted up
- so that they are at the top of the matrix.
- */
- Matrix Msub = null;
- int i, iL, iI, m, msub;
- int d;
- int[] roworder;
-
- newpos = null;
-
- m = M.rowsize;
- d = M.colsize;
- msub = m;
- if (m >= 0 && d >= 0)
- {
- roworder = new int[m + 1];
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i)) msub -= 1;
- }
- Msub = CreateMatrix(msub, d);
- iL = 1; iI = M.linset.Count + 1; /* starting positions */
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i))
- {
- roworder[i] = 0; /* zero means the row i is removed */
- }
- else
- {
- if (M.linset.Get(i))
- {
- CopyArow(Msub.matrix[iL - 1], M.matrix[i - 1], d);
- Msub.linset.Set(i, false);
- Msub.linset.Set(iL, true);
- roworder[i] = iL;
- iL += 1;
- }
- else
- {
- CopyArow(Msub.matrix[iI - 1], M.matrix[i - 1], d);
- roworder[i] = iI;
- iI += 1;
- }
- }
- }
- newpos = roworder;
- CopyArow(Msub.rowvec, M.rowvec, d);
- Msub.numbtype = M.numbtype;
- Msub.representation = M.representation;
- Msub.objective = M.objective;
- }
- return Msub;
- }
-
- public int MatrixRowsRemove(ref Matrix M, BitArray delset) /* 094 */
- {
- Matrix Msub = null;
- int success;
-
- Msub = MatrixSubmatrix(M, delset);
-
- M = Msub;
- success = 1;
- return success;
- }
-
- public int MatrixRowsRemove2(ref Matrix M, BitArray delset, out int[] newpos) /* 094 */
- {
- Matrix Msub = null;
- int success;
-
- Msub = MatrixSubmatrix2(M, delset, out newpos);
-
- M = Msub;
- success = 1;
- return success;
- }
-
- public int MatrixShiftupLinearity(ref Matrix M, out int[] newpos) /* 094 */
- {
- Matrix Msub = null;
- int success;
- BitArray delset;
-
- delset = new BitArray(M.rowsize); /* emptyset */
- Msub = MatrixSubmatrix2L(M, delset, out newpos);
-
- M = Msub;
-
- success = 1;
- return success;
- }
-
- public Polyhedra CreatePolyhedraData(int m, int d)
- {
- int i;
- Polyhedra poly = null;
-
- poly = new Polyhedra();
- poly.child = null; /* this links the homogenized cone data */
- poly.m = m;
- poly.d = d;
- poly.n = -1; /* the size of output is not known */
- poly.m_alloc = m + 2; /* the allocated row size of matrix A */
- poly.d_alloc = d; /* the allocated col size of matrix A */
- poly.numbtype = NumberType.Real;
- InitializeAmatrix(poly.m_alloc, poly.d_alloc, out poly.A);
- InitializeArow(d, out poly.c); /* cost vector */
- poly.representation = RepresentationType.Inequality;
- poly.homogeneous = false;
-
- poly.EqualityIndex = new int[m + 2];
- /* size increased to m+2 in 092b because it is used by the child cone,
- This is a bug fix suggested by Thao Dang. */
- /* ith component is 1 if it is equality, -1 if it is strict inequality, 0 otherwise. */
- for (i = 0; i <= m + 1; i++) poly.EqualityIndex[i] = 0;
-
- poly.IsEmpty = -1; /* initially set to -1, neither true nor false, meaning unknown */
-
- poly.NondegAssumed = false;
- poly.InitBasisAtBottom = false;
- poly.RestrictedEnumeration = false;
- poly.RelaxedEnumeration = false;
-
- poly.AincGenerated = false; /* Ainc is a set array to store the input incidence. */
-
- return poly;
- }
-
- public bool InitializeConeData(int m, int d, out Cone cone)
- {
- bool success = true;
- int j;
-
- cone = new Cone();
-
- /* INPUT: A given representation of a cone: inequality */
- cone.m = m;
- cone.d = d;
- cone.m_alloc = m + 2; /* allocated row size of matrix A */
- cone.d_alloc = d; /* allocated col size of matrix A, B and Bsave */
- cone.numbtype = NumberType.Real;
- cone.parent = null;
-
- /* CONTROL: variables to control computation */
- cone.Iteration = 0;
-
- cone.HalfspaceOrder = RowOrderType.LexMin;
-
- cone.ArtificialRay = null;
- cone.FirstRay = null;
- cone.LastRay = null; /* The second description: Generator */
- cone.PosHead = null;
- cone.ZeroHead = null;
- cone.NegHead = null;
- cone.PosLast = null;
- cone.ZeroLast = null;
- cone.NegLast = null;
- cone.RecomputeRowOrder = true;
- cone.PreOrderedRun = false;
- cone.GroundSet = new BitArray(cone.m_alloc);
- cone.EqualitySet = new BitArray(cone.m_alloc);
- cone.NonequalitySet = new BitArray(cone.m_alloc);
- cone.AddedHalfspaces = new BitArray(cone.m_alloc);
- cone.WeaklyAddedHalfspaces = new BitArray(cone.m_alloc);
- cone.InitialHalfspaces = new BitArray(cone.m_alloc);
- cone.RayCount = 0;
- cone.FeasibleRayCount = 0;
- cone.WeaklyFeasibleRayCount = 0;
- cone.TotalRayCount = 0;
- cone.ZeroRayCount = 0;
- cone.EdgeCount = 0;
- cone.TotalEdgeCount = 0;
- cone.count_int = 0;
- cone.count_int_good = 0;
- cone.count_int_bad = 0;
- cone.rseed = 1; /* random seed for random row permutation */
-
- InitializeBmatrix(cone.d_alloc, out cone.B);
- InitializeBmatrix(cone.d_alloc, out cone.Bsave);
- InitializeAmatrix(cone.m_alloc, cone.d_alloc, out cone.A);
-
- cone.Edges = new Adjacency[cone.m_alloc];
- cone.InitialRayIndex = new int[d + 1];
- cone.OrderVector = new int[cone.m_alloc + 1];
-
- cone.newcol = new int[cone.d + 1];
- for (j = 0; j <= cone.d; j++) cone.newcol[j] = j; /* identity map, initially */
- cone.LinearityDim = -2; /* -2 if it is not computed */
- cone.ColReduced = false;
- cone.d_orig = d;
-
- /* STATES: variables to represent current state. */
- /*cone.Error;
- cone.CompStatus;
- cone.starttime;
- cone.endtime;
- */
-
- return success;
- }
-
- public Cone ConeDataLoad(Polyhedra poly)
- {
- Cone cone = null;
- int d, j;
- int m, i;
-
- m = poly.m;
- d = poly.d;
- if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)
- {
- m = poly.m + 1;
- }
- poly.m1 = m;
-
- InitializeConeData(m, d, out cone);
- cone.representation = poly.representation;
-
- /* Points to the original polyhedra data, and reversely */
- cone.parent = poly;
- poly.child = cone;
-
- for (i = 1; i <= poly.m; i++)
- for (j = 1; j <= cone.d; j++)
- cone.A[i - 1][j - 1] = poly.A[i - 1][j - 1];
-
- if (poly.representation == RepresentationType.Inequality && !(poly.homogeneous))
- {
- cone.A[m - 1][0] = 1.0;
- for (j = 2; j <= d; j++) cone.A[m - 1][j - 1] = 0.0;
- }
-
- return cone;
- }
-
- #if nemkell
- void SetLinearity(Matrix M, char *line)
- {
- int i=0;
- int eqsize,var;
- char *next;
- const char ct[]=", "; /* allows separators "," and " ". */
-
- next=strtok(line,ct);
- eqsize=atol(next);
- while (i < eqsize && (next=strtok(null,ct))!=null) {
- var=atol(next);
- set_addelem(M.linset,var); i++;
- }
- if (i!=eqsize) {
- fprintf(stderr,"* Warning: there are inconsistencies in linearity setting.\n");
- }
- return;
- }
-
- Matrix PolyFile2Matrix (FILE *f, ErrorType *Error)
- {
- Matrix M=null;
- int m_input,i;
- int d_input,j;
- RepresentationType rep=Inequality;
- double value;
- bool found=false, newformat=false, successful=false, linearity=false;
- char command[linelenmax], comsave[linelenmax], numbtype[wordlenmax];
- NumberType NT;
- //#if !defined(GMPRATIONAL)
- double rvalue;
- //#endif
-
- init(value);
- (*Error)=NoError;
- while (!found)
- {
- if (fscanf(f,"%s",command)==EOF) {
- (*Error)=ImproperInputFormat;
- goto _L99;
- }
- else {
- if (strncmp(command, "V-representation", 16)==0) {
- rep=Generator; newformat=true;
- }
- if (strncmp(command, "H-representation", 16)==0){
- rep=Inequality; newformat=true;
- }
- if (strncmp(command, "partial_enum", 12)==0 ||
- strncmp(command, "equality", 8)==0 ||
- strncmp(command, "linearity", 9)==0 ) {
- linearity=true;
- fgets(comsave,linelenmax,f);
- }
- if (strncmp(command, "begin", 5)==0) found=true;
- }
- }
- fscanf(f, "%ld %ld %s", &m_input, &d_input, numbtype);
- fprintf(stderr,"size = %ld x %ld\nNumber Type = %s\n", m_input, d_input, numbtype);
- NT=GetNumberType(numbtype);
- if (NT==Unknown) {
- (*Error)=ImproperInputFormat;
- goto _L99;
- }
- M=CreateMatrix(m_input, d_input);
- M.representation=rep;
- M.numbtype=NT;
-
- for (i = 1; i <= m_input; i++) {
- for (j = 1; j <= d_input; j++) {
- if (NT==Real) {
- //#if defined GMPRATIONAL
- *Error=NoRealNumberSupport;
- goto _L99;
- //#else
- fscanf(f, "%lf", &rvalue);
- set_d(value, rvalue);
- //#endif
- } else {
- fread_rational_value (f, value);
- }
- set(M.matrix[i-1][j - 1],value);
- if (debug) {fprintf(stderr,"a(%3ld,%5ld) = ",i,j); WriteNumber(stderr,value);}
- } /*of j*/
- } /*of i*/
- if (fscanf(f,"%s",command)==EOF) {
- (*Error)=ImproperInputFormat;
- goto _L99;
- }
- else if (strncmp(command, "end", 3)!=0) {
- if (debug) fprintf(stderr,"'end' missing or illegal extra data: %s\n",command);
- (*Error)=ImproperInputFormat;
- goto _L99;
- }
-
- successful=true;
- if (linearity) {
- SetLinearity(M,comsave);
- }
- while (!feof(f)) {
- fscanf(f,"%s", command);
- ProcessCommandLine(f, M, command);
- fgets(command,linelenmax,f); /* skip the CR/LF */
- }
-
- _L99: ;
- clear(value);
- /* if (f!=null) fclose(f); */
- return M;
- }
-
- #endif
-
- public Polyhedra DDMatrix2Poly(Matrix M, out ErrorType err)
- {
- int i;
- int j;
- Polyhedra poly = null;
-
- err = ErrorType.NoError;
- if (M.rowsize < 0 || M.colsize < 0)
- {
- err = ErrorType.NegativeMatrixSize;
- return null;
- }
- poly = CreatePolyhedraData(M.rowsize, M.colsize);
- poly.representation = M.representation;
- poly.homogeneous = true;
-
- for (i = 1; i <= M.rowsize; i++)
- {
- if (M.linset.Get(i))
- {
- poly.EqualityIndex[i] = 1;
- }
- for (j = 1; j <= M.colsize; j++)
- {
- poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
- if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;
- } /*of j*/
- } /*of i*/
- DoubleDescription(poly, out err);
-
- return poly;
- }
-
- public Polyhedra DDMatrix2Poly2(Matrix M, RowOrderType horder, out ErrorType err)
- {
- int i;
- int j;
- Polyhedra poly = null;
-
- err = ErrorType.NoError;
- if (M.rowsize < 0 || M.colsize < 0)
- {
- err = ErrorType.NegativeMatrixSize;
- return null;
- }
- poly = CreatePolyhedraData(M.rowsize, M.colsize);
- poly.representation = M.representation;
- poly.homogeneous = true;
-
- for (i = 1; i <= M.rowsize; i++)
- {
- if (M.linset.Get(i))
- {
- poly.EqualityIndex[i] = 1;
- }
- for (j = 1; j <= M.colsize; j++)
- {
- poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
- if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;
- } /*of j*/
- } /*of i*/
- DoubleDescription2(poly, horder, out err);
-
- return poly;
- }
-
- void MatrixIntegerFilter(Matrix M)
- { /* setting an almost integer to the integer. */
- int i;
- int j;
- double x;
-
- x = 0.0;
- for (i = 0; i < M.rowsize; i++)
- for (j = 0; j < M.colsize; j++)
- {
- SnapToInteger(out x, M.matrix[i][j]);
- M.matrix[i][j] = x;
- }
-
- }
-
- void CopyRay(double[] a, int d_origsize, Ray RR,
- RepresentationType rep, int[] reducedcol)
- {
- int j, j1;
- double b;
-
- b = 0.0;
- for (j = 1; j <= d_origsize; j++)
- {
- j1 = reducedcol[j];
- if (j1 > 0)
- {
- a[j - 1] = RR.RRay[j1 - 1];
- /* the original column j is mapped to j1, and thus
- copy the corresponding component */
- }
- else
- {
- a[j - 1] = 0.0;
- /* original column is redundant and removed for computation */
- }
- }
-
- b = a[0];
- if (rep == RepresentationType.Generator && b != 0)
- {
- a[0] = 1.0;
- for (j = 2; j <= d_origsize; j++)
- a[j - 1] /= b; /* normalization for generators */
- }
- }
-
-
- void WriteRay(TextWriter f, int d_origsize, Ray RR, RepresentationType rep, int[] reducedcol)
- {
- int j;
- double[] a;
-
-
- InitializeArow(d_origsize + 1, out a);
-
- CopyRay(a, d_origsize, RR, rep, reducedcol);
- for (j = 0; j < d_origsize; j++) WriteNumber(f, a[j]);
- f.WriteLine();
- }
-
- void WriteArow(TextWriter f, double[] a, int d)
- {
- int j;
-
- for (j = 0; j < d; j++) WriteNumber(f, a[j]);
- f.WriteLine();
- }
-
- void WriteAmatrix(TextWriter f, double[][] A, int rowmax, int colmax)
- {
- int i, j;
-
- if (A == null)
- {
- f.WriteLine("WriteAmatrix: The requested matrix is empty");
- return;
- }
- f.WriteLine("begin");
- #if GMPRATIONAL
- fprintf(f, " %ld %ld rational\n",rowmax, colmax);
- #else
- f.WriteLine(" {0} {1} real", rowmax, colmax);
- #endif
- for (i = 1; i <= rowmax; i++)
- {
- for (j = 1; j <= colmax; j++)
- {
- WriteNumber(f, A[i - 1][j - 1]);
- }
- f.WriteLine();
- }
- f.WriteLine("end");
-
- }
-
- void WriteBmatrix(TextWriter f, int d_size, double[][] B)
- {
- int j1, j2;
-
- if (B == null)
- {
- f.WriteLine("WriteBmatrix: The requested matrix is empty");
- return;
- }
- for (j1 = 0; j1 < d_size; j1++)
- {
- for (j2 = 0; j2 < d_size; j2++)
- {
- WriteNumber(f, B[j1][j2]);
- } /*of j2*/
- f.WriteLine();
- } /*of j1*/
- f.WriteLine();
-
- }
- #if nemkell
- void WriteSetFamily(TextWriter f, SetFamily F)
- {
- int i;
-
- if (F==null){
- f.WriteLine("WriteSetFamily: The requested family is empty");
- return;
- }
- f.WriteLine("begin");
- f.WriteLine(" {0} {1}", F.famsize, F.setsize);
- for (i=0; i<F.famsize; i++) {
- f.Write(" {0} {1} : ", i+1, F.set[i].Count);
- //*****set_fwrite(f, F.set[i]);
- }
- f.WriteLine("end");
-
- }
-
- void WriteSetFamilyCompressed(FILE *f, SetFamily F)
- {
- int i,card;
-
- if (F==null){
- fprintf(f, "WriteSetFamily: The requested family is empty\n");
- goto _L99;
- }
- fprintf(f,"begin\n");
- fprintf(f," %ld %ld\n", F.famsize, F.setsize);
- for (i=0; i<F.famsize; i++) {
- card=set_card(F.set[i]);
- if (F.setsize - card >= card){
- fprintf(f, " %ld %ld : ", i+1, card);
- set_fwrite(f, F.set[i]);
- } else {
- fprintf(f, " %ld %ld : ", i+1, -card);
- set_fwrite_compl(f, F.set[i]);
- }
- }
- fprintf(f,"end\n");
- _L99:;
- }
- #endif
- void WriteMatrix(TextWriter f, Matrix M)
- {
- int i, linsize;
-
- if (M == null)
- {
- f.WriteLine("WriteMatrix: The requested matrix is empty");
- return;
- }
- switch (M.representation)
- {
- case RepresentationType.Inequality:
- f.WriteLine("H-representation"); break;
- case RepresentationType.Generator:
- f.WriteLine("V-representation"); break;
- case RepresentationType.Unspecified:
- break;
- }
- linsize = M.linset.Count;
- if (linsize > 0)
- {
- f.Write("linearity {0} ", linsize);
- for (i = 1; i <= M.rowsize; i++)
- if (M.linset.Get(i)) f.Write(" {0}", i);
- f.WriteLine();
- }
- WriteAmatrix(f, M.matrix, M.rowsize, M.colsize);
- if (M.objective != LPObjectiveType.LPnone)
- {
- if (M.objective == LPObjectiveType.LPmax)
- f.WriteLine("maximize");
- else
- f.WriteLine("minimize");
- WriteArow(f, M.rowvec, M.colsize);
- }
-
- }
-
- void WriteLP(TextWriter f, LP lp)
- {
- if (lp == null)
- {
- f.WriteLine("WriteLP: The requested lp is empty");
- return;
- }
- f.WriteLine("H-representation");
-
- WriteAmatrix(f, lp.A, (lp.m) - 1, lp.d);
- if (lp.objective != LPObjectiveType.LPnone)
- {
- if (lp.objective == LPObjectiveType.LPmax)
- f.WriteLine("maximize");
- else
- f.WriteLine("minimize");
- WriteArow(f, lp.A[lp.objrow - 1], lp.d);
- }
-
- }
-
-
- void SnapToInteger(out double y, double x)
- {
- /* this is broken. It does nothing. */
- y = x;
- }
-
-
- void WriteReal(TextWriter f, double x)
- {
- /*
- int ix1,ix2,ix;
- double ax;
-
- ax=x;
- ix1= (int) (Math.Abs(ax) * 10000.0 + 0.5);
- ix2= (int) (Math.Abs(ax) + 0.5);
- ix2= ix2*10000;
- if ( ix1 == ix2) {
- if (x > 0) {
- ix = (int)(ax + 0.5);
- } else {
- ix = (int)(-ax + 0.5);
- ix = -ix;
- }
- f.Write(" {0}", ix); //%2ld
- } else
- f.Write(" {0}.9E",ax); //% .9E
- * */
-
- f.Write(x);
- }
-
- void WriteNumber(TextWriter f, double x)
- {
- #if GMPRATIONAL
- mpz_t zn,zd;
-
- mpz_init(zn); mpz_init(zd);
- mpq_canonicalize(x);
- mpq_get_num(zn,x);
- mpq_get_den(zd,x);
- fprintf(f," ");
- if (mpz_sgn(zn)==0){
- fprintf(f,"0");
- } else if (mpz_cmp_ui(zd,1U)==0){
- mpz_out_str(f,10,zn);
- } else {
- mpz_out_str(f,10,zn);fprintf(f,"/");mpz_out_str(f,10,zd);
- }
- mpz_clear(zn); mpz_clear(zd);
- #else
- WriteReal(f, x);
- #endif
- }
-
- #if nemkell
- void WriteIncidence(FILE *f, Polyhedra poly)
- {
- SetFamily I;
-
- switch (poly.representation) {
- case Inequality:
- fprintf(f, "ecd_file: Incidence of generators and inequalities\n");
- break;
- case Generator:
- fprintf(f, "icd_file: Incidence of inequalities and generators\n");
- break;
-
- default:
- break;
- }
- I=CopyIncidence(poly);
- WriteSetFamilyCompressed(f,I);
- FreeSetFamily(I);
- }
-
- void WriteAdjacency(FILE *f, Polyhedra poly)
- {
- SetFamily A;
-
- switch (poly.representation) {
- case Inequality:
- fprintf(f, "ead_file: Adjacency of generators\n");
- break;
- case Generator:
- fprintf(f, "iad_file: Adjacency of inequalities\n");
- break;
-
- default:
- break;
- }
- A=CopyAdjacency(poly);
- WriteSetFamilyCompressed(f,A);
- FreeSetFamily(A);
- }
- #endif
-
- void ComputeAinc(Polyhedra poly)
- {
- /* This generates the input incidence array poly.Ainc, and
- two sets: poly.Ared, poly.Adom.
- */
- int k;
- int i, m1;
- int j;
- bool redundant;
- Matrix M = null;
- double sum, temp;
-
- sum = 0.0; temp = 0.0;
- if (poly.AincGenerated == true) return;
-
- M = CopyOutput(poly);
- poly.n = M.rowsize;
- m1 = poly.m1;
- /* this number is same as poly.m, except when
- poly is given by nonhomogeneous inequalty:
- !(poly.homogeneous) && poly.representation==Inequality,
- it is poly.m+1. See ConeDataLoad.
- */
- poly.Ainc = new BitArray[m1];
- for (i = 1; i <= m1; i++) poly.Ainc[i - 1] = new BitArray(poly.n);
- poly.Ared = new BitArray(m1);
- poly.Adom = new BitArray(m1);
-
- for (k = 1; k <= poly.n; k++)
- {
- for (i = 1; i <= poly.m; i++)
- {
- sum = 0.0;
- for (j = 1; j <= poly.d; j++)
- {
- temp = poly.A[i - 1][j - 1] * M.matrix[k - 1][j - 1];
- sum += temp;
- }
- if (sum == 0.0)
- {
- poly.Ainc[i - 1].Set(k, true);
- }
- }
- if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)
- {
- if (M.matrix[k - 1][0] == 0.0)
- {
- poly.Ainc[m1 - 1].Set(k, true); /* added infinity inequality (1,0,0,...,0) */
- }
- }
- }
-
- for (i = 1; i <= m1; i++)
- {
- if (poly.Ainc[i - 1].Count == M.rowsize)
- {
- poly.Adom.Set(i, true);
- }
- }
- for (i = m1; i >= 1; i--)
- {
- if (poly.Ainc[i - 1].Count == 0)
- {
- redundant = true;
- poly.Ared.Set(i, true);
- }
- else
- {
- redundant = false;
- for (k = 1; k <= m1; k++)
- {
- if (k != i && !poly.Ared.Get(k) && !poly.Adom.Get(k) &&
- poly.Ainc[i - 1].IsSubsetOf(poly.Ainc[k - 1]))
- { //****
- if (!redundant)
- {
- redundant = true;
- }
- poly.Ared.Set(i, true);
- }
- }
- }
- }
-
- poly.AincGenerated = true;
-
-
- }
-
- bool InputAdjacentQ(Polyhedra poly,
- int i1, int i2)
- /* Before calling this function, RedundantSet must be
- a set of row indices whose removal results in a minimal
- nonredundant system to represent the input polyhedron,
- DominantSet must be the set of row indices which are
- active at every extreme points/rays.
- */
- {
- bool adj = true;
- int i;
- BitArray common;
- int lastn = 0;
-
- if (poly.AincGenerated == false) ComputeAinc(poly);
- if (lastn != poly.n)
- {
- //if (lastn >0) set_free(common);
- common = new BitArray(poly.n);
- lastn = poly.n;
- }
- if (poly.Ared.Get(i1) || poly.Ared.Get(i2))
- {
- adj = false;
- return adj;
- }
- if (poly.Adom.Get(i1) || poly.Adom.Get(i2))
- {
- // dominant inequality is considered adjacencent to all others.
- adj = true;
- return adj;
- }
- common = poly.Ainc[i1 - 1].And(poly.Ainc[i2 - 1]); //**** intersection
- i = 0;
- while (i < poly.m1 && adj == true)
- {
- i++;
- if (i != i1 && i != i2 && !poly.Ared.Get(i) &&
- !poly.Adom.Get(i) && common.IsSubsetOf(poly.Ainc[i - 1]))
- {
- adj = false;
- }
- }
-
- return adj;
- }
-
- #if nemkell
- void WriteInputIncidence(FILE *f, Polyhedra poly)
- {
- SetFamily I;
-
- if (poly.AincGenerated==false) ComputeAinc(poly);
- switch (poly.representation) {
- case Inequality:
- fprintf(f,"icd_file: Incidence of inequalities and generators\n");
- break;
-
- case Generator:
- fprintf(f,"ecd_file: Incidence of generators and inequalities\n");
- break;
-
- default:
- break;
- }
-
- I=CopyInputIncidence(poly);
- WriteSetFamilyCompressed(f,I);
- FreeSetFamily(I);
- }
-
- void WriteInputAdjacency(FILE *f, Polyhedra poly)
- {
- SetFamily A;
-
- if (poly.AincGenerated==false){
- ComputeAinc(poly);
- }
- switch (poly.representation) {
- case Inequality:
- fprintf(f, "iad_file: Adjacency of inequalities\n");
- break;
-
- case Generator:
- fprintf(f, "ead_file: Adjacency of generators\n");
- break;
-
- default:
- break;
- }
- A=CopyInputAdjacency(poly);
- WriteSetFamilyCompressed(f,A);
- FreeSetFamily(A);
- }
-
-
- void WriteProgramDescription(FILE *f)
- {
- fprintf(f, "* cddlib: a double description library:%s\n", DDVERSION);
- fprintf(f, "* compiled for %s arithmetic.\n", ARITHMETIC);
- fprintf(f,"* %s\n",COPYRIGHT);
- }
-
- void WriteTimes(FILE *f, time_t starttime, time_t endtime)
- {
- int ptime,ptime_sec,ptime_minu, ptime_hour;
-
- /*
- ptime=difftime(endtime,starttime);
- This function is ANSI standard, but not available sometime
- */
- ptime=(endtime - starttime);
- /* This is to replace the line above, but it may not give
- correct time in seconds */
- ptime_hour=ptime/3600;
- ptime_minu=(ptime-ptime_hour*3600)/60;
- ptime_sec=ptime%60;
- fprintf(f,"* Computation started at %s",asctime(localtime(&starttime)));
- fprintf(f,"* ended at %s",asctime(localtime(&endtime)));
- fprintf(f,"* Total processor time = %ld seconds\n",ptime);
- fprintf(f,"* = %ld h %ld m %ld s\n",ptime_hour, ptime_minu, ptime_sec);
- }
-
- void WriteDDTimes(FILE *f, Polyhedra poly)
- {
- WriteTimes(f,poly.child.starttime,poly.child.endtime);
- }
-
- void WriteLPTimes(FILE *f, LPPtr lp)
- {
- WriteTimes(f,lp.starttime,lp.endtime);
- }
-
- void WriteLPStats(FILE *f)
- {
- time_t currenttime;
-
- time(¤ttime);
-
- fprintf(f, "\n*--- Statistics of pivots ---\n");
- //#if defined GMPRATIONAL
- fprintf(f, "* f0 = %ld (float basis finding pivots)\n",ddf_statBApivots);
- fprintf(f, "* fc = %ld (float CC pivots)\n",ddf_statCCpivots);
- fprintf(f, "* f1 = %ld (float dual simplex phase I pivots)\n",ddf_statDS1pivots);
- fprintf(f, "* f2 = %ld (float dual simplex phase II pivots)\n",ddf_statDS2pivots);
- fprintf(f, "* f3 = %ld (float anticycling CC pivots)\n",ddf_statACpivots);
- fprintf(f, "* e0 = %ld (exact basis finding pivots)\n",statBApivots);
- fprintf(f, "* ec = %ld (exact CC pivots)\n",statCCpivots);
- fprintf(f, "* e1 = %ld (exact dual simplex phase I pivots)\n",statDS1pivots);
- fprintf(f, "* e2 = %ld (exact dual simplex phase II pivots)\n",statDS2pivots);
- fprintf(f, "* e3 = %ld (exact anticycling CC pivots)\n",statACpivots);
- fprintf(f, "* e4 = %ld (exact basis verification pivots)\n",statBSpivots);
- //#else
- fprintf(f, "f0 = %ld (float basis finding pivots)\n",statBApivots);
- fprintf(f, "fc = %ld (float CC pivots)\n",statCCpivots);
- fprintf(f, "f1 = %ld (float dual simplex phase I pivots)\n",statDS1pivots);
- fprintf(f, "f2 = %ld (float dual simplex phase II pivots)\n",statDS2pivots);
- fprintf(f, "f3 = %ld (float anticycling CC pivots)\n",statACpivots);
- //#endif
- WriteLPMode(f);
- WriteTimes(f,statStartTime, currenttime);
- }
-
- void WriteLPMode(FILE *f)
- {
- fprintf(f, "\n* LP solver: ");
- switch (choiceLPSolverDefault) {
- case DualSimplex:
- fprintf(f, "DualSimplex\n");
- break;
- case CrissCross:
- fprintf(f, "Criss-Cross\n");
- break;
- default: break;
- }
-
- fprintf(f, "* Redundancy cheking solver: ");
- switch (choiceRedcheckAlgorithm) {
- case DualSimplex:
- fprintf(f, "DualSimplex\n");
- break;
- case CrissCross:
- fprintf(f, "Criss-Cross\n");
- break;
- default: break;
- }
-
- fprintf(f, "* Lexicographic pivot: ");
- if (choiceLexicoPivotQ) fprintf(f, " on\n");
- else fprintf(f, " off\n");
-
- }
-
-
- void WriteRunningMode(FILE *f, Polyhedra poly)
- {
- if (poly.child!=null){
- fprintf(f,"* roworder: ");
- switch (poly.child.HalfspaceOrder) {
-
- case MinIndex:
- fprintf(f, "minindex\n");
- break;
-
- case MaxIndex:
- fprintf(f, "maxindex\n");
- break;
-
- case MinCutoff:
- fprintf(f, "mincutoff\n");
- break;
-
- case MaxCutoff:
- fprintf(f, "maxcutoff\n");
- break;
-
- case MixCutoff:
- fprintf(f, "mixcutoff\n");
- break;
-
- case LexMin:
- fprintf(f, "lexmin\n");
- break;
-
- case LexMax:
- fprintf(f, "lexmax\n");
- break;
-
- case RandomRow:
- fprintf(f, "random %d\n",poly.child.rseed);
- break;
-
- default: break;
- }
- }
- }
-
-
- void WriteCompletionStatus(FILE *f, Cone cone)
- {
- if (cone.Iteration<cone.m && cone.CompStatus==AllFound) {
- fprintf(f,"*Computation completed at Iteration %4ld.\n", cone.Iteration);
- }
- if (cone.CompStatus == RegionEmpty) {
- fprintf(f,"*Computation completed at Iteration %4ld because the region found empty.\n",cone.Iteration);
- }
- }
-
- void WritePolyFile(FILE *f, Polyhedra poly)
- {
- WriteAmatrix(f,poly.A,poly.m,poly.d);
- }
-
-
- void WriteErrorMessages(FILE *f, ErrorType Error)
- {
- switch (Error) {
-
- case DimensionTooLarge:
- fprintf(f, "*Input Error: Input matrix is too large:\n");
- fprintf(f, "*Please increase MMAX and/or NMAX in the source code and recompile.\n");
- break;
-
- case IFileNotFound:
- fprintf(f, "*Input Error: Specified input file does not exist.\n");
- break;
-
- case OFileNotOpen:
- fprintf(f, "*Output Error: Specified output file cannot be opened.\n");
- break;
-
- case NegativeMatrixSize:
- fprintf(f, "*Input Error: Input matrix has a negative size:\n");
- fprintf(f, "*Please check rowsize or colsize.\n");
- break;
-
- case ImproperInputFormat:
- fprintf(f,"*Input Error: Input format is not correct.\n");
- fprintf(f,"*Format:\n");
- fprintf(f," begin\n");
- fprintf(f," m n NumberType(real, rational or integer)\n");
- fprintf(f," b -A\n");
- fprintf(f," end\n");
- break;
-
- case EmptyVrepresentation:
- fprintf(f, "*Input Error: V-representation is empty:\n");
- fprintf(f, "*cddlib does not accept this trivial case for which output can be any inconsistent system.\n");
- break;
-
- case EmptyHrepresentation:
- fprintf(f, "*Input Error: H-representation is empty.\n");
- break;
-
- case EmptyRepresentation:
- fprintf(f, "*Input Error: Representation is empty.\n");
- break;
-
- case NoLPObjective:
- fprintf(f, "*LP Error: No LP objective (max or min) is set.\n");
- break;
-
- case NoRealNumberSupport:
- fprintf(f, "*LP Error: The binary (with GMP Rational) does not support Real number input.\n");
- fprintf(f, " : Use a binary compiled without -DGMPRATIONAL option.\n");
- break;
-
- case NotAvailForH:
- fprintf(f, "*Error: A function is called with H-rep which does not support an H-representation.\n");
- break;
-
- case NotAvailForV:
- fprintf(f, "*Error: A function is called with V-rep which does not support an V-representation.\n");
- break;
-
- case CannotHandleLinearity:
- fprintf(f, "*Error: The function called cannot handle linearity.\n");
- break;
-
- case RowIndexOutOfRange:
- fprintf(f, "*Error: Specified row index is out of range\n");
- break;
-
- case ColIndexOutOfRange:
- fprintf(f, "*Error: Specified column index is out of range\n");
- break;
-
- case LPCycling:
- fprintf(f, "*Error: Possibly an LP cycling occurs. Use the Criss-Cross method.\n");
- break;
-
- case NumericallyInconsistent:
- fprintf(f, "*Error: Numerical inconsistency is found. Use the GMP exact arithmetic.\n");
- break;
-
- case NoError:
- fprintf(f,"*No Error found.\n");
- break;
- }
- }
- #endif
-
- SetFamily CopyIncidence(Polyhedra poly)
- {
- SetFamily F = null;
- int k;
- int i;
-
- if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;
- if (poly.AincGenerated == false) ComputeAinc(poly);
- F = CreateSetFamily(poly.n, poly.m1);
- for (i = 1; i <= poly.m1; i++)
- for (k = 1; k <= poly.n; k++)
- if (poly.Ainc[i - 1].Get(k)) F.set[k - 1].Set(i, true);
-
- return F;
- }
-
- SetFamily CopyInputIncidence(Polyhedra poly)
- {
- int i;
- SetFamily F = null;
-
- if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;
- if (poly.AincGenerated == false) ComputeAinc(poly);
- F = CreateSetFamily(poly.m1, poly.n);
- for (i = 0; i < poly.m1; i++)
- {
- poly.Ainc[i].CopyTo(F.set[i]);
- }
-
- return F;
- }
-
- SetFamily CopyAdjacency(Polyhedra poly)
- {
- Ray RayPtr1, RayPtr2;
- SetFamily F = null;
- int pos1, pos2;
- int lstart, k, n;
- BitArray linset, allset;
- bool adj;
-
- if (poly.n == 0 && poly.homogeneous && poly.representation == RepresentationType.Inequality)
- {
- n = 1; /* the origin (the unique vertex) should be output. */
- }
- else n = poly.n;
- linset = new BitArray(n);
- allset = new BitArray(n);
- if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;
- F = CreateSetFamily(n, n);
- if (n <= 0) return F;
- poly.child.LastRay.Next = null;
- for (RayPtr1 = poly.child.FirstRay, pos1 = 1; RayPtr1 != null;
- RayPtr1 = RayPtr1.Next, pos1++)
- {
- for (RayPtr2 = poly.child.FirstRay, pos2 = 1; RayPtr2 != null;
- RayPtr2 = RayPtr2.Next, pos2++)
- {
- if (RayPtr1 != RayPtr2)
- {
- CheckAdjacency(poly.child, RayPtr1, RayPtr2, out adj);
- if (adj)
- {
- F.set[pos1 - 1].Set(pos2, true);
- }
- }
- }
- }
- lstart = poly.n - poly.ldim + 1;
- allset = allset.Not(); /* allset is set to the ground set. */
- for (k = lstart; k <= poly.n; k++)
- {
- linset.Set(k, true); /* linearity set */
- allset.CopyTo(F.set[k - 1]); /* linearity generator is adjacent to all */
- }
- for (k = 1; k < lstart; k++)
- {
- F.set[k - 1] = F.set[k - 1].Or(linset); // union
- /* every generator is adjacent to all linearity generators */
- }
-
- return F;
- }
-
- SetFamily CopyInputAdjacency(Polyhedra poly)
- {
- int i, j;
- SetFamily F = null;
-
- if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;
- if (poly.AincGenerated == false) ComputeAinc(poly);
- F = CreateSetFamily(poly.m1, poly.m1);
- for (i = 1; i <= poly.m1; i++)
- {
- for (j = 1; j <= poly.m1; j++)
- {
- if (i != j && InputAdjacentQ(poly, i, j))
- {
- F.set[i - 1].Set(j, true);
- }
- }
- }
-
- return F;
- }
-
- Matrix CopyOutput(Polyhedra poly)
- {
- Ray Ray;
- Matrix M = null;
- int i = 0, total;
- int j, j1;
- double b;
- RepresentationType outputrep = RepresentationType.Inequality;
- bool outputorigin = false;
-
- b = 0.0;
- total = poly.child.LinearityDim + poly.child.FeasibleRayCount;
-
- if (poly.child.d <= 0 || poly.child.newcol[1] == 0) total = total - 1;
- if (poly.representation == RepresentationType.Inequality) outputrep = RepresentationType.Generator;
- if (total == 0 && poly.homogeneous && poly.representation == RepresentationType.Inequality)
- {
- total = 1;
- outputorigin = true;
- /* the origin (the unique vertex) should be output. */
- }
- if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return M;
-
- M = CreateMatrix(total, poly.d);
- Ray = poly.child.FirstRay;
- while (Ray != null)
- {
- if (Ray.feasible)
- {
- CopyRay(M.matrix[i], poly.d, Ray, outputrep, poly.child.newcol);
- i++; /* 086 */
- }
- Ray = Ray.Next;
- }
- for (j = 2; j <= poly.d; j++)
- {
- if (poly.child.newcol[j] == 0)
- {
- /* original column j is dependent on others and removed for the cone */
- b = poly.child.Bsave[0][j - 1];
- if (outputrep == RepresentationType.Generator && b > 0.0)
- {
- M.matrix[i][0] = 1.0; /* Normalize */
- for (j1 = 1; j1 < poly.d; j1++)
- M.matrix[i][j1] = poly.child.Bsave[j1][j - 1] / b;
- }
- else
- {
- for (j1 = 0; j1 < poly.d; j1++)
- M.matrix[i][j1] = poly.child.Bsave[j1][j - 1];
- }
- M.linset.Set(i + 1, true);
- i++;
- }
- }
- if (outputorigin)
- {
- /* output the origin for homogeneous H-polyhedron with no rays. */
- M.matrix[0][0] = 1.0;
- for (j = 1; j < poly.d; j++)
- {
- M.matrix[0][j] = 0.0;
- }
- }
- MatrixIntegerFilter(M);
- if (poly.representation == RepresentationType.Inequality)
- M.representation = RepresentationType.Generator;
- else
- M.representation = RepresentationType.Inequality;
-
-
- return M;
- }
-
- Matrix CopyInput(Polyhedra poly)
- {
- Matrix M = null;
- int i;
-
- M = CreateMatrix(poly.m, poly.d);
- CopyAmatrix(M.matrix, poly.A, poly.m, poly.d);
-
- for (i = 1; i <= poly.m; i++)
- if (poly.EqualityIndex[i] == 1) M.linset.Set(i, true);
- MatrixIntegerFilter(M);
- if (poly.representation == RepresentationType.Generator)
- M.representation = RepresentationType.Generator;
- else
- M.representation = RepresentationType.Inequality;
- return M;
- }
-
- public Matrix CopyGenerators(Polyhedra poly)
- {
- Matrix M = null;
-
- if (poly.representation == RepresentationType.Generator)
- {
- M = CopyInput(poly);
- }
- else
- {
- M = CopyOutput(poly);
- }
- return M;
- }
-
- Matrix CopyInequalities(Polyhedra poly)
- {
- Matrix M = null;
-
- if (poly.representation == RepresentationType.Inequality)
- {
- M = CopyInput(poly);
- }
- else
- {
- M = CopyOutput(poly);
- }
- return M;
- }
-
- #if nemkell
- /****************************************************************************************/
- /* rational number (a/b) read is taken from Vinci by Benno Bueeler and Andreas Enge */
- /****************************************************************************************/
- void sread_rational_value (char *s, double value)
- /* reads a rational value from the specified string "s" and assigns it to "value" */
-
- {
- char *numerator_s=null, *denominator_s=null, *position;
- int sign = 1;
- double numerator, denominator;
- //#if defined GMPRATIONAL
- mpz_t znum, zden;
- //#else
- double rvalue;
- //#endif
-
- /* determine the sign of the number */
- numerator_s = s;
- if (s [0] == '-')
- { sign = -1;
- numerator_s++;
- }
- else if (s [0] == '+')
- numerator_s++;
-
- /* look for a sign '/' and eventually split the number in numerator and denominator */
- position = strchr (numerator_s, '/');
- if (position != null)
- { *position = '\0'; /* terminates the numerator */
- denominator_s = position + 1;
- };
-
- /* determine the floating point values of numerator and denominator */
- numerator=atol (numerator_s);
-
- if (position != null)
- {
- denominator=atol (denominator_s);
- }
- else denominator = 1;
-
- /*
- fprintf(stderr,"\nrational_read: numerator %f\n",numerator);
- fprintf(stderr,"rational_read: denominator %f\n",denominator);
- fprintf(stderr,"rational_read: sign %d\n",sign);
- */
-
- //#if defined GMPRATIONAL
- mpz_init_set_str(znum,numerator_s,10);
- if (sign<0) mpz_neg(znum,znum);
- mpz_init(zden); mpz_set_ui(zden,1);
- if (denominator_s!=null) mpz_init_set_str(zden,denominator_s,10);
- mpq_set_num(value,znum); mpq_set_den(value,zden);
- mpq_canonicalize(value);
- mpz_clear(znum); mpz_clear(zden);
- /* num=(int)sign * numerator; */
- /* den=(unsigned int) denominator; */
- /* mpq_set_si(value, num, den); */
- //#elif defined GMPFLOAT
- rvalue=sign * numerator/ (signed int) denominator;
- mpf_set_d(value, rvalue);
- //#else
- rvalue=sign * numerator/ (signed int) denominator;
- dset_d(value, rvalue);
- //#endif
- if (debug) {
- fprintf(stderr,"rational_read: ");
- WriteNumber(stderr,value); fprintf(stderr,"\n");
- }
- }
-
-
- void fread_rational_value (FILE *f, double value)
- /* reads a rational value from the specified file "f" and assigns it to "value" */
-
- {
- char number_s [255];
- double rational_value;
-
- init(rational_value);
- fscanf(f, "%s ", number_s);
- sread_rational_value (number_s, rational_value);
- set(value,rational_value);
- clear(rational_value);
- }
-
- /****************************************************************************************/
-
-
- /* end of cddio.c */
- #endif
-
- }
- }