/* 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.*/#include "setoper.h"  /* set operation library header (Ver. June 1, 2000 or later) */#include "cdd.h"#include <stdio.h>#include <stdlib.h>#include <time.h>#include <math.h>#include <string.h>/* void dd_fread_rational_value (FILE *, mytype *); */void dd_SetLinearity(dd_MatrixPtr, char *);void dd_SetInputFile(FILE **f,dd_DataFileType inputfile,dd_ErrorType *Error){  int opened=0,stop,quit=0;  int i,dotpos=0,trial=0;  char ch;  char *tempname;      *Error=dd_NoError;  while (!opened && !quit) {    fprintf(stderr,"\n>> Input file: ");    scanf("%s",inputfile);    ch=getchar();    stop=dd_FALSE;    for (i=0; i<dd_filenamelen && !stop; i++){      ch=inputfile[i];      switch (ch) {        case '.':           dotpos=i+1;          break;        case ';':  case ' ':  case '\0':  case '\n':  case '\t':               stop=dd_TRUE;          tempname=(char*)calloc(dd_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=dd_NoError;    }    else{      fprintf(stderr,"The file %s not found\n",inputfile);      trial++;      if (trial>5) {        *Error=dd_IFileNotFound;        quit=1;      }    }  }}void dd_SetWriteFileName(dd_DataFileType inputfile, dd_DataFileType outfile, char cflag, dd_RepresentationType rep){  char *extension;  dd_DataFileType ifilehead="";  unsigned int i,dotpos;    switch (cflag) {    case 'o':      switch (rep) {        case dd_Generator:          extension=".ine"; break;     /* output file for ine data */        case dd_Inequality:          extension=".ext"; break;     /* output file for ext data */        default:           extension=".xxx";break;      }      break;    case 'a':         /* decide for output adjacence */      if (rep==dd_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==dd_Inequality)        extension=".ecd";       /* ext incidence file */      else        extension=".icd";       /* ine incidence file */      break;    case 'n':         /* decide for input incidence */      if (rep==dd_Inequality)        extension=".icd";       /* ine incidence file */      else        extension=".ecd";       /* ext incidence file */      break;    case 'j':        /* decide for input adjacency */      if (rep==dd_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);  */}dd_NumberType dd_GetNumberType(char *line){  dd_NumberType nt;  if (strncmp(line, "integer", 7)==0) {    nt = dd_Integer;  }  else if (strncmp(line, "rational", 8)==0) {    nt = dd_Rational;  }  else if (strncmp(line, "real", 4)==0) {    nt = dd_Real;  }  else {     nt=dd_Unknown;  }  return nt;}void dd_ProcessCommandLine(FILE *f, dd_MatrixPtr M, char *line){  char newline[dd_linelenmax];  dd_colrange j;  mytype value;  dd_init(value);  if (strncmp(line, "hull", 4)==0) {    M->representation = dd_Generator;  }  if (strncmp(line, "debug", 5)==0) {    dd_debug = dd_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,dd_linelenmax,f);        dd_SetLinearity(M,newline);  }  if (strncmp(line, "maximize", 8)==0 ||      strncmp(line, "minimize", 8)==0) {    if (strncmp(line, "maximize", 8)==0) M->objective=dd_LPmax;    else M->objective=dd_LPmin;    for (j = 1; j <= M->colsize; j++) {    if (M->numbtype==dd_Real) {#if !defined(GMPRATIONAL)        double rvalue;        fscanf(f, "%lf", &rvalue);        dd_set_d(value, rvalue);#endif      } else {        dd_fread_rational_value (f, value);      }      dd_set(M->rowvec[j - 1],value);      if (dd_debug) {fprintf(stderr,"cost(%5ld) =",j); dd_WriteNumber(stderr,value);}    }  /*of j*/  }  dd_clear(value);}dd_boolean dd_AppendMatrix2Poly(dd_PolyhedraPtr *poly, dd_MatrixPtr M){  dd_boolean success=dd_FALSE;  dd_MatrixPtr Mpoly,Mnew=NULL;  dd_ErrorType err;  if ((*poly)!=NULL && (*poly)->m >=0 && (*poly)->d>=0 &&      (*poly)->d==M->colsize && M->rowsize>0){    Mpoly=dd_CopyInput(*poly);    Mnew=dd_AppendMatrix(Mpoly, M);    dd_FreePolyhedra(*poly);    *poly=dd_DDMatrix2Poly(Mnew,&err);    dd_FreeMatrix(Mpoly);    dd_FreeMatrix(Mnew);    if (err==dd_NoError) success=dd_TRUE;  }  return success;}dd_MatrixPtr dd_MatrixCopy(dd_MatrixPtr M){  dd_MatrixPtr Mcopy=NULL;  dd_rowrange m;  dd_colrange d;  m= M->rowsize;  d= M->colsize;  if (m >=0 && d >=0){    Mcopy=dd_CreateMatrix(m, d);    dd_CopyAmatrix(Mcopy->matrix, M->matrix, m, d);    dd_CopyArow(Mcopy->rowvec, M->rowvec, d);    set_copy(Mcopy->linset,M->linset);    Mcopy->numbtype=M->numbtype;    Mcopy->representation=M->representation;    Mcopy->objective=M->objective;  }  return Mcopy;}dd_MatrixPtr dd_CopyMatrix(dd_MatrixPtr M){  return dd_MatrixCopy(M);}dd_MatrixPtr dd_MatrixNormalizedCopy(dd_MatrixPtr M){  dd_MatrixPtr Mcopy=NULL;  dd_rowrange m;  dd_colrange d;  m= M->rowsize;  d= M->colsize;  if (m >=0 && d >=0){    Mcopy=dd_CreateMatrix(m, d);    dd_CopyNormalizedAmatrix(Mcopy->matrix, M->matrix, m, d);    dd_CopyArow(Mcopy->rowvec, M->rowvec, d);    set_copy(Mcopy->linset,M->linset);    Mcopy->numbtype=M->numbtype;    Mcopy->representation=M->representation;    Mcopy->objective=M->objective;  }  return Mcopy;}dd_MatrixPtr dd_MatrixAppend(dd_MatrixPtr M1, dd_MatrixPtr M2){  dd_MatrixPtr M=NULL;  dd_rowrange i, m,m1,m2;  dd_colrange 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=dd_CreateMatrix(m, d);    dd_CopyAmatrix(M->matrix, M1->matrix, m1, d);    dd_CopyArow(M->rowvec, M1->rowvec, d);    for (i=0; i<m1; i++){      if (set_member(i+1,M1->linset)) set_addelem(M->linset,i+1);    }    for (i=0; i<m2; i++){       for (j=0; j<d; j++)         dd_set(M->matrix[m1+i][j],M2->matrix[i][j]);         /* append the second matrix */       if (set_member(i+1,M2->linset)) set_addelem(M->linset,m1+i+1);    }    M->numbtype=M1->numbtype;  }  return M;}dd_MatrixPtr dd_MatrixNormalizedSortedCopy(dd_MatrixPtr M,dd_rowindex *newpos)  /* 094 */{   /* Sort the rows of Amatrix 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). */  dd_MatrixPtr Mcopy=NULL,Mnorm=NULL;  dd_rowrange m,i;  dd_colrange d;  dd_rowindex roworder;  /* if (newpos!=NULL) free(newpos); */  m= M->rowsize;  d= M->colsize;  roworder=(long *)calloc(m+1,sizeof(long*));  *newpos=(long *)calloc(m+1,sizeof(long*));  if (m >=0 && d >=0){    Mnorm=dd_MatrixNormalizedCopy(M);    Mcopy=dd_CreateMatrix(m, d);    for(i=1; i<=m; i++) roworder[i]=i;        dd_RandomPermutation(roworder, m, 123);    dd_QuickSort(roworder,1,m,Mnorm->matrix,d);    dd_PermuteCopyAmatrix(Mcopy->matrix, Mnorm->matrix, m, d, roworder);    dd_CopyArow(Mcopy->rowvec, M->rowvec, d);    for(i=1; i<=m; i++) {      if (set_member(roworder[i],M->linset)) set_addelem(Mcopy->linset, i);      (*newpos)[roworder[i]]=i;    }    Mcopy->numbtype=M->numbtype;    Mcopy->representation=M->representation;    Mcopy->objective=M->objective;    dd_FreeMatrix(Mnorm);  }  free(roworder);  return Mcopy;}dd_MatrixPtr dd_MatrixUniqueCopy(dd_MatrixPtr M,dd_rowindex *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 dd_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.  */  dd_MatrixPtr Mcopy=NULL;  dd_rowrange m,i,uniqrows;  dd_rowset preferredrows;  dd_colrange d;  dd_rowindex roworder;  /* if (newpos!=NULL) free(newpos); */  m= M->rowsize;  d= M->colsize;  preferredrows=M->linset;  roworder=(long *)calloc(m+1,sizeof(long*));  if (m >=0 && d >=0){    for(i=1; i<=m; i++) roworder[i]=i;    dd_UniqueRows(roworder, 1, m, M->matrix, d,preferredrows, &uniqrows);    Mcopy=dd_CreateMatrix(uniqrows, d);    dd_PermutePartialCopyAmatrix(Mcopy->matrix, M->matrix, m, d, roworder,1,m);    dd_CopyArow(Mcopy->rowvec, M->rowvec, d);    for(i=1; i<=m; i++) {      if (roworder[i]>0 && set_member(i,M->linset)) set_addelem(Mcopy->linset, roworder[i]);    }    Mcopy->numbtype=M->numbtype;    Mcopy->representation=M->representation;    Mcopy->objective=M->objective;  }  *newpos=roworder;  return Mcopy;}dd_MatrixPtr dd_MatrixNormalizedSortedUniqueCopy(dd_MatrixPtr M,dd_rowindex *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 dd_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.  */  dd_MatrixPtr M1=NULL,M2=NULL;  dd_rowrange m,i;  dd_colrange d;  dd_rowindex newpos1=NULL,newpos1r=NULL,newpos2=NULL;  /* if (newpos!=NULL) free(newpos); */  m= M->rowsize;  d= M->colsize;  *newpos=(long *)calloc(m+1,sizeof(long*));    newpos1r=(long *)calloc(m+1,sizeof(long*));    if (m>=0 && d>=0){    M1=dd_MatrixNormalizedSortedCopy(M,&newpos1);    for (i=1; i<=m;i++) newpos1r[newpos1[i]]=i;  /* reverse of newpos1 */    M2=dd_MatrixUniqueCopy(M1,&newpos2);    set_emptyset(M2->linset);    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 (set_member(i,M->linset)) set_addelem(M2->linset, newpos2[newpos1[i]]);         (*newpos)[i]=newpos2[newpos1[i]];      } else {         (*newpos)[i]=-newpos1r[-newpos2[newpos1[i]]];      }    }  dd_FreeMatrix(M1);free(newpos1);free(newpos2);free(newpos1r);  }    return M2;}dd_MatrixPtr dd_MatrixSortedUniqueCopy(dd_MatrixPtr M,dd_rowindex *newpos)  /* 094 */{   /* Same as dd_MatrixNormalizedSortedUniqueCopy except that it returns a unnormalized origial data     with original ordering.  */  dd_MatrixPtr M1=NULL,M2=NULL;  dd_rowrange m,i,k,ii;  dd_colrange d;  dd_rowindex newpos1=NULL,newpos1r=NULL,newpos2=NULL;  /* if (newpos!=NULL) free(newpos); */  m= M->rowsize;  d= M->colsize;  *newpos=(long *)calloc(m+1,sizeof(long*));    newpos1r=(long *)calloc(m+1,sizeof(long*));    if (m>=0 && d>=0){    M1=dd_MatrixNormalizedSortedCopy(M,&newpos1);    for (i=1; i<=m;i++) newpos1r[newpos1[i]]=i;  /* reverse of newpos1 */    M2=dd_MatrixUniqueCopy(M1,&newpos2);    dd_FreeMatrix(M1);    set_emptyset(M2->linset);    for(i=1; i<=m; i++) {      if (newpos2[newpos1[i]]>0){         if (set_member(i,M->linset)) set_addelem(M2->linset, newpos2[newpos1[i]]);         (*newpos)[i]=newpos2[newpos1[i]];      } else {         (*newpos)[i]=-newpos1r[-newpos2[newpos1[i]]];      }    }    ii=0;    set_emptyset(M2->linset);    for (i = 1; i<=m ; i++) {      k=(*newpos)[i];      if (k>0) {        ii+=1;        (*newpos)[i]=ii;        dd_CopyArow(M2->matrix[ii-1],M->matrix[i-1],d);        if (set_member(i,M->linset)) set_addelem(M2->linset, ii);      }    }    free(newpos1);free(newpos2);free(newpos1r);  }    return M2;}dd_MatrixPtr dd_AppendMatrix(dd_MatrixPtr M1, dd_MatrixPtr M2){  return dd_MatrixAppend(M1,M2);}int dd_MatrixAppendTo(dd_MatrixPtr *M1, dd_MatrixPtr M2){  dd_MatrixPtr M=NULL;  dd_rowrange i, m,m1,m2;  dd_colrange j, d,d1,d2;  dd_boolean success=0;    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=dd_CreateMatrix(m, d);    dd_CopyAmatrix(M->matrix, (*M1)->matrix, m1, d);    dd_CopyArow(M->rowvec, (*M1)->rowvec, d);    for (i=0; i<m1; i++){      if (set_member(i+1,(*M1)->linset)) set_addelem(M->linset,i+1);    }    for (i=0; i<m2; i++){       for (j=0; j<d; j++)         dd_set(M->matrix[m1+i][j],M2->matrix[i][j]);         /* append the second matrix */       if (set_member(i+1,M2->linset)) set_addelem(M->linset,m1+i+1);    }    M->numbtype=(*M1)->numbtype;    dd_FreeMatrix(*M1);    *M1=M;    success=1;  }  return success;}int dd_MatrixRowRemove(dd_MatrixPtr *M, dd_rowrange r) /* 092 */{  dd_rowrange i,m;  dd_colrange d;  dd_boolean success=0;