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

cddlib.c

Go to the documentation of this file.
00001 
00043 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00044 
00045 /* cddlib.c: cdd library  (library version of cdd)
00046    written by Komei Fukuda, fukuda@ifor.math.ethz.ch
00047    Version 0.94b, Aug. 25, 2005
00048    Standard ftp site: ftp.ifor.math.ethz.ch, Directory: pub/fukuda/cdd
00049 */
00050 
00051 /* cdd : C-Implementation of the double description method for
00052    computing all vertices and extreme rays of the polyhedron 
00053    P= {x :  b - A x >= 0}.
00054    Please read COPYING (GNU General Public Licence) and
00055    the manual cddlibman.tex for detail.
00056 */
00057 
00058 /*  This program is free software; you can redistribute it and/or modify
00059     it under the terms of the GNU General Public License as published by
00060     the Free Software Foundation; either version 2 of the License, or
00061     (at your option) any later version.
00062 
00063     This program is distributed in the hope that it will be useful,
00064     but WITHOUT ANY WARRANTY; without even the implied warranty of
00065     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00066     GNU General Public License for more details.
00067 
00068     You should have received a copy of the GNU General Public License
00069     along with this program; if not, write to the Free Software
00070     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00071 */
00072 
00073 /* The first version C0.21 was created on November 10,1993 
00074    with Dave Gillespie's p2c translator 
00075    from the Pascal program pdd.p written by Komei Fukuda. 
00076 */
00077 
00078 #include "setoper.h" 
00079   /* set operation library header (Ver. June 1, 2000 or later) */
00080 #include "cdd.h"
00081 #include <stdio.h>
00082 #include <stdlib.h>
00083 #include <time.h>
00084 #include <math.h>
00085 #include <string.h>
00086 
00087 /* Global Variables */
00088 dd_boolean dd_debug               =dd_FALSE;
00089 dd_boolean dd_log                 =dd_FALSE;
00090 /* GLOBAL CONSTANTS and STATICS VARIABLES (to be set by dd_set_global_constants() */
00091 mytype dd_zero;
00092 mytype dd_one;
00093 mytype dd_purezero;
00094 mytype dd_minuszero;
00095 mytype dd_minusone;
00096 
00097 time_t dd_statStartTime; /* cddlib starting time */
00098 long dd_statBApivots;  /* basis finding pivots */
00099 long dd_statCCpivots;  /* criss-cross pivots */
00100 long dd_statDS1pivots; /* phase 1 pivots */
00101 long dd_statDS2pivots; /* phase 2 pivots */
00102 long dd_statACpivots;  /* anticycling (cc) pivots */
00103 #ifdef GMPRATIONAL
00104 long dd_statBSpivots;  /* basis status checking pivots */
00105 #endif
00106 dd_LPSolverType dd_choiceLPSolverDefault;  /* Default LP solver Algorithm */
00107 dd_LPSolverType dd_choiceRedcheckAlgorithm;  /* Redundancy Checking Algorithm */
00108 dd_boolean dd_choiceLexicoPivotQ;    /* whether to use the lexicographic pivot */
00109 
00110 /* #include <profile.h>    THINK C PROFILER */
00111 /* #include <console.h>    THINK C PROFILER */
00112 
00113 void dd_DDInit(dd_ConePtr cone)
00114 {
00115   cone->Error=dd_NoError;
00116   cone->CompStatus=dd_InProgress;
00117   cone->RayCount = 0;
00118   cone->TotalRayCount = 0;
00119   cone->FeasibleRayCount = 0;
00120   cone->WeaklyFeasibleRayCount = 0;
00121   cone->EdgeCount=0; /* active edge count */
00122   cone->TotalEdgeCount=0; /* active edge count */
00123   dd_SetInequalitySets(cone);
00124   dd_ComputeRowOrderVector(cone);
00125   cone->RecomputeRowOrder=dd_FALSE;
00126 }
00127 
00128 void dd_DDMain(dd_ConePtr cone)
00129 {
00130   dd_rowrange hh, itemp, otemp;
00131   dd_boolean locallog=dd_log; /* if dd_log=dd_FALSE, no log will be written.  */
00132 
00133   if (cone->d<=0){
00134     cone->Iteration=cone->m;
00135     cone->FeasibleRayCount=0;
00136     cone->CompStatus=dd_AllFound;
00137     goto _L99;
00138   }
00139   if (locallog) {
00140      fprintf(stderr,"(Initially added rows ) = ");
00141      set_fwrite(stderr,cone->InitialHalfspaces);
00142   }
00143   while (cone->Iteration <= cone->m) {
00144     dd_SelectNextHalfspace(cone, cone->WeaklyAddedHalfspaces, &hh);
00145     if (set_member(hh,cone->NonequalitySet)){  /* Skip the row hh */
00146       if (dd_debug) {
00147         fprintf(stderr,"*The row # %3ld should be inactive and thus skipped.\n", hh);
00148       }
00149       set_addelem(cone->WeaklyAddedHalfspaces, hh);
00150     } else {
00151       if (cone->PreOrderedRun)
00152         dd_AddNewHalfspace2(cone, hh);
00153       else{
00154         dd_AddNewHalfspace1(cone, hh);
00155       }
00156       set_addelem(cone->AddedHalfspaces, hh);
00157       set_addelem(cone->WeaklyAddedHalfspaces, hh);
00158     }
00159     if (!cone->PreOrderedRun){
00160       for (itemp=1; cone->OrderVector[itemp]!=hh; itemp++);
00161         otemp=cone->OrderVector[cone->Iteration];
00162       cone->OrderVector[cone->Iteration]=hh;
00163         /* store the dynamic ordering in ordervec */
00164       cone->OrderVector[itemp]=otemp;
00165         /* store the dynamic ordering in ordervec */
00166     }
00167     if (locallog){
00168       fprintf(stderr,"(Iter, Row, #Total, #Curr, #Feas)= %5ld %5ld %9ld %6ld %6ld\n",
00169         cone->Iteration, hh, cone->TotalRayCount, cone->RayCount,
00170         cone->FeasibleRayCount);
00171     }
00172     if (cone->CompStatus==dd_AllFound||cone->CompStatus==dd_RegionEmpty) {
00173       set_addelem(cone->AddedHalfspaces, hh);
00174       goto _L99;
00175     }
00176     (cone->Iteration)++;
00177   }
00178   _L99:;
00179   if (cone->d<=0 || cone->newcol[1]==0){ /* fixing the number of output */
00180      cone->parent->n=cone->LinearityDim + cone->FeasibleRayCount -1;
00181      cone->parent->ldim=cone->LinearityDim - 1;
00182   } else {
00183     cone->parent->n=cone->LinearityDim + cone->FeasibleRayCount;
00184     cone->parent->ldim=cone->LinearityDim;
00185   }
00186 }
00187 
00188 
00189 void dd_InitialDataSetup(dd_ConePtr cone)
00190 {
00191   long j, r;
00192   dd_rowset ZSet;
00193   static dd_Arow Vector1,Vector2;
00194   static dd_colrange last_d=0;
00195 
00196   if (last_d < cone->d){
00197     if (last_d>0) {
00198     for (j=0; j<last_d; j++){
00199       dd_init(Vector1[j]);
00200       dd_init(Vector2[j]);
00201     }
00202     free(Vector1); free(Vector2);
00203     }
00204     Vector1=(mytype*)calloc(cone->d,sizeof(mytype));
00205     Vector2=(mytype*)calloc(cone->d,sizeof(mytype));
00206     for (j=0; j<cone->d; j++){
00207       dd_init(Vector1[j]);
00208       dd_init(Vector2[j]);
00209     }
00210     last_d=cone->d;
00211   }
00212 
00213   cone->RecomputeRowOrder=dd_FALSE;
00214   cone->ArtificialRay = NULL;
00215   cone->FirstRay = NULL;
00216   cone->LastRay = NULL;
00217   set_initialize(&ZSet,cone->m);
00218   dd_AddArtificialRay(cone);
00219   set_copy(cone->AddedHalfspaces, cone->InitialHalfspaces);
00220   set_copy(cone->WeaklyAddedHalfspaces, cone->InitialHalfspaces);
00221   dd_UpdateRowOrderVector(cone, cone->InitialHalfspaces);
00222   for (r = 1; r <= cone->d; r++) {
00223     for (j = 0; j < cone->d; j++){
00224       dd_set(Vector1[j], cone->B[j][r-1]);
00225       dd_neg(Vector2[j], cone->B[j][r-1]);
00226     }
00227     dd_Normalize(cone->d, Vector1);
00228     dd_Normalize(cone->d, Vector2);
00229     dd_ZeroIndexSet(cone->m, cone->d, cone->A, Vector1, ZSet);
00230     if (set_subset(cone->EqualitySet, ZSet)){
00231       if (dd_debug) {
00232         fprintf(stderr,"add an initial ray with zero set:");
00233         set_fwrite(stderr,ZSet);
00234       }
00235       dd_AddRay(cone, Vector1);
00236       if (cone->InitialRayIndex[r]==0) {
00237         dd_AddRay(cone, Vector2);
00238         if (dd_debug) {
00239           fprintf(stderr,"and add its negative also.\n");
00240         }
00241       }
00242     }
00243   }
00244   dd_CreateInitialEdges(cone);
00245   cone->Iteration = cone->d + 1;
00246   if (cone->Iteration > cone->m) cone->CompStatus=dd_AllFound; /* 0.94b  */
00247   set_free(ZSet);
00248 }
00249 
00250 dd_boolean dd_CheckEmptiness(dd_PolyhedraPtr poly, dd_ErrorType *err)
00251 {
00252   dd_rowset R, S;
00253   dd_MatrixPtr M=NULL;
00254   dd_boolean answer=dd_FALSE;
00255 
00256   *err=dd_NoError;
00257 
00258   if (poly->representation==dd_Inequality){
00259         M=dd_CopyInequalities(poly);
00260         set_initialize(&R, M->rowsize);
00261         set_initialize(&S, M->rowsize);
00262         if (!dd_ExistsRestrictedFace(M, R, S, err)){
00263           poly->child->CompStatus=dd_AllFound;
00264           poly->IsEmpty=dd_TRUE;
00265           poly->n=0;
00266           answer=dd_TRUE;
00267         }
00268         set_free(R);
00269         set_free(S);
00270         dd_FreeMatrix(M);
00271   } else if (poly->representation==dd_Generator && poly->m<=0){
00272         *err=dd_EmptyVrepresentation;
00273         poly->IsEmpty=dd_TRUE;
00274         poly->child->CompStatus=dd_AllFound;
00275         answer=dd_TRUE;
00276         poly->child->Error=*err;  
00277   }
00278   
00279   return answer;
00280 }
00281 
00282 
00283 dd_boolean dd_DoubleDescription(dd_PolyhedraPtr poly, dd_ErrorType *err)
00284 {
00285   dd_ConePtr cone=NULL;
00286   dd_boolean found=dd_FALSE;
00287 
00288   *err=dd_NoError;
00289 
00290   if (poly!=NULL && (poly->child==NULL || poly->child->CompStatus!=dd_AllFound)){
00291     cone=dd_ConeDataLoad(poly);
00292     /* create a cone associated with poly by homogenization */
00293     time(&cone->starttime);
00294     dd_DDInit(cone);
00295     if (poly->representation==dd_Generator && poly->m<=0){
00296        *err=dd_EmptyVrepresentation;
00297        cone->Error=*err;
00298            goto _L99;
00299         }
00300         /* Check emptiness of the polyhedron */
00301         dd_CheckEmptiness(poly,err);
00302 
00303     if (cone->CompStatus!=dd_AllFound){
00304           dd_FindInitialRays(cone, &found);
00305           if (found) {
00306             dd_InitialDataSetup(cone);
00307             if (cone->CompStatus==dd_AllFound) goto _L99;
00308             dd_DDMain(cone);
00309             if (cone->FeasibleRayCount!=cone->RayCount) *err=dd_NumericallyInconsistent; /* cddlib-093d */
00310           }
00311         }
00312     time(&cone->endtime);
00313   }
00314     
00315 _L99: ;
00316 
00317   return found;
00318 }
00319 
00320 dd_boolean dd_DoubleDescription2(dd_PolyhedraPtr poly, dd_RowOrderType horder, dd_ErrorType *err)
00321 {
00322   dd_ConePtr cone=NULL;
00323   dd_boolean found=dd_FALSE;
00324 
00325   *err=dd_NoError;
00326 
00327   if (poly!=NULL && (poly->child==NULL || poly->child->CompStatus!=dd_AllFound)){
00328     cone=dd_ConeDataLoad(poly);
00329     /* create a cone associated with poly by homogenization */
00330         cone->HalfspaceOrder=horder;  /* set the row order */
00331     time(&cone->starttime);
00332     dd_DDInit(cone);
00333     if (poly->representation==dd_Generator && poly->m<=0){
00334        *err=dd_EmptyVrepresentation;
00335        cone->Error=*err;
00336            goto _L99;
00337         }
00338         /* Check emptiness of the polyhedron */
00339         dd_CheckEmptiness(poly,err);
00340 
00341     if (cone->CompStatus!=dd_AllFound){
00342           dd_FindInitialRays(cone, &found);
00343           if (found) {
00344             dd_InitialDataSetup(cone);
00345             if (cone->CompStatus==dd_AllFound) goto _L99;
00346             dd_DDMain(cone);
00347             if (cone->FeasibleRayCount!=cone->RayCount) *err=dd_NumericallyInconsistent; /* cddlib-093d */
00348           }
00349         }
00350     time(&cone->endtime);
00351   }
00352     
00353 _L99: ;
00354 
00355   return found;
00356 }
00357 
00358 dd_boolean dd_DDInputAppend(dd_PolyhedraPtr *poly, dd_MatrixPtr M,
00359   dd_ErrorType *err)
00360 {
00361   /* This is imcomplete.  It simply solves the problem from scratch.  */
00362   dd_boolean found;
00363   dd_ErrorType error;
00364 
00365   if ((*poly)->child!=NULL) dd_FreeDDMemory(*poly);
00366   dd_AppendMatrix2Poly(poly, M);
00367   (*poly)->representation=dd_Inequality;
00368   found=dd_DoubleDescription(*poly, &error);
00369   *err=error;
00370   return found;
00371 }
00372 
00373 dd_boolean dd_DDFile2File(char *ifile, char *ofile, dd_ErrorType *err)
00374 {
00375   /* The representation conversion from an input file to an outfile.  */
00376   /* modified by D. Avis to allow stdin/stdout */
00377   dd_boolean found=dd_TRUE;
00378   FILE *reading=NULL,*writing=NULL;
00379   dd_PolyhedraPtr poly;
00380   dd_MatrixPtr M, A, G;
00381 
00382   if (strcmp(ifile,"**stdin") == 0 )
00383      reading = stdin;
00384   else if ( ( reading = fopen(ifile, "r") )!= NULL) {
00385     fprintf(stderr,"input file %s is open\n", ifile);
00386    }
00387   else{
00388     fprintf(stderr,"The input file %s not found\n",ifile);
00389     found=dd_FALSE;
00390     *err=dd_IFileNotFound;
00391     goto _L99;
00392   }
00393 
00394   if (found){
00395     if (strcmp(ofile,"**stdout") == 0 )
00396       writing = stdout;
00397     else if ( (writing = fopen(ofile, "w") ) != NULL){
00398       fprintf(stderr,"output file %s is open\n",ofile);
00399       found=dd_TRUE;
00400     } else {
00401       fprintf(stderr,"The output file %s cannot be opened\n",ofile);
00402       found=dd_FALSE;
00403       *err=dd_OFileNotOpen;
00404       goto _L99;
00405     }
00406   }
00407 
00408   M=dd_PolyFile2Matrix(reading, err);
00409   if (*err!=dd_NoError){
00410     goto _L99;
00411   } poly=dd_DDMatrix2Poly(M, err);  /* compute the second representation */ dd_FreeMatrix(M);
00412 
00413   if (*err==dd_NoError) {
00414     dd_WriteRunningMode(writing, poly);
00415     A=dd_CopyInequalities(poly);
00416     G=dd_CopyGenerators(poly);
00417 
00418     if (poly->representation==dd_Inequality) {
00419       dd_WriteMatrix(writing,G);
00420      } else {
00421       dd_WriteMatrix(writing,A);
00422     }
00423 
00424     dd_FreePolyhedra(poly);
00425     dd_FreeMatrix(A);
00426     dd_FreeMatrix(G);
00427   } 
00428 
00429 _L99: ;
00430   if (*err!=dd_NoError) dd_WriteErrorMessages(stderr,*err);
00431   if (reading!=NULL) fclose(reading);
00432   if (writing!=NULL) fclose(writing);
00433   return found;
00434 }
00435 
00436 /* end of cddlib.c */
00437 
00438 #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