Sandia Home Sandia Home
Main Page | Publications | Downloads | Configuration | Running the Code | Solver Parameters | FAQ | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members | Related Pages

cddio.c

Go to the documentation of this file.
00001 
00043 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00044 
00045 /* cddio.c:  Basic Input and Output Procedures for cddlib
00046    written by Komei Fukuda, fukuda@ifor.math.ethz.ch
00047    Version 0.94, Aug. 4, 2005
00048 */
00049 
00050 /* cddlib : C-library of the double description method for
00051    computing all vertices and extreme rays of the polyhedron 
00052    P= {x :  b - A x >= 0}.  
00053    Please read COPYING (GNU General Public Licence) and
00054    the manual cddlibman.tex for detail.
00055 */
00056 
00057 #include "setoper.h"  /* set operation library header (Ver. June 1, 2000 or later) */
00058 #include "cdd.h"
00059 #include <stdio.h>
00060 #include <stdlib.h>
00061 #include <time.h>
00062 #include <math.h>
00063 #include <string.h>
00064 
00065 /* void dd_fread_rational_value (FILE *, mytype *); */
00066 void dd_SetLinearity(dd_MatrixPtr, char *);
00067 
00068 void dd_SetInputFile(FILE **f,dd_DataFileType inputfile,dd_ErrorType *Error)
00069 {
00070   int opened=0,stop,quit=0;
00071   int i,dotpos=0,trial=0;
00072   char ch;
00073   char *tempname;
00074   
00075   
00076   *Error=dd_NoError;
00077   while (!opened && !quit) {
00078     fprintf(stderr,"\n>> Input file: ");
00079     scanf("%s",inputfile);
00080     ch=getchar();
00081     stop=dd_FALSE;
00082     for (i=0; i<dd_filenamelen && !stop; i++){
00083       ch=inputfile[i];
00084       switch (ch) {
00085         case '.': 
00086           dotpos=i+1;
00087           break;
00088         case ';':  case ' ':  case '\0':  case '\n':  case '\t':     
00089           stop=dd_TRUE;
00090           tempname=(char*)calloc(dd_filenamelen,sizeof(ch));
00091           strncpy(tempname,inputfile,i);
00092           strcpy(inputfile,tempname);
00093           free(tempname);
00094           break;
00095       }
00096     }
00097     if ( ( *f = fopen(inputfile,"r") )!= NULL) {
00098       fprintf(stderr,"input file %s is open\n",inputfile);
00099       opened=1;
00100       *Error=dd_NoError;
00101     }
00102     else{
00103       fprintf(stderr,"The file %s not found\n",inputfile);
00104       trial++;
00105       if (trial>5) {
00106         *Error=dd_IFileNotFound;
00107         quit=1;
00108       }
00109     }
00110   }
00111 }
00112 
00113 void dd_SetWriteFileName(dd_DataFileType inputfile, dd_DataFileType outfile, char cflag, dd_RepresentationType rep)
00114 {
00115   char *extension;
00116   dd_DataFileType ifilehead="";
00117   int i,dotpos;
00118   
00119   switch (cflag) {
00120     case 'o':
00121       switch (rep) {
00122         case dd_Generator:
00123           extension=".ine"; break;     /* output file for ine data */
00124         case dd_Inequality:
00125           extension=".ext"; break;     /* output file for ext data */
00126         default: 
00127           extension=".xxx";break;
00128       }
00129       break;
00130 
00131     case 'a':         /* decide for output adjacence */
00132       if (rep==dd_Inequality)
00133         extension=".ead";       /* adjacency file for ext data */
00134       else
00135         extension=".iad";       /* adjacency file for ine data */
00136       break;
00137     case 'i':         /* decide for output incidence */
00138       if (rep==dd_Inequality)
00139         extension=".ecd";       /* ext incidence file */
00140       else
00141         extension=".icd";       /* ine incidence file */
00142       break;
00143     case 'n':         /* decide for input incidence */
00144       if (rep==dd_Inequality)
00145         extension=".icd";       /* ine incidence file */
00146       else
00147         extension=".ecd";       /* ext incidence file */
00148       break;
00149     case 'j':        /* decide for input adjacency */
00150       if (rep==dd_Inequality)
00151         extension=".iad";       /* ine adjacency file */
00152       else
00153         extension=".ead";       /* ext adjacency file */
00154       break;
00155     case 'l':
00156       extension=".ddl";break;   /* log file */
00157     case 'd':
00158       extension=".dex";break;   /* decomposition output */
00159     case 'p':
00160       extension="sub.ine";break;  /* preprojection sub inequality file */
00161     case 'v':
00162       extension=".solved";break;  /* verify_input file */
00163     case 's':
00164       extension=".lps";break;   /* LP solution file */
00165     default:
00166       extension=".xxx";break;
00167   }
00168   dotpos=-1;
00169   for (i=0; i< strlen(inputfile); i++){
00170     if (inputfile[i]=='.') dotpos=i;
00171   }
00172   if (dotpos>1) strncpy(ifilehead, inputfile, dotpos);
00173   else strcpy(ifilehead,inputfile);
00174   if (strlen(inputfile)<=0) strcpy(ifilehead,"tempcdd");
00175   strcpy(outfile,ifilehead); 
00176   strcat(outfile,extension); 
00177   if (strcmp(inputfile, outfile)==0) {
00178     strcpy(outfile,inputfile); 
00179     strcat(outfile,extension); 
00180   }
00181 /*  fprintf(stderr,"outfile name = %s\n",outfile);  */
00182 }
00183 
00184 
00185 dd_NumberType dd_GetNumberType(char *line)
00186 {
00187   dd_NumberType nt;
00188 
00189   if (strncmp(line, "integer", 7)==0) {
00190     nt = dd_Integer;
00191   }
00192   else if (strncmp(line, "rational", 8)==0) {
00193     nt = dd_Rational;
00194   }
00195   else if (strncmp(line, "real", 4)==0) {
00196     nt = dd_Real;
00197   }
00198   else { 
00199     nt=dd_Unknown;
00200   }
00201   return nt;
00202 }
00203 
00204 void dd_ProcessCommandLine(FILE *f, dd_MatrixPtr M, char *line)
00205 {
00206   char newline[dd_linelenmax];
00207   dd_colrange j;
00208   mytype value;
00209 
00210   dd_init(value);
00211   if (strncmp(line, "hull", 4)==0) {
00212     M->representation = dd_Generator;
00213   }
00214   if (strncmp(line, "debug", 5)==0) {
00215     dd_debug = dd_TRUE;
00216 #ifdef GMPRATIONAL
00217     ddf_debug = ddf_TRUE;
00218 #endif
00219   }
00220   if (strncmp(line, "partial_enum", 12)==0 ||
00221        strncmp(line, "equality", 8)==0  ||
00222        strncmp(line, "linearity", 9)==0 ) {
00223     fgets(newline,dd_linelenmax,f);    
00224     dd_SetLinearity(M,newline);
00225   }
00226   if (strncmp(line, "maximize", 8)==0 ||
00227       strncmp(line, "minimize", 8)==0) {
00228     if (strncmp(line, "maximize", 8)==0) M->objective=dd_LPmax;
00229     else M->objective=dd_LPmin;
00230     for (j = 1; j <= M->colsize; j++) {
00231     if (M->numbtype==dd_Real) {
00232 #if !defined(GMPRATIONAL)
00233         double rvalue;
00234         fscanf(f, "%lf", &rvalue);
00235         dd_set_d(value, rvalue);
00236 #endif
00237       } else {
00238         dd_fread_rational_value (f, value);
00239       }
00240       dd_set(M->rowvec[j - 1],value);
00241       if (dd_debug) {fprintf(stderr,"cost(%5ld) =",j); dd_WriteNumber(stderr,value);}
00242     }  /*of j*/
00243   }
00244   dd_clear(value);
00245 }
00246 
00247 dd_boolean dd_AppendMatrix2Poly(dd_PolyhedraPtr *poly, dd_MatrixPtr M)
00248 {
00249   dd_boolean success=dd_FALSE;
00250   dd_MatrixPtr Mpoly,Mnew=NULL;
00251   dd_ErrorType err;
00252 
00253   if ((*poly)!=NULL && (*poly)->m >=0 && (*poly)->d>=0 &&
00254       (*poly)->d==M->colsize && M->rowsize>0){
00255     Mpoly=dd_CopyInput(*poly);
00256     Mnew=dd_AppendMatrix(Mpoly, M);
00257     dd_FreePolyhedra(*poly);
00258     *poly=dd_DDMatrix2Poly(Mnew,&err);
00259     dd_FreeMatrix(Mpoly);
00260     dd_FreeMatrix(Mnew);
00261     if (err==dd_NoError) success=dd_TRUE;
00262   }
00263   return success;
00264 }
00265 
00266 dd_MatrixPtr dd_MatrixCopy(dd_MatrixPtr M)
00267 {
00268   dd_MatrixPtr Mcopy=NULL;
00269   dd_rowrange m;
00270   dd_colrange d;
00271 
00272   m= M->rowsize;
00273   d= M->colsize;
00274   if (m >=0 && d >=0){
00275     Mcopy=dd_CreateMatrix(m, d);
00276     dd_CopyAmatrix(Mcopy->matrix, M->matrix, m, d);
00277     dd_CopyArow(Mcopy->rowvec, M->rowvec, d);
00278     set_copy(Mcopy->linset,M->linset);
00279     Mcopy->numbtype=M->numbtype;
00280     Mcopy->representation=M->representation;
00281     Mcopy->objective=M->objective;
00282   }
00283   return Mcopy;
00284 }
00285 
00286 dd_MatrixPtr dd_CopyMatrix(dd_MatrixPtr M)
00287 {
00288   return dd_MatrixCopy(M);
00289 }
00290 
00291 dd_MatrixPtr dd_MatrixNormalizedCopy(dd_MatrixPtr M)
00292 {
00293   dd_MatrixPtr Mcopy=NULL;
00294   dd_rowrange m;
00295   dd_colrange d;
00296 
00297   m= M->rowsize;
00298   d= M->colsize;
00299   if (m >=0 && d >=0){
00300     Mcopy=dd_CreateMatrix(m, d);
00301     dd_CopyNormalizedAmatrix(Mcopy->matrix, M->matrix, m, d);
00302     dd_CopyArow(Mcopy->rowvec, M->rowvec, d);
00303     set_copy(Mcopy->linset,M->linset);
00304     Mcopy->numbtype=M->numbtype;
00305     Mcopy->representation=M->representation;
00306     Mcopy->objective=M->objective;
00307   }
00308   return Mcopy;
00309 }
00310 
00311 
00312 dd_MatrixPtr dd_MatrixAppend(dd_MatrixPtr M1, dd_MatrixPtr M2)
00313 {
00314   dd_MatrixPtr M=NULL;
00315   dd_rowrange i, m,m1,m2;
00316   dd_colrange j, d,d1,d2;
00317   
00318   m1=M1->rowsize;
00319   d1=M1->colsize;
00320   m2=M2->rowsize;
00321   d2=M2->colsize;
00322 
00323   m=m1+m2;
00324   d=d1;
00325 
00326   if (d1>=0 && d1==d2 && m1>=0 && m2>=0){
00327     M=dd_CreateMatrix(m, d);
00328     dd_CopyAmatrix(M->matrix, M1->matrix, m1, d);
00329     dd_CopyArow(M->rowvec, M1->rowvec, d);
00330     for (i=0; i<m1; i++){
00331       if (set_member(i+1,M1->linset)) set_addelem(M->linset,i+1);
00332     }
00333     for (i=0; i<m2; i++){
00334        for (j=0; j<d; j++)
00335          dd_set(M->matrix[m1+i][j],M2->matrix[i][j]);
00336          /* append the second matrix */
00337        if (set_member(i+1,M2->linset)) set_addelem(M->linset,m1+i+1);
00338     }
00339     M->numbtype=M1->numbtype;
00340   }
00341   return M;
00342 }
00343 
00344 dd_MatrixPtr dd_MatrixNormalizedSortedCopy(dd_MatrixPtr M,dd_rowindex *newpos)  /* 094 */
00345 { 
00346   /* Sort the rows of Amatrix lexicographically, and return a link to this sorted copy. 
00347   The vector newpos is allocated, where newpos[i] returns the new row index
00348   of the original row i (i=1,...,M->rowsize). */
00349   dd_MatrixPtr Mcopy=NULL,Mnorm=NULL;
00350   dd_rowrange m,i;
00351   dd_colrange d;
00352   dd_rowindex roworder;
00353 
00354   /* if (newpos!=NULL) free(newpos); */
00355   m= M->rowsize;
00356   d= M->colsize;
00357   roworder=(long *)calloc(m+1,sizeof(long*));
00358   *newpos=(long *)calloc(m+1,sizeof(long*));
00359   if (m >=0 && d >=0){
00360     Mnorm=dd_MatrixNormalizedCopy(M);
00361     Mcopy=dd_CreateMatrix(m, d);
00362     for(i=1; i<=m; i++) roworder[i]=i;
00363     
00364     dd_RandomPermutation(roworder, m, 123);
00365     dd_QuickSort(roworder,1,m,Mnorm->matrix,d);
00366 
00367     dd_PermuteCopyAmatrix(Mcopy->matrix, Mnorm->matrix, m, d, roworder);
00368     dd_CopyArow(Mcopy->rowvec, M->rowvec, d);
00369     for(i=1; i<=m; i++) {
00370       if (set_member(roworder[i],M->linset)) set_addelem(Mcopy->linset, i);
00371       (*newpos)[roworder[i]]=i;
00372     }
00373     Mcopy->numbtype=M->numbtype;
00374     Mcopy->representation=M->representation;
00375     Mcopy->objective=M->objective;
00376     dd_FreeMatrix(Mnorm);
00377   }
00378   free(roworder);
00379   return Mcopy;
00380 }
00381 
00382 dd_MatrixPtr dd_MatrixUniqueCopy(dd_MatrixPtr M,dd_rowindex *newpos)
00383 { 
00384   /* Remove row duplicates, and return a link to this sorted copy. 
00385      Linearity rows have priority over the other rows.
00386      It is better to call this after sorting with dd_MatrixNormalizedSortedCopy.
00387      The vector newpos is allocated, where *newpos[i] returns the new row index
00388      of the original row i (i=1,...,M->rowsize).  *newpos[i] is negative if the original
00389      row is dominated by -*newpos[i] and eliminated in the new copy.
00390   */
00391   dd_MatrixPtr Mcopy=NULL;
00392   dd_rowrange m,i,uniqrows;
00393   dd_rowset preferredrows;
00394   dd_colrange d;
00395   dd_rowindex roworder;
00396 
00397   /* if (newpos!=NULL) free(newpos); */
00398   m= M->rowsize;
00399   d= M->colsize;
00400   preferredrows=M->linset;
00401   roworder=(long *)calloc(m+1,sizeof(long*));
00402   if (m >=0 && d >=0){
00403     for(i=1; i<=m; i++) roworder[i]=i;
00404     dd_UniqueRows(roworder, 1, m, M->matrix, d,preferredrows, &uniqrows);
00405 
00406     Mcopy=dd_CreateMatrix(uniqrows, d);
00407     dd_PermutePartialCopyAmatrix(Mcopy->matrix, M->matrix, m, d, roworder,1,m);
00408     dd_CopyArow(Mcopy->rowvec, M->rowvec, d);
00409     for(i=1; i<=m; i++) {
00410       if (roworder[i]>0 && set_member(i,M->linset)) set_addelem(Mcopy->linset, roworder[i]);
00411     }
00412     Mcopy->numbtype=M->numbtype;
00413     Mcopy->representation=M->representation;
00414     Mcopy->objective=M->objective;
00415   }
00416   *newpos=roworder;
00417   return Mcopy;
00418 }
00419 
00420 
00421 dd_MatrixPtr dd_MatrixNormalizedSortedUniqueCopy(dd_MatrixPtr M,dd_rowindex *newpos) /* 094 */
00422 { 
00423   /* Sort and remove row duplicates, and return a link to this sorted copy. 
00424      Linearity rows have priority over the other rows.
00425      It is better to call this after sorting with dd_MatrixNormalizedSortedCopy.
00426      The vector newpos is allocated, where *newpos[i] returns the new row index
00427      of the original row i (i=1,...,M->rowsize).  *newpos[i] is negative if the original
00428      row is dominated by -*newpos[i] and eliminated in the new copy.
00429   */
00430   dd_MatrixPtr M1=NULL,M2=NULL;
00431   dd_rowrange m,i;
00432   dd_colrange d;
00433   dd_rowindex newpos1=NULL,newpos1r=NULL,newpos2=NULL;
00434 
00435   /* if (newpos!=NULL) free(newpos); */
00436   m= M->rowsize;
00437   d= M->colsize;
00438   *newpos=(long *)calloc(m+1,sizeof(long*));  
00439   newpos1r=(long *)calloc(m+1,sizeof(long*));  
00440   if (m>=0 && d>=0){
00441     M1=dd_MatrixNormalizedSortedCopy(M,&newpos1);
00442     for (i=1; i<=m;i++) newpos1r[newpos1[i]]=i;  /* reverse of newpos1 */
00443     M2=dd_MatrixUniqueCopy(M1,&newpos2);
00444     set_emptyset(M2->linset);
00445     for(i=1; i<=m; i++) {
00446       if (newpos2[newpos1[i]]>0){
00447          printf("newpos1[%ld]=%ld, newpos2[newpos1[%ld]]=%ld\n",i,newpos1[i], i,newpos2[newpos1[i]]);
00448          if (set_member(i,M->linset)) set_addelem(M2->linset, newpos2[newpos1[i]]);
00449          (*newpos)[i]=newpos2[newpos1[i]];
00450       } else {
00451          (*newpos)[i]=-newpos1r[-newpos2[newpos1[i]]];
00452       }
00453     }
00454   dd_FreeMatrix(M1);free(newpos1);free(newpos2);free(newpos1r);
00455   }
00456   
00457   return M2;
00458 }
00459 
00460 dd_MatrixPtr dd_MatrixSortedUniqueCopy(dd_MatrixPtr M,dd_rowindex *newpos)  /* 094 */
00461 { 
00462   /* Same as dd_MatrixNormalizedSortedUniqueCopy except that it returns a unnormalized origial data
00463      with original ordering.
00464   */
00465   dd_MatrixPtr M1=NULL,M2=NULL;
00466   dd_rowrange m,i,k,ii;
00467   dd_colrange d;
00468   dd_rowindex newpos1=NULL,newpos1r=NULL,newpos2=NULL;
00469 
00470   /* if (newpos!=NULL) free(newpos); */
00471   m= M->rowsize;
00472   d= M->colsize;
00473   *newpos=(long *)calloc(m+1,sizeof(long*));  
00474   newpos1r=(long *)calloc(m+1,sizeof(long*));  
00475   if (m>=0 && d>=0){
00476     M1=dd_MatrixNormalizedSortedCopy(M,&newpos1);
00477     for (i=1; i<=m;i++) newpos1r[newpos1[i]]=i;  /* reverse of newpos1 */
00478     M2=dd_MatrixUniqueCopy(M1,&newpos2);
00479     dd_FreeMatrix(M1);
00480     set_emptyset(M2->linset);
00481     for(i=1; i<=m; i++) {
00482       if (newpos2[newpos1[i]]>0){
00483          if (set_member(i,M->linset)) set_addelem(M2->linset, newpos2[newpos1[i]]);
00484          (*newpos)[i]=newpos2[newpos1[i]];
00485       } else {
00486          (*newpos)[i]=-newpos1r[-newpos2[newpos1[i]]];
00487       }
00488     }
00489 
00490     ii=0;
00491     set_emptyset(M2->linset);
00492     for (i = 1; i<=m ; i++) {
00493       k=(*newpos)[i];
00494       if (k>0) {
00495         ii+=1;
00496         (*newpos)[i]=ii;
00497         dd_CopyArow(M2->matrix[ii-1],M->matrix[i-1],d);
00498         if (set_member(i,M->linset)) set_addelem(M2->linset, ii);
00499       }
00500     }
00501 
00502     free(newpos1);free(newpos2);free(newpos1r);
00503   }
00504   
00505   return M2;
00506 }
00507 
00508 dd_MatrixPtr dd_AppendMatrix(dd_MatrixPtr M1, dd_MatrixPtr M2)
00509 {
00510   return dd_MatrixAppend(M1,M2);
00511 }
00512 
00513 int dd_MatrixAppendTo(dd_MatrixPtr *M1, dd_MatrixPtr M2)
00514 {
00515   dd_MatrixPtr M=NULL;
00516   dd_rowrange i, m,m1,m2;
00517   dd_colrange j, d,d1,d2;
00518   dd_boolean success=0;
00519   
00520   m1=(*M1)->rowsize;
00521   d1=(*M1)->colsize;
00522   m2=M2->rowsize;
00523   d2=M2->colsize;
00524 
00525   m=m1+m2;
00526   d=d1;
00527 
00528   if (d1>=0 && d1==d2 && m1>=0 && m2>=0){
00529     M=dd_CreateMatrix(m, d);
00530     dd_CopyAmatrix(M->matrix, (*M1)->matrix, m1, d);
00531     dd_CopyArow(M->rowvec, (*M1)->rowvec, d);
00532     for (i=0; i<m1; i++){
00533       if (set_member(i+1,(*M1)->linset)) set_addelem(M->linset,i+1);
00534     }
00535     for (i=0; i<m2; i++){
00536        for (j=0; j<d; j++)
00537          dd_set(M->matrix[m1+i][j],M2->matrix[i][j]);
00538          /* append the second matrix */
00539        if (set_member(i+1,M2->linset)) set_addelem(M->linset,m1+i+1);
00540     }
00541     M->numbtype=(*M1)->numbtype;
00542     dd_FreeMatrix(*M1);
00543     *M1=M;
00544     success=1;
00545   }
00546   return success;
00547 }
00548 
00549 int dd_MatrixRowRemove(dd_MatrixPtr *M, dd_rowrange r) /* 092 */
00550 {
00551   dd_rowrange i,m;
00552   dd_colrange d;
00553   dd_boolean success=0;
00554   
00555   m=(*M)->rowsize;
00556   d=(*M)->colsize;
00557 
00558   if (r >= 1 && r <=m){
00559     (*M)->rowsize=m-1;
00560     dd_FreeArow(d, (*M)->matrix[r-1]);
00561     set_delelem((*M)->linset,r);
00562     /* slide the row headers */
00563     for (i=r; i<m; i++){
00564       (*M)->matrix[i-1]=(*M)->matrix[i];
00565       if (set_member(i+1, (*M)->linset)){
00566         set_delelem((*M)->linset,i+1);
00567         set_addelem((*M)->linset,i);
00568       }
00569     }
00570     success=1;
00571   }
00572   return success;
00573 }
00574 
00575 int dd_MatrixRowRemove2(dd_MatrixPtr *M, dd_rowrange r, dd_rowindex *newpos) /* 094 */
00576 {
00577   dd_rowrange i,m;
00578   dd_colrange d;
00579   dd_boolean success=0;
00580   dd_rowindex roworder;
00581   
00582   m=(*M)->rowsize;
00583   d=(*M)->colsize;
00584 
00585   if (r >= 1 && r <=m){
00586     roworder=(long *)calloc(m+1,sizeof(long*));
00587     (*M)->rowsize=m-1;
00588     dd_FreeArow(d, (*M)->matrix[r-1]);
00589     set_delelem((*M)->linset,r);
00590     /* slide the row headers */
00591     for (i=1; i<r; i++) roworder[i]=i;
00592     roworder[r]=0; /* meaning it is removed */
00593     for (i=r; i<m; i++){
00594       (*M)->matrix[i-1]=(*M)->matrix[i];
00595       roworder[i+1]=i;
00596       if (set_member(i+1, (*M)->linset)){
00597         set_delelem((*M)->linset,i+1);
00598         set_addelem((*M)->linset,i);
00599       }
00600     }
00601     success=1;
00602   }
00603   return success;
00604 }
00605 
00606 dd_MatrixPtr dd_MatrixSubmatrix(dd_MatrixPtr M, dd_rowset delset) /* 092 */
00607 {
00608   dd_MatrixPtr Msub=NULL;
00609   dd_rowrange i,isub=1, m,msub;
00610   dd_colrange d;
00611  
00612   m= M->rowsize;
00613   d= M->colsize;
00614   msub=m;
00615   if (m >=0 && d >=0){
00616     for (i=1; i<=m; i++) {
00617        if (set_member(i,delset)) msub-=1;
00618     }
00619     Msub=dd_CreateMatrix(msub, d);
00620     for (i=1; i<=m; i++){
00621       if (!set_member(i,delset)){
00622         dd_CopyArow(Msub->matrix[isub-1], M->matrix[i-1], d);
00623         if (set_member(i, M->linset)){
00624           set_addelem(Msub->linset,isub);
00625         }
00626         isub++;
00627       }
00628     }
00629     dd_CopyArow(Msub->rowvec, M->rowvec, d);
00630     Msub->numbtype=M->numbtype;
00631     Msub->representation=M->representation;
00632     Msub->objective=M->objective;
00633   }
00634   return Msub;
00635 }
00636 
00637 dd_MatrixPtr dd_MatrixSubmatrix2(dd_MatrixPtr M, dd_rowset delset,dd_rowindex *newpos) /* 092 */
00638 { /* returns a pointer to a new matrix which is a submatrix of M with rows in delset
00639   removed.  *newpos[i] returns the position of the original row i in the new matrix.
00640   It is -1 if and only if it is deleted.
00641   */
00642   
00643   dd_MatrixPtr Msub=NULL;
00644   dd_rowrange i,isub=1, m,msub;
00645   dd_colrange d;
00646   dd_rowindex roworder;
00647 
00648   m= M->rowsize;
00649   d= M->colsize;
00650   msub=m;
00651   if (m >=0 && d >=0){
00652     roworder=(long *)calloc(m+1,sizeof(long*));
00653     for (i=1; i<=m; i++) {
00654        if (set_member(i,delset)) msub-=1;
00655     }
00656     Msub=dd_CreateMatrix(msub, d);
00657     for (i=1; i<=m; i++){
00658       if (set_member(i,delset)){
00659         roworder[i]=0; /* zero means the row i is removed */
00660       } else {
00661         dd_CopyArow(Msub->matrix[isub-1], M->matrix[i-1], d);
00662         if (set_member(i, M->linset)){
00663           set_addelem(Msub->linset,isub);
00664         }
00665         roworder[i]=isub;
00666         isub++;
00667       }
00668     }
00669     *newpos=roworder;
00670     dd_CopyArow(Msub->rowvec, M->rowvec, d);
00671     Msub->numbtype=M->numbtype;
00672     Msub->representation=M->representation;
00673     Msub->objective=M->objective;
00674   }
00675   return Msub;
00676 }
00677 
00678 
00679 dd_MatrixPtr dd_MatrixSubmatrix2L(dd_MatrixPtr M, dd_rowset delset,dd_rowindex *newpos) /* 094 */
00680 { /* This is same as dd_MatrixSubmatrix2 except that the linearity rows will be shifted up
00681      so that they are at the top of the matrix.
00682   */
00683   dd_MatrixPtr Msub=NULL;
00684   dd_rowrange i,iL, iI, m,msub;
00685   dd_colrange d;
00686   dd_rowindex roworder;
00687 
00688   m= M->rowsize;
00689   d= M->colsize;
00690   msub=m;
00691   if (m >=0 && d >=0){
00692     roworder=(long *)calloc(m+1,sizeof(long*));
00693     for (i=1; i<=m; i++) {
00694        if (set_member(i,delset)) msub-=1;
00695     }
00696     Msub=dd_CreateMatrix(msub, d);
00697     iL=1; iI=set_card(M->linset)+1;  /* starting positions */
00698     for (i=1; i<=m; i++){
00699       if (set_member(i,delset)){
00700         roworder[i]=0; /* zero means the row i is removed */
00701       } else {
00702         if (set_member(i,M->linset)){
00703           dd_CopyArow(Msub->matrix[iL-1], M->matrix[i-1], d);
00704           set_delelem(Msub->linset,i);
00705           set_addelem(Msub->linset,iL);
00706           roworder[i]=iL;
00707           iL+=1;
00708         } else {
00709           dd_CopyArow(Msub->matrix[iI-1], M->matrix[i-1], d);
00710           roworder[i]=iI;
00711           iI+=1;
00712         }
00713        }
00714     }
00715     *newpos=roworder;
00716     dd_CopyArow(Msub->rowvec, M->rowvec, d);
00717     Msub->numbtype=M->numbtype;
00718     Msub->representation=M->representation;
00719     Msub->objective=M->objective;
00720   }
00721   return Msub;
00722 }
00723 
00724 int dd_MatrixRowsRemove(dd_MatrixPtr *M, dd_rowset delset) /* 094 */
00725 {
00726   dd_MatrixPtr Msub=NULL;
00727   int success;
00728   
00729   Msub=dd_MatrixSubmatrix(*M, delset);
00730   dd_FreeMatrix(*M);
00731   *M=Msub;
00732   success=1;
00733   return success;
00734 }
00735 
00736 int dd_MatrixRowsRemove2(dd_MatrixPtr *M, dd_rowset delset,dd_rowindex *newpos) /* 094 */
00737 {
00738   dd_MatrixPtr Msub=NULL;
00739   int success;
00740   
00741   Msub=dd_MatrixSubmatrix2(*M, delset,newpos);
00742   dd_FreeMatrix(*M);
00743   *M=Msub;
00744   success=1;
00745   return success;
00746 }
00747 
00748 int dd_MatrixShiftupLinearity(dd_MatrixPtr *M,dd_rowindex *newpos) /* 094 */
00749 {
00750   dd_MatrixPtr Msub=NULL;
00751   int success;
00752   dd_rowset delset;
00753   
00754   set_initialize(&delset,(*M)->rowsize);  /* emptyset */
00755   Msub=dd_MatrixSubmatrix2L(*M, delset,newpos);
00756   dd_FreeMatrix(*M);
00757   *M=Msub;
00758   
00759   set_free(delset);
00760   success=1;
00761   return success;
00762 }
00763 
00764 dd_PolyhedraPtr dd_CreatePolyhedraData(dd_rowrange m, dd_colrange d)
00765 {
00766   dd_rowrange i;
00767   dd_PolyhedraPtr poly=NULL;
00768 
00769   poly=(dd_PolyhedraPtr) malloc (sizeof(dd_PolyhedraType));
00770   poly->child       =NULL; /* this links the homogenized cone data */
00771   poly->m           =m;
00772   poly->d           =d;  
00773   poly->n           =-1;  /* the size of output is not known */
00774   poly->m_alloc     =m+2; /* the allocated row size of matrix A */
00775   poly->d_alloc     =d;   /* the allocated col size of matrix A */
00776   poly->numbtype=dd_Real;
00777   dd_InitializeAmatrix(poly->m_alloc,poly->d_alloc,&(poly->A));
00778   dd_InitializeArow(d,&(poly->c));           /* cost vector */
00779   poly->representation       =dd_Inequality;
00780   poly->homogeneous =dd_FALSE;
00781 
00782   poly->EqualityIndex=(int *)calloc(m+2, sizeof(int));  
00783     /* size increased to m+2 in 092b because it is used by the child cone, 
00784        This is a bug fix suggested by Thao Dang. */
00785     /* ith component is 1 if it is equality, -1 if it is strict inequality, 0 otherwise. */
00786   for (i = 0; i <= m+1; i++) poly->EqualityIndex[i]=0;
00787 
00788   poly->IsEmpty                 = -1;  /* initially set to -1, neither TRUE nor FALSE, meaning unknown */
00789   
00790   poly->NondegAssumed           = dd_FALSE;
00791   poly->InitBasisAtBottom       = dd_FALSE;
00792   poly->RestrictedEnumeration   = dd_FALSE;
00793   poly->RelaxedEnumeration      = dd_FALSE;
00794 
00795   poly->AincGenerated=dd_FALSE;  /* Ainc is a set array to store the input incidence. */
00796 
00797   return poly;
00798 }
00799 
00800 dd_boolean dd_InitializeConeData(dd_rowrange m, dd_colrange d, dd_ConePtr *cone)
00801 {
00802   dd_boolean success=dd_TRUE;
00803   dd_colrange j;
00804 
00805   (*cone)=(dd_ConePtr)calloc(1, sizeof(dd_ConeType));
00806 
00807 /* INPUT: A given representation of a cone: inequality */
00808   (*cone)->m=m;
00809   (*cone)->d=d;
00810   (*cone)->m_alloc=m+2; /* allocated row size of matrix A */
00811   (*cone)->d_alloc=d;   /* allocated col size of matrix A, B and Bsave */
00812   (*cone)->numbtype=dd_Real;
00813   (*cone)->parent=NULL;
00814 
00815 /* CONTROL: variables to control computation */
00816   (*cone)->Iteration=0;
00817 
00818   (*cone)->HalfspaceOrder=dd_LexMin;
00819 
00820   (*cone)->ArtificialRay=NULL;
00821   (*cone)->FirstRay=NULL;
00822   (*cone)->LastRay=NULL; /* The second description: Generator */
00823   (*cone)->PosHead=NULL;
00824   (*cone)->ZeroHead=NULL;
00825   (*cone)->NegHead=NULL;
00826   (*cone)->PosLast=NULL;
00827   (*cone)->ZeroLast=NULL;
00828   (*cone)->NegLast=NULL;
00829   (*cone)->RecomputeRowOrder  = dd_TRUE;
00830   (*cone)->PreOrderedRun      = dd_FALSE;
00831   set_initialize(&((*cone)->GroundSet),(*cone)->m_alloc);
00832   set_initialize(&((*cone)->EqualitySet),(*cone)->m_alloc);
00833   set_initialize(&((*cone)->NonequalitySet),(*cone)->m_alloc);
00834   set_initialize(&((*cone)->AddedHalfspaces),(*cone)->m_alloc);
00835   set_initialize(&((*cone)->WeaklyAddedHalfspaces),(*cone)->m_alloc);
00836   set_initialize(&((*cone)->InitialHalfspaces),(*cone)->m_alloc);
00837   (*cone)->RayCount=0;
00838   (*cone)->FeasibleRayCount=0;
00839   (*cone)->WeaklyFeasibleRayCount=0;
00840   (*cone)->TotalRayCount=0;
00841   (*cone)->ZeroRayCount=0;
00842   (*cone)->EdgeCount=0;
00843   (*cone)->TotalEdgeCount=0;
00844   (*cone)->count_int=0;
00845   (*cone)->count_int_good=0;
00846   (*cone)->count_int_bad=0;
00847   (*cone)->rseed=1;  /* random seed for random row permutation */
00848  
00849   dd_InitializeBmatrix((*cone)->d_alloc, &((*cone)->B));
00850   dd_InitializeBmatrix((*cone)->d_alloc, &((*cone)->Bsave));
00851   dd_InitializeAmatrix((*cone)->m_alloc,(*cone)->d_alloc,&((*cone)->A));
00852 
00853   (*cone)->Edges
00854      =(dd_AdjacencyType**) calloc((*cone)->m_alloc,sizeof(dd_AdjacencyType*));
00855   (*cone)->InitialRayIndex=(long*)calloc(d+1,sizeof(long));
00856   (*cone)->OrderVector=(long*)calloc((*cone)->m_alloc+1,sizeof(long));
00857 
00858   (*cone)->newcol=(long*)calloc(((*cone)->d)+1,sizeof(long));
00859   for (j=0; j<=(*cone)->d; j++) (*cone)->newcol[j]=j;  /* identity map, initially */
00860   (*cone)->LinearityDim = -2; /* -2 if it is not computed */
00861   (*cone)->ColReduced   = dd_FALSE;
00862   (*cone)->d_orig = d;
00863 
00864 /* STATES: variables to represent current state. */
00865 /*(*cone)->Error;
00866   (*cone)->CompStatus;
00867   (*cone)->starttime;
00868   (*cone)->endtime;
00869 */
00870     
00871   return success;
00872 }
00873 
00874 dd_ConePtr dd_ConeDataLoad(dd_PolyhedraPtr poly)
00875 {
00876   dd_ConePtr cone=NULL;
00877   dd_colrange d,j;
00878   dd_rowrange m,i;
00879 
00880   m=poly->m;
00881   d=poly->d;
00882   if (!(poly->homogeneous) && poly->representation==dd_Inequality){
00883     m=poly->m+1;
00884   }
00885   poly->m1=m;
00886 
00887   dd_InitializeConeData(m, d, &cone);
00888   cone->representation=poly->representation;
00889 
00890 /* Points to the original polyhedra data, and reversely */
00891   cone->parent=poly;
00892   poly->child=cone;
00893 
00894   for (i=1; i<=poly->m; i++)
00895     for (j=1; j<=cone->d; j++)
00896       dd_set(cone->A[i-1][j-1],poly->A[i-1][j-1]);  
00897   
00898   if (poly->representation==dd_Inequality && !(poly->homogeneous)){
00899     dd_set(cone->A[m-1][0],dd_one);
00900     for (j=2; j<=d; j++) dd_set(cone->A[m-1][j-1],dd_purezero);
00901   }
00902 
00903   return cone;
00904 }
00905 
00906 void dd_SetLinearity(dd_MatrixPtr M, char *line)
00907 {
00908   int i=0;
00909   dd_rowrange eqsize,var;
00910   char *next;
00911   const char ct[]=", ";  /* allows separators "," and " ". */
00912 
00913   next=strtok(line,ct);
00914   eqsize=atol(next); 
00915   while (i < eqsize && (next=strtok(NULL,ct))!=NULL) {
00916      var=atol(next);
00917      set_addelem(M->linset,var); i++;
00918   }
00919   if (i!=eqsize) {
00920     fprintf(stderr,"* Warning: there are inconsistencies in linearity setting.\n");
00921   }
00922   return;
00923 }
00924 
00925 dd_MatrixPtr dd_PolyFile2Matrix (FILE *f, dd_ErrorType *Error)
00926 {
00927   dd_MatrixPtr M=NULL;
00928   dd_rowrange m_input,i;
00929   dd_colrange d_input,j;
00930   dd_RepresentationType rep=dd_Inequality;
00931   mytype value;
00932   dd_boolean found=dd_FALSE, newformat=dd_FALSE, successful=dd_FALSE, linearity=dd_FALSE;
00933   char command[dd_linelenmax], comsave[dd_linelenmax], numbtype[dd_wordlenmax];
00934   dd_NumberType NT;
00935 #if !defined(GMPRATIONAL)
00936   double rvalue;
00937 #endif
00938 
00939   dd_init(value);
00940   (*Error)=dd_NoError;
00941   while (!found)
00942   {
00943     if (fscanf(f,"%s",command)==EOF) {
00944       (*Error)=dd_ImproperInputFormat;
00945       goto _L99;
00946     }
00947     else {
00948       if (strncmp(command, "V-representation", 16)==0) {
00949         rep=dd_Generator; newformat=dd_TRUE;
00950       }
00951       if (strncmp(command, "H-representation", 16)==0){
00952         rep=dd_Inequality; newformat=dd_TRUE;
00953       }
00954       if (strncmp(command, "partial_enum", 12)==0 || 
00955           strncmp(command, "equality", 8)==0  ||
00956           strncmp(command, "linearity", 9)==0 ) {
00957         linearity=dd_TRUE;
00958         fgets(comsave,dd_linelenmax,f);
00959       }
00960       if (strncmp(command, "begin", 5)==0) found=dd_TRUE;
00961     }
00962   }
00963   fscanf(f, "%ld %ld %s", &m_input, &d_input, numbtype);
00964   fprintf(stderr,"size = %ld x %ld\nNumber Type = %s\n", m_input, d_input, numbtype);
00965   NT=dd_GetNumberType(numbtype);
00966   if (NT==dd_Unknown) {
00967       (*Error)=dd_ImproperInputFormat;
00968       goto _L99;
00969     } 
00970   M=dd_CreateMatrix(m_input, d_input);
00971   M->representation=rep;
00972   M->numbtype=NT;
00973 
00974   for (i = 1; i <= m_input; i++) {
00975     for (j = 1; j <= d_input; j++) {
00976       if (NT==dd_Real) {
00977 #if defined GMPRATIONAL
00978         *Error=dd_NoRealNumberSupport;
00979         goto _L99;
00980 #else
00981         fscanf(f, "%lf", &rvalue);
00982         dd_set_d(value, rvalue);
00983 #endif
00984       } else {
00985         dd_fread_rational_value (f, value);
00986       }
00987       dd_set(M->matrix[i-1][j - 1],value);
00988       if (dd_debug) {fprintf(stderr,"a(%3ld,%5ld) = ",i,j); dd_WriteNumber(stderr,value);}
00989     }  /*of j*/
00990   }  /*of i*/
00991   if (fscanf(f,"%s",command)==EOF) {
00992          (*Error)=dd_ImproperInputFormat;
00993          goto _L99;
00994   }
00995   else if (strncmp(command, "end", 3)!=0) {
00996      if (dd_debug) fprintf(stderr,"'end' missing or illegal extra data: %s\n",command);
00997      (*Error)=dd_ImproperInputFormat;
00998      goto _L99;
00999   }
01000   
01001   successful=dd_TRUE;
01002   if (linearity) {
01003     dd_SetLinearity(M,comsave);
01004   }
01005   while (!feof(f)) {
01006     fscanf(f,"%s", command);
01007     dd_ProcessCommandLine(f, M, command);
01008     fgets(command,dd_linelenmax,f); /* skip the CR/LF */
01009   } 
01010 
01011 _L99: ;
01012   dd_clear(value);
01013   /* if (f!=NULL) fclose(f); */
01014   return M;
01015 }
01016 
01017 
01018 dd_PolyhedraPtr dd_DDMatrix2Poly(dd_MatrixPtr M, dd_ErrorType *err)
01019 {
01020   dd_rowrange i;
01021   dd_colrange j;
01022   dd_PolyhedraPtr poly=NULL;
01023 
01024   *err=dd_NoError;
01025   if (M->rowsize<0 || M->colsize<0){
01026     *err=dd_NegativeMatrixSize;
01027     goto _L99;
01028   }
01029   poly=dd_CreatePolyhedraData(M->rowsize, M->colsize);
01030   poly->representation=M->representation;
01031   poly->homogeneous=dd_TRUE;
01032 
01033   for (i = 1; i <= M->rowsize; i++) {
01034     if (set_member(i, M->linset)) {
01035       poly->EqualityIndex[i]=1;
01036     }
01037     for (j = 1; j <= M->colsize; j++) {
01038       dd_set(poly->A[i-1][j-1], M->matrix[i-1][j-1]);
01039       if (j==1 && dd_Nonzero(M->matrix[i-1][j-1])) poly->homogeneous = dd_FALSE;
01040     }  /*of j*/
01041   }  /*of i*/
01042   dd_DoubleDescription(poly,err);
01043 _L99:
01044   return poly; 
01045 }
01046 
01047 dd_PolyhedraPtr dd_DDMatrix2Poly2(dd_MatrixPtr M, dd_RowOrderType horder, dd_ErrorType *err)
01048 {
01049   dd_rowrange i;
01050   dd_colrange j;
01051   dd_PolyhedraPtr poly=NULL;
01052 
01053   *err=dd_NoError;
01054   if (M->rowsize<0 || M->colsize<0){
01055     *err=dd_NegativeMatrixSize;
01056     goto _L99;
01057   }
01058   poly=dd_CreatePolyhedraData(M->rowsize, M->colsize);
01059   poly->representation=M->representation;
01060   poly->homogeneous=dd_TRUE;
01061 
01062   for (i = 1; i <= M->rowsize; i++) {
01063     if (set_member(i, M->linset)) {
01064       poly->EqualityIndex[i]=1;
01065     }
01066     for (j = 1; j <= M->colsize; j++) {
01067       dd_set(poly->A[i-1][j-1], M->matrix[i-1][j-1]);
01068       if (j==1 && dd_Nonzero(M->matrix[i-1][j-1])) poly->homogeneous = dd_FALSE;
01069     }  /*of j*/
01070   }  /*of i*/
01071   dd_DoubleDescription2(poly, horder, err);
01072 _L99:
01073   return poly; 
01074 }
01075 
01076 void dd_MatrixIntegerFilter(dd_MatrixPtr M)
01077 {   /* setting an almost integer to the integer. */
01078   dd_rowrange i;
01079   dd_colrange j;
01080   mytype x;
01081 
01082   dd_init(x);
01083   for (i=0; i< M->rowsize; i++)
01084     for (j=0; j< M->colsize; j++){
01085        dd_SnapToInteger(x, M->matrix[i][j]);
01086        dd_set(M->matrix[i][j],x);
01087     }
01088   dd_clear(x);
01089 }
01090 
01091 void dd_CopyRay(mytype *a, dd_colrange d_origsize, dd_RayPtr RR, 
01092   dd_RepresentationType rep, dd_colindex reducedcol)
01093 {
01094   long j,j1;
01095   mytype b;
01096 
01097   dd_init(b);
01098   for (j = 1; j <= d_origsize; j++){
01099     j1=reducedcol[j];
01100     if (j1>0){
01101       dd_set(a[j-1],RR->Ray[j1-1]); 
01102         /* the original column j is mapped to j1, and thus
01103            copy the corresponding component */
01104     } else {
01105       dd_set(a[j-1],dd_purezero);  
01106         /* original column is redundant and removed for computation */
01107     }
01108   }
01109 
01110   dd_set(b,a[0]);
01111   if (rep==dd_Generator && dd_Nonzero(b)){
01112     dd_set(a[0],dd_one);
01113     for (j = 2; j <= d_origsize; j++)
01114        dd_div(a[j-1],a[j-1],b);    /* normalization for generators */
01115   }
01116   dd_clear(b);
01117 }
01118 
01119 void dd_WriteRay(FILE *f, dd_colrange d_origsize, dd_RayPtr RR, dd_RepresentationType rep, dd_colindex reducedcol)
01120 {
01121   dd_colrange j;
01122   static dd_colrange d_last=0;
01123   static dd_Arow a;
01124 
01125   if (d_last< d_origsize){
01126     if (d_last>0) free(a);
01127     dd_InitializeArow(d_origsize+1, &a);
01128     d_last=d_origsize+1;
01129   }
01130 
01131  dd_CopyRay(a, d_origsize, RR, rep, reducedcol);
01132   for (j = 0; j < d_origsize; j++) dd_WriteNumber(f, a[j]);
01133   fprintf(f, "\n");
01134 }
01135 
01136 void dd_WriteArow(FILE *f, dd_Arow a, dd_colrange d)
01137 {
01138   dd_colrange j;
01139 
01140   for (j = 0; j < d; j++) dd_WriteNumber(f, a[j]);
01141   fprintf(f, "\n");
01142 }
01143 
01144 void dd_WriteAmatrix(FILE *f, dd_Amatrix A, long rowmax, long colmax)
01145 {
01146   long i,j;
01147 
01148   if (A==NULL){
01149     fprintf(f, "WriteAmatrix: The requested matrix is empty\n");
01150     goto _L99;
01151   }
01152   fprintf(f, "begin\n");
01153 #if defined GMPRATIONAL
01154   fprintf(f, " %ld %ld rational\n",rowmax, colmax);
01155 #else
01156   fprintf(f, " %ld %ld real\n",rowmax, colmax);
01157 #endif
01158   for (i=1; i <= rowmax; i++) {
01159     for (j=1; j <= colmax; j++) {
01160       dd_WriteNumber(f, A[i-1][j-1]);
01161     }
01162     fprintf(f,"\n");
01163   }
01164   fprintf(f, "end\n");
01165 _L99:;
01166 }
01167 
01168 void dd_WriteBmatrix(FILE *f, dd_colrange d_size, dd_Bmatrix B)
01169 {
01170   dd_colrange j1, j2;
01171 
01172   if (B==NULL){
01173     fprintf(f, "WriteBmatrix: The requested matrix is empty\n");
01174     goto _L99;
01175   }
01176   for (j1 = 0; j1 < d_size; j1++) {
01177     for (j2 = 0; j2 < d_size; j2++) {
01178       dd_WriteNumber(f, B[j1][j2]);
01179     }  /*of j2*/
01180     putc('\n', f);
01181   }  /*of j1*/
01182   putc('\n', f);
01183 _L99:;
01184 }
01185 
01186 void dd_WriteSetFamily(FILE *f, dd_SetFamilyPtr F)
01187 {
01188   dd_bigrange i;
01189 
01190   if (F==NULL){
01191     fprintf(f, "WriteSetFamily: The requested family is empty\n");
01192     goto _L99;
01193   }
01194   fprintf(f,"begin\n");
01195   fprintf(f,"  %ld    %ld\n", F->famsize, F->setsize);
01196   for (i=0; i<F->famsize; i++) {
01197     fprintf(f, " %ld %ld : ", i+1, set_card(F->set[i]));
01198     set_fwrite(f, F->set[i]);
01199   }
01200   fprintf(f,"end\n");
01201 _L99:;
01202 }
01203 
01204 void dd_WriteSetFamilyCompressed(FILE *f, dd_SetFamilyPtr F)
01205 {
01206   dd_bigrange i,card;
01207 
01208   if (F==NULL){
01209     fprintf(f, "WriteSetFamily: The requested family is empty\n");
01210     goto _L99;
01211   }
01212   fprintf(f,"begin\n");
01213   fprintf(f,"  %ld    %ld\n", F->famsize, F->setsize);
01214   for (i=0; i<F->famsize; i++) {
01215     card=set_card(F->set[i]);
01216     if (F->setsize - card >= card){
01217       fprintf(f, " %ld %ld : ", i+1, card);
01218       set_fwrite(f, F->set[i]);
01219     } else {
01220       fprintf(f, " %ld %ld : ", i+1, -card);
01221       set_fwrite_compl(f, F->set[i]);
01222     }
01223   }
01224   fprintf(f,"end\n");
01225 _L99:;
01226 }
01227 
01228 void dd_WriteMatrix(FILE *f, dd_MatrixPtr M)
01229 {
01230   dd_rowrange i, linsize;
01231 
01232   if (M==NULL){
01233     fprintf(f, "WriteMatrix: The requested matrix is empty\n");
01234     goto _L99;
01235   }
01236   switch (M->representation) {
01237     case dd_Inequality:
01238       fprintf(f, "H-representation\n"); break;
01239     case dd_Generator:
01240       fprintf(f, "V-representation\n"); break;
01241     case dd_Unspecified:
01242       break;
01243   }
01244   linsize=set_card(M->linset);
01245   if (linsize>0) {
01246     fprintf(f, "linearity %ld ", linsize);
01247     for (i=1; i<=M->rowsize; i++) 
01248       if (set_member(i, M->linset)) fprintf(f, " %ld", i);
01249     fprintf(f, "\n");
01250   }
01251   dd_WriteAmatrix(f, M->matrix, M->rowsize, M->colsize);
01252   if (M->objective!=dd_LPnone){
01253     if (M->objective==dd_LPmax)
01254       fprintf(f, "maximize\n");
01255     else
01256       fprintf(f, "minimize\n");      
01257     dd_WriteArow(f, M->rowvec, M->colsize);
01258   }
01259 _L99:;
01260 }
01261 
01262 void dd_WriteLP(FILE *f, dd_LPPtr lp)
01263 {
01264   if (lp==NULL){
01265     fprintf(f, "WriteLP: The requested lp is empty\n");
01266     goto _L99;
01267   }
01268   fprintf(f, "H-representation\n");
01269 
01270   dd_WriteAmatrix(f, lp->A, (lp->m)-1, lp->d);
01271   if (lp->objective!=dd_LPnone){
01272     if (lp->objective==dd_LPmax)
01273       fprintf(f, "maximize\n");
01274     else
01275       fprintf(f, "minimize\n");      
01276     dd_WriteArow(f, lp->A[lp->objrow-1], lp->d);
01277   }
01278 _L99:;
01279 }
01280 
01281 
01282 void dd_SnapToInteger(mytype y, mytype x)
01283 {
01284  /* this is broken.  It does nothing. */
01285    dd_set(y,x);
01286 }
01287 
01288 
01289 void dd_WriteReal(FILE *f, mytype x)
01290 {
01291   long ix1,ix2,ix;
01292   double ax;
01293 
01294   ax=dd_get_d(x);  
01295   ix1= (long) (fabs(ax) * 10000. + 0.5);
01296   ix2= (long) (fabs(ax) + 0.5);
01297   ix2= ix2*10000;
01298   if ( ix1 == ix2) {
01299     if (dd_Positive(x)) {
01300       ix = (long)(ax + 0.5);
01301     } else {
01302       ix = (long)(-ax + 0.5);
01303       ix = -ix;
01304     }
01305     fprintf(f, " %2ld", ix);
01306   } else
01307     fprintf(f, " % .9E",ax);
01308 }
01309 
01310 void dd_WriteNumber(FILE *f, mytype x)
01311 {
01312 #if defined GMPRATIONAL
01313   mpz_t zn,zd;
01314 
01315   mpz_init(zn); mpz_init(zd);
01316   mpq_canonicalize(x);
01317   mpq_get_num(zn,x);
01318   mpq_get_den(zd,x);
01319   fprintf(f," ");
01320   if (mpz_sgn(zn)==0){
01321     fprintf(f,"0");
01322   } else if (mpz_cmp_ui(zd,1U)==0){
01323     mpz_out_str(f,10,zn);
01324   } else {
01325     mpz_out_str(f,10,zn);fprintf(f,"/");mpz_out_str(f,10,zd);
01326   }
01327   mpz_clear(zn); mpz_clear(zd);
01328 #else
01329   dd_WriteReal(f, x);
01330 #endif
01331 }
01332 
01333 
01334 void dd_WriteIncidence(FILE *f, dd_PolyhedraPtr poly)
01335 {
01336   dd_SetFamilyPtr I;
01337 
01338   switch (poly->representation) {
01339     case dd_Inequality:
01340        fprintf(f, "ecd_file: Incidence of generators and inequalities\n");
01341       break;
01342     case dd_Generator:
01343        fprintf(f, "icd_file: Incidence of inequalities and generators\n");
01344       break;
01345 
01346     default:
01347       break;
01348   }
01349   I=dd_CopyIncidence(poly);
01350   dd_WriteSetFamilyCompressed(f,I);
01351   dd_FreeSetFamily(I);
01352 }
01353 
01354 void dd_WriteAdjacency(FILE *f, dd_PolyhedraPtr poly)
01355 {
01356   dd_SetFamilyPtr A;
01357 
01358   switch (poly->representation) {
01359     case dd_Inequality:
01360        fprintf(f, "ead_file: Adjacency of generators\n");
01361       break;
01362     case dd_Generator:
01363        fprintf(f, "iad_file: Adjacency of inequalities\n");
01364       break;
01365 
01366     default:
01367       break;
01368   }
01369   A=dd_CopyAdjacency(poly);
01370   dd_WriteSetFamilyCompressed(f,A);
01371   dd_FreeSetFamily(A);
01372 }
01373 
01374 
01375 void dd_ComputeAinc(dd_PolyhedraPtr poly)
01376 {
01377 /* This generates the input incidence array poly->Ainc, and
01378    two sets: poly->Ared, poly->Adom. 
01379 */
01380   dd_bigrange k;
01381   dd_rowrange i,m1;
01382   dd_colrange j;
01383   dd_boolean redundant;
01384   dd_MatrixPtr M=NULL;
01385   mytype sum,temp;
01386 
01387   dd_init(sum); dd_init(temp);
01388   if (poly->AincGenerated==dd_TRUE) goto _L99;
01389 
01390   M=dd_CopyOutput(poly);
01391   poly->n=M->rowsize;
01392   m1=poly->m1;  
01393    /* this number is same as poly->m, except when
01394       poly is given by nonhomogeneous inequalty:
01395       !(poly->homogeneous) && poly->representation==Inequality,
01396       it is poly->m+1.   See dd_ConeDataLoad.
01397    */
01398   poly->Ainc=(set_type*)calloc(m1, sizeof(set_type));
01399   for(i=1; i<=m1; i++) set_initialize(&(poly->Ainc[i-1]),poly->n);
01400   set_initialize(&(poly->Ared), m1); 
01401   set_initialize(&(poly->Adom), m1); 
01402 
01403   for (k=1; k<=poly->n; k++){
01404     for (i=1; i<=poly->m; i++){
01405       dd_set(sum,dd_purezero);
01406       for (j=1; j<=poly->d; j++){
01407         dd_mul(temp,poly->A[i-1][j-1],M->matrix[k-1][j-1]);
01408         dd_add(sum,sum,temp);
01409       }
01410       if (dd_EqualToZero(sum)) {
01411         set_addelem(poly->Ainc[i-1], k);
01412       }
01413     }
01414     if (!(poly->homogeneous) && poly->representation==dd_Inequality){
01415       if (dd_EqualToZero(M->matrix[k-1][0])) {
01416         set_addelem(poly->Ainc[m1-1], k);  /* added infinity inequality (1,0,0,...,0) */
01417       }
01418     }
01419   }
01420 
01421   for (i=1; i<=m1; i++){
01422     if (set_card(poly->Ainc[i-1])==M->rowsize){
01423       set_addelem(poly->Adom, i);
01424     }  
01425   }
01426   for (i=m1; i>=1; i--){
01427     if (set_card(poly->Ainc[i-1])==0){
01428       redundant=dd_TRUE;
01429       set_addelem(poly->Ared, i);
01430     }else {
01431       redundant=dd_FALSE;
01432       for (k=1; k<=m1; k++) {
01433         if (k!=i && !set_member(k, poly->Ared)  && !set_member(k, poly->Adom) && 
01434             set_subset(poly->Ainc[i-1], poly->Ainc[k-1])){
01435           if (!redundant){
01436             redundant=dd_TRUE;
01437           }
01438           set_addelem(poly->Ared, i);
01439         }
01440       }
01441     }
01442   }
01443   dd_FreeMatrix(M);
01444   poly->AincGenerated=dd_TRUE;
01445 _L99:;
01446   dd_clear(sum);  dd_clear(temp);
01447 }
01448 
01449 dd_boolean dd_InputAdjacentQ(dd_PolyhedraPtr poly, 
01450   dd_rowrange i1, dd_rowrange i2)
01451 /* Before calling this function, RedundantSet must be 
01452    a set of row indices whose removal results in a minimal
01453    nonredundant system to represent the input polyhedron,
01454    DominantSet must be the set of row indices which are
01455    active at every extreme points/rays.
01456 */
01457 {
01458   dd_boolean adj=dd_TRUE;
01459   dd_rowrange i;
01460   static set_type common;
01461   static long lastn=0;
01462 
01463   if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
01464   if (lastn!=poly->n){
01465     if (lastn >0) set_free(common);
01466     set_initialize(&common, poly->n);
01467     lastn=poly->n;
01468   }
01469   if (set_member(i1, poly->Ared) || set_member(i2, poly->Ared)){
01470     adj=dd_FALSE;
01471     goto _L99;
01472   }
01473   if (set_member(i1, poly->Adom) || set_member(i2, poly->Adom)){
01474     /* dominant inequality is considered adjacencent to all others. */
01475     adj=dd_TRUE;
01476     goto _L99;
01477   }
01478   set_int(common, poly->Ainc[i1-1], poly->Ainc[i2-1]);
01479   i=0;
01480   while (i<poly->m1 && adj==dd_TRUE){ 
01481     i++; 
01482     if (i!=i1 && i!=i2 && !set_member(i, poly->Ared) &&
01483         !set_member(i, poly->Adom) && set_subset(common,poly->Ainc[i-1])){
01484       adj=dd_FALSE;
01485     }
01486   }
01487  _L99:;
01488   return adj;
01489 } 
01490 
01491 
01492 void dd_WriteInputIncidence(FILE *f, dd_PolyhedraPtr poly)
01493 {
01494   dd_SetFamilyPtr I;
01495 
01496   if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
01497   switch (poly->representation) {
01498   case dd_Inequality:
01499     fprintf(f,"icd_file: Incidence of inequalities and generators\n");
01500     break;
01501 
01502   case dd_Generator:
01503    fprintf(f,"ecd_file: Incidence of generators and inequalities\n");
01504     break;
01505 
01506   default:
01507     break;
01508   }
01509 
01510   I=dd_CopyInputIncidence(poly);
01511   dd_WriteSetFamilyCompressed(f,I);
01512   dd_FreeSetFamily(I);
01513 }
01514 
01515 void dd_WriteInputAdjacency(FILE *f, dd_PolyhedraPtr poly)
01516 {
01517   dd_SetFamilyPtr A;
01518 
01519   if (poly->AincGenerated==dd_FALSE){
01520     dd_ComputeAinc(poly);
01521   }
01522   switch (poly->representation) {
01523   case dd_Inequality:
01524     fprintf(f, "iad_file: Adjacency of inequalities\n");
01525     break;
01526 
01527   case dd_Generator:
01528     fprintf(f, "ead_file: Adjacency of generators\n");
01529     break;
01530 
01531   default:
01532     break;
01533   }
01534   A=dd_CopyInputAdjacency(poly);
01535   dd_WriteSetFamilyCompressed(f,A);
01536   dd_FreeSetFamily(A);
01537 }
01538 
01539 
01540 void dd_WriteProgramDescription(FILE *f)
01541 {
01542   fprintf(f, "* cddlib: a double description library:%s\n", dd_DDVERSION);
01543   fprintf(f, "* compiled for %s arithmetic.\n", dd_ARITHMETIC);
01544   fprintf(f,"* %s\n",dd_COPYRIGHT);
01545 }
01546 
01547 void dd_WriteTimes(FILE *f, time_t starttime, time_t endtime)
01548 {
01549   long ptime,ptime_sec,ptime_minu, ptime_hour;
01550   
01551 /* 
01552    ptime=difftime(endtime,starttime); 
01553    This function is ANSI standard, but not available sometime 
01554 */
01555   ptime=(endtime - starttime);     
01556  /* This is to replace the line above, but it may not give 
01557     correct time in seconds */ 
01558   ptime_hour=ptime/3600;
01559   ptime_minu=(ptime-ptime_hour*3600)/60;
01560   ptime_sec=ptime%60;
01561   fprintf(f,"* Computation started at %s",asctime(localtime(&starttime)));
01562   fprintf(f,"*             ended   at %s",asctime(localtime(&endtime)));
01563   fprintf(f,"* Total processor time = %ld seconds\n",ptime);
01564   fprintf(f,"*                      = %ld h %ld m %ld s\n",ptime_hour, ptime_minu, ptime_sec);
01565 }
01566 
01567 void dd_WriteDDTimes(FILE *f, dd_PolyhedraPtr poly)
01568 {
01569   dd_WriteTimes(f,poly->child->starttime,poly->child->endtime);
01570 }
01571 
01572 void dd_WriteLPTimes(FILE *f, dd_LPPtr lp)
01573 {
01574   dd_WriteTimes(f,lp->starttime,lp->endtime);
01575 }
01576 
01577 void dd_WriteLPStats(FILE *f)
01578 {
01579   time_t currenttime;
01580   
01581   time(&currenttime);
01582   
01583   fprintf(f, "\n*--- Statistics of pivots ---\n");
01584 #if defined GMPRATIONAL
01585   fprintf(f, "* f0 = %ld (float basis finding pivots)\n",ddf_statBApivots);
01586   fprintf(f, "* fc = %ld (float CC pivots)\n",ddf_statCCpivots);
01587   fprintf(f, "* f1 = %ld (float dual simplex phase I pivots)\n",ddf_statDS1pivots);
01588   fprintf(f, "* f2 = %ld (float dual simplex phase II pivots)\n",ddf_statDS2pivots);
01589   fprintf(f, "* f3 = %ld (float anticycling CC pivots)\n",ddf_statACpivots);
01590   fprintf(f, "* e0 = %ld (exact basis finding pivots)\n",dd_statBApivots);
01591   fprintf(f, "* ec = %ld (exact CC pivots)\n",dd_statCCpivots);
01592   fprintf(f, "* e1 = %ld (exact dual simplex phase I pivots)\n",dd_statDS1pivots);
01593   fprintf(f, "* e2 = %ld (exact dual simplex phase II pivots)\n",dd_statDS2pivots);
01594   fprintf(f, "* e3 = %ld (exact anticycling CC pivots)\n",dd_statACpivots);
01595   fprintf(f, "* e4 = %ld (exact basis verification pivots)\n",dd_statBSpivots);
01596 #else
01597   fprintf(f, "f0 = %ld (float basis finding pivots)\n",dd_statBApivots);
01598   fprintf(f, "fc = %ld (float CC pivots)\n",dd_statCCpivots);
01599   fprintf(f, "f1 = %ld (float dual simplex phase I pivots)\n",dd_statDS1pivots);
01600   fprintf(f, "f2 = %ld (float dual simplex phase II pivots)\n",dd_statDS2pivots);
01601   fprintf(f, "f3 = %ld (float anticycling CC pivots)\n",dd_statACpivots);
01602 #endif
01603  dd_WriteLPMode(f);
01604   dd_WriteTimes(f,dd_statStartTime, currenttime);
01605 }
01606 
01607 void dd_WriteLPMode(FILE *f)
01608 {
01609   fprintf(f, "\n* LP solver: ");
01610   switch (dd_choiceLPSolverDefault) {
01611     case dd_DualSimplex:
01612       fprintf(f, "DualSimplex\n");
01613       break;
01614     case dd_CrissCross:
01615       fprintf(f, "Criss-Cross\n");
01616       break;
01617     default: break;
01618   }
01619   
01620   fprintf(f, "* Redundancy cheking solver: ");
01621   switch (dd_choiceRedcheckAlgorithm) {
01622     case dd_DualSimplex:
01623       fprintf(f, "DualSimplex\n");
01624       break;
01625     case dd_CrissCross:
01626       fprintf(f, "Criss-Cross\n");
01627       break;
01628     default: break;
01629   }
01630   
01631   fprintf(f, "* Lexicographic pivot: ");
01632   if (dd_choiceLexicoPivotQ)  fprintf(f, " on\n"); 
01633   else fprintf(f, " off\n"); 
01634 
01635 }
01636 
01637 
01638 void dd_WriteRunningMode(FILE *f, dd_PolyhedraPtr poly)
01639 {
01640   if (poly->child!=NULL){
01641     fprintf(f,"* roworder: ");
01642     switch (poly->child->HalfspaceOrder) {
01643 
01644     case dd_MinIndex:
01645       fprintf(f, "minindex\n");
01646       break;
01647 
01648     case dd_MaxIndex:
01649       fprintf(f, "maxindex\n");
01650       break;
01651 
01652     case dd_MinCutoff:
01653       fprintf(f, "mincutoff\n");
01654       break;
01655 
01656     case dd_MaxCutoff:
01657       fprintf(f, "maxcutoff\n");
01658     break;
01659 
01660     case dd_MixCutoff:
01661       fprintf(f, "mixcutoff\n");
01662       break;
01663 
01664     case dd_LexMin:
01665       fprintf(f, "lexmin\n");
01666       break;
01667 
01668     case dd_LexMax:
01669       fprintf(f, "lexmax\n");
01670       break;
01671 
01672     case dd_RandomRow:
01673       fprintf(f, "random  %d\n",poly->child->rseed);
01674       break;
01675 
01676     default: break;
01677     }
01678   }
01679 }
01680 
01681 
01682 void dd_WriteCompletionStatus(FILE *f, dd_ConePtr cone)
01683 {
01684   if (cone->Iteration<cone->m && cone->CompStatus==dd_AllFound) {
01685     fprintf(f,"*Computation completed at Iteration %4ld.\n", cone->Iteration);
01686   } 
01687   if (cone->CompStatus == dd_RegionEmpty) {
01688     fprintf(f,"*Computation completed at Iteration %4ld because the region found empty.\n",cone->Iteration);
01689   }   
01690 }
01691 
01692 void dd_WritePolyFile(FILE *f, dd_PolyhedraPtr poly)
01693 {
01694   dd_WriteAmatrix(f,poly->A,poly->m,poly->d);
01695 }
01696 
01697 
01698 void dd_WriteErrorMessages(FILE *f, dd_ErrorType Error)
01699 {
01700   switch (Error) {
01701 
01702   case dd_DimensionTooLarge:
01703     fprintf(f, "*Input Error: Input matrix is too large:\n");
01704     fprintf(f, "*Please increase MMAX and/or NMAX in the source code and recompile.\n");
01705     break;
01706 
01707   case dd_IFileNotFound:
01708     fprintf(f, "*Input Error: Specified input file does not exist.\n");
01709     break;
01710 
01711   case dd_OFileNotOpen:
01712     fprintf(f, "*Output Error: Specified output file cannot be opened.\n");
01713     break;
01714 
01715   case dd_NegativeMatrixSize:
01716     fprintf(f, "*Input Error: Input matrix has a negative size:\n");
01717     fprintf(f, "*Please check rowsize or colsize.\n");
01718     break;
01719 
01720   case dd_ImproperInputFormat:
01721     fprintf(f,"*Input Error: Input format is not correct.\n");
01722     fprintf(f,"*Format:\n");
01723     fprintf(f," begin\n");
01724     fprintf(f,"   m   n  NumberType(real, rational or integer)\n");
01725     fprintf(f,"   b  -A\n");
01726     fprintf(f," end\n");
01727     break;
01728 
01729   case dd_EmptyVrepresentation:
01730     fprintf(f, "*Input Error: V-representation is empty:\n");
01731     fprintf(f, "*cddlib does not accept this trivial case for which output can be any inconsistent system.\n");
01732     break;
01733 
01734   case dd_EmptyHrepresentation:
01735     fprintf(f, "*Input Error: H-representation is empty.\n");
01736     break;
01737 
01738   case dd_EmptyRepresentation:
01739     fprintf(f, "*Input Error: Representation is empty.\n");
01740     break;
01741 
01742   case dd_NoLPObjective:
01743     fprintf(f, "*LP Error: No LP objective (max or min) is set.\n");
01744     break;
01745 
01746   case dd_NoRealNumberSupport:
01747     fprintf(f, "*LP Error: The binary (with GMP Rational) does not support Real number input.\n");
01748     fprintf(f, "         : Use a binary compiled without -DGMPRATIONAL option.\n");
01749     break;
01750 
01751  case dd_NotAvailForH:
01752     fprintf(f, "*Error: A function is called with H-rep which does not support an H-representation.\n");
01753     break;
01754 
01755  case dd_NotAvailForV:
01756     fprintf(f, "*Error: A function is called with V-rep which does not support an V-representation.\n");
01757     break;
01758 
01759  case dd_CannotHandleLinearity:
01760     fprintf(f, "*Error: The function called cannot handle linearity.\n");
01761     break;
01762 
01763  case dd_RowIndexOutOfRange:
01764     fprintf(f, "*Error: Specified row index is out of range\n");
01765     break;
01766 
01767  case dd_ColIndexOutOfRange:
01768     fprintf(f, "*Error: Specified column index is out of range\n");
01769     break;
01770 
01771  case dd_LPCycling:
01772     fprintf(f, "*Error: Possibly an LP cycling occurs.  Use the Criss-Cross method.\n");
01773     break;
01774     
01775  case dd_NumericallyInconsistent:
01776     fprintf(f, "*Error: Numerical inconsistency is found.  Use the GMP exact arithmetic.\n");
01777     break;
01778     
01779   case dd_NoError:
01780     fprintf(f,"*No Error found.\n");
01781     break;
01782   }
01783 }
01784 
01785 
01786 dd_SetFamilyPtr dd_CopyIncidence(dd_PolyhedraPtr poly)
01787 {
01788   dd_SetFamilyPtr F=NULL;
01789   dd_bigrange k;
01790   dd_rowrange i;
01791 
01792   if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
01793   if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
01794   F=dd_CreateSetFamily(poly->n, poly->m1);
01795   for (i=1; i<=poly->m1; i++)
01796     for (k=1; k<=poly->n; k++)
01797       if (set_member(k,poly->Ainc[i-1])) set_addelem(F->set[k-1],i);
01798 _L99:;
01799   return F;
01800 }
01801 
01802 dd_SetFamilyPtr dd_CopyInputIncidence(dd_PolyhedraPtr poly)
01803 {
01804   dd_rowrange i;
01805   dd_SetFamilyPtr F=NULL;
01806 
01807   if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
01808   if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
01809   F=dd_CreateSetFamily(poly->m1, poly->n);
01810   for(i=0; i< poly->m1; i++){
01811     set_copy(F->set[i], poly->Ainc[i]);
01812   }
01813 _L99:;
01814   return F;
01815 }
01816 
01817 dd_SetFamilyPtr dd_CopyAdjacency(dd_PolyhedraPtr poly)
01818 {
01819   dd_RayPtr RayPtr1,RayPtr2;
01820   dd_SetFamilyPtr F=NULL;
01821   long pos1, pos2;
01822   dd_bigrange lstart,k,n;
01823   set_type linset,allset;
01824   dd_boolean adj;
01825 
01826   if (poly->n==0 && poly->homogeneous && poly->representation==dd_Inequality){
01827     n=1; /* the origin (the unique vertex) should be output. */
01828   } else n=poly->n;
01829   set_initialize(&linset, n);
01830   set_initialize(&allset, n);
01831   if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
01832   F=dd_CreateSetFamily(n, n);
01833   if (n<=0) goto _L99;
01834   poly->child->LastRay->Next=NULL;
01835   for (RayPtr1=poly->child->FirstRay, pos1=1;RayPtr1 != NULL; 
01836                                 RayPtr1 = RayPtr1->Next, pos1++){
01837     for (RayPtr2=poly->child->FirstRay, pos2=1; RayPtr2 != NULL; 
01838                                         RayPtr2 = RayPtr2->Next, pos2++){
01839       if (RayPtr1!=RayPtr2){
01840         dd_CheckAdjacency(poly->child, &RayPtr1, &RayPtr2, &adj);
01841         if (adj){
01842           set_addelem(F->set[pos1-1], pos2);
01843         }
01844       }
01845     }
01846   }
01847   lstart=poly->n - poly->ldim + 1;
01848   set_compl(allset,allset);  /* allset is set to the ground set. */
01849   for (k=lstart; k<=poly->n; k++){
01850     set_addelem(linset,k);     /* linearity set */
01851     set_copy(F->set[k-1],allset);  /* linearity generator is adjacent to all */
01852   }
01853   for (k=1; k<lstart; k++){
01854     set_uni(F->set[k-1],F->set[k-1],linset);
01855      /* every generator is adjacent to all linearity generators */
01856   }
01857 _L99:;
01858   set_free(allset); set_free(linset);
01859   return F;
01860 }
01861 
01862 dd_SetFamilyPtr dd_CopyInputAdjacency(dd_PolyhedraPtr poly)
01863 {
01864   dd_rowrange i,j;
01865   dd_SetFamilyPtr F=NULL;
01866 
01867   if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
01868   if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
01869   F=dd_CreateSetFamily(poly->m1, poly->m1);
01870   for (i=1; i<=poly->m1; i++){
01871     for (j=1; j<=poly->m1; j++){
01872       if (i!=j && dd_InputAdjacentQ(poly, i, j)) {
01873         set_addelem(F->set[i-1],j);
01874       }
01875     }
01876   }
01877 _L99:;
01878   return F;
01879 }
01880 
01881 dd_MatrixPtr dd_CopyOutput(dd_PolyhedraPtr poly)
01882 {
01883   dd_RayPtr RayPtr;
01884   dd_MatrixPtr M=NULL;
01885   dd_rowrange i=0,total;
01886   dd_colrange j,j1;
01887   mytype b;
01888   dd_RepresentationType outputrep=dd_Inequality;
01889   dd_boolean outputorigin=dd_FALSE;
01890 
01891   dd_init(b);
01892   total=poly->child->LinearityDim + poly->child->FeasibleRayCount;
01893   
01894   if (poly->child->d<=0 || poly->child->newcol[1]==0) total=total-1;
01895   if (poly->representation==dd_Inequality) outputrep=dd_Generator;
01896   if (total==0 && poly->homogeneous && poly->representation==dd_Inequality){
01897     total=1;
01898     outputorigin=dd_TRUE;
01899     /* the origin (the unique vertex) should be output. */
01900   }
01901   if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
01902 
01903   M=dd_CreateMatrix(total, poly->d);
01904   RayPtr = poly->child->FirstRay;
01905   while (RayPtr != NULL) {
01906     if (RayPtr->feasible) {
01907      dd_CopyRay(M->matrix[i], poly->d, RayPtr, outputrep, poly->child->newcol);
01908       i++;  /* 086 */
01909     }
01910     RayPtr = RayPtr->Next;
01911   }
01912   for (j=2; j<=poly->d; j++){
01913     if (poly->child->newcol[j]==0){ 
01914        /* original column j is dependent on others and removed for the cone */
01915       dd_set(b,poly->child->Bsave[0][j-1]);
01916       if (outputrep==dd_Generator && dd_Positive(b)){
01917         dd_set(M->matrix[i][0],dd_one);  /* dd_Normalize */
01918         for (j1=1; j1<poly->d; j1++) 
01919           dd_div(M->matrix[i][j1],(poly->child->Bsave[j1][j-1]),b);
01920       } else {
01921         for (j1=0; j1<poly->d; j1++)
01922           dd_set(M->matrix[i][j1],poly->child->Bsave[j1][j-1]);
01923       }
01924       set_addelem(M->linset, i+1);
01925       i++;
01926     }     
01927   }
01928   if (outputorigin){ 
01929     /* output the origin for homogeneous H-polyhedron with no rays. */
01930     dd_set(M->matrix[0][0],dd_one);
01931     for (j=1; j<poly->d; j++){
01932       dd_set(M->matrix[0][j],dd_purezero);
01933     }
01934   }
01935   dd_MatrixIntegerFilter(M);
01936   if (poly->representation==dd_Inequality)
01937     M->representation=dd_Generator;
01938   else
01939     M->representation=dd_Inequality;
01940 _L99:;
01941   dd_clear(b);
01942   return M;
01943 }
01944 
01945 dd_MatrixPtr dd_CopyInput(dd_PolyhedraPtr poly)
01946 {
01947   dd_MatrixPtr M=NULL;
01948   dd_rowrange i;
01949 
01950   M=dd_CreateMatrix(poly->m, poly->d);
01951   dd_CopyAmatrix(M->matrix, poly->A, poly->m, poly->d);
01952   for (i=1; i<=poly->m; i++) 
01953     if (poly->EqualityIndex[i]==1) set_addelem(M->linset,i);
01954   dd_MatrixIntegerFilter(M);
01955   if (poly->representation==dd_Generator)
01956     M->representation=dd_Generator;
01957   else
01958     M->representation=dd_Inequality;
01959   return M;
01960 }
01961 
01962 dd_MatrixPtr dd_CopyGenerators(dd_PolyhedraPtr poly)
01963 {
01964   dd_MatrixPtr M=NULL;
01965 
01966   if (poly->representation==dd_Generator){
01967     M=dd_CopyInput(poly);
01968   } else {
01969     M=dd_CopyOutput(poly);
01970   }
01971   return M;
01972 }
01973 
01974 dd_MatrixPtr dd_CopyInequalities(dd_PolyhedraPtr poly)
01975 {
01976   dd_MatrixPtr M=NULL;
01977 
01978   if (poly->representation==dd_Inequality){
01979     M=dd_CopyInput(poly);
01980   } else {
01981     M=dd_CopyOutput(poly);
01982   }
01983   return M;
01984 }
01985 
01986 /****************************************************************************************/
01987 /*  rational number (a/b) read is taken from Vinci by Benno Bueeler and Andreas Enge    */
01988 /****************************************************************************************/
01989 void dd_sread_rational_value (char *s, mytype value)
01990    /* reads a rational value from the specified string "s" and assigns it to "value"    */
01991    
01992 {
01993    char     *numerator_s=NULL, *denominator_s=NULL, *position;
01994    int      sign = 1;
01995    double   numerator, denominator;
01996 #if defined GMPRATIONAL
01997    mpz_t znum, zden;
01998 #else
01999    double rvalue;
02000 #endif
02001   
02002    /* determine the sign of the number */
02003    numerator_s = s;
02004    if (s [0] == '-')
02005    {  sign = -1;
02006       numerator_s++;
02007    }
02008    else if (s [0] == '+')
02009       numerator_s++;
02010       
02011    /* look for a sign '/' and eventually split the number in numerator and denominator */
02012    position = strchr (numerator_s, '/');
02013    if (position != NULL)
02014    {  *position = '\0'; /* terminates the numerator */
02015       denominator_s = position + 1;
02016    };
02017 
02018    /* determine the floating point values of numerator and denominator */
02019    numerator=atol (numerator_s);
02020   
02021    if (position != NULL)
02022    { 
02023      denominator=atol (denominator_s);  
02024    }
02025    else denominator = 1;
02026 
02027 /* 
02028    fprintf(stderr,"\nrational_read: numerator %f\n",numerator);
02029    fprintf(stderr,"rational_read: denominator %f\n",denominator);
02030    fprintf(stderr,"rational_read: sign %d\n",sign); 
02031 */
02032 
02033 #if defined GMPRATIONAL
02034    mpz_init_set_str(znum,numerator_s,10);
02035    if (sign<0) mpz_neg(znum,znum);
02036    mpz_init(zden); mpz_set_ui(zden,1);
02037    if (denominator_s!=NULL) mpz_init_set_str(zden,denominator_s,10);
02038    mpq_set_num(value,znum); mpq_set_den(value,zden);
02039    mpq_canonicalize(value);
02040    mpz_clear(znum); mpz_clear(zden);
02041    /* num=(long)sign * numerator; */
02042    /* den=(unsigned long) denominator; */
02043    /* mpq_set_si(value, num, den); */
02044 #elif defined GMPFLOAT
02045    rvalue=sign * numerator/ (signed long) denominator;
02046    mpf_set_d(value, rvalue);
02047 #else
02048    rvalue=sign * numerator/ (signed long) denominator;
02049    ddd_set_d(value, rvalue);
02050 #endif
02051    if (dd_debug) {
02052      fprintf(stderr,"rational_read: "); 
02053      dd_WriteNumber(stderr,value); fprintf(stderr,"\n");
02054    }
02055 }
02056    
02057 
02058 void dd_fread_rational_value (FILE *f, mytype value)
02059    /* reads a rational value from the specified file "f" and assigns it to "value"      */
02060    
02061 {
02062    char     number_s [255];
02063    mytype rational_value;
02064    
02065    dd_init(rational_value);
02066    fscanf(f, "%s ", number_s);
02067    dd_sread_rational_value (number_s, rational_value);
02068    dd_set(value,rational_value);
02069    dd_clear(rational_value);
02070 }
02071    
02072 /****************************************************************************************/
02073 
02074 
02075 /* end of cddio.c */
02076 
02077 #endif

 

© Sandia Corporation | Site Contact | Privacy and Security

Generated on Fri Feb 16 10:33:35 2007 for APPSPACK 5.0.1 by doxygen 1.3.9.1 written by Dimitri van Heesch, © 1997-2002