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

cddcore.c

Go to the documentation of this file.
00001 
00043 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00044 
00045 /* cddcore.c:  Core 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_CheckAdjacency(dd_ConePtr cone,
00066     dd_RayPtr *RP1, dd_RayPtr *RP2, dd_boolean *adjacent)
00067 {
00068   dd_RayPtr TempRay;
00069   dd_boolean localdebug=dd_FALSE;
00070   static dd_rowset Face, Face1;
00071   static dd_rowrange last_m=0;
00072   
00073   if (last_m!=cone->m) {
00074     if (last_m>0){
00075       set_free(Face); set_free(Face1);
00076     }
00077     set_initialize(&Face, cone->m); 
00078     set_initialize(&Face1, cone->m); 
00079     last_m=cone->m;
00080   }
00081 
00082   if (dd_debug) localdebug=dd_TRUE;
00083   *adjacent = dd_TRUE;
00084   set_int(Face1, (*RP1)->ZeroSet, (*RP2)->ZeroSet);
00085   set_int(Face, Face1, cone->AddedHalfspaces);
00086   if (set_card(Face)< cone->d - 2) {
00087     *adjacent = dd_FALSE;
00088     if (localdebug) {
00089       fprintf(stderr,"non adjacent: set_card(face) %ld < %ld = cone->d.\n",
00090         set_card(Face),cone->d);
00091     }
00092     return;
00093   }
00094   else if (cone->parent->NondegAssumed) {
00095         *adjacent = dd_TRUE;
00096         return;
00097   }
00098   TempRay = cone->FirstRay;
00099   while (TempRay != NULL && *adjacent) {
00100     if (TempRay != *RP1 && TempRay != *RP2) {
00101         set_int(Face1, TempRay->ZeroSet, cone->AddedHalfspaces);
00102         if (set_subset(Face, Face1)) *adjacent = dd_FALSE;
00103     }
00104     TempRay = TempRay->Next;
00105   }
00106 }
00107 
00108 void dd_Eliminate(dd_ConePtr cone, dd_RayPtr*Ptr)
00109 {
00110   /*eliminate the record pointed by Ptr^.Next*/
00111   dd_RayPtr TempPtr;
00112   dd_colrange j;
00113 
00114   TempPtr = (*Ptr)->Next;
00115   (*Ptr)->Next = (*Ptr)->Next->Next;
00116   if (TempPtr == cone->FirstRay)   /*Update the first pointer*/
00117     cone->FirstRay = (*Ptr)->Next;
00118   if (TempPtr == cone->LastRay)   /*Update the last pointer*/
00119     cone->LastRay = *Ptr;
00120 
00121   /* Added, Marc Pfetsch 010219 */
00122   for (j=0;j < cone->d;j++)
00123 dd_clear(TempPtr->Ray[j]);
00124   dd_clear(TempPtr->ARay);
00125 
00126   free(TempPtr->Ray);          /* free the ray vector memory */
00127   set_free(TempPtr->ZeroSet);  /* free the ZeroSet memory */
00128   free(TempPtr);   /* free the dd_Ray structure memory */
00129   cone->RayCount--; 
00130 }
00131 
00132 void dd_SetInequalitySets(dd_ConePtr cone)
00133 {
00134   dd_rowrange i;
00135   
00136   set_emptyset(cone->GroundSet);
00137   set_emptyset(cone->EqualitySet);
00138   set_emptyset(cone->NonequalitySet);  
00139   for (i = 1; i <= (cone->parent->m); i++){
00140     set_addelem(cone->GroundSet, i);
00141     if (cone->parent->EqualityIndex[i]==1) set_addelem(cone->EqualitySet,i);
00142     if (cone->parent->EqualityIndex[i]==-1) set_addelem(cone->NonequalitySet,i);
00143   }
00144 }
00145 
00146 
00147 void dd_AValue(mytype *val, dd_colrange d_size, dd_Amatrix A, mytype *p, dd_rowrange i)
00148 {
00149   /*return the ith component of the vector  A x p */
00150   dd_colrange j;
00151   mytype x;
00152 
00153   dd_init(x); 
00154   dd_set(*val,dd_purezero);
00155  /* Changed by Marc Pfetsch 010219 */
00156 
00157   for (j = 0; j < d_size; j++){
00158     dd_mul(x,A[i - 1][j], p[j]);
00159     dd_add(*val, *val, x);
00160   }
00161   dd_clear(x);
00162 }
00163 
00164 void dd_StoreRay1(dd_ConePtr cone, mytype *p, dd_boolean *feasible)
00165 {  /* Original ray storing routine when RelaxedEnumeration is dd_FALSE */
00166   dd_rowrange i,k,fii=cone->m+1;
00167   dd_colrange j;
00168   mytype temp;
00169   dd_RayPtr RR;
00170   dd_boolean localdebug=dd_debug;
00171 
00172   dd_init(temp);
00173   RR=cone->LastRay;
00174   *feasible = dd_TRUE;
00175   set_initialize(&(RR->ZeroSet),cone->m);
00176   for (j = 0; j < cone->d; j++){
00177     dd_set(RR->Ray[j],p[j]);
00178   }
00179   for (i = 1; i <= cone->m; i++) {
00180     k=cone->OrderVector[i];
00181     dd_AValue(&temp, cone->d, cone->A, p, k);
00182     if (localdebug) {
00183       fprintf(stderr,"dd_StoreRay1: dd_AValue at row %ld =",k);
00184       dd_WriteNumber(stderr, temp);
00185       fprintf(stderr,"\n");
00186     }
00187     if (dd_EqualToZero(temp)) {
00188       set_addelem(RR->ZeroSet, k);
00189       if (localdebug) {
00190         fprintf(stderr,"recognized zero!\n");
00191       }
00192     }
00193     if (dd_Negative(temp)){
00194       if (localdebug) {
00195         fprintf(stderr,"recognized negative!\n");
00196       }
00197       *feasible = dd_FALSE;
00198       if (fii>cone->m) fii=i;  /* the first violating inequality index */
00199       if (localdebug) {
00200         fprintf(stderr,"this ray is not feasible, neg comp = %ld\n", fii);
00201         dd_WriteNumber(stderr, temp);  fprintf(stderr,"\n");
00202       }
00203     }
00204   }
00205   RR->FirstInfeasIndex=fii;
00206   RR->feasible = *feasible;
00207   dd_clear(temp);
00208 }
00209 
00210 void dd_StoreRay2(dd_ConePtr cone, mytype *p, 
00211     dd_boolean *feasible, dd_boolean *weaklyfeasible)
00212    /* Ray storing routine when RelaxedEnumeration is dd_TRUE.
00213        weaklyfeasible is true iff it is feasible with
00214        the strict_inequality conditions deleted. */
00215 {
00216   dd_RayPtr RR;
00217   dd_rowrange i,k,fii=cone->m+1;
00218   dd_colrange j;
00219   mytype temp;
00220   dd_boolean localdebug=dd_debug;
00221 
00222   dd_init(temp);
00223   RR=cone->LastRay;
00224   if (dd_debug) localdebug=dd_TRUE;
00225   *feasible = dd_TRUE;
00226   *weaklyfeasible = dd_TRUE;
00227   set_initialize(&(RR->ZeroSet),cone->m);
00228   for (j = 0; j < cone->d; j++){
00229     dd_set(RR->Ray[j],p[j]);
00230   }
00231   for (i = 1; i <= cone->m; i++) {
00232     k=cone->OrderVector[i];
00233     dd_AValue(&temp, cone->d, cone->A, p, k);
00234     if (dd_EqualToZero(temp)){
00235       set_addelem(RR->ZeroSet, k);
00236       if (cone->parent->EqualityIndex[k]==-1) 
00237         *feasible=dd_FALSE;  /* strict inequality required */
00238     }
00239 /*    if (temp < -zero){ */
00240     if (dd_Negative(temp)){
00241       *feasible = dd_FALSE;
00242       if (fii>cone->m && cone->parent->EqualityIndex[k]>=0) {
00243         fii=i;  /* the first violating inequality index */
00244         *weaklyfeasible=dd_FALSE;
00245       }
00246     }
00247   }
00248   RR->FirstInfeasIndex=fii;
00249   RR->feasible = *feasible;
00250   dd_clear(temp);
00251 }
00252 
00253 
00254 void dd_AddRay(dd_ConePtr cone, mytype *p)
00255 {  
00256   dd_boolean feasible, weaklyfeasible;
00257   dd_colrange j;
00258 
00259   if (cone->FirstRay == NULL) {
00260     cone->FirstRay = (dd_RayPtr) malloc(sizeof(dd_RayType));
00261     cone->FirstRay->Ray = (mytype *) calloc(cone->d, sizeof(mytype));
00262     for (j=0; j<cone->d; j++) dd_init(cone->FirstRay->Ray[j]);
00263     dd_init(cone->FirstRay->ARay);
00264     if (dd_debug)
00265       fprintf(stderr,"Create the first ray pointer\n");
00266     cone->LastRay = cone->FirstRay;
00267     cone->ArtificialRay->Next = cone->FirstRay;
00268   } else {
00269     cone->LastRay->Next = (dd_RayPtr) malloc(sizeof(dd_RayType));
00270     cone->LastRay->Next->Ray = (mytype *) calloc(cone->d, sizeof(mytype));
00271     for (j=0; j<cone->d; j++) dd_init(cone->LastRay->Next->Ray[j]);
00272     dd_init(cone->LastRay->Next->ARay);
00273     if (dd_debug) fprintf(stderr,"Create a new ray pointer\n");
00274     cone->LastRay = cone->LastRay->Next;
00275   }
00276   cone->LastRay->Next = NULL;
00277   cone->RayCount++;
00278   cone->TotalRayCount++;
00279   if (dd_debug) {
00280     if (cone->TotalRayCount % 100 == 0) {
00281       fprintf(stderr,"*Rays (Total, Currently Active, Feasible) =%8ld%8ld%8ld\n",
00282          cone->TotalRayCount, cone->RayCount, cone->FeasibleRayCount);
00283     }
00284   }
00285   if (cone->parent->RelaxedEnumeration){
00286     dd_StoreRay2(cone, p, &feasible, &weaklyfeasible);
00287     if (weaklyfeasible) (cone->WeaklyFeasibleRayCount)++;
00288   } else {
00289     dd_StoreRay1(cone, p, &feasible);
00290     if (feasible) (cone->WeaklyFeasibleRayCount)++;
00291     /* weaklyfeasible is equiv. to feasible in this case. */
00292   }
00293   if (!feasible) return;
00294   else {
00295     (cone->FeasibleRayCount)++;
00296   }
00297 }
00298 
00299 void dd_AddArtificialRay(dd_ConePtr cone)
00300 {  
00301   dd_Arow zerovector;
00302   dd_colrange j,d1;
00303   dd_boolean feasible;
00304 
00305   if (cone->d<=0) d1=1; else d1=cone->d;
00306   dd_InitializeArow(d1, &zerovector);
00307   if (cone->ArtificialRay != NULL) {
00308     fprintf(stderr,"Warning !!!  FirstRay in not nil.  Illegal Call\n");
00309     free(zerovector); /* 086 */
00310     return;
00311   }
00312   cone->ArtificialRay = (dd_RayPtr) malloc(sizeof(dd_RayType));
00313   cone->ArtificialRay->Ray = (mytype *) calloc(d1, sizeof(mytype));
00314   for (j=0; j<d1; j++) dd_init(cone->ArtificialRay->Ray[j]);
00315   dd_init(cone->ArtificialRay->ARay);
00316 
00317   if (dd_debug) fprintf(stderr,"Create the artificial ray pointer\n");
00318 
00319   cone->LastRay=cone->ArtificialRay;
00320   dd_StoreRay1(cone, zerovector, &feasible);  
00321     /* This stores a vector to the record pointed by cone->LastRay */
00322   cone->ArtificialRay->Next = NULL;
00323   for (j = 0; j < d1; j++){
00324     dd_clear(zerovector[j]);
00325   }
00326   free(zerovector); /* 086 */
00327 }
00328 
00329 void dd_ConditionalAddEdge(dd_ConePtr cone, 
00330     dd_RayPtr Ray1, dd_RayPtr Ray2, dd_RayPtr ValidFirstRay)
00331 {
00332   long it,it_row,fii1,fii2,fmin,fmax;
00333   dd_boolean adjacent,lastchance;
00334   dd_RayPtr TempRay,Rmin,Rmax;
00335   dd_AdjacencyType *NewEdge;
00336   dd_boolean localdebug=dd_FALSE;
00337   dd_rowset ZSmin, ZSmax;
00338   static dd_rowset Face, Face1;
00339   static dd_rowrange last_m=0;
00340   
00341   if (last_m!=cone->m) {
00342     if (last_m>0){
00343       set_free(Face); set_free(Face1);
00344     }
00345     set_initialize(&Face, cone->m);
00346     set_initialize(&Face1, cone->m);
00347     last_m=cone->m;
00348   }
00349   
00350   fii1=Ray1->FirstInfeasIndex;
00351   fii2=Ray2->FirstInfeasIndex;
00352   if (fii1<fii2){
00353     fmin=fii1; fmax=fii2;
00354     Rmin=Ray1;
00355     Rmax=Ray2;
00356   }
00357   else{
00358     fmin=fii2; fmax=fii1;
00359     Rmin=Ray2;
00360     Rmax=Ray1;
00361   }
00362   ZSmin = Rmin->ZeroSet;
00363   ZSmax = Rmax->ZeroSet;
00364   if (localdebug) {
00365     fprintf(stderr,"dd_ConditionalAddEdge: FMIN = %ld (row%ld)   FMAX=%ld\n",
00366       fmin, cone->OrderVector[fmin], fmax);
00367   }
00368   if (fmin==fmax){
00369     if (localdebug) fprintf(stderr,"dd_ConditionalAddEdge: equal FII value-> No edge added\n");
00370   }
00371   else if (set_member(cone->OrderVector[fmin],ZSmax)){
00372     if (localdebug) fprintf(stderr,"dd_ConditionalAddEdge: No strong separation -> No edge added\n");
00373   }
00374   else {  /* the pair will be separated at the iteration fmin */
00375     lastchance=dd_TRUE;
00376     /* flag to check it will be the last chance to store the edge candidate */
00377     set_int(Face1, ZSmax, ZSmin);
00378     (cone->count_int)++;
00379     if (localdebug){
00380       fprintf(stderr,"Face: ");
00381       for (it=1; it<=cone->m; it++) {
00382         it_row=cone->OrderVector[it];
00383         if (set_member(it_row, Face1)) fprintf(stderr,"%ld ",it_row);
00384       }
00385       fprintf(stderr,"\n");
00386     }
00387     for (it=cone->Iteration+1; it < fmin && lastchance; it++){
00388       it_row=cone->OrderVector[it];
00389       if (cone->parent->EqualityIndex[it_row]>=0 && set_member(it_row, Face1)){
00390         lastchance=dd_FALSE;
00391         (cone->count_int_bad)++;
00392         if (localdebug){
00393           fprintf(stderr,"There will be another chance iteration %ld (row %ld) to store the pair\n", it, it_row);
00394         }
00395       }
00396     }
00397     if (lastchance){
00398       adjacent = dd_TRUE;
00399       (cone->count_int_good)++;
00400       /* adjacent checking */
00401       set_int(Face, Face1, cone->AddedHalfspaces);
00402       if (localdebug){
00403         fprintf(stderr,"Check adjacency\n");
00404         fprintf(stderr,"AddedHalfspaces: "); set_fwrite(stderr,cone->AddedHalfspaces);
00405         fprintf(stderr,"Face: ");
00406         for (it=1; it<=cone->m; it++) {
00407           it_row=cone->OrderVector[it];
00408           if (set_member(it_row, Face)) fprintf(stderr,"%ld ",it_row);
00409         }
00410         fprintf(stderr,"\n");
00411       }
00412       if (set_card(Face)< cone->d - 2) {
00413         adjacent = dd_FALSE;
00414       }
00415       else if (cone->parent->NondegAssumed) {
00416         adjacent = dd_TRUE;
00417       }
00418       else{
00419         TempRay = ValidFirstRay;  /* the first ray for adjacency checking */
00420         while (TempRay != NULL && adjacent) {
00421           if (TempRay != Ray1 && TempRay != Ray2) {
00422             set_int(Face1, TempRay->ZeroSet, cone->AddedHalfspaces);
00423             if (set_subset(Face, Face1)) {
00424               if (localdebug) set_fwrite(stderr,Face1);
00425               adjacent = dd_FALSE;
00426             }
00427           }
00428           TempRay = TempRay->Next;
00429         }
00430       }
00431       if (adjacent){
00432         if (localdebug) fprintf(stderr,"The pair is adjacent and the pair must be stored for iteration %ld (row%ld)\n",
00433           fmin, cone->OrderVector[fmin]);
00434         NewEdge=(dd_AdjacencyPtr) malloc(sizeof *NewEdge);
00435         NewEdge->Ray1=Rmax;  /* save the one remains in iteration fmin in the first */
00436         NewEdge->Ray2=Rmin;  /* save the one deleted in iteration fmin in the second */
00437         NewEdge->Next=NULL;
00438         (cone->EdgeCount)++; 
00439         (cone->TotalEdgeCount)++;
00440         if (cone->Edges[fmin]==NULL){
00441           cone->Edges[fmin]=NewEdge;
00442           if (localdebug) fprintf(stderr,"Create a new edge list of %ld\n", fmin);
00443         }else{
00444           NewEdge->Next=cone->Edges[fmin];
00445           cone->Edges[fmin]=NewEdge;
00446         }
00447       }
00448     }
00449   }
00450 }
00451 
00452 void dd_CreateInitialEdges(dd_ConePtr cone)
00453 {
00454   dd_RayPtr Ptr1, Ptr2;
00455   dd_rowrange fii1,fii2;
00456   long count=0;
00457   dd_boolean adj,localdebug=dd_FALSE;
00458 
00459   cone->Iteration=cone->d;  /* CHECK */
00460   if (cone->FirstRay ==NULL || cone->LastRay==NULL){
00461     /* fprintf(stderr,"Warning: dd_ CreateInitialEdges called with NULL pointer(s)\n"); */
00462     goto _L99;
00463   }
00464   Ptr1=cone->FirstRay;
00465   while(Ptr1!=cone->LastRay && Ptr1!=NULL){
00466     fii1=Ptr1->FirstInfeasIndex;
00467     Ptr2=Ptr1->Next;
00468     while(Ptr2!=NULL){
00469       fii2=Ptr2->FirstInfeasIndex;
00470       count++;
00471       if (localdebug) fprintf(stderr,"dd_ CreateInitialEdges: edge %ld \n",count);
00472       dd_CheckAdjacency(cone, &Ptr1, &Ptr2, &adj);
00473       if (fii1!=fii2 && adj) 
00474         dd_ConditionalAddEdge(cone, Ptr1, Ptr2, cone->FirstRay);
00475       Ptr2=Ptr2->Next;
00476     }
00477     Ptr1=Ptr1->Next;
00478   }
00479 _L99:;  
00480 }
00481 
00482 
00483 void dd_UpdateEdges(dd_ConePtr cone, dd_RayPtr RRbegin, dd_RayPtr RRend)
00484 /* This procedure must be called after the ray list is sorted
00485    by dd_EvaluateARay2 so that FirstInfeasIndex's are monotonically
00486    increasing.
00487 */
00488 {
00489   dd_RayPtr Ptr1, Ptr2begin, Ptr2;
00490   dd_rowrange fii1;
00491   dd_boolean ptr2found,quit,localdebug=dd_FALSE;
00492   long count=0,pos1, pos2;
00493   float workleft,prevworkleft=110.0,totalpairs;
00494 
00495   totalpairs=(cone->ZeroRayCount-1.0)*(cone->ZeroRayCount-2.0)+1.0;
00496   Ptr2begin = NULL; 
00497   if (RRbegin ==NULL || RRend==NULL){
00498     if (1) fprintf(stderr,"Warning: dd_UpdateEdges called with NULL pointer(s)\n");
00499     goto _L99;
00500   }
00501   Ptr1=RRbegin;
00502   pos1=1;
00503   do{
00504     ptr2found=dd_FALSE;
00505     quit=dd_FALSE;
00506     fii1=Ptr1->FirstInfeasIndex;
00507     pos2=2;
00508     for (Ptr2=Ptr1->Next; !ptr2found && !quit; Ptr2=Ptr2->Next,pos2++){
00509       if  (Ptr2->FirstInfeasIndex > fii1){
00510         Ptr2begin=Ptr2;
00511         ptr2found=dd_TRUE;
00512       }
00513       else if (Ptr2==RRend) quit=dd_TRUE;
00514     }
00515     if (ptr2found){
00516       quit=dd_FALSE;
00517       for (Ptr2=Ptr2begin; !quit ; Ptr2=Ptr2->Next){
00518         count++;
00519         if (localdebug) fprintf(stderr,"dd_UpdateEdges: edge %ld \n",count);
00520         dd_ConditionalAddEdge(cone, Ptr1,Ptr2,RRbegin);
00521         if (Ptr2==RRend || Ptr2->Next==NULL) quit=dd_TRUE;
00522       }
00523     }
00524     Ptr1=Ptr1->Next;
00525     pos1++;
00526     workleft = 100.0 * (cone->ZeroRayCount-pos1) * (cone->ZeroRayCount - pos1-1.0) / totalpairs;
00527     if (cone->ZeroRayCount>=500 && dd_debug && pos1%10 ==0 && prevworkleft-workleft>=10 ) {
00528       fprintf(stderr,"*Work of iteration %5ld(/%ld): %4ld/%4ld => %4.1f%% left\n",
00529              cone->Iteration, cone->m, pos1, cone->ZeroRayCount, workleft);
00530       prevworkleft=workleft;
00531     }    
00532   }while(Ptr1!=RRend && Ptr1!=NULL);
00533 _L99:;  
00534 }
00535 
00536 void dd_FreeDDMemory0(dd_ConePtr cone)
00537 {
00538   dd_RayPtr Ptr, PrevPtr;
00539   long count;
00540   dd_colrange j;
00541   dd_boolean localdebug=dd_FALSE;
00542   
00543   /* THIS SHOULD BE REWRITTEN carefully */
00544   PrevPtr=cone->ArtificialRay;
00545   if (PrevPtr!=NULL){
00546     count=0;
00547     for (Ptr=cone->ArtificialRay->Next; Ptr!=NULL; Ptr=Ptr->Next){
00548       /* Added Marc Pfetsch 2/19/01 */
00549       for (j=0;j < cone->d;j++)
00550 dd_clear(PrevPtr->Ray[j]);
00551       dd_clear(PrevPtr->ARay);
00552 
00553       free(PrevPtr->Ray);
00554       free(PrevPtr->ZeroSet);
00555       free(PrevPtr);
00556       count++;
00557       PrevPtr=Ptr;
00558     };
00559     cone->FirstRay=NULL;
00560     /* Added Marc Pfetsch 010219 */
00561     for (j=0;j < cone->d;j++)
00562 dd_clear(cone->LastRay->Ray[j]);
00563     dd_clear(cone->LastRay->ARay);
00564 
00565     free(cone->LastRay->Ray);
00566     cone->LastRay->Ray = NULL;
00567     set_free(cone->LastRay->ZeroSet);
00568     cone->LastRay->ZeroSet = NULL;
00569     free(cone->LastRay);
00570     cone->LastRay = NULL;
00571     cone->ArtificialRay=NULL;
00572     if (localdebug) fprintf(stderr,"%ld ray storage spaces freed\n",count);
00573   }
00574 /* must add (by Sato) */
00575   free(cone->Edges);
00576   
00577   set_free(cone->GroundSet); 
00578   set_free(cone->EqualitySet); 
00579   set_free(cone->NonequalitySet); 
00580   set_free(cone->AddedHalfspaces); 
00581   set_free(cone->WeaklyAddedHalfspaces); 
00582   set_free(cone->InitialHalfspaces);
00583   free(cone->InitialRayIndex);
00584   free(cone->OrderVector);
00585   free(cone->newcol);
00586 
00587 /* Fixed by Shawn Rusaw.  Originally it was cone->d instead of cone->d_alloc */
00588   dd_FreeBmatrix(cone->d_alloc,cone->B);
00589   dd_FreeBmatrix(cone->d_alloc,cone->Bsave);
00590 
00591 /* Fixed by Marc Pfetsch 010219*/
00592   dd_FreeAmatrix(cone->m_alloc,cone->d_alloc,cone->A);
00593   cone->A = NULL;
00594 
00595   free(cone);
00596 }
00597 
00598 void dd_FreeDDMemory(dd_PolyhedraPtr poly)
00599 {
00600   dd_FreeDDMemory0(poly->child);
00601   poly->child=NULL;
00602 }
00603 
00604 void dd_FreePolyhedra(dd_PolyhedraPtr poly)
00605 {
00606   dd_bigrange i;
00607 
00608   if ((poly)->child != NULL) dd_FreeDDMemory(poly);
00609   dd_FreeAmatrix((poly)->m_alloc,poly->d_alloc, poly->A);
00610   dd_FreeArow((poly)->d_alloc,(poly)->c);
00611   free((poly)->EqualityIndex);
00612   if (poly->AincGenerated){
00613     for (i=1; i<=poly->m1; i++){
00614       set_free(poly->Ainc[i-1]);
00615     }
00616     free(poly->Ainc);
00617     set_free(poly->Ared);
00618     set_free(poly->Adom);
00619     poly->Ainc=NULL;
00620   }
00621 
00622   free(poly);
00623 }
00624 
00625 void dd_Normalize(dd_colrange d_size, mytype *V)
00626 {
00627   long j,jmin=0;
00628   mytype temp,min;
00629   dd_boolean nonzerofound=dd_FALSE;
00630 
00631   if (d_size>0){
00632     dd_init(min);  dd_init(temp);
00633     dd_abs(min,V[0]);  jmin=0; /* set the minmizer to 0 */
00634     if (dd_Positive(min)) nonzerofound=dd_TRUE;
00635     for (j = 1; j < d_size; j++) {
00636       dd_abs(temp,V[j]);
00637       if (dd_Positive(temp)){
00638         if (!nonzerofound || dd_Smaller(temp,min)){
00639           nonzerofound=dd_TRUE;
00640           dd_set(min, temp);  jmin=j;
00641         }
00642       }
00643     }
00644     if (dd_Positive(min)){
00645       for (j = 0; j < d_size; j++) dd_div(V[j], V[j], min);
00646     }
00647     dd_clear(min); dd_clear(temp);
00648   }
00649 }
00650 
00651 
00652 void dd_ZeroIndexSet(dd_rowrange m_size, dd_colrange d_size, dd_Amatrix A, mytype *x, dd_rowset ZS)
00653 {
00654   dd_rowrange i;
00655   mytype temp;
00656 
00657   /* Changed by Marc Pfetsch 010219 */
00658   dd_init(temp);
00659   set_emptyset(ZS);
00660   for (i = 1; i <= m_size; i++) {
00661     dd_AValue(&temp, d_size, A, x, i);
00662     if (dd_EqualToZero(temp)) set_addelem(ZS, i);
00663   }
00664 
00665   /* Changed by Marc Pfetsch 010219 */
00666   dd_clear(temp);
00667 }
00668 
00669 void dd_CopyBmatrix(dd_colrange d_size, dd_Bmatrix T, dd_Bmatrix TCOPY)
00670 {
00671   dd_rowrange i;
00672   dd_colrange j;
00673 
00674   for (i=0; i < d_size; i++) {
00675     for (j=0; j < d_size; j++) {
00676       dd_set(TCOPY[i][j],T[i][j]);
00677     }
00678   }
00679 }
00680 
00681 
00682 void dd_CopyArow(mytype *acopy, mytype *a, dd_colrange d)
00683 {
00684   dd_colrange j;
00685 
00686   for (j = 0; j < d; j++) {
00687     dd_set(acopy[j],a[j]);
00688   }
00689 }
00690 
00691 void dd_CopyNormalizedArow(mytype *acopy, mytype *a, dd_colrange d)
00692 {
00693   dd_CopyArow(acopy, a, d);
00694   dd_Normalize(d,acopy);
00695 }
00696 
00697 void dd_CopyAmatrix(mytype **Acopy, mytype **A, dd_rowrange m, dd_colrange d)
00698 {
00699   dd_rowrange i;
00700 
00701   for (i = 0; i< m; i++) {
00702     dd_CopyArow(Acopy[i],A[i],d);
00703   }
00704 }
00705 
00706 void dd_CopyNormalizedAmatrix(mytype **Acopy, mytype **A, dd_rowrange m, dd_colrange d)
00707 {
00708   dd_rowrange i;
00709 
00710   for (i = 0; i< m; i++) {
00711     dd_CopyNormalizedArow(Acopy[i],A[i],d);
00712   }
00713 }
00714 
00715 void dd_PermuteCopyAmatrix(mytype **Acopy, mytype **A, dd_rowrange m, dd_colrange d, dd_rowindex roworder)
00716 {
00717   dd_rowrange i;
00718 
00719   for (i = 1; i<= m; i++) {
00720     dd_CopyArow(Acopy[i-1],A[roworder[i]-1],d);
00721   }
00722 }
00723 
00724 void dd_PermutePartialCopyAmatrix(mytype **Acopy, mytype **A, dd_rowrange m, dd_colrange d, dd_rowindex roworder,dd_rowrange p, dd_rowrange q)
00725 {
00726  /* copy the rows of A whose roworder is positive.  roworder[i] is the row index of the copied row. */
00727   dd_rowrange i,k;
00728 
00729   k=0;
00730   for (i = 1; i<= m; i++) {
00731     if (roworder[i]>0) dd_CopyArow(Acopy[roworder[i]-1],A[i-1],d);
00732   }
00733 }
00734 
00735 void dd_InitializeArow(dd_colrange d,dd_Arow *a)
00736 {
00737   dd_colrange j;
00738 
00739   if (d>0) *a=(mytype*) calloc(d,sizeof(mytype));
00740   for (j = 0; j < d; j++) {
00741       dd_init((*a)[j]);
00742   }
00743 }
00744 
00745 void dd_InitializeAmatrix(dd_rowrange m,dd_colrange d,dd_Amatrix *A)
00746 {
00747   dd_rowrange i;
00748 
00749   if (m>0) (*A)=(mytype**) calloc(m,sizeof(mytype*));
00750   for (i = 0; i < m; i++) {
00751     dd_InitializeArow(d,&((*A)[i]));
00752   }
00753 }
00754 
00755 void dd_FreeAmatrix(dd_rowrange m,dd_colrange d,dd_Amatrix A)
00756 {
00757   dd_rowrange i;
00758   dd_colrange j;
00759 
00760   for (i = 0; i < m; i++) {
00761     for (j = 0; j < d; j++) {
00762       dd_clear(A[i][j]);
00763     }
00764   }
00765   if (A!=NULL) {
00766     for (i = 0; i < m; i++) {
00767       free(A[i]);
00768     }
00769     free(A);
00770   }
00771 }
00772 
00773 void dd_FreeArow(dd_colrange d, dd_Arow a)
00774 {
00775   dd_colrange j;
00776 
00777   for (j = 0; j < d; j++) {
00778     dd_clear(a[j]);
00779   }
00780   free(a);
00781 }
00782 
00783 
00784 void dd_InitializeBmatrix(dd_colrange d,dd_Bmatrix *B)
00785 {
00786   dd_colrange i,j;
00787 
00788   (*B)=(mytype**) calloc(d,sizeof(mytype*));
00789   for (j = 0; j < d; j++) {
00790     (*B)[j]=(mytype*) calloc(d,sizeof(mytype));
00791   }
00792   for (i = 0; i < d; i++) {
00793     for (j = 0; j < d; j++) {
00794       dd_init((*B)[i][j]);
00795     }
00796   }
00797 }
00798 
00799 void dd_FreeBmatrix(dd_colrange d,dd_Bmatrix B)
00800 {
00801   dd_colrange i,j;
00802 
00803   for (i = 0; i < d; i++) {
00804     for (j = 0; j < d; j++) {
00805       dd_clear(B[i][j]);
00806     }
00807   }
00808   if (B!=NULL) {
00809     for (j = 0; j < d; j++) {
00810       free(B[j]);
00811     }
00812     free(B);
00813   }
00814 }
00815 
00816 dd_SetFamilyPtr dd_CreateSetFamily(dd_bigrange fsize, dd_bigrange ssize)
00817 {
00818   dd_SetFamilyPtr F;
00819   dd_bigrange i,f0,f1,s0,s1;
00820 
00821   if (fsize<=0) {
00822     f0=0; f1=1;  
00823     /* if fsize<=0, the fsize is set to zero and the created size is one */
00824   } else {
00825     f0=fsize; f1=fsize;
00826   }
00827   if (ssize<=0) {
00828     s0=0; s1=1;  
00829     /* if ssize<=0, the ssize is set to zero and the created size is one */
00830   } else {
00831     s0=ssize; s1=ssize;
00832   }
00833 
00834   F=(dd_SetFamilyPtr) malloc (sizeof(dd_SetFamilyType));
00835   F->set=(set_type*) calloc(f1,sizeof(set_type));
00836   for (i=0; i<f1; i++) {
00837     set_initialize(&(F->set[i]), s1);
00838   }
00839   F->famsize=f0;
00840   F->setsize=s0;
00841   return F;
00842 }
00843 
00844 
00845 void dd_FreeSetFamily(dd_SetFamilyPtr F)
00846 {
00847   dd_bigrange i,f1;
00848 
00849   if (F!=NULL){
00850     if (F->famsize<=0) f1=1; else f1=F->famsize; 
00851       /* the smallest created size is one */
00852     for (i=0; i<f1; i++) {
00853       set_free(F->set[i]);
00854     }
00855     free(F->set);
00856     free(F);
00857   }
00858 }
00859 
00860 dd_MatrixPtr dd_CreateMatrix(dd_rowrange m_size,dd_colrange d_size)
00861 {
00862   dd_MatrixPtr M;
00863   dd_rowrange m0,m1;
00864   dd_colrange d0,d1;
00865 
00866   if (m_size<=0){ 
00867     m0=0; m1=1;  
00868     /* if m_size <=0, the number of rows is set to zero, the actual size is 1 */
00869   } else {
00870     m0=m_size; m1=m_size;
00871   }
00872   if (d_size<=0){ 
00873     d0=0; d1=1;  
00874     /* if d_size <=0, the number of cols is set to zero, the actual size is 1 */
00875   } else {
00876     d0=d_size; d1=d_size;
00877   }
00878   M=(dd_MatrixPtr) malloc (sizeof(dd_MatrixType));
00879   dd_InitializeAmatrix(m1,d1,&(M->matrix));
00880   dd_InitializeArow(d1,&(M->rowvec));
00881   M->rowsize=m0;
00882   set_initialize(&(M->linset), m1);
00883   M->colsize=d0;
00884   M->objective=dd_LPnone;
00885   M->numbtype=dd_Unknown;
00886   M->representation=dd_Unspecified;
00887   return M;
00888 }
00889 
00890 void dd_FreeMatrix(dd_MatrixPtr M)
00891 {
00892   dd_rowrange m1;
00893   dd_colrange d1;
00894 
00895   if (M!=NULL) {
00896     if (M->rowsize<=0) m1=1; else m1=M->rowsize;
00897     if (M->colsize<=0) d1=1; else d1=M->colsize;
00898     if (M!=NULL) {
00899       dd_FreeAmatrix(m1,d1,M->matrix);
00900       dd_FreeArow(d1,M->rowvec);
00901       set_free(M->linset);
00902       free(M);
00903     }
00904   }
00905 }
00906 
00907 void dd_SetToIdentity(dd_colrange d_size, dd_Bmatrix T)
00908 {
00909   dd_colrange j1, j2;
00910 
00911   for (j1 = 1; j1 <= d_size; j1++) {
00912     for (j2 = 1; j2 <= d_size; j2++) {
00913       if (j1 == j2)
00914         dd_set(T[j1 - 1][j2 - 1],dd_one);
00915       else
00916         dd_set(T[j1 - 1][j2 - 1],dd_purezero);
00917     }
00918   }
00919 }
00920 
00921 void dd_ColumnReduce(dd_ConePtr cone)
00922 {
00923   dd_colrange j,j1=0;
00924   dd_rowrange i;
00925   dd_boolean localdebug=dd_FALSE;
00926 
00927   for (j=1;j<=cone->d;j++) {
00928     if (cone->InitialRayIndex[j]>0){
00929       j1=j1+1;
00930       if (j1<j) {
00931         for (i=1; i<=cone->m; i++) dd_set(cone->A[i-1][j1-1],cone->A[i-1][j-1]);
00932         cone->newcol[j]=j1;
00933         if (localdebug){
00934           fprintf(stderr,"shifting the column %ld to column %ld\n", j, j1);
00935         }
00936           /* shifting forward */
00937       }
00938     } else {
00939       cone->newcol[j]=0;
00940       if (localdebug) {
00941         fprintf(stderr,"a generator (or an equation) of the linearity space: ");
00942         for (i=1; i<=cone->d; i++) dd_WriteNumber(stderr, cone->B[i-1][j-1]);
00943         fprintf(stderr,"\n");
00944       }
00945     }
00946   }
00947   cone->d=j1;  /* update the dimension. cone->d_orig remembers the old. */
00948   dd_CopyBmatrix(cone->d_orig, cone->B, cone->Bsave);  
00949     /* save the dual basis inverse as Bsave.  This matrix contains the linearity space generators. */
00950   cone->ColReduced=dd_TRUE;
00951 }
00952 
00953 long dd_MatrixRank(dd_MatrixPtr M, dd_rowset ignoredrows, dd_colset ignoredcols, dd_rowset *rowbasis, dd_colset *colbasis)
00954 {
00955   dd_boolean stop, chosen, localdebug=dd_debug;
00956   dd_rowset NopivotRow,PriorityRow;
00957   dd_colset ColSelected;
00958   dd_Bmatrix B;
00959   dd_rowindex roworder;
00960   dd_rowrange r;
00961   dd_colrange s;
00962   long rank;
00963 
00964   rank = 0;
00965   stop = dd_FALSE;
00966   set_initialize(&ColSelected, M->colsize);
00967   set_initialize(&NopivotRow, M->rowsize);
00968   set_initialize(rowbasis, M->rowsize);
00969   set_initialize(colbasis, M->colsize);
00970   set_initialize(&PriorityRow, M->rowsize);
00971   set_copy(NopivotRow,ignoredrows);
00972   set_copy(ColSelected,ignoredcols);
00973   dd_InitializeBmatrix(M->colsize, &B);
00974   dd_SetToIdentity(M->colsize, B);
00975   roworder=(long *)calloc(M->rowsize+1,sizeof(long*));
00976   for (r=0; r<M->rowsize; r++){roworder[r+1]=r+1;
00977   }
00978 
00979   do {   /* Find a set of rows for a basis */
00980       dd_SelectPivot2(M->rowsize, M->colsize,M->matrix,B,dd_MinIndex,roworder,
00981        PriorityRow,M->rowsize, NopivotRow, ColSelected, &r, &s, &chosen);
00982       if (dd_debug && chosen) 
00983         fprintf(stderr,"Procedure dd_MatrixRank: pivot on (r,s) =(%ld, %ld).\n", r, s);
00984       if (chosen) {
00985         set_addelem(NopivotRow, r);
00986         set_addelem(*rowbasis, r);
00987         set_addelem(ColSelected, s);
00988         set_addelem(*colbasis, s);
00989         rank++;
00990         dd_GaussianColumnPivot(M->rowsize, M->colsize, M->matrix, B, r, s);
00991         if (localdebug) dd_WriteBmatrix(stderr,M->colsize,B);
00992       } else {
00993         stop=dd_TRUE;
00994       }
00995       if (rank==M->colsize) stop = dd_TRUE;
00996   } while (!stop);
00997   dd_FreeBmatrix(M->colsize,B);
00998   free(roworder);
00999   set_free(ColSelected);
01000   set_free(NopivotRow);
01001   set_free(PriorityRow);
01002   return rank;
01003 }
01004 
01005 
01006 void dd_FindBasis(dd_ConePtr cone, long *rank)
01007 {
01008   dd_boolean stop, chosen, localdebug=dd_debug;
01009   dd_rowset NopivotRow;
01010   dd_colset ColSelected;
01011   dd_rowrange r;
01012   dd_colrange j,s;
01013 
01014   *rank = 0;
01015   stop = dd_FALSE;
01016   for (j=0;j<=cone->d;j++) cone->InitialRayIndex[j]=0;
01017   set_emptyset(cone->InitialHalfspaces);
01018   set_initialize(&ColSelected, cone->d);
01019   set_initialize(&NopivotRow, cone->m);
01020   set_copy(NopivotRow,cone->NonequalitySet);
01021   dd_SetToIdentity(cone->d, cone->B);
01022   do {   /* Find a set of rows for a basis */
01023       dd_SelectPivot2(cone->m, cone->d,cone->A,cone->B,cone->HalfspaceOrder,cone->OrderVector,
01024        cone->EqualitySet,cone->m, NopivotRow, ColSelected, &r, &s, &chosen);
01025       if (dd_debug && chosen) 
01026         fprintf(stderr,"Procedure dd_FindBasis: pivot on (r,s) =(%ld, %ld).\n", r, s);
01027       if (chosen) {
01028         set_addelem(cone->InitialHalfspaces, r);
01029         set_addelem(NopivotRow, r);
01030         set_addelem(ColSelected, s);
01031         cone->InitialRayIndex[s]=r;    /* cone->InitialRayIndex[s] stores the corr. row index */
01032         (*rank)++;
01033         dd_GaussianColumnPivot(cone->m, cone->d, cone->A, cone->B, r, s);
01034         if (localdebug) dd_WriteBmatrix(stderr,cone->d,cone->B);
01035       } else {
01036         stop=dd_TRUE;
01037       }
01038       if (*rank==cone->d) stop = dd_TRUE;
01039   } while (!stop);
01040   set_free(ColSelected);
01041   set_free(NopivotRow);
01042 }
01043 
01044 
01045 void dd_FindInitialRays(dd_ConePtr cone, dd_boolean *found)
01046 {
01047   dd_rowset CandidateRows;
01048   dd_rowrange i;
01049   long rank;
01050   dd_RowOrderType roworder_save=dd_LexMin;
01051 
01052   *found = dd_FALSE;
01053   set_initialize(&CandidateRows, cone->m);
01054   if (cone->parent->InitBasisAtBottom==dd_TRUE) {
01055     roworder_save=cone->HalfspaceOrder;
01056     cone->HalfspaceOrder=dd_MaxIndex;
01057     cone->PreOrderedRun=dd_FALSE;
01058   }
01059   else cone->PreOrderedRun=dd_TRUE;
01060   if (dd_debug) dd_WriteBmatrix(stderr, cone->d, cone->B);
01061   for (i = 1; i <= cone->m; i++)
01062     if (!set_member(i,cone->NonequalitySet)) set_addelem(CandidateRows, i);
01063     /*all rows not in NonequalitySet are candidates for initial cone*/
01064   dd_FindBasis(cone, &rank);
01065   if (dd_debug) dd_WriteBmatrix(stderr, cone->d, cone->B);
01066   if (dd_debug) fprintf(stderr,"dd_FindInitialRays: rank of Amatrix = %ld\n", rank);
01067   cone->LinearityDim=cone->d - rank;
01068   if (dd_debug) fprintf(stderr,"Linearity Dimension = %ld\n", cone->LinearityDim);
01069   if (cone->LinearityDim > 0) {
01070      dd_ColumnReduce(cone);
01071      dd_FindBasis(cone, &rank);
01072   }
01073   if (!set_subset(cone->EqualitySet,cone->InitialHalfspaces)) {
01074     if (dd_debug) {
01075       fprintf(stderr,"Equality set is dependent. Equality Set and an initial basis:\n");
01076       set_fwrite(stderr,cone->EqualitySet);
01077       set_fwrite(stderr,cone->InitialHalfspaces);
01078     };
01079   }
01080   *found = dd_TRUE;
01081   set_free(CandidateRows);
01082   if (cone->parent->InitBasisAtBottom==dd_TRUE) {
01083     cone->HalfspaceOrder=roworder_save;
01084   }
01085   if (cone->HalfspaceOrder==dd_MaxCutoff||
01086       cone->HalfspaceOrder==dd_MinCutoff||
01087       cone->HalfspaceOrder==dd_MixCutoff){
01088     cone->PreOrderedRun=dd_FALSE;
01089   } else cone->PreOrderedRun=dd_TRUE;
01090 }
01091 
01092 void dd_CheckEquality(dd_colrange d_size, dd_RayPtr*RP1, dd_RayPtr*RP2, dd_boolean *equal)
01093 {
01094   long j;
01095 
01096   if (dd_debug)
01097     fprintf(stderr,"Check equality of two rays\n");
01098   *equal = dd_TRUE;
01099   j = 1;
01100   while (j <= d_size && *equal) {
01101     if (!dd_Equal((*RP1)->Ray[j - 1],(*RP2)->Ray[j - 1]))
01102       *equal = dd_FALSE;
01103     j++;
01104   }
01105   if (*equal)
01106     fprintf(stderr,"Equal records found !!!!\n");
01107 }
01108 
01109 void dd_CreateNewRay(dd_ConePtr cone, 
01110     dd_RayPtr Ptr1, dd_RayPtr Ptr2, dd_rowrange ii)
01111 {
01112   /*Create a new ray by taking a linear combination of two rays*/
01113   dd_colrange j;
01114   mytype a1, a2, v1, v2;
01115   static dd_Arow NewRay;
01116   static dd_colrange last_d=0;
01117   dd_boolean localdebug=dd_debug;
01118 
01119   dd_init(a1); dd_init(a2); dd_init(v1); dd_init(v2);
01120   if (last_d!=cone->d){
01121     if (last_d>0) {
01122       for (j=0; j<last_d; j++) dd_clear(NewRay[j]);
01123       free(NewRay);
01124     }
01125     NewRay=(mytype*)calloc(cone->d,sizeof(mytype));
01126     for (j=0; j<cone->d; j++) dd_init(NewRay[j]);
01127     last_d=cone->d;
01128   }
01129 
01130   dd_AValue(&a1, cone->d, cone->A, Ptr1->Ray, ii);
01131   dd_AValue(&a2, cone->d, cone->A, Ptr2->Ray, ii);
01132   if (localdebug) {
01133     fprintf(stderr,"CreatNewRay: Ray1 ="); dd_WriteArow(stderr, Ptr1->Ray, cone->d);
01134     fprintf(stderr,"CreatNewRay: Ray2 ="); dd_WriteArow(stderr, Ptr2->Ray, cone->d);
01135   }
01136   dd_abs(v1,a1);
01137   dd_abs(v2,a2);
01138   if (localdebug){
01139     fprintf(stderr,"dd_AValue1 and ABS");  dd_WriteNumber(stderr,a1); dd_WriteNumber(stderr,v1); fprintf(stderr,"\n");
01140     fprintf(stderr,"dd_AValue2 and ABS");  dd_WriteNumber(stderr,a2); dd_WriteNumber(stderr,v2); fprintf(stderr,"\n");
01141   }
01142   for (j = 0; j < cone->d; j++){
01143     dd_LinearComb(NewRay[j], Ptr1->Ray[j],v2,Ptr2->Ray[j],v1);
01144   }
01145   if (localdebug) {
01146     fprintf(stderr,"CreatNewRay: New ray ="); dd_WriteArow(stderr, NewRay, cone->d);
01147   }
01148   dd_Normalize(cone->d, NewRay);
01149   if (localdebug) {
01150     fprintf(stderr,"CreatNewRay: dd_Normalized ray ="); dd_WriteArow(stderr, NewRay, cone->d);
01151   }
01152   dd_AddRay(cone, NewRay);
01153   dd_clear(a1); dd_clear(a2); dd_clear(v1); dd_clear(v2);
01154 }
01155 
01156 void dd_EvaluateARay1(dd_rowrange i, dd_ConePtr cone)
01157 /* Evaluate the ith component of the vector  A x RD.Ray 
01158     and rearrange the linked list so that
01159     the infeasible rays with respect to  i  will be
01160     placed consecutively from First 
01161  */
01162 {
01163   dd_colrange j;
01164   mytype temp,tnext;
01165   dd_RayPtr Ptr, PrevPtr, TempPtr;
01166 
01167   dd_init(temp); dd_init(tnext);
01168   Ptr = cone->FirstRay;
01169   PrevPtr = cone->ArtificialRay;
01170   if (PrevPtr->Next != Ptr) {
01171     fprintf(stderr,"Error.  Artificial Ray does not point to FirstRay!!!\n");
01172   }
01173   while (Ptr != NULL) {
01174     dd_set(temp,dd_purezero);
01175     for (j = 0; j < cone->d; j++){
01176       dd_mul(tnext,cone->A[i - 1][j],Ptr->Ray[j]);
01177       dd_add(temp,temp,tnext);
01178     }
01179     dd_set(Ptr->ARay,temp);
01180 /*    if ( temp <= -zero && Ptr != cone->FirstRay) {*/
01181     if ( dd_Negative(temp) && Ptr != cone->FirstRay) {
01182       /* fprintf(stderr,"Moving an infeasible record w.r.t. %ld to FirstRay\n",i); */
01183       if (Ptr==cone->LastRay) cone->LastRay=PrevPtr;
01184       TempPtr=Ptr;
01185       Ptr = Ptr->Next;
01186       PrevPtr->Next = Ptr;
01187       cone->ArtificialRay->Next = TempPtr;
01188       TempPtr->Next = cone->FirstRay;
01189       cone->FirstRay = TempPtr;
01190     }
01191     else {
01192       PrevPtr = Ptr;
01193       Ptr = Ptr->Next;
01194     }
01195   }
01196   dd_clear(temp); dd_clear(tnext);
01197 }
01198 
01199 void dd_EvaluateARay2(dd_rowrange i, dd_ConePtr cone)
01200 /* Evaluate the ith component of the vector  A x RD.Ray 
01201    and rearrange the linked list so that
01202    the infeasible rays with respect to  i  will be
01203    placed consecutively from First. Also for all feasible rays,
01204    "positive" rays and "zero" rays will be placed consecutively.
01205  */
01206 {
01207   dd_colrange j;
01208   mytype temp,tnext;
01209   dd_RayPtr Ptr, NextPtr;
01210   dd_boolean zerofound=dd_FALSE,negfound=dd_FALSE,posfound=dd_FALSE;
01211 
01212   if (cone==NULL || cone->TotalRayCount<=0) goto _L99;  
01213   dd_init(temp); dd_init(tnext);
01214   cone->PosHead=NULL;cone->ZeroHead=NULL;cone->NegHead=NULL;
01215   cone->PosLast=NULL;cone->ZeroLast=NULL;cone->NegLast=NULL;
01216   Ptr = cone->FirstRay;
01217   while (Ptr != NULL) {
01218     NextPtr=Ptr->Next;  /* remember the Next record */
01219     Ptr->Next=NULL;     /* then clear the Next pointer */
01220     dd_set(temp,dd_purezero);
01221     for (j = 0; j < cone->d; j++){
01222       dd_mul(tnext,cone->A[i - 1][j],Ptr->Ray[j]);
01223       dd_add(temp,temp,tnext);
01224     }
01225     dd_set(Ptr->ARay,temp);
01226 /*    if ( temp < -zero) {*/
01227     if ( dd_Negative(temp)) {
01228       if (!negfound){
01229         negfound=dd_TRUE;
01230         cone->NegHead=Ptr;
01231         cone->NegLast=Ptr;
01232       }
01233       else{
01234         Ptr->Next=cone->NegHead;
01235         cone->NegHead=Ptr;
01236       }
01237     }
01238 /*    else if (temp > zero){*/
01239     else if (dd_Positive(temp)){
01240       if (!posfound){
01241         posfound=dd_TRUE;
01242         cone->PosHead=Ptr;
01243         cone->PosLast=Ptr;
01244       }
01245       else{  
01246         Ptr->Next=cone->PosHead;
01247         cone->PosHead=Ptr;
01248        }
01249     }
01250     else {
01251       if (!zerofound){
01252         zerofound=dd_TRUE;
01253         cone->ZeroHead=Ptr;
01254         cone->ZeroLast=Ptr;
01255       }
01256       else{
01257         Ptr->Next=cone->ZeroHead;
01258         cone->ZeroHead=Ptr;
01259       }
01260     }
01261     Ptr=NextPtr;
01262   }
01263   /* joining three neg, pos and zero lists */
01264   if (negfound){                 /* -list nonempty */
01265     cone->FirstRay=cone->NegHead;
01266     if (posfound){               /* -list & +list nonempty */
01267       cone->NegLast->Next=cone->PosHead;
01268       if (zerofound){            /* -list, +list, 0list all nonempty */
01269         cone->PosLast->Next=cone->ZeroHead;
01270         cone->LastRay=cone->ZeroLast;
01271       } 
01272       else{                      /* -list, +list nonempty but  0list empty */
01273         cone->LastRay=cone->PosLast;      
01274       }
01275     }
01276     else{                        /* -list nonempty & +list empty */
01277       if (zerofound){            /* -list,0list nonempty & +list empty */
01278         cone->NegLast->Next=cone->ZeroHead;
01279         cone->LastRay=cone->ZeroLast;
01280       } 
01281       else {                      /* -list nonempty & +list,0list empty */
01282         cone->LastRay=cone->NegLast;
01283       }
01284     }
01285   }
01286   else if (posfound){            /* -list empty & +list nonempty */
01287     cone->FirstRay=cone->PosHead;
01288     if (zerofound){              /* -list empty & +list,0list nonempty */
01289       cone->PosLast->Next=cone->ZeroHead;
01290       cone->LastRay=cone->ZeroLast;
01291     } 
01292     else{                        /* -list,0list empty & +list nonempty */
01293       cone->LastRay=cone->PosLast;
01294     }
01295   }
01296   else{                          /* -list,+list empty & 0list nonempty */
01297     cone->FirstRay=cone->ZeroHead;
01298     cone->LastRay=cone->ZeroLast;
01299   }
01300   cone->ArtificialRay->Next=cone->FirstRay;
01301   cone->LastRay->Next=NULL;
01302   dd_clear(temp); dd_clear(tnext);
01303   _L99:;  
01304 }
01305 
01306 void dd_DeleteNegativeRays(dd_ConePtr cone)
01307 /* Eliminate the infeasible rays with respect to  i  which
01308    are supposed to be consecutive from the head of the dd_Ray list,
01309    and sort the zero list assumed to be consecutive at the
01310    end of the list.
01311  */
01312 {
01313   dd_rowrange fii,fiitest;
01314   mytype temp;
01315   dd_RayPtr Ptr, PrevPtr,NextPtr,ZeroPtr1,ZeroPtr0;
01316   dd_boolean found, completed, zerofound=dd_FALSE,negfound=dd_FALSE,posfound=dd_FALSE;
01317   dd_boolean localdebug=dd_FALSE;
01318   
01319   dd_init(temp);
01320   cone->PosHead=NULL;cone->ZeroHead=NULL;cone->NegHead=NULL;
01321   cone->PosLast=NULL;cone->ZeroLast=NULL;cone->NegLast=NULL;
01322 
01323   /* Delete the infeasible rays  */
01324   PrevPtr= cone->ArtificialRay;
01325   Ptr = cone->FirstRay;
01326   if (PrevPtr->Next != Ptr) 
01327     fprintf(stderr,"Error at dd_DeleteNegativeRays: ArtificialRay does not point the FirstRay.\n");
01328   completed=dd_FALSE;
01329   while (Ptr != NULL && !completed){
01330 /*    if ( (Ptr->ARay) < -zero ){ */
01331     if ( dd_Negative(Ptr->ARay)){
01332       dd_Eliminate(cone, &PrevPtr);
01333       Ptr=PrevPtr->Next;
01334     }
01335     else{
01336       completed=dd_TRUE;
01337     }
01338   }
01339   
01340   /* Sort the zero rays */
01341   Ptr = cone->FirstRay;
01342   cone->ZeroRayCount=0;
01343   while (Ptr != NULL) {
01344     NextPtr=Ptr->Next;  /* remember the Next record */
01345     dd_set(temp,Ptr->ARay);
01346     if (localdebug) {fprintf(stderr,"Ptr->ARay :"); dd_WriteNumber(stderr, temp);}
01347 /*    if ( temp < -zero) {*/
01348     if ( dd_Negative(temp)) {
01349       if (!negfound){
01350         fprintf(stderr,"Error: An infeasible ray found after their removal\n");
01351         negfound=dd_TRUE;
01352       }
01353     }
01354 /*    else if (temp > zero){*/
01355     else if (dd_Positive(temp)){
01356       if (!posfound){
01357         posfound=dd_TRUE;
01358         cone->PosHead=Ptr;
01359         cone->PosLast=Ptr;
01360       }
01361       else{  
01362         cone->PosLast=Ptr;
01363        }
01364     }
01365     else {
01366       (cone->ZeroRayCount)++;
01367       if (!zerofound){
01368         zerofound=dd_TRUE;
01369         cone->ZeroHead=Ptr;
01370         cone->ZeroLast=Ptr;
01371         cone->ZeroLast->Next=NULL;
01372       }
01373       else{/* Find a right position to store the record sorted w.r.t. FirstInfeasIndex */
01374         fii=Ptr->FirstInfeasIndex; 
01375         found=dd_FALSE;
01376         ZeroPtr1=NULL;
01377         for (ZeroPtr0=cone->ZeroHead; !found && ZeroPtr0!=NULL ; ZeroPtr0=ZeroPtr0->Next){
01378           fiitest=ZeroPtr0->FirstInfeasIndex;
01379           if (fiitest >= fii){
01380             found=dd_TRUE;
01381           }
01382           else ZeroPtr1=ZeroPtr0;
01383         }
01384         /* fprintf(stderr,"insert position found \n %d  index %ld\n",found, fiitest); */
01385         if (!found){           /* the new record must be stored at the end of list */
01386           cone->ZeroLast->Next=Ptr;
01387           cone->ZeroLast=Ptr;
01388           cone->ZeroLast->Next=NULL;
01389         }
01390         else{
01391           if (ZeroPtr1==NULL){ /* store the new one at the head, and update the head ptr */
01392             /* fprintf(stderr,"Insert at the head\n"); */
01393             Ptr->Next=cone->ZeroHead;
01394             cone->ZeroHead=Ptr;
01395           }
01396           else{                /* store the new one inbetween ZeroPtr1 and 0 */
01397             /* fprintf(stderr,"Insert inbetween\n");  */
01398             Ptr->Next=ZeroPtr1->Next;
01399             ZeroPtr1->Next=Ptr;
01400           }
01401         }
01402         /*
01403         Ptr->Next=cone->ZeroHead;
01404         cone->ZeroHead=Ptr;
01405         */
01406       }
01407     }
01408     Ptr=NextPtr;
01409   }
01410   /* joining the pos and zero lists */
01411   if (posfound){            /* -list empty & +list nonempty */
01412     cone->FirstRay=cone->PosHead;
01413     if (zerofound){              /* +list,0list nonempty */
01414       cone->PosLast->Next=cone->ZeroHead;
01415       cone->LastRay=cone->ZeroLast;
01416     } 
01417     else{                        /* 0list empty & +list nonempty */
01418       cone->LastRay=cone->PosLast;
01419     }
01420   }
01421   else{                          /* +list empty & 0list nonempty */
01422     cone->FirstRay=cone->ZeroHead;
01423     cone->LastRay=cone->ZeroLast;
01424   }
01425   cone->ArtificialRay->Next=cone->FirstRay;
01426   cone->LastRay->Next=NULL;
01427   dd_clear(temp);
01428 }
01429 
01430 void dd_FeasibilityIndices(long *fnum, long *infnum, dd_rowrange i, dd_ConePtr cone)
01431 {
01432   /*Evaluate the number of feasible rays and infeasible rays*/
01433   /*  w.r.t the hyperplane  i*/
01434   dd_colrange j;
01435   mytype temp, tnext;
01436   dd_RayPtr Ptr;
01437 
01438   dd_init(temp); dd_init(tnext);
01439   *fnum = 0;
01440   *infnum = 0;
01441   Ptr = cone->FirstRay;
01442   while (Ptr != NULL) {
01443     dd_set(temp,dd_purezero);
01444     for (j = 0; j < cone->d; j++){
01445       dd_mul(tnext, cone->A[i - 1][j],Ptr->Ray[j]);
01446       dd_add(temp, temp, tnext);
01447     }
01448     if (dd_Nonnegative(temp))
01449       (*fnum)++;
01450     else
01451       (*infnum)++;
01452     Ptr = Ptr->Next;
01453   }
01454   dd_clear(temp); dd_clear(tnext);
01455 }
01456 
01457 dd_boolean dd_LexSmaller(mytype *v1, mytype *v2, long dmax)
01458 { /* dmax is the size of vectors v1,v2 */
01459   dd_boolean determined, smaller;
01460   dd_colrange j;
01461 
01462   smaller = dd_FALSE;
01463   determined = dd_FALSE;
01464   j = 1;
01465   do {
01466     if (!dd_Equal(v1[j - 1],v2[j - 1])) {  /* 086 */
01467       if (dd_Smaller(v1[j - 1],v2[j - 1])) {  /*086 */
01468             smaller = dd_TRUE;
01469           }
01470       determined = dd_TRUE;
01471     } else
01472       j++;
01473   } while (!(determined) && (j <= dmax));
01474   return smaller;
01475 }
01476 
01477 
01478 dd_boolean dd_LexLarger(mytype *v1, mytype *v2, long dmax)
01479 {
01480   return dd_LexSmaller(v2, v1, dmax);
01481 }
01482 
01483 dd_boolean dd_LexEqual(mytype *v1, mytype *v2, long dmax)
01484 { /* dmax is the size of vectors v1,v2 */
01485   dd_boolean determined, equal;
01486   dd_colrange j;
01487 
01488   equal = dd_TRUE;
01489   determined = dd_FALSE;
01490   j = 1;
01491   do {
01492     if (!dd_Equal(v1[j - 1],v2[j - 1])) {  /* 093c */
01493         equal = dd_FALSE;
01494         determined = dd_TRUE;
01495     } else {
01496       j++;
01497     }
01498   } while (!(determined) && (j <= dmax));
01499   return equal;
01500 }
01501 
01502 void dd_AddNewHalfspace1(dd_ConePtr cone, dd_rowrange hnew)
01503 /* This procedure 1 must be used with PreorderedRun=dd_FALSE 
01504    This procedure is the most elementary implementation of
01505    DD and can be used with any type of ordering, including
01506    dynamic ordering of rows, e.g. MaxCutoff, MinCutoff.
01507    The memory requirement is minimum because it does not
01508    store any adjacency among the rays.
01509 */
01510 {
01511   dd_RayPtr RayPtr0,RayPtr1,RayPtr2,RayPtr2s,RayPtr3;
01512   long pos1, pos2;
01513   double prevprogress, progress;
01514   mytype value1, value2;
01515   dd_boolean adj, equal, completed;
01516 
01517   dd_init(value1); dd_init(value2);
01518   dd_EvaluateARay1(hnew, cone);  
01519    /*Check feasibility of rays w.r.t. hnew 
01520      and put all infeasible ones consecutively */
01521 
01522   RayPtr0 = cone->ArtificialRay;   /*Pointer pointing RayPrt1*/
01523   RayPtr1 = cone->FirstRay;        /*1st hnew-infeasible ray to scan and compare with feasible rays*/
01524   dd_set(value1,cone->FirstRay->ARay);
01525   if (dd_Nonnegative(value1)) {
01526     if (cone->RayCount==cone->WeaklyFeasibleRayCount) cone->CompStatus=dd_AllFound;
01527     goto _L99;        /* Sicne there is no hnew-infeasible ray and nothing to do */
01528   }
01529   else {
01530     RayPtr2s = RayPtr1->Next;/* RayPtr2s must point the first feasible ray */
01531     pos2=1;
01532     while (RayPtr2s!=NULL && dd_Negative(RayPtr2s->ARay)) {
01533       RayPtr2s = RayPtr2s->Next;
01534       pos2++;
01535     }
01536   }
01537   if (RayPtr2s==NULL) {
01538     cone->FirstRay=NULL;
01539     cone->ArtificialRay->Next=cone->FirstRay;
01540     cone->RayCount=0;
01541     cone->CompStatus=dd_AllFound;
01542     goto _L99;   /* All rays are infeasible, and the computation must stop */
01543   }
01544   RayPtr2 = RayPtr2s;   /*2nd feasible ray to scan and compare with 1st*/
01545   RayPtr3 = cone->LastRay;    /*Last feasible for scanning*/
01546   prevprogress=-10.0;
01547   pos1 = 1;
01548   completed=dd_FALSE;
01549   while ((RayPtr1 != RayPtr2s) && !completed) {
01550     dd_set(value1,RayPtr1->ARay);
01551     dd_set(value2,RayPtr2->ARay);
01552     dd_CheckEquality(cone->d, &RayPtr1, &RayPtr2, &equal);
01553     if ((dd_Positive(value1) && dd_Negative(value2)) || (dd_Negative(value1) && dd_Positive(value2))){
01554       dd_CheckAdjacency(cone, &RayPtr1, &RayPtr2, &adj);
01555       if (adj) dd_CreateNewRay(cone, RayPtr1, RayPtr2, hnew);
01556     }
01557     if (RayPtr2 != RayPtr3) {
01558       RayPtr2 = RayPtr2->Next;
01559       continue;
01560     }
01561     if (dd_Negative(value1) || equal) {
01562       dd_Eliminate(cone, &RayPtr0);
01563       RayPtr1 = RayPtr0->Next;
01564       RayPtr2 = RayPtr2s;
01565     } else {
01566       completed=dd_TRUE;
01567     }
01568     pos1++;
01569     progress = 100.0 * ((double)pos1 / pos2) * (2.0 * pos2 - pos1) / pos2;
01570     if (progress-prevprogress>=10 && pos1%10==0 && dd_debug) {
01571       fprintf(stderr,"*Progress of iteration %5ld(/%ld):   %4ld/%4ld => %4.1f%% done\n",
01572              cone->Iteration, cone->m, pos1, pos2, progress);
01573       prevprogress=progress;
01574     }
01575   }
01576   if (cone->RayCount==cone->WeaklyFeasibleRayCount) cone->CompStatus=dd_AllFound;
01577   _L99:;
01578   dd_clear(value1); dd_clear(value2);
01579 }
01580 
01581 void dd_AddNewHalfspace2(dd_ConePtr cone, dd_rowrange hnew)
01582 /* This procedure must be used under PreOrderedRun mode */
01583 {
01584   dd_RayPtr RayPtr0,RayPtr1,RayPtr2;
01585   dd_AdjacencyType *EdgePtr, *EdgePtr0;
01586   long pos1;
01587   dd_rowrange fii1, fii2;
01588   dd_boolean localdebug=dd_FALSE;
01589 
01590   dd_EvaluateARay2(hnew, cone);
01591    /* Check feasibility of rays w.r.t. hnew 
01592       and sort them. ( -rays, +rays, 0rays)*/
01593 
01594   if (cone->PosHead==NULL && cone->ZeroHead==NULL) {
01595     cone->FirstRay=NULL;
01596     cone->ArtificialRay->Next=cone->FirstRay;
01597     cone->RayCount=0;
01598     cone->CompStatus=dd_AllFound;
01599     goto _L99;   /* All rays are infeasible, and the computation must stop */
01600   }
01601   
01602   if (localdebug){
01603     pos1=0;
01604     fprintf(stderr,"(pos, FirstInfeasIndex, A Ray)=\n");
01605     for (RayPtr0=cone->FirstRay; RayPtr0!=NULL; RayPtr0=RayPtr0->Next){
01606       pos1++;
01607       fprintf(stderr,"(%ld,%ld,",pos1,RayPtr0->FirstInfeasIndex);
01608       dd_WriteNumber(stderr,RayPtr0->ARay); 
01609       fprintf(stderr,") ");
01610    }
01611     fprintf(stderr,"\n");
01612   }
01613   
01614   if (cone->ZeroHead==NULL) cone->ZeroHead=cone->LastRay;
01615 
01616   EdgePtr=cone->Edges[cone->Iteration];
01617   while (EdgePtr!=NULL){
01618     RayPtr1=EdgePtr->Ray1;
01619     RayPtr2=EdgePtr->Ray2;
01620     fii1=RayPtr1->FirstInfeasIndex;   
01621     dd_CreateNewRay(cone, RayPtr1, RayPtr2, hnew);
01622     fii2=cone->LastRay->FirstInfeasIndex;
01623     if (fii1 != fii2) 
01624       dd_ConditionalAddEdge(cone,RayPtr1,cone->LastRay,cone->PosHead);
01625     EdgePtr0=EdgePtr;
01626     EdgePtr=EdgePtr->Next;
01627     free(EdgePtr0);
01628     (cone->EdgeCount)--;
01629   }
01630   cone->Edges[cone->Iteration]=NULL;
01631   
01632   dd_DeleteNegativeRays(cone);
01633     
01634   set_addelem(cone->AddedHalfspaces, hnew);
01635 
01636   if (cone->Iteration<cone->m){
01637     if (cone->ZeroHead!=NULL && cone->ZeroHead!=cone->LastRay){
01638       if (cone->ZeroRayCount>200 && dd_debug) fprintf(stderr,"*New edges being scanned...\n");
01639       dd_UpdateEdges(cone, cone->ZeroHead, cone->LastRay);
01640     }
01641   }
01642 
01643   if (cone->RayCount==cone->WeaklyFeasibleRayCount) cone->CompStatus=dd_AllFound;
01644 _L99:;
01645 }
01646 
01647 
01648 void dd_SelectNextHalfspace0(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
01649 {
01650   /*A natural way to choose the next hyperplane.  Simply the largest index*/
01651   long i;
01652   dd_boolean determined;
01653 
01654   i = cone->m;
01655   determined = dd_FALSE;
01656   do {
01657     if (set_member(i, excluded))
01658       i--;
01659     else
01660       determined = dd_TRUE;
01661   } while (!determined && i>=1);
01662   if (determined) 
01663     *hnext = i;
01664   else
01665     *hnext = 0;
01666 }
01667 
01668 void dd_SelectNextHalfspace1(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
01669 {
01670   /*Natural way to choose the next hyperplane.  Simply the least index*/
01671   long i;
01672   dd_boolean determined;
01673 
01674   i = 1;
01675   determined = dd_FALSE;
01676   do {
01677     if (set_member(i, excluded))
01678       i++;
01679     else
01680       determined = dd_TRUE;
01681   } while (!determined && i<=cone->m);
01682   if (determined) 
01683     *hnext = i;
01684   else 
01685     *hnext=0;
01686 }
01687 
01688 void dd_SelectNextHalfspace2(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
01689 {
01690   /*Choose the next hyperplane with maximum infeasibility*/
01691   long i, fea, inf, infmin, fi=0;   /*feasibility and infeasibility numbers*/
01692 
01693   infmin = cone->RayCount + 1;
01694   for (i = 1; i <= cone->m; i++) {
01695     if (!set_member(i, excluded)) {
01696       dd_FeasibilityIndices(&fea, &inf, i, cone);
01697       if (inf < infmin) {
01698         infmin = inf;
01699         fi = fea;
01700         *hnext = i;
01701       }
01702     }
01703   }
01704   if (dd_debug) {
01705     fprintf(stderr,"*infeasible rays (min) =%5ld, #feas rays =%5ld\n", infmin, fi);
01706   }
01707 }
01708 
01709 void dd_SelectNextHalfspace3(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
01710 {
01711   /*Choose the next hyperplane with maximum infeasibility*/
01712   long i, fea, inf, infmax, fi=0;   /*feasibility and infeasibility numbers*/
01713   dd_boolean localdebug=dd_debug;
01714 
01715   infmax = -1;
01716   for (i = 1; i <= cone->m; i++) {
01717     if (!set_member(i, excluded)) {
01718       dd_FeasibilityIndices(&fea, &inf, i, cone);
01719       if (inf > infmax) {
01720         infmax = inf;
01721         fi = fea;
01722         *hnext = i;
01723       }
01724     }
01725   }
01726   if (localdebug) {
01727     fprintf(stderr,"*infeasible rays (max) =%5ld, #feas rays =%5ld\n", infmax, fi);
01728   }
01729 }
01730 
01731 void dd_SelectNextHalfspace4(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
01732 {
01733   /*Choose the next hyperplane with the most unbalanced cut*/
01734   long i, fea, inf, max, tmax, fi=0, infi=0;
01735       /*feasibility and infeasibility numbers*/
01736 
01737   max = -1;
01738   for (i = 1; i <= cone->m; i++) {
01739     if (!set_member(i, excluded)) {
01740       dd_FeasibilityIndices(&fea, &inf, i, cone);
01741       if (fea <= inf)
01742         tmax = inf;
01743       else
01744         tmax = fea;
01745       if (tmax > max) {
01746         max = tmax;
01747         fi = fea;
01748         infi = inf;
01749         *hnext = i;
01750       }
01751     }
01752   }
01753   if (!dd_debug)
01754     return;
01755   if (max == fi) {
01756     fprintf(stderr,"*infeasible rays (min) =%5ld, #feas rays =%5ld\n", infi, fi);
01757   } else {
01758     fprintf(stderr,"*infeasible rays (max) =%5ld, #feas rays =%5ld\n", infi, fi);
01759   }
01760 }
01761 
01762 void dd_SelectNextHalfspace5(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
01763 {
01764   /*Choose the next hyperplane which is lexico-min*/
01765   long i, minindex;
01766   mytype *v1, *v2;
01767 
01768   minindex = 0;
01769   v1 = NULL;
01770   for (i = 1; i <= cone->m; i++) {
01771     if (!set_member(i, excluded)) {
01772           v2 = cone->A[i - 1];
01773       if (minindex == 0) {
01774             minindex = i;
01775             v1=v2;
01776       } else if (dd_LexSmaller(v2,v1,cone->d)) {
01777         minindex = i;
01778             v1=v2;
01779       }
01780     }
01781   }
01782   *hnext = minindex;
01783 }
01784 
01785 
01786 void dd_SelectNextHalfspace6(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
01787 {
01788   /*Choose the next hyperplane which is lexico-max*/
01789   long i, maxindex;
01790   mytype *v1, *v2;
01791 
01792   maxindex = 0;
01793   v1 = NULL;
01794   for (i = 1; i <= cone->m; i++) {
01795     if (!set_member(i, excluded)) {
01796       v2= cone->A[i - 1];
01797       if (maxindex == 0) {
01798         maxindex = i;
01799         v1=v2;
01800       } else if (dd_LexLarger(v2, v1, cone->d)) {
01801         maxindex = i;
01802         v1=v2;
01803      }
01804     }
01805   }
01806   *hnext = maxindex;
01807 }
01808 
01809 void dd_UniqueRows(dd_rowindex OV, long p, long r, dd_Amatrix A, long dmax, dd_rowset preferred, long *uniqrows)
01810 {
01811  /* Select a subset of rows of A (with range [p, q] up to dimension dmax) by removing duplicates.
01812     When a row subset preferred is nonempty, those row indices in the set have priority.  If
01813     two priority rows define the same row vector, one is chosen.
01814     For a selected unique row i, OV[i] returns a new position of the unique row i. 
01815     For other nonuniqu row i, OV[i] returns a negative of the original row j dominating i.
01816     Thus the original contents of OV[p..r] will be rewritten.  Other components remain the same.
01817     *uniqrows returns the number of unique rows.
01818 */
01819   long i,iuniq,j;
01820   mytype *x;
01821   dd_boolean localdebug=dd_FALSE;
01822   
01823   if (p<=0 || p > r) goto _L99;
01824   iuniq=p; j=1;  /* the first unique row candidate */
01825   x=A[p-1];
01826   OV[p]=j;  /* tentative row index of the candidate */
01827   for (i=p+1;i<=r; i++){
01828     if (!dd_LexEqual(x,A[i-1],dmax)) {
01829       /* a new row vector found. */
01830       iuniq=i;
01831       j=j+1;
01832       OV[i]=j;    /* Tentatively register the row i.  */
01833       x=A[i-1];
01834     } else {
01835       /* rows are equal */
01836       if (set_member(i, preferred) && !set_member(iuniq, preferred)){
01837         OV[iuniq]=-i;  /* the row iuniq is dominated by the row i */
01838         iuniq=i;  /* the row i is preferred.  Change the candidate. */
01839         OV[i]=j;  /* the row i is tentatively registered. */
01840         x=A[i-1];
01841       } else {
01842         OV[i]=-iuniq;  /* the row iuniq is dominated by the row i */
01843       }
01844     }
01845   }
01846   *uniqrows=j;
01847   if (localdebug){
01848     printf("The number of unique rows are %ld\n with order vector",*uniqrows);
01849     for (i=p;i<=r; i++){
01850       printf(" %ld:%ld ",i,OV[i]);
01851     }
01852     printf("\n");
01853   }
01854   _L99:;
01855 }
01856 
01857 long dd_Partition(dd_rowindex OV, long p, long r, dd_Amatrix A, long dmax)
01858 {
01859   mytype *x;
01860   long i,j,ovi;
01861   
01862   x=A[OV[p]-1];
01863 
01864   i=p-1;
01865   j=r+1;
01866   while (dd_TRUE){
01867     do{
01868       j--;
01869     } while (dd_LexLarger(A[OV[j]-1],x,dmax));
01870     do{
01871       i++;
01872     } while (dd_LexSmaller(A[OV[i]-1],x,dmax));
01873     if (i<j){
01874       ovi=OV[i];
01875       OV[i]=OV[j];
01876       OV[j]=ovi;
01877     }
01878     else{
01879       return j;
01880     }
01881   }
01882 }
01883 
01884 void dd_QuickSort(dd_rowindex OV, long p, long r, dd_Amatrix A, long dmax)
01885 {
01886   long q;
01887   
01888   if (p < r){
01889     q = dd_Partition(OV, p, r, A, dmax);
01890     dd_QuickSort(OV, p, q, A, dmax);
01891     dd_QuickSort(OV, q+1, r, A, dmax);
01892   }
01893 }
01894 
01895 
01896 #ifndef RAND_MAX 
01897 #define RAND_MAX 32767 
01898 #endif
01899 
01900 void dd_RandomPermutation(dd_rowindex OV, long t, unsigned int seed)
01901 {
01902   long k,j,ovj;
01903   double u,xk,r,rand_max=(double) RAND_MAX;
01904   dd_boolean localdebug=dd_FALSE;
01905 
01906   srand(seed);
01907   for (j=t; j>1 ; j--) {
01908     r=rand();
01909     u=r/rand_max;
01910     xk=(double)(j*u +1);
01911     k=(long)xk;
01912     if (localdebug) fprintf(stderr,"u=%g, k=%ld, r=%g, randmax= %g\n",u,k,r,rand_max);
01913     ovj=OV[j];
01914     OV[j]=OV[k];
01915     OV[k]=ovj;
01916     if (localdebug) fprintf(stderr,"row %ld is exchanged with %ld\n",j,k); 
01917   }
01918 }
01919 
01920 void dd_ComputeRowOrderVector(dd_ConePtr cone)
01921 {
01922   long i,itemp;
01923 
01924   cone->OrderVector[0]=0;
01925   switch (cone->HalfspaceOrder){
01926   case dd_MaxIndex:
01927     for(i=1; i<=cone->m; i++) cone->OrderVector[i]=cone->m-i+1;
01928     break;
01929 
01930   case dd_MinIndex: 
01931     for(i=1; i<=cone->m; i++) cone->OrderVector[i]=i;    
01932     break;
01933 
01934   case dd_LexMin: case dd_MinCutoff: case dd_MixCutoff: case dd_MaxCutoff:
01935     for(i=1; i<=cone->m; i++) cone->OrderVector[i]=i;
01936     dd_RandomPermutation(cone->OrderVector, cone->m, cone->rseed);
01937     dd_QuickSort(cone->OrderVector, 1, cone->m, cone->A, cone->d);
01938     break;
01939 
01940   case dd_LexMax:
01941     for(i=1; i<=cone->m; i++) cone->OrderVector[i]=i;
01942     dd_RandomPermutation(cone->OrderVector, cone->m, cone->rseed);
01943     dd_QuickSort(cone->OrderVector, 1, cone->m, cone->A, cone->d);
01944     for(i=1; i<=cone->m/2;i++){   /* just reverse the order */
01945       itemp=cone->OrderVector[i];
01946       cone->OrderVector[i]=cone->OrderVector[cone->m-i+1];
01947       cone->OrderVector[cone->m-i+1]=itemp;
01948     }
01949     break;
01950 
01951   case dd_RandomRow:
01952     for(i=1; i<=cone->m; i++) cone->OrderVector[i]=i;
01953     dd_RandomPermutation(cone->OrderVector, cone->m, cone->rseed);
01954     break;
01955 
01956   }
01957 }
01958 
01959 void dd_UpdateRowOrderVector(dd_ConePtr cone, dd_rowset PriorityRows)
01960 /* Update the RowOrder vector to shift selected rows
01961 in highest order.
01962 */
01963 {
01964   dd_rowrange i,j,k,j1=0,oj=0;
01965   long rr;
01966   dd_boolean found, localdebug=dd_FALSE;
01967   
01968   if (dd_debug) localdebug=dd_TRUE;
01969   found=dd_TRUE;
01970   rr=set_card(PriorityRows);
01971   if (localdebug) set_fwrite(stderr,PriorityRows);
01972   for (i=1; i<=rr; i++){
01973     found=dd_FALSE;
01974     for (j=i; j<=cone->m && !found; j++){
01975       oj=cone->OrderVector[j];
01976       if (set_member(oj, PriorityRows)){
01977         found=dd_TRUE;
01978         if (localdebug) fprintf(stderr,"%ldth in sorted list (row %ld) is in PriorityRows\n", j, oj);
01979         j1=j;
01980       }
01981     }
01982     if (found){
01983       if (j1>i) {
01984         /* shift everything lower: ov[i]->cone->ov[i+1]..ov[j1-1]->cone->ov[j1] */
01985         for (k=j1; k>=i; k--) cone->OrderVector[k]=cone->OrderVector[k-1];
01986         cone->OrderVector[i]=oj;
01987         if (localdebug){
01988           fprintf(stderr,"OrderVector updated to:\n");
01989           for (j = 1; j <= cone->m; j++) fprintf(stderr," %2ld", cone->OrderVector[j]);
01990           fprintf(stderr,"\n");
01991         }
01992       }
01993     } else {
01994       fprintf(stderr,"UpdateRowOrder: Error.\n");
01995       goto _L99;
01996     }
01997   }
01998 _L99:;
01999 }
02000 
02001 void dd_SelectPreorderedNext(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hh)
02002 {
02003   dd_rowrange i,k;
02004   
02005   *hh=0;
02006   for (i=1; i<=cone->m && *hh==0; i++){
02007     k=cone->OrderVector[i];
02008     if (!set_member(k, excluded)) *hh=k ;
02009   }
02010 }
02011 
02012 void dd_SelectNextHalfspace(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hh)
02013 {
02014   if (cone->PreOrderedRun){
02015     if (dd_debug) {
02016       fprintf(stderr,"debug dd_SelectNextHalfspace: Use PreorderNext\n");
02017     }
02018     dd_SelectPreorderedNext(cone, excluded, hh);
02019   }
02020   else {
02021     if (dd_debug) {
02022       fprintf(stderr,"debug dd_SelectNextHalfspace: Use DynamicOrderedNext\n");
02023     }
02024 
02025     switch (cone->HalfspaceOrder) {
02026 
02027     case dd_MaxIndex:
02028       dd_SelectNextHalfspace0(cone, excluded, hh);
02029       break;
02030 
02031     case dd_MinIndex:
02032       dd_SelectNextHalfspace1(cone, excluded, hh);
02033       break;
02034 
02035     case dd_MinCutoff:
02036       dd_SelectNextHalfspace2(cone, excluded, hh);
02037       break;
02038 
02039     case dd_MaxCutoff:
02040       dd_SelectNextHalfspace3(cone, excluded, hh);
02041       break;
02042 
02043     case dd_MixCutoff:
02044       dd_SelectNextHalfspace4(cone, excluded, hh);
02045       break;
02046 
02047     default:
02048       dd_SelectNextHalfspace0(cone, excluded, hh);
02049       break;
02050     }
02051   }
02052 }
02053 
02054 dd_boolean dd_Nonnegative(mytype val)
02055 {
02056 /*  if (val>=-dd_zero) return dd_TRUE;  */
02057   if (dd_cmp(val,dd_minuszero)>=0) return dd_TRUE;
02058   else return dd_FALSE;
02059 }
02060 
02061 dd_boolean dd_Nonpositive(mytype val)
02062 {
02063 /*  if (val<=dd_zero) return dd_TRUE;  */
02064   if (dd_cmp(val,dd_zero)<=0) return dd_TRUE;
02065   else return dd_FALSE;
02066 }
02067 
02068 dd_boolean dd_Positive(mytype val)
02069 {
02070   return !dd_Nonpositive(val);
02071 }
02072 
02073 dd_boolean dd_Negative(mytype val)
02074 {
02075   return !dd_Nonnegative(val);
02076 }
02077 
02078 dd_boolean dd_EqualToZero(mytype val)
02079 {
02080   return (dd_Nonnegative(val) && dd_Nonpositive(val));
02081 }
02082 
02083 dd_boolean dd_Nonzero(mytype val)
02084 {
02085   return (dd_Positive(val) || dd_Negative(val));
02086 }
02087 
02088 dd_boolean dd_Equal(mytype val1,mytype val2)
02089 {
02090   return (!dd_Larger(val1,val2) && !dd_Smaller(val1,val2));
02091 }
02092 
02093 dd_boolean dd_Larger(mytype val1,mytype val2)
02094 {
02095   mytype temp;
02096   dd_boolean answer;
02097 
02098   dd_init(temp);
02099   dd_sub(temp,val1, val2);
02100   answer=dd_Positive(temp);
02101   dd_clear(temp);
02102   return answer;
02103 }
02104 
02105 dd_boolean dd_Smaller(mytype val1,mytype val2)
02106 {
02107   return dd_Larger(val2,val1);
02108 }
02109 
02110 void dd_abs(mytype absval, mytype val)
02111 {
02112   if (dd_Negative(val)) dd_neg(absval,val);
02113   else dd_set(absval,val); 
02114 }
02115 
02116 void dd_LinearComb(mytype lc, mytype v1, mytype c1, mytype v2, mytype c2)
02117 /*  lc := v1 * c1 + v2 * c2   */
02118 {
02119   mytype temp;
02120 
02121   dd_init(temp);
02122   dd_mul(lc,v1,c1);
02123   dd_mul(temp,v2,c2); 
02124   dd_add(lc,lc,temp);
02125   dd_clear(temp);
02126 }
02127 
02128 void dd_InnerProduct(mytype prod, dd_colrange d, dd_Arow v1, dd_Arow v2)
02129 {
02130   mytype temp;
02131   dd_colrange j;
02132   dd_boolean localdebug=dd_FALSE;
02133 
02134   dd_init(temp);
02135   dd_set_si(prod, 0);
02136   for (j = 0; j < d; j++){
02137     dd_mul(temp,v1[j],v2[j]); 
02138     dd_add(prod,prod,temp);
02139   }
02140   if (localdebug) {
02141     fprintf(stderr,"InnerProduct:\n"); 
02142     dd_WriteArow(stderr, v1, d);
02143     dd_WriteArow(stderr, v2, d);
02144     fprintf(stderr,"prod ="); 
02145     dd_WriteNumber(stderr, prod);
02146     fprintf(stderr,"\n");
02147   }
02148   
02149   dd_clear(temp);
02150 }
02151 
02152 /* end of cddcore.c */
02153 
02154 #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