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

cddlp.c

Go to the documentation of this file.
00001 
00043 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00044 
00045 /* cddlp.c:  dual simplex method c-code
00046    written by Komei Fukuda, fukuda@ifor.math.ethz.ch
00047    Version 0.94a, Aug. 24, 2005
00048 */
00049 
00050 /* cddlp.c : C-Implementation of the dual simplex method for
00051    solving an LP: max/min  A_(m-1).x subject to  x in P, where
00052    P= {x :  A_i.x >= 0, i=0,...,m-2, and  x_0=1}, and
00053    A_i is the i-th row of an m x n matrix A.
00054    Please read COPYING (GNU General Public Licence) and
00055    the manual cddlibman.tex for detail.
00056 */
00057 
00058 #include "setoper.h"  /* set operation library header (Ver. May 18, 2000 or later) */
00059 #include "cdd.h"
00060 #include <stdio.h>
00061 #include <stdlib.h>
00062 #include <time.h>
00063 #include <math.h>
00064 #include <string.h>
00065 
00066 #if defined GMPRATIONAL
00067 #include "cdd_f.h"
00068 #endif
00069 
00070 #define dd_CDDLPVERSION  "Version 0.94b (August 25, 2005)"
00071 
00072 #define dd_FALSE 0
00073 #define dd_TRUE 1
00074 
00075 typedef set_type rowset;  /* set_type defined in setoper.h */
00076 typedef set_type colset;
00077 
00078 void dd_CrissCrossSolve(dd_LPPtr lp,dd_ErrorType *);
00079 void dd_DualSimplexSolve(dd_LPPtr lp,dd_ErrorType *);
00080 void dd_CrissCrossMinimize(dd_LPPtr,dd_ErrorType *);
00081 void dd_CrissCrossMaximize(dd_LPPtr,dd_ErrorType *);
00082 void dd_DualSimplexMinimize(dd_LPPtr,dd_ErrorType *);
00083 void dd_DualSimplexMaximize(dd_LPPtr,dd_ErrorType *);
00084 void dd_FindLPBasis(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowindex,dd_rowset,
00085     dd_colindex,dd_rowindex,dd_rowrange,dd_colrange,
00086     dd_colrange *,int *,dd_LPStatusType *,long *);
00087 void dd_FindDualFeasibleBasis(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowindex,
00088     dd_colindex,long *,dd_rowrange,dd_colrange,dd_boolean,
00089     dd_colrange *,dd_ErrorType *,dd_LPStatusType *,long *, long maxpivots);
00090 
00091 
00092 #ifdef GMPRATIONAL
00093 void dd_BasisStatus(ddf_LPPtr lpf, dd_LPPtr lp, dd_boolean*);
00094 void dd_BasisStatusMinimize(dd_rowrange,dd_colrange, dd_Amatrix,dd_Bmatrix,dd_rowset,
00095     dd_rowrange,dd_colrange,ddf_LPStatusType,mytype *,dd_Arow,dd_Arow,dd_rowset,ddf_colindex,
00096     ddf_rowrange,ddf_colrange,long *, int *, int *);
00097 void dd_BasisStatusMaximize(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowset,
00098     dd_rowrange,dd_colrange,ddf_LPStatusType,mytype *,dd_Arow,dd_Arow,dd_rowset,ddf_colindex,
00099     ddf_rowrange,ddf_colrange,long *, int *, int *);
00100 #endif
00101 
00102 void dd_WriteBmatrix(FILE *f,dd_colrange d_size,dd_Bmatrix T);
00103 void dd_SetNumberType(char *line,dd_NumberType *number,dd_ErrorType *Error);
00104 void dd_ComputeRowOrderVector2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,
00105     dd_rowindex OV,dd_RowOrderType ho,unsigned int rseed);
00106 void dd_SelectPreorderedNext2(dd_rowrange m_size,dd_colrange d_size,
00107     rowset excluded,dd_rowindex OV,dd_rowrange *hnext);
00108 void dd_SetSolutions(dd_rowrange,dd_colrange,
00109    dd_Amatrix,dd_Bmatrix,dd_rowrange,dd_colrange,dd_LPStatusType,
00110    mytype *,dd_Arow,dd_Arow,dd_rowset,dd_colindex,dd_rowrange,dd_colrange,dd_rowindex);
00111    
00112 void dd_WriteTableau(FILE *,dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,
00113   dd_colindex,dd_rowindex);
00114 
00115 void dd_WriteSignTableau(FILE *,dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,
00116   dd_colindex,dd_rowindex);
00117 
00118 
00119 dd_LPSolutionPtr dd_CopyLPSolution(dd_LPPtr lp)
00120 {
00121   dd_LPSolutionPtr lps;
00122   dd_colrange j;
00123   long i;
00124 
00125   lps=(dd_LPSolutionPtr) calloc(1,sizeof(dd_LPSolutionType));
00126   for (i=1; i<=dd_filenamelen; i++) lps->filename[i-1]=lp->filename[i-1];
00127   lps->objective=lp->objective;
00128   lps->solver=lp->solver; 
00129   lps->m=lp->m;
00130   lps->d=lp->d;
00131   lps->numbtype=lp->numbtype;
00132 
00133   lps->LPS=lp->LPS;  /* the current solution status */
00134   dd_init(lps->optvalue);
00135   dd_set(lps->optvalue,lp->optvalue);  /* optimal value */
00136   dd_InitializeArow(lp->d+1,&(lps->sol));
00137   dd_InitializeArow(lp->d+1,&(lps->dsol));
00138   lps->nbindex=(long*) calloc((lp->d)+1,sizeof(long));  /* dual solution */
00139   for (j=0; j<=lp->d; j++){
00140     dd_set(lps->sol[j],lp->sol[j]);
00141     dd_set(lps->dsol[j],lp->dsol[j]);
00142     lps->nbindex[j]=lp->nbindex[j];
00143   }
00144   lps->pivots[0]=lp->pivots[0];
00145   lps->pivots[1]=lp->pivots[1];
00146   lps->pivots[2]=lp->pivots[2];
00147   lps->pivots[3]=lp->pivots[3];
00148   lps->pivots[4]=lp->pivots[4];
00149   lps->total_pivots=lp->total_pivots;
00150 
00151   return lps;
00152 }
00153 
00154 
00155 dd_LPPtr dd_CreateLPData(dd_LPObjectiveType obj,
00156    dd_NumberType nt,dd_rowrange m,dd_colrange d)
00157 {
00158   dd_LPType *lp;
00159 
00160   lp=(dd_LPPtr) calloc(1,sizeof(dd_LPType));
00161   lp->solver=dd_choiceLPSolverDefault;  /* set the default lp solver */
00162   lp->d=d;
00163   lp->m=m;
00164   lp->numbtype=nt;
00165   lp->objrow=m;
00166   lp->rhscol=1L;
00167   lp->objective=dd_LPnone;
00168   lp->LPS=dd_LPSundecided;
00169   lp->eqnumber=0;  /* the number of equalities */
00170 
00171   lp->nbindex=(long*) calloc(d+1,sizeof(long));
00172   lp->given_nbindex=(long*) calloc(d+1,sizeof(long));
00173   set_initialize(&(lp->equalityset),m);  
00174     /* i must be in the set iff i-th row is equality . */
00175 
00176   lp->redcheck_extensive=dd_FALSE; /* this is on only for RedundantExtensive */
00177   lp->ired=0;
00178   set_initialize(&(lp->redset_extra),m);  
00179     /* i is in the set if i-th row is newly recognized redundant (during the checking the row ired). */
00180   set_initialize(&(lp->redset_accum),m);  
00181     /* i is in the set if i-th row is recognized redundant (during the checking the row ired). */
00182    set_initialize(&(lp->posset_extra),m);  
00183     /* i is in the set if i-th row is recognized non-linearity (during the course of computation). */
00184  lp->lexicopivot=dd_choiceLexicoPivotQ;  /* dd_choice... is set in dd_set_global_constants() */
00185 
00186   lp->m_alloc=lp->m+2;
00187   lp->d_alloc=lp->d+2;
00188   lp->objective=obj;
00189   dd_InitializeBmatrix(lp->d_alloc,&(lp->B));
00190   dd_InitializeAmatrix(lp->m_alloc,lp->d_alloc,&(lp->A));
00191   dd_InitializeArow(lp->d_alloc,&(lp->sol));
00192   dd_InitializeArow(lp->d_alloc,&(lp->dsol));
00193   dd_init(lp->optvalue);
00194   return lp;
00195 }
00196 
00197 
00198 dd_LPPtr dd_Matrix2LP(dd_MatrixPtr M, dd_ErrorType *err)
00199 {
00200   dd_rowrange m, i, irev, linc;
00201   dd_colrange d, j;
00202   dd_LPType *lp;
00203   dd_boolean localdebug=dd_FALSE;
00204 
00205   *err=dd_NoError;
00206   linc=set_card(M->linset);
00207   m=M->rowsize+1+linc; 
00208      /* We represent each equation by two inequalities.
00209         This is not the best way but makes the code simple. */
00210   d=M->colsize;
00211   if (localdebug) fprintf(stderr,"number of equalities = %ld\n", linc);
00212   
00213   lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
00214   lp->Homogeneous = dd_TRUE;
00215   lp->eqnumber=linc;  /* this records the number of equations */
00216 
00217   irev=M->rowsize; /* the first row of the linc reversed inequalities. */
00218   for (i = 1; i <= M->rowsize; i++) {
00219     if (set_member(i, M->linset)) {
00220       irev=irev+1;
00221       set_addelem(lp->equalityset,i);    /* it is equality. */
00222                                          /* the reversed row irev is not in the equality set. */
00223       for (j = 1; j <= M->colsize; j++) {
00224         dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
00225       }  /*of j*/
00226       if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
00227     }
00228     for (j = 1; j <= M->colsize; j++) {
00229       dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
00230       if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
00231     }  /*of j*/
00232   }  /*of i*/
00233   for (j = 1; j <= M->colsize; j++) {
00234     dd_set(lp->A[m-1][j-1],M->rowvec[j-1]);  /* objective row */
00235   }  /*of j*/
00236 
00237   return lp;
00238 }
00239 
00240 dd_LPPtr dd_Matrix2Feasibility(dd_MatrixPtr M, dd_ErrorType *err)
00241 /* Load a matrix to create an LP object for feasibility.   It is 
00242    essentially the dd_Matrix2LP except that the objject function
00243    is set to identically ZERO (maximization).
00244    
00245 */  
00246          /*  094 */
00247 {
00248   dd_rowrange m, linc;
00249   dd_colrange j;
00250   dd_LPType *lp;
00251 
00252   *err=dd_NoError;
00253   linc=set_card(M->linset);
00254   m=M->rowsize+1+linc; 
00255      /* We represent each equation by two inequalities.
00256         This is not the best way but makes the code simple. */
00257   
00258   lp=dd_Matrix2LP(M, err);
00259   lp->objective = dd_LPmax;   /* since the objective is zero, this is not important */
00260   for (j = 1; j <= M->colsize; j++) {
00261     dd_set(lp->A[m-1][j-1],dd_purezero);  /* set the objective to zero. */
00262   }  /*of j*/
00263 
00264   return lp;
00265 }
00266 
00267 dd_LPPtr dd_Matrix2Feasibility2(dd_MatrixPtr M, dd_rowset R, dd_rowset S, dd_ErrorType *err)
00268 /* Load a matrix to create an LP object for feasibility with additional equality and
00269    strict inequality constraints given by R and S.  There are three types of inequalities:
00270    
00271    b_r + A_r x =  0     Linearity (Equations) specified by M
00272    b_s + A_s x >  0     Strict Inequalities specified by row index set S
00273    b_t + A_t x >= 0     The rest inequalities in M
00274    
00275    Where the linearity is considered here as the union of linearity specified by
00276    M and the additional set R.  When S contains any linearity rows, those
00277    rows are considered linearity (equation).  Thus S does not overlide linearity.
00278    To find a feasible solution, we set an LP
00279    
00280    maximize  z
00281    subject to
00282    b_r + A_r x     =  0      all r in Linearity
00283    b_s + A_s x - z >= 0      for all s in S
00284    b_t + A_t x     >= 0      for all the rest rows t
00285    1           - z >= 0      to make the LP bounded.
00286    
00287    Clearly, the feasibility problem has a solution iff the LP has a positive optimal value. 
00288    The variable z will be the last variable x_{d+1}.
00289    
00290 */  
00291 /*  094 */
00292 {
00293   dd_rowrange m, i, irev, linc;
00294   dd_colrange d, j;
00295   dd_LPType *lp;
00296   dd_rowset L;
00297   dd_boolean localdebug=dd_FALSE;
00298 
00299   *err=dd_NoError;
00300   set_initialize(&L, M->rowsize);
00301   set_uni(L,M->linset,R);
00302   linc=set_card(L);
00303   m=M->rowsize+1+linc+1; 
00304      /* We represent each equation by two inequalities.
00305         This is not the best way but makes the code simple. */
00306   d=M->colsize+1;
00307   if (localdebug) fprintf(stderr,"number of equalities = %ld\n", linc);
00308   
00309   lp=dd_CreateLPData(dd_LPmax, M->numbtype, m, d);
00310   lp->Homogeneous = dd_TRUE;
00311   lp->eqnumber=linc;  /* this records the number of equations */
00312 
00313   irev=M->rowsize; /* the first row of the linc reversed inequalities. */
00314   for (i = 1; i <= M->rowsize; i++) {
00315     if (set_member(i, L)) {
00316       irev=irev+1;
00317       set_addelem(lp->equalityset,i);    /* it is equality. */
00318                                          /* the reversed row irev is not in the equality set. */
00319       for (j = 1; j <= M->colsize; j++) {
00320         dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
00321       }  /*of j*/
00322       if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
00323     } else if (set_member(i, S)) {
00324           dd_set(lp->A[i-1][M->colsize],dd_minusone);
00325     }
00326     for (j = 1; j <= M->colsize; j++) {
00327       dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
00328       if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
00329     }  /*of j*/
00330   }  /*of i*/
00331   for (j = 1; j <= d; j++) {
00332     dd_set(lp->A[m-2][j-1],dd_purezero);  /* initialize */
00333   }  /*of j*/
00334   dd_set(lp->A[m-2][0],dd_one);  /* the bounding constraint. */
00335   dd_set(lp->A[m-2][M->colsize],dd_minusone);  /* the bounding constraint. */
00336   for (j = 1; j <= d; j++) {
00337     dd_set(lp->A[m-1][j-1],dd_purezero);  /* initialize */
00338   }  /*of j*/
00339   dd_set(lp->A[m-1][M->colsize],dd_one);  /* maximize  z */
00340   
00341   set_free(L);
00342   return lp;
00343 }
00344 
00345 
00346 
00347 void dd_FreeLPData(dd_LPPtr lp)
00348 {
00349   if ((lp)!=NULL){
00350     dd_clear(lp->optvalue);
00351     dd_FreeArow(lp->d_alloc,lp->dsol);
00352     dd_FreeArow(lp->d_alloc,lp->sol);
00353     dd_FreeBmatrix(lp->d_alloc,lp->B);
00354     dd_FreeAmatrix(lp->m_alloc,lp->d_alloc,lp->A);
00355     set_free(lp->equalityset);
00356     set_free(lp->redset_extra);
00357     set_free(lp->redset_accum);
00358     set_free(lp->posset_extra);
00359     free(lp->nbindex);
00360     free(lp->given_nbindex);
00361     free(lp);
00362   }
00363 }
00364 
00365 void dd_FreeLPSolution(dd_LPSolutionPtr lps)
00366 {
00367   if (lps!=NULL){
00368     free(lps->nbindex);
00369     dd_FreeArow(lps->d+1,lps->dsol);
00370     dd_FreeArow(lps->d+1,lps->sol);
00371     dd_clear(lps->optvalue);
00372     
00373     free(lps);
00374   }
00375 }
00376 
00377 int dd_LPReverseRow(dd_LPPtr lp, dd_rowrange i)
00378 {
00379   dd_colrange j;
00380   int success=0;
00381 
00382   if (i>=1 && i<=lp->m){
00383     lp->LPS=dd_LPSundecided;
00384     for (j=1; j<=lp->d; j++) {
00385       dd_neg(lp->A[i-1][j-1],lp->A[i-1][j-1]);
00386       /* negating the i-th constraint of A */
00387     }
00388     success=1;
00389   }
00390   return success;
00391 }
00392 
00393 int dd_LPReplaceRow(dd_LPPtr lp, dd_rowrange i, dd_Arow a)
00394 {
00395   dd_colrange j;
00396   int success=0;
00397 
00398   if (i>=1 && i<=lp->m){
00399     lp->LPS=dd_LPSundecided;
00400     for (j=1; j<=lp->d; j++) {
00401       dd_set(lp->A[i-1][j-1],a[j-1]);
00402       /* replacing the i-th constraint by a */
00403     }
00404     success=1;
00405   }
00406   return success;
00407 }
00408 
00409 dd_Arow dd_LPCopyRow(dd_LPPtr lp, dd_rowrange i)
00410 {
00411   dd_colrange j;
00412   dd_Arow a;
00413 
00414   if (i>=1 && i<=lp->m){
00415     dd_InitializeArow(lp->d, &a);
00416     for (j=1; j<=lp->d; j++) {
00417       dd_set(a[j-1],lp->A[i-1][j-1]);
00418       /* copying the i-th row to a */
00419     }
00420   }
00421   return a;
00422 }
00423 
00424 
00425 void dd_SetNumberType(char *line,dd_NumberType *number,dd_ErrorType *Error)
00426 {
00427   if (strncmp(line,"integer",7)==0) {
00428     *number = dd_Integer;
00429     return;
00430   }
00431   else if (strncmp(line,"rational",8)==0) {
00432     *number = dd_Rational;
00433     return;
00434   }
00435   else if (strncmp(line,"real",4)==0) {
00436     *number = dd_Real;
00437     return;
00438   }
00439   else { 
00440     *number=dd_Unknown;
00441     *Error=dd_ImproperInputFormat;
00442   }
00443 }
00444 
00445 
00446 void dd_WriteTableau(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
00447   dd_colindex nbindex,dd_rowindex bflag)
00448 /* Write the tableau  A.T   */
00449 {
00450   dd_colrange j;
00451   dd_rowrange i;
00452   mytype x;
00453   
00454   dd_init(x);
00455   fprintf(f," %ld  %ld  real\n",m_size,d_size);
00456   fprintf(f,"          |");
00457   for (j=1; j<= d_size; j++) {
00458     fprintf(f," %ld",nbindex[j]);
00459   } fprintf(f,"\n");
00460   for (j=1; j<= d_size+1; j++) {
00461     fprintf(f," ----");
00462   } fprintf(f,"\n");
00463   for (i=1; i<= m_size; i++) {
00464     fprintf(f," %3ld(%3ld) |",i,bflag[i]);  
00465     for (j=1; j<= d_size; j++) {
00466       dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
00467       dd_WriteNumber(f,x);
00468     }
00469     fprintf(f,"\n");
00470   }
00471   fprintf(f,"end\n");
00472   dd_clear(x);
00473 }
00474 
00475 void dd_WriteSignTableau(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
00476   dd_colindex nbindex,dd_rowindex bflag)
00477 /* Write the sign tableau  A.T   */
00478 {
00479   dd_colrange j;
00480   dd_rowrange i;
00481   mytype x;
00482   
00483   dd_init(x);
00484   fprintf(f," %ld  %ld  real\n",m_size,d_size);
00485   fprintf(f,"          |");
00486   for (j=1; j<= d_size; j++) {
00487     fprintf(f,"%3ld",nbindex[j]);
00488   } fprintf(f,"\n  ------- | ");
00489   for (j=1; j<= d_size; j++) {
00490     fprintf(f,"---");
00491   } fprintf(f,"\n");
00492   for (i=1; i<= m_size; i++) {
00493     fprintf(f," %3ld(%3ld) |",i,bflag[i]);
00494     for (j=1; j<= d_size; j++) {
00495       dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
00496       if (dd_Positive(x)) fprintf(f, "  +");
00497       else if (dd_Negative(x)) fprintf(f, "  -");
00498         else  fprintf(f, "  0");
00499     }
00500     fprintf(f,"\n");
00501   }
00502   fprintf(f,"end\n");
00503   dd_clear(x);
00504 }
00505 
00506 void dd_WriteSignTableau2(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
00507   dd_colindex nbindex_ref, dd_colindex nbindex,dd_rowindex bflag)
00508 /* Write the sign tableau  A.T  and the reference basis */
00509 {
00510   dd_colrange j;
00511   dd_rowrange i;
00512   mytype x;
00513   
00514   dd_init(x);
00515   fprintf(f," %ld  %ld  real\n",m_size,d_size);
00516   fprintf(f,"          |");
00517   for (j=1; j<= d_size; j++) fprintf(f,"%3ld",nbindex_ref[j]);
00518   fprintf(f,"\n          |");
00519   for (j=1; j<= d_size; j++) {
00520     fprintf(f,"%3ld",nbindex[j]);
00521   } fprintf(f,"\n  ------- | ");
00522   for (j=1; j<= d_size; j++) {
00523     fprintf(f,"---");
00524   } fprintf(f,"\n");
00525   for (i=1; i<= m_size; i++) {
00526     fprintf(f," %3ld(%3ld) |",i,bflag[i]);
00527     for (j=1; j<= d_size; j++) {
00528       dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
00529       if (dd_Positive(x)) fprintf(f, "  +");
00530       else if (dd_Negative(x)) fprintf(f, "  -");
00531         else  fprintf(f, "  0");
00532     }
00533     fprintf(f,"\n");
00534   }
00535   fprintf(f,"end\n");
00536   dd_clear(x);
00537 }
00538 
00539 
00540 void dd_GetRedundancyInformation(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
00541   dd_colindex nbindex,dd_rowindex bflag, dd_rowset redset)
00542 /* Some basic variables that are forced to be nonnegative will be output.  These are
00543    variables whose dictionary row components are all nonnegative.   */
00544 {
00545   dd_colrange j;
00546   dd_rowrange i;
00547   mytype x;
00548   dd_boolean red=dd_FALSE,localdebug=dd_FALSE;
00549   long numbred=0;
00550   
00551   dd_init(x);
00552   for (i=1; i<= m_size; i++) {
00553     red=dd_TRUE;
00554     for (j=1; j<= d_size; j++) {
00555       dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
00556       if (red && dd_Negative(x)) red=dd_FALSE;
00557     }
00558     if (bflag[i]<0 && red) {
00559       numbred+=1;
00560       set_addelem(redset,i);
00561     }
00562   }
00563   if (localdebug) fprintf(stderr,"\ndd_GetRedundancyInformation: %ld redundant rows over %ld\n",numbred, m_size);  
00564   dd_clear(x);
00565 }
00566 
00567 
00568 void dd_SelectDualSimplexPivot(dd_rowrange m_size,dd_colrange d_size,
00569     int Phase1,dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,
00570     dd_colindex nbindex_ref, dd_colindex nbindex,dd_rowindex bflag,
00571     dd_rowrange objrow,dd_colrange rhscol, dd_boolean lexicopivot,
00572     dd_rowrange *r,dd_colrange *s,int *selected,dd_LPStatusType *lps)
00573 { 
00574   /* selects a dual simplex pivot (*r,*s) if the current
00575      basis is dual feasible and not optimal. If not dual feasible,
00576      the procedure returns *selected=dd_FALSE and *lps=LPSundecided.
00577      If Phase1=dd_TRUE, the RHS column will be considered as the negative
00578      of the column of the largest variable (==m_size).  For this case, it is assumed
00579      that the caller used the auxiliary row (with variable m_size) to make the current
00580      dictionary dual feasible before calling this routine so that the nonbasic
00581      column for m_size corresponds to the auxiliary variable.
00582   */
00583   dd_boolean colselected=dd_FALSE,rowselected=dd_FALSE,
00584     dualfeasible=dd_TRUE,localdebug=dd_FALSE;
00585   dd_rowrange i,iref;
00586   dd_colrange j,k;
00587   mytype val,valn, minval,rat,minrat;
00588   static dd_Arow rcost;
00589   static dd_colrange d_last=0;
00590   static dd_colset tieset,stieset;  /* store the column indices with tie */
00591 
00592   dd_init(val); dd_init(valn); dd_init(minval); dd_init(rat); dd_init(minrat);
00593   if (d_last<d_size) {
00594     if (d_last>0) {
00595       for (j=1; j<=d_last; j++){ dd_clear(rcost[j-1]);}
00596       free(rcost);
00597       set_free(tieset);
00598       set_free(stieset);
00599     }
00600     rcost=(mytype*) calloc(d_size,sizeof(mytype));
00601     for (j=1; j<=d_size; j++){ dd_init(rcost[j-1]);}
00602     set_initialize(&tieset,d_size);
00603     set_initialize(&stieset,d_size);
00604   }
00605   d_last=d_size;
00606 
00607   *r=0; *s=0;
00608   *selected=dd_FALSE;
00609   *lps=dd_LPSundecided;
00610   for (j=1; j<=d_size; j++){
00611     if (j!=rhscol){
00612       dd_TableauEntry(&(rcost[j-1]),m_size,d_size,A,T,objrow,j);
00613       if (dd_Positive(rcost[j-1])) { 
00614         dualfeasible=dd_FALSE;
00615       }
00616     }
00617   }
00618   if (dualfeasible){
00619     while ((*lps==dd_LPSundecided) && (!rowselected) && (!colselected)) {
00620       for (i=1; i<=m_size; i++) {
00621         if (i!=objrow && bflag[i]==-1) {  /* i is a basic variable */
00622           if (Phase1){
00623             dd_TableauEntry(&val, m_size,d_size,A,T,i,bflag[m_size]);
00624             dd_neg(val,val);
00625             /* for dual Phase I.  The RHS (dual objective) is the negative of the auxiliary variable column. */
00626           } 
00627           else {dd_TableauEntry(&val,m_size,d_size,A,T,i,rhscol);}
00628           if (dd_Smaller(val,minval)) {
00629             *r=i;
00630             dd_set(minval,val);
00631           }
00632         }
00633       }
00634       if (dd_Nonnegative(minval)) {
00635         *lps=dd_Optimal;
00636       }
00637       else {
00638         rowselected=dd_TRUE;
00639         set_emptyset(tieset);
00640         for (j=1; j<=d_size; j++){
00641           dd_TableauEntry(&val,m_size,d_size,A,T,*r,j);
00642           if (j!=rhscol && dd_Positive(val)) {
00643             dd_div(rat,rcost[j-1],val);
00644             dd_neg(rat,rat);
00645             if (*s==0 || dd_Smaller(rat,minrat)){
00646               dd_set(minrat,rat);
00647               *s=j;
00648               set_emptyset(tieset);
00649               set_addelem(tieset, j);
00650             } else if (dd_Equal(rat,minrat)){
00651               set_addelem(tieset,j);
00652             }
00653           }
00654         }
00655         if (*s>0) {
00656           if (!lexicopivot || set_card(tieset)==1){
00657             colselected=dd_TRUE; *selected=dd_TRUE;
00658           } else { /* lexicographic rule with respect to the given reference cobasis.  */
00659             if (localdebug) {printf("Tie occurred at:"); set_write(tieset); printf("\n");
00660               dd_WriteTableau(stderr,m_size,d_size,A,T,nbindex,bflag);
00661             }
00662             *s=0;
00663             k=2; /* k runs through the column indices except RHS. */
00664             do {
00665               iref=nbindex_ref[k];  /* iref runs though the reference basic indices */
00666               if (iref>0) {
00667                 j=bflag[iref];
00668                 if (j>0) {
00669                   if (set_member(j,tieset) && set_card(tieset)==1) {
00670                     *s=j;
00671                      colselected=dd_TRUE;
00672                   } else {
00673                     set_delelem(tieset, j);
00674                     /* iref is cobasic, and the corresponding col is not the pivot column except it is the last one. */
00675                   }
00676                 } else {
00677                   *s=0;
00678                   for (j=1; j<=d_size; j++){
00679                     if (set_member(j,tieset)) {
00680                       dd_TableauEntry(&val,m_size,d_size,A,T,*r,j);
00681                       dd_TableauEntry(&valn,m_size,d_size,A,T,iref,j);
00682                       if (j!=rhscol && dd_Positive(val)) {
00683                         dd_div(rat,valn,val);
00684                         if (*s==0 || dd_Smaller(rat,minrat)){
00685                           dd_set(minrat,rat);
00686                           *s=j;
00687                           set_emptyset(stieset);
00688                           set_addelem(stieset, j);
00689                         } else if (dd_Equal(rat,minrat)){
00690                           set_addelem(stieset,j);
00691                         }
00692                       }
00693                     }
00694                   }
00695                   set_copy(tieset,stieset);              
00696                   if (set_card(tieset)==1) colselected=dd_TRUE;
00697                 }
00698               }
00699               k+=1;
00700             } while (!colselected && k<=d_size);
00701             *selected=dd_TRUE;
00702           }
00703         } else *lps=dd_Inconsistent;
00704       }
00705     } /* end of while */
00706   }
00707   if (localdebug) {
00708      if (Phase1) fprintf(stderr,"Phase 1 : select %ld,%ld\n",*r,*s);
00709      else fprintf(stderr,"Phase 2 : select %ld,%ld\n",*r,*s);
00710   }
00711   dd_clear(val); dd_clear(valn); dd_clear(minval); dd_clear(rat); dd_clear(minrat);
00712 }
00713 
00714 void dd_TableauEntry(mytype *x,dd_rowrange m_size, dd_colrange d_size, dd_Amatrix X, dd_Bmatrix T,
00715                                 dd_rowrange r, dd_colrange s)
00716 /* Compute the (r,s) entry of X.T   */
00717 {
00718   dd_colrange j;
00719   mytype temp;
00720 
00721   dd_init(temp);
00722   dd_set(*x,dd_purezero);
00723   for (j=0; j< d_size; j++) {
00724     dd_mul(temp,X[r-1][j], T[j][s-1]);
00725     dd_add(*x, *x, temp);
00726   }
00727   dd_clear(temp);
00728 }
00729 
00730 void dd_SelectPivot2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
00731             dd_RowOrderType roworder,dd_rowindex ordervec, rowset equalityset,
00732             dd_rowrange rowmax,rowset NopivotRow,
00733             colset NopivotCol,dd_rowrange *r,dd_colrange *s,
00734             dd_boolean *selected)
00735 /* Select a position (*r,*s) in the matrix A.T such that (A.T)[*r][*s] is nonzero
00736    The choice is feasible, i.e., not on NopivotRow and NopivotCol, and
00737    best with respect to the specified roworder 
00738  */
00739 {
00740   int stop;
00741   dd_rowrange i,rtemp;
00742   rowset rowexcluded;
00743   mytype Xtemp;
00744   dd_boolean localdebug=dd_FALSE;
00745 
00746   stop = dd_FALSE;
00747   localdebug=dd_debug;
00748   dd_init(Xtemp);
00749   set_initialize(&rowexcluded,m_size);
00750   set_copy(rowexcluded,NopivotRow);
00751   for (i=rowmax+1;i<=m_size;i++) {
00752     set_addelem(rowexcluded,i);   /* cannot pivot on any row > rmax */
00753   }
00754   *selected = dd_FALSE;
00755   do {
00756     rtemp=0; i=1;
00757     while (i<=m_size && rtemp==0) {  /* equalityset vars have highest priorities */
00758       if (set_member(i,equalityset) && !set_member(i,rowexcluded)){
00759         if (localdebug) fprintf(stderr,"marked set %ld chosen as a candidate\n",i);
00760         rtemp=i;
00761       }
00762       i++;
00763     }
00764     if (rtemp==0) dd_SelectPreorderedNext2(m_size,d_size,rowexcluded,ordervec,&rtemp);;
00765     if (rtemp>=1) {
00766       *r=rtemp;
00767       *s=1;
00768       while (*s <= d_size && !*selected) {
00769         dd_TableauEntry(&Xtemp,m_size,d_size,A,T,*r,*s);
00770         if (!set_member(*s,NopivotCol) && dd_Nonzero(Xtemp)) {
00771           *selected = dd_TRUE;
00772           stop = dd_TRUE;
00773         } else {
00774           (*s)++;
00775         }
00776       }
00777       if (!*selected) {
00778         set_addelem(rowexcluded,rtemp);
00779       }
00780     }
00781     else {
00782       *r = 0;
00783       *s = 0;
00784       stop = dd_TRUE;
00785     }
00786   } while (!stop);
00787   set_free(rowexcluded); dd_clear(Xtemp);
00788 }
00789 
00790 void dd_GaussianColumnPivot(dd_rowrange m_size, dd_colrange d_size, 
00791     dd_Amatrix X, dd_Bmatrix T, dd_rowrange r, dd_colrange s)
00792 /* Update the Transformation matrix T with the pivot operation on (r,s) 
00793    This procedure performs a implicit pivot operation on the matrix X by
00794    updating the dual basis inverse  T.
00795  */
00796 {
00797   dd_colrange j, j1;
00798   mytype Xtemp0, Xtemp1, Xtemp;
00799   static dd_Arow Rtemp;
00800   static dd_colrange last_d=0;
00801 
00802   dd_init(Xtemp0); dd_init(Xtemp1); dd_init(Xtemp);
00803   if (last_d!=d_size){
00804     if (last_d>0) {
00805       for (j=1; j<=last_d; j++) dd_clear(Rtemp[j-1]);
00806       free(Rtemp);
00807     }
00808     Rtemp=(mytype*)calloc(d_size,sizeof(mytype));
00809     for (j=1; j<=d_size; j++) dd_init(Rtemp[j-1]);
00810     last_d=d_size;
00811   }
00812 
00813   for (j=1; j<=d_size; j++) {
00814     dd_TableauEntry(&(Rtemp[j-1]), m_size, d_size, X, T, r,j);
00815   }
00816   dd_set(Xtemp0,Rtemp[s-1]);
00817   for (j = 1; j <= d_size; j++) {
00818     if (j != s) {
00819       dd_div(Xtemp,Rtemp[j-1],Xtemp0);
00820       dd_set(Xtemp1,dd_purezero);
00821       for (j1 = 1; j1 <= d_size; j1++){
00822         dd_mul(Xtemp1,Xtemp,T[j1-1][s - 1]);
00823         dd_sub(T[j1-1][j-1],T[j1-1][j-1],Xtemp1);
00824  /*     T[j1-1][j-1] -= T[j1-1][s - 1] * Xtemp / Xtemp0;  */
00825       }
00826     }
00827   }
00828   for (j = 1; j <= d_size; j++)
00829     dd_div(T[j-1][s - 1],T[j-1][s - 1],Xtemp0);
00830 
00831   dd_clear(Xtemp0); dd_clear(Xtemp1); dd_clear(Xtemp);
00832 }
00833 
00834 void dd_GaussianColumnPivot2(dd_rowrange m_size,dd_colrange d_size,
00835     dd_Amatrix A,dd_Bmatrix T,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange r,dd_colrange s)
00836 /* Update the Transformation matrix T with the pivot operation on (r,s) 
00837    This procedure performs a implicit pivot operation on the matrix A by
00838    updating the dual basis inverse  T.
00839  */
00840 {
00841   int localdebug=dd_FALSE;
00842   long entering;
00843 
00844   if (dd_debug) localdebug=dd_debug;
00845   dd_GaussianColumnPivot(m_size,d_size,A,T,r,s);
00846   entering=nbindex[s];
00847   bflag[r]=s;     /* the nonbasic variable r corresponds to column s */
00848   nbindex[s]=r;   /* the nonbasic variable on s column is r */
00849 
00850   if (entering>0) bflag[entering]=-1;
00851      /* original variables have negative index and should not affect the row index */
00852 
00853   if (localdebug) {
00854     fprintf(stderr,"dd_GaussianColumnPivot2\n");
00855     fprintf(stderr," pivot: (leaving, entering) = (%ld, %ld)\n", r,entering);
00856     fprintf(stderr, " bflag[%ld] is set to %ld\n", r, s);
00857   }
00858 }
00859 
00860 
00861 void dd_ResetTableau(dd_rowrange m_size,dd_colrange d_size,dd_Bmatrix T,
00862     dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol)
00863 {
00864   dd_rowrange i;
00865   dd_colrange j;
00866   
00867   /* Initialize T and nbindex */
00868   for (j=1; j<=d_size; j++) nbindex[j]=-j;
00869   nbindex[rhscol]=0; 
00870     /* RHS is already in nonbasis and is considered to be associated
00871        with the zero-th row of input. */
00872   dd_SetToIdentity(d_size,T);
00873   
00874   /* Set the bflag according to nbindex */
00875   for (i=1; i<=m_size; i++) bflag[i]=-1;  
00876     /* all basic variables have index -1 */
00877   bflag[objrow]= 0; 
00878     /* bflag of the objective variable is 0,
00879        different from other basic variables which have -1 */
00880   for (j=1; j<=d_size; j++) if (nbindex[j]>0) bflag[nbindex[j]]=j;
00881     /* bflag of a nonbasic variable is its column number */
00882 
00883 }
00884 
00885 void dd_SelectCrissCrossPivot(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
00886     dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
00887     dd_rowrange *r,dd_colrange *s,
00888     int *selected,dd_LPStatusType *lps)
00889 {
00890   int colselected=dd_FALSE,rowselected=dd_FALSE;
00891   dd_rowrange i;
00892   mytype val;
00893   
00894   dd_init(val);
00895   *selected=dd_FALSE;
00896   *lps=dd_LPSundecided;
00897   while ((*lps==dd_LPSundecided) && (!rowselected) && (!colselected)) {
00898     for (i=1; i<=m_size; i++) {
00899       if (i!=objrow && bflag[i]==-1) {  /* i is a basic variable */
00900         dd_TableauEntry(&val,m_size,d_size,A,T,i,rhscol);
00901         if (dd_Negative(val)) {
00902           rowselected=dd_TRUE;
00903           *r=i;
00904           break;
00905         }
00906       }
00907       else if (bflag[i] >0) { /* i is nonbasic variable */
00908         dd_TableauEntry(&val,m_size,d_size,A,T,objrow,bflag[i]);
00909         if (dd_Positive(val)) {
00910           colselected=dd_TRUE;
00911           *s=bflag[i];
00912           break;
00913         }
00914       }
00915     }
00916     if  ((!rowselected) && (!colselected)) {
00917       *lps=dd_Optimal;
00918       return;
00919     }
00920     else if (rowselected) {
00921      for (i=1; i<=m_size; i++) {
00922        if (bflag[i] >0) { /* i is nonbasic variable */
00923           dd_TableauEntry(&val,m_size,d_size,A,T,*r,bflag[i]);
00924           if (dd_Positive(val)) {
00925             colselected=dd_TRUE;
00926             *s=bflag[i];
00927             *selected=dd_TRUE;
00928             break;
00929           }
00930         }
00931       }
00932     }
00933     else if (colselected) {
00934       for (i=1; i<=m_size; i++) {
00935         if (i!=objrow && bflag[i]==-1) {  /* i is a basic variable */
00936           dd_TableauEntry(&val,m_size,d_size,A,T,i,*s);
00937           if (dd_Negative(val)) {
00938             rowselected=dd_TRUE;
00939             *r=i;
00940             *selected=dd_TRUE;
00941             break;
00942           }
00943         }
00944       }
00945     }
00946     if (!rowselected) {
00947       *lps=dd_DualInconsistent;
00948     }
00949     else if (!colselected) {
00950       *lps=dd_Inconsistent;
00951     }
00952   }
00953   dd_clear(val);
00954 }
00955 
00956 void dd_CrissCrossSolve(dd_LPPtr lp, dd_ErrorType *err)
00957 {
00958   switch (lp->objective) {
00959     case dd_LPmax:
00960          dd_CrissCrossMaximize(lp,err);
00961       break;
00962       
00963     case dd_LPmin:
00964          dd_CrissCrossMinimize(lp,err);
00965       break;
00966 
00967     case dd_LPnone: *err=dd_NoLPObjective; break;
00968   }
00969 
00970 }
00971 
00972 void dd_DualSimplexSolve(dd_LPPtr lp, dd_ErrorType *err)
00973 {
00974   switch (lp->objective) {
00975     case dd_LPmax:
00976          dd_DualSimplexMaximize(lp,err);
00977       break;
00978       
00979     case dd_LPmin:
00980          dd_DualSimplexMinimize(lp,err);
00981       break;
00982 
00983     case dd_LPnone: *err=dd_NoLPObjective; break;
00984   }
00985 }
00986 
00987 #ifdef GMPRATIONAL
00988 
00989 dd_LPStatusType LPSf2LPS(ddf_LPStatusType lpsf)
00990 {
00991    dd_LPStatusType lps=dd_LPSundecided;
00992 
00993    switch (lpsf) {
00994    case ddf_LPSundecided: lps=dd_LPSundecided; break;
00995    case ddf_Optimal: lps=dd_Optimal; break;
00996    case ddf_Inconsistent: lps=dd_Inconsistent; break; 
00997    case ddf_DualInconsistent: lps=dd_DualInconsistent; break;
00998    case ddf_StrucInconsistent: lps=dd_StrucInconsistent; break; 
00999    case ddf_StrucDualInconsistent: lps=dd_StrucDualInconsistent; break;
01000    case ddf_Unbounded: lps=dd_Unbounded; break;
01001    case ddf_DualUnbounded: lps=dd_DualUnbounded; break;
01002    }
01003    return lps;
01004 }
01005 
01006 
01007 void dd_BasisStatus(ddf_LPPtr lpf, dd_LPPtr lp, dd_boolean *LPScorrect)
01008 {
01009   int i;
01010   dd_colrange j;
01011   dd_boolean basisfound; 
01012  
01013   switch (lp->objective) {
01014     case dd_LPmax:
01015       dd_BasisStatusMaximize(lp->m,lp->d,lp->A,lp->B,lp->equalityset,lp->objrow,lp->rhscol,
01016            lpf->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lpf->nbindex,lpf->re,lpf->se,lp->pivots, 
01017            &basisfound, LPScorrect);
01018       if (*LPScorrect) {
01019          /* printf("BasisStatus Check: the current basis is verified with GMP\n"); */
01020          lp->LPS=LPSf2LPS(lpf->LPS);
01021          lp->re=lpf->re;
01022          lp->se=lpf->se;
01023          for (j=1; j<=lp->d; j++) lp->nbindex[j]=lpf->nbindex[j]; 
01024       }
01025       for (i=1; i<=5; i++) lp->pivots[i-1]+=lpf->pivots[i-1]; 
01026       break;
01027     case dd_LPmin:
01028       dd_BasisStatusMinimize(lp->m,lp->d,lp->A,lp->B,lp->equalityset,lp->objrow,lp->rhscol,
01029            lpf->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lpf->nbindex,lpf->re,lpf->se,lp->pivots, 
01030            &basisfound, LPScorrect);
01031       if (*LPScorrect) {
01032          /* printf("BasisStatus Check: the current basis is verified with GMP\n"); */
01033          lp->LPS=LPSf2LPS(lpf->LPS);
01034          lp->re=lpf->re;
01035          lp->se=lpf->se;
01036          for (j=1; j<=lp->d; j++) lp->nbindex[j]=lpf->nbindex[j]; 
01037       }
01038       for (i=1; i<=5; i++) lp->pivots[i-1]+=lpf->pivots[i-1]; 
01039       break;
01040     case dd_LPnone:  break;
01041    }      
01042 }
01043 #endif
01044 
01045 void dd_FindLPBasis(dd_rowrange m_size,dd_colrange d_size,
01046     dd_Amatrix A, dd_Bmatrix T,dd_rowindex OV,dd_rowset equalityset, dd_colindex nbindex,
01047     dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
01048     dd_colrange *cs,int *found,dd_LPStatusType *lps,long *pivot_no)
01049 { 
01050   /* Find a LP basis using Gaussian pivots.
01051      If the problem has an LP basis,
01052      the procedure returns *found=dd_TRUE,*lps=LPSundecided and an LP basis.
01053      If the constraint matrix A (excluding the rhs and objective) is not
01054      column indepent, there are two cases.  If the dependency gives a dual
01055      inconsistency, this returns *found=dd_FALSE, *lps=dd_StrucDualInconsistent and 
01056      the evidence column *s.  Otherwise, this returns *found=dd_TRUE, 
01057      *lps=LPSundecided and an LP basis of size less than d_size.  Columns j
01058      that do not belong to the basis (i.e. cannot be chosen as pivot because
01059      they are all zero) will be indicated in nbindex vector: nbindex[j] will
01060      be negative and set to -j.
01061   */
01062   int chosen,stop;
01063   long pivots_p0=0,rank;
01064   colset ColSelected;
01065   rowset RowSelected;
01066   mytype val;
01067 
01068   dd_rowrange r;
01069   dd_colrange j,s;
01070 
01071   dd_init(val);
01072   *found=dd_FALSE; *cs=0; rank=0;
01073   *lps=dd_LPSundecided;
01074 
01075   set_initialize(&RowSelected,m_size);
01076   set_initialize(&ColSelected,d_size);
01077   set_addelem(RowSelected,objrow);
01078   set_addelem(ColSelected,rhscol);
01079 
01080   stop=dd_FALSE;
01081   do {   /* Find a LP basis */
01082     dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset,
01083       m_size,RowSelected,ColSelected,&r,&s,&chosen);
01084     if (chosen) {
01085       set_addelem(RowSelected,r);
01086       set_addelem(ColSelected,s);
01087       dd_GaussianColumnPivot2(m_size,d_size,A,T,nbindex,bflag,r,s);
01088       pivots_p0++;
01089       rank++;
01090     } else {
01091       for (j=1;j<=d_size  && *lps==dd_LPSundecided; j++) {
01092         if (j!=rhscol && nbindex[j]<0){
01093           dd_TableauEntry(&val,m_size,d_size,A,T,objrow,j);
01094           if (dd_Nonzero(val)){  /* dual inconsistent */
01095             *lps=dd_StrucDualInconsistent;
01096             *cs=j;
01097             /* dual inconsistent because the nonzero reduced cost */
01098           }
01099         }
01100       }
01101       if (*lps==dd_LPSundecided) *found=dd_TRUE;  
01102          /* dependent columns but not dual inconsistent. */
01103       stop=dd_TRUE;
01104     }
01105     if (rank==d_size-1) {
01106       stop = dd_TRUE;
01107       *found=dd_TRUE;
01108     }
01109   } while (!stop);
01110 
01111   *pivot_no=pivots_p0;
01112   dd_statBApivots+=pivots_p0;
01113   set_free(RowSelected);
01114   set_free(ColSelected);
01115   dd_clear(val);
01116 }
01117 
01118 
01119 void dd_FindLPBasis2(dd_rowrange m_size,dd_colrange d_size,
01120     dd_Amatrix A, dd_Bmatrix T,dd_rowindex OV,dd_rowset equalityset, dd_colindex nbindex,
01121     dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
01122     dd_colrange *cs,int *found,long *pivot_no)
01123 { 
01124   /* Similar to dd_FindLPBasis but it is much simpler.  This tries to recompute T for
01125   the specified basis given by nbindex.  It will return *found=dd_FALSE if the specified
01126   basis is not a basis.
01127   */
01128   int chosen,stop;
01129   long pivots_p0=0,rank;
01130   dd_colset ColSelected,DependentCols;
01131   dd_rowset RowSelected, NopivotRow;
01132   mytype val;
01133   dd_boolean localdebug=dd_FALSE;
01134 
01135   dd_rowrange r,negcount=0;
01136   dd_colrange j,s;
01137 
01138   dd_init(val);
01139   *found=dd_FALSE; *cs=0; rank=0;
01140 
01141   set_initialize(&RowSelected,m_size);
01142   set_initialize(&DependentCols,d_size);
01143   set_initialize(&ColSelected,d_size);
01144   set_initialize(&NopivotRow,m_size);
01145   set_addelem(RowSelected,objrow);
01146   set_addelem(ColSelected,rhscol);
01147   set_compl(NopivotRow, NopivotRow);  /* set NopivotRow to be the groundset */
01148   
01149   for (j=2; j<=d_size; j++) 
01150     if (nbindex[j]>0) 
01151        set_delelem(NopivotRow, nbindex[j]);
01152     else if (nbindex[j]<0){ 
01153        negcount++;       
01154        set_addelem(DependentCols, -nbindex[j]); 
01155        set_addelem(ColSelected, -nbindex[j]); 
01156     }
01157  
01158   set_uni(RowSelected, RowSelected, NopivotRow);  /* RowSelected is the set of rows not allowed to poviot on */
01159 
01160   stop=dd_FALSE;
01161   do {   /* Find a LP basis */
01162     dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset, m_size,RowSelected,ColSelected,&r,&s,&chosen);
01163     if (chosen) {
01164       set_addelem(RowSelected,r);
01165       set_addelem(ColSelected,s);
01166 
01167       dd_GaussianColumnPivot2(m_size,d_size,A,T,nbindex,bflag,r,s);
01168       if (localdebug && m_size <=10){
01169         dd_WriteBmatrix(stderr,d_size,T);
01170         dd_WriteTableau(stderr,m_size,d_size,A,T,nbindex,bflag);
01171       }
01172       pivots_p0++;
01173       rank++;
01174     } else{
01175       *found=dd_FALSE;   /* cannot pivot on any of the spacified positions. */
01176       stop=dd_TRUE;
01177     }
01178     if (rank==d_size-1-negcount) {
01179       if (negcount){
01180         /* Now it tries to pivot on rows that are supposed to be dependent. */ 
01181         set_diff(ColSelected, ColSelected, DependentCols); 
01182         dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset, m_size,RowSelected,ColSelected,&r,&s,&chosen);
01183         if (chosen) *found=dd_FALSE;  /* not supposed to be independent */
01184         else *found=dd_TRUE;
01185         if (localdebug){
01186           printf("Try to check the dependent cols:");
01187           set_write(DependentCols);
01188           if (chosen) printf("They are not dependent.  Can still pivot on (%ld, %ld)\n",r, s);
01189           else printf("They are indeed dependent.\n");
01190         }
01191       } else {
01192         *found=dd_TRUE;
01193      }   
01194      stop = dd_TRUE;
01195     }
01196   } while (!stop);
01197 
01198   for (j=1; j<=d_size; j++) if (nbindex[j]>0) bflag[nbindex[j]]=j;
01199   *pivot_no=pivots_p0;
01200   set_free(RowSelected);
01201   set_free(ColSelected);
01202   set_free(NopivotRow);
01203   set_free(DependentCols);
01204   dd_clear(val);
01205 }
01206 
01207 void dd_FindDualFeasibleBasis(dd_rowrange m_size,dd_colrange d_size,
01208     dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,
01209     dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol, dd_boolean lexicopivot,
01210     dd_colrange *s,dd_ErrorType *err,dd_LPStatusType *lps,long *pivot_no, long maxpivots)
01211 { 
01212   /* Find a dual feasible basis using Phase I of Dual Simplex method.
01213      If the problem is dual feasible,
01214      the procedure returns *err=NoError, *lps=LPSundecided and a dual feasible
01215      basis.   If the problem is dual infeasible, this returns
01216      *err=NoError, *lps=DualInconsistent and the evidence column *s.
01217      Caution: matrix A must have at least one extra row:  the row space A[m_size] must
01218      have been allocated.
01219   */
01220   dd_boolean phase1,dualfeasible=dd_TRUE;
01221   dd_boolean localdebug=dd_FALSE,chosen,stop;
01222   dd_LPStatusType LPSphase1;
01223   long pivots_p1=0;
01224   dd_rowrange i,r_val;
01225   dd_colrange j,l,ms=0,s_val,local_m_size;
01226   mytype x,val,maxcost,axvalue,maxratio;
01227   static dd_colrange d_last=0;
01228   static dd_Arow rcost;
01229   static dd_colindex nbindex_ref; /* to be used to store the initial feasible basis for lexico rule */
01230 
01231   mytype scaling,svalue;  /* random scaling mytype value */
01232   mytype minval;
01233 
01234   if (dd_debug) localdebug=dd_debug;
01235   dd_init(x); dd_init(val); dd_init(scaling); dd_init(svalue);  dd_init(axvalue);
01236   dd_init(maxcost);  dd_set(maxcost,dd_minuszero);
01237   dd_init(maxratio);  dd_set(maxratio,dd_minuszero);
01238   if (d_last<d_size) {
01239     if (d_last>0) {
01240       for (j=1; j<=d_last; j++){ dd_clear(rcost[j-1]);}
01241       free(rcost);
01242       free(nbindex_ref);
01243     }
01244     rcost=(mytype*) calloc(d_size,sizeof(mytype));
01245     nbindex_ref=(long*) calloc(d_size+1,sizeof(long));
01246     for (j=1; j<=d_size; j++){ dd_init(rcost[j-1]);}
01247   }
01248   d_last=d_size;
01249 
01250   *err=dd_NoError; *lps=dd_LPSundecided; *s=0;
01251   local_m_size=m_size+1;  /* increase m_size by 1 */
01252 
01253   ms=0;  /* ms will be the index of column which has the largest reduced cost */
01254   for (j=1; j<=d_size; j++){
01255     if (j!=rhscol){
01256       dd_TableauEntry(&(rcost[j-1]),local_m_size,d_size,A,T,objrow,j);
01257       if (dd_Larger(rcost[j-1],maxcost)) {dd_set(maxcost,rcost[j-1]); ms = j;}
01258     }
01259   }
01260   if (dd_Positive(maxcost)) dualfeasible=dd_FALSE;
01261 
01262   if (!dualfeasible){
01263     for (j=1; j<=d_size; j++){
01264       dd_set(A[local_m_size-1][j-1], dd_purezero);
01265       for (l=1; l<=d_size; l++){
01266         if (nbindex[l]>0) {
01267           dd_set_si(scaling,l+10);
01268           dd_mul(svalue,A[nbindex[l]-1][j-1],scaling); 
01269           dd_sub(A[local_m_size-1][j-1],A[local_m_size-1][j-1],svalue); 
01270           /* To make the auxiliary row (0,-11,-12,...,-d-10).
01271              It is likely to be better than  (0, -1, -1, ..., -1)
01272              to avoid a degenerate LP.  Version 093c. */
01273         }
01274       }
01275     }
01276     
01277     if (localdebug){
01278       fprintf(stderr,"\ndd_FindDualFeasibleBasis: curruent basis is not dual feasible.\n");
01279       fprintf(stderr,"because of the column %ld assoc. with var %ld   dual cost =",
01280        ms,nbindex[ms]);
01281       dd_WriteNumber(stderr, maxcost);
01282       if (localdebug) {
01283         if (m_size <=100 && d_size <=30){
01284           printf("\ndd_FindDualFeasibleBasis: the starting dictionary.\n");
01285           dd_WriteSignTableau(stdout,m_size+1,d_size,A,T,nbindex,bflag);
01286         }
01287       }
01288     }
01289     
01290     ms=0; 
01291      /* Ratio Test: ms will be now the index of column which has the largest reduced cost 
01292         over the auxiliary row entry */
01293     for (j=1; j<=d_size; j++){
01294       if ((j!=rhscol) && dd_Positive(rcost[j-1])){
01295         dd_TableauEntry(&axvalue,local_m_size,d_size,A,T,local_m_size,j);
01296         if (dd_Nonnegative(axvalue)) {
01297           *err=dd_NumericallyInconsistent; 
01298            /* This should not happen as they are set negative above.  Quit the phase I.*/
01299           if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected.\n");
01300           goto _L99;
01301         }
01302         dd_neg(axvalue,axvalue);
01303         dd_div(axvalue,rcost[j-1],axvalue);  /* axvalue is the negative of ratio that is to be maximized. */
01304         if (dd_Larger(axvalue,maxratio)) {
01305           dd_set(maxratio,axvalue); 
01306           ms = j;
01307         }
01308       }
01309     }
01310 
01311     if (ms==0) {
01312       *err=dd_NumericallyInconsistent; /* This should not happen. Quit the phase I.*/
01313       if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected.\n");
01314       goto _L99;
01315     }
01316 
01317     /* Pivot on (local_m_size,ms) so that the dual basic solution becomes feasible */
01318     dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,local_m_size,ms);
01319     pivots_p1=pivots_p1+1;
01320     if (localdebug) {
01321       printf("\ndd_FindDualFeasibleBasis: Pivot on %ld %ld.\n",local_m_size,ms);
01322     }
01323 
01324   for (j=1; j<=d_size; j++) nbindex_ref[j]=nbindex[j];
01325      /* set the reference basis to be the current feasible basis. */
01326   if (localdebug){
01327     fprintf(stderr, "Store the current feasible basis:");
01328     for (j=1; j<=d_size; j++) fprintf(stderr, " %ld", nbindex_ref[j]);
01329     fprintf(stderr, "\n");   
01330     if (m_size <=100 && d_size <=30)
01331       dd_WriteSignTableau2(stdout,m_size+1,d_size,A,T,nbindex_ref,nbindex,bflag);
01332   }
01333 
01334     phase1=dd_TRUE; stop=dd_FALSE;
01335     do {   /* Dual Simplex Phase I */
01336       chosen=dd_FALSE; LPSphase1=dd_LPSundecided;
01337       if (pivots_p1>maxpivots) {
01338         *err=dd_LPCycling;
01339         fprintf(stderr,"max number %ld of pivots performed in Phase I. Switch to the anticycling phase.\n", maxpivots);
01340         goto _L99;  /* failure due to max no. of pivots performed */
01341       }
01342       dd_SelectDualSimplexPivot(local_m_size,d_size,phase1,A,T,OV,nbindex_ref,nbindex,bflag,
01343         objrow,rhscol,lexicopivot,&r_val,&s_val,&chosen,&LPSphase1);
01344       if (!chosen) {
01345         /* The current dictionary is terminal.  There are two cases:
01346            dd_TableauEntry(local_m_size,d_size,A,T,objrow,ms) is negative or zero.
01347            The first case implies dual infeasible,
01348            and the latter implies dual feasible but local_m_size is still in nonbasis.
01349            We must pivot in the auxiliary variable local_m_size. 
01350         */
01351         dd_TableauEntry(&x,local_m_size,d_size,A,T,objrow,ms);
01352         if (dd_Negative(x)){
01353           *err=dd_NoError; *lps=dd_DualInconsistent;  *s=ms;
01354         }
01355         if (localdebug) {
01356           fprintf(stderr,"\ndd_FindDualFeasibleBasis: the auxiliary variable was forced to enter the basis (# pivots = %ld).\n",pivots_p1);
01357           fprintf(stderr," -- objrow %ld, ms %ld entry: ",objrow,ms);
01358           dd_WriteNumber(stderr, x); fprintf(stderr,"\n");
01359           if (dd_Negative(x)){
01360             fprintf(stderr,"->The basis is dual inconsistent. Terminate.\n");
01361           } else {
01362             fprintf(stderr,"->The basis is feasible. Go to phase II.\n");
01363           }
01364         }
01365 
01366         dd_init(minval);
01367         r_val=0;
01368         for (i=1; i<=local_m_size; i++){
01369           if (bflag[i]<0) { 
01370              /* i is basic and not the objective variable */
01371             dd_TableauEntry(&val,local_m_size,d_size,A,T,i,ms);  /* auxiliary column*/
01372             if (dd_Smaller(val, minval)) {
01373               r_val=i;
01374               dd_set(minval,val);
01375             }
01376           }
01377         }
01378         dd_clear(minval);
01379         
01380         if (r_val==0) {
01381           *err=dd_NumericallyInconsistent; /* This should not happen. Quit the phase I.*/
01382           if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected (r_val is 0).\n");
01383           goto _L99;
01384         }
01385 
01386         dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,r_val,ms);
01387         pivots_p1=pivots_p1+1;
01388         if (localdebug) {
01389           printf("\ndd_FindDualFeasibleBasis: make the %ld-th pivot on %ld  %ld to force the auxiliary variable to enter the basis.\n",pivots_p1,r_val,ms);
01390           if (m_size <=100 && d_size <=30)
01391             dd_WriteSignTableau2(stdout,m_size+1,d_size,A,T,nbindex_ref,nbindex,bflag);
01392         }
01393 
01394         stop=dd_TRUE;
01395 
01396       } else {
01397         dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,r_val,s_val);  
01398         pivots_p1=pivots_p1+1;
01399         if (localdebug) {
01400           printf("\ndd_FindDualFeasibleBasis: make a %ld-th pivot on %ld  %ld\n",pivots_p1,r_val,s_val);
01401           if (m_size <=100 && d_size <=30)
01402             dd_WriteSignTableau2(stdout,local_m_size,d_size,A,T,nbindex_ref,nbindex,bflag);
01403         }
01404 
01405 
01406         if (bflag[local_m_size]<0) {
01407           stop=dd_TRUE; 
01408           if (localdebug) 
01409             fprintf(stderr,"\nDualSimplex Phase I: the auxiliary variable entered the basis (# pivots = %ld).\nGo to phase II\n",pivots_p1);
01410         }
01411       }
01412     } while(!stop);
01413   }
01414 _L99:
01415   *pivot_no=pivots_p1;
01416   dd_statDS1pivots+=pivots_p1;
01417   dd_clear(x); dd_clear(val); dd_clear(maxcost); dd_clear(maxratio);
01418   dd_clear(scaling); dd_clear(svalue); dd_clear(axvalue);
01419 }
01420 
01421 void dd_DualSimplexMinimize(dd_LPPtr lp,dd_ErrorType *err)
01422 {
01423    dd_colrange j;
01424 
01425    *err=dd_NoError;
01426    for (j=1; j<=lp->d; j++)
01427      dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);
01428    dd_DualSimplexMaximize(lp,err);
01429    dd_neg(lp->optvalue,lp->optvalue);
01430    for (j=1; j<=lp->d; j++){
01431      dd_neg(lp->dsol[j-1],lp->dsol[j-1]);
01432      dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);
01433    }
01434 }
01435 
01436 void dd_DualSimplexMaximize(dd_LPPtr lp,dd_ErrorType *err)
01437 /* 
01438 When LP is inconsistent then lp->re returns the evidence row.
01439 When LP is dual-inconsistent then lp->se returns the evidence column.
01440 */
01441 {
01442   int stop,chosen,phase1,found;
01443   long pivots_ds=0,pivots_p0=0,pivots_p1=0,pivots_pc=0,maxpivots,maxpivfactor=20;
01444   dd_boolean localdebug=dd_FALSE,localdebug1=dd_FALSE;
01445 
01446 #if !defined GMPRATIONAL
01447   long maxccpivots,maxccpivfactor=100; 
01448     /* criss-cross should not cycle, but with floating-point arithmetics, it happens
01449        (very rarely).  Jorg Rambau reported such an LP, in August 2003.  Thanks Jorg!
01450     */
01451 #endif
01452 
01453   dd_rowrange i,r;
01454   dd_colrange j,s;
01455   static dd_rowindex bflag;
01456   static long mlast=0,nlast=0;
01457   static dd_rowindex OrderVector;  /* the permutation vector to store a preordered row indeces */
01458   static dd_colindex nbindex_ref; /* to be used to store the initial feasible basis for lexico rule */
01459 
01460   double redpercent=0,redpercent_prev=0,redgain=0;
01461   unsigned int rseed=1;
01462   
01463   /* *err=dd_NoError; */
01464   if (dd_debug) localdebug=dd_debug;
01465   set_emptyset(lp->redset_extra);
01466   for (i=0; i<= 4; i++) lp->pivots[i]=0;
01467   maxpivots=maxpivfactor*lp->d;  /* maximum pivots to be performed before cc pivot is applied. */
01468 #if !defined GMPRATIONAL
01469   maxccpivots=maxccpivfactor*lp->d;  /* maximum pivots to be performed with emergency cc pivots. */
01470 #endif
01471   if (mlast!=lp->m || nlast!=lp->d){
01472      if (mlast>0) { /* called previously with different lp->m */
01473        free(OrderVector);
01474        free(bflag);
01475        free(nbindex_ref);
01476      }
01477      OrderVector=(long *)calloc(lp->m+1,sizeof(*OrderVector));
01478      bflag=(long *) calloc(lp->m+2,sizeof(*bflag));  /* one more element for an auxiliary variable  */
01479      nbindex_ref=(long*) calloc(lp->d+1,sizeof(long));
01480      mlast=lp->m;nlast=lp->d;
01481   }
01482   /* Initializing control variables. */
01483   dd_ComputeRowOrderVector2(lp->m,lp->d,lp->A,OrderVector,dd_MinIndex,rseed);
01484 
01485   lp->re=0; lp->se=0;
01486   
01487   dd_ResetTableau(lp->m,lp->d,lp->B,lp->nbindex,bflag,lp->objrow,lp->rhscol);
01488    
01489   dd_FindLPBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->equalityset,lp->nbindex,bflag,
01490       lp->objrow,lp->rhscol,&s,&found,&(lp->LPS),&pivots_p0);
01491   lp->pivots[0]=pivots_p0;
01492 
01493   if (!found){
01494      lp->se=s;
01495      goto _L99;
01496      /* No LP basis is found, and thus Inconsistent.  
01497      Output the evidence column. */
01498   }
01499 
01500   dd_FindDualFeasibleBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->nbindex,bflag,
01501       lp->objrow,lp->rhscol,lp->lexicopivot,&s, err,&(lp->LPS),&pivots_p1, maxpivots);
01502   lp->pivots[1]=pivots_p1;
01503 
01504   for (j=1; j<=lp->d; j++) nbindex_ref[j]=lp->nbindex[j];
01505      /* set the reference basis to be the current feasible basis. */
01506   if (localdebug){
01507     fprintf(stderr, "dd_DualSimplexMaximize: Store the current feasible basis:");
01508     for (j=1; j<=lp->d; j++) fprintf(stderr, " %ld", nbindex_ref[j]);
01509     fprintf(stderr, "\n");
01510     if (lp->m <=100 && lp->d <=30)
01511       dd_WriteSignTableau2(stdout,lp->m+1,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag); 
01512   }
01513   
01514   if (*err==dd_LPCycling || *err==dd_NumericallyInconsistent){
01515     if (localdebug) fprintf(stderr, "Phase I failed and thus switch to the Criss-Cross method\n");
01516     dd_CrissCrossMaximize(lp,err);
01517     return;
01518   }
01519 
01520   if (lp->LPS==dd_DualInconsistent){
01521      lp->se=s;
01522      goto _L99;
01523      /* No dual feasible basis is found, and thus DualInconsistent.  
01524      Output the evidence column. */
01525   }
01526 
01527   /* Dual Simplex Method */
01528   stop=dd_FALSE;
01529   do {
01530     chosen=dd_FALSE; lp->LPS=dd_LPSundecided; phase1=dd_FALSE;
01531     if (pivots_ds<maxpivots) {
01532       dd_SelectDualSimplexPivot(lp->m,lp->d,phase1,lp->A,lp->B,OrderVector,nbindex_ref,lp->nbindex,bflag,
01533         lp->objrow,lp->rhscol,lp->lexicopivot,&r,&s,&chosen,&(lp->LPS));
01534     }
01535     if (chosen) {
01536       pivots_ds=pivots_ds+1;
01537       if (lp->redcheck_extensive) {
01538         dd_GetRedundancyInformation(lp->m,lp->d,lp->A,lp->B,lp->nbindex, bflag, lp->redset_extra);
01539         set_uni(lp->redset_accum, lp->redset_accum,lp->redset_extra);
01540         redpercent=100*(double)set_card(lp->redset_extra)/(double)lp->m;
01541         redgain=redpercent-redpercent_prev;
01542         redpercent_prev=redpercent;
01543         if (localdebug1){
01544           fprintf(stderr,"\ndd_DualSimplexMaximize: Phase II pivot %ld on (%ld, %ld).\n",pivots_ds,r,s);
01545           fprintf(stderr,"  redundancy %f percent: redset size = %ld\n",redpercent,set_card(lp->redset_extra));
01546         }
01547       }
01548     }
01549     if (!chosen && lp->LPS==dd_LPSundecided) {  
01550       if (localdebug1){
01551          fprintf(stderr,"Warning: an emergency CC pivot in Phase II is performed\n");
01552          /* In principle this should not be executed because we already have dual feasibility
01553             attained and dual simplex pivot should have been chosen.  This might occur
01554             under floating point computation, or the case of cycling.
01555          */
01556       if (localdebug && lp->m <=100 && lp->d <=30){
01557           fprintf(stderr,"\ndd_DualSimplexMaximize: The current dictionary.\n");
01558           dd_WriteSignTableau2(stdout,lp->m,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag);
01559       }
01560     }
01561 
01562 #if !defined GMPRATIONAL
01563       if (pivots_pc>maxccpivots) {
01564         *err=dd_LPCycling;
01565         stop=dd_TRUE;
01566         goto _L99;
01567       }
01568 #endif
01569       
01570       dd_SelectCrissCrossPivot(lp->m,lp->d,lp->A,lp->B,bflag,
01571         lp->objrow,lp->rhscol,&r,&s,&chosen,&(lp->LPS));
01572       if (chosen) pivots_pc=pivots_pc+1;
01573     }
01574     if (chosen) {
01575       dd_GaussianColumnPivot2(lp->m,lp->d,lp->A,lp->B,lp->nbindex,bflag,r,s);
01576       if (localdebug && lp->m <=100 && lp->d <=30){
01577           fprintf(stderr,"\ndd_DualSimplexMaximize: The current dictionary.\n");
01578           dd_WriteSignTableau2(stdout,lp->m,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag);
01579       }
01580     } else {
01581       switch (lp->LPS){
01582         case dd_Inconsistent: lp->re=r;
01583         case dd_DualInconsistent: lp->se=s;
01584 
01585         default: break;
01586       }
01587       stop=dd_TRUE;
01588     }
01589   } while(!stop);
01590 _L99: 
01591   lp->pivots[2]=pivots_ds;
01592   lp->pivots[3]=pivots_pc;
01593   dd_statDS2pivots+=pivots_ds;
01594   dd_statACpivots+=pivots_pc;
01595 
01596   dd_SetSolutions(lp->m,lp->d,lp->A,lp->B,lp->objrow,lp->rhscol,lp->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lp->nbindex,lp->re,lp->se,bflag);
01597 }
01598 
01599 
01600 
01601 void dd_CrissCrossMinimize(dd_LPPtr lp,dd_ErrorType *err)
01602 {
01603    dd_colrange j;
01604 
01605    *err=dd_NoError;
01606    for (j=1; j<=lp->d; j++)
01607      dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);
01608    dd_CrissCrossMaximize(lp,err);
01609    dd_neg(lp->optvalue,lp->optvalue);
01610    for (j=1; j<=lp->d; j++){
01611      dd_neg(lp->dsol[j-1],lp->dsol[j-1]);
01612      dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);
01613    }
01614 }
01615 
01616 void dd_CrissCrossMaximize(dd_LPPtr lp,dd_ErrorType *err)
01617 /* 
01618 When LP is inconsistent then lp->re returns the evidence row.
01619 When LP is dual-inconsistent then lp->se returns the evidence column.
01620 */
01621 {
01622   int stop,chosen,found;
01623   long pivots0,pivots1;
01624 #if !defined GMPRATIONAL
01625   long maxpivots,maxpivfactor=1000; 
01626     /* criss-cross should not cycle, but with floating-point arithmetics, it happens
01627        (very rarely).  Jorg Rambau reported such an LP, in August 2003.  Thanks Jorg!
01628     */
01629 #endif
01630 
01631   dd_rowrange i,r;
01632   dd_colrange s;
01633   static dd_rowindex bflag;
01634   static long mlast=0;
01635   static dd_rowindex OrderVector;  /* the permutation vector to store a preordered row indeces */
01636   unsigned int rseed=1;
01637   dd_colindex nbtemp;
01638 
01639   *err=dd_NoError;
01640 #if !defined GMPRATIONAL
01641   maxpivots=maxpivfactor*lp->d;  /* maximum pivots to be performed when floating-point arithmetics is used. */
01642 #endif
01643   nbtemp=(long *) calloc(lp->d+1,sizeof(long*));
01644   for (i=0; i<= 4; i++) lp->pivots[i]=0;
01645   if (bflag==NULL || mlast!=lp->m){
01646      if (mlast!=lp->m && mlast>0) {
01647        free(bflag);   /* called previously with different lp->m */
01648        free(OrderVector);
01649      }
01650      bflag=(long *) calloc(lp->m+1,sizeof(long*));
01651      OrderVector=(long *)calloc(lp->m+1,sizeof(long*)); 
01652      /* initialize only for the first time or when a larger space is needed */
01653      
01654      mlast=lp->m;
01655   }
01656   /* Initializing control variables. */
01657   dd_ComputeRowOrderVector2(lp->m,lp->d,lp->A,OrderVector,dd_MinIndex,rseed);
01658 
01659   lp->re=0; lp->se=0; pivots1=0;
01660 
01661   dd_ResetTableau(lp->m,lp->d,lp->B,lp->nbindex,bflag,lp->objrow,lp->rhscol);
01662 
01663   dd_FindLPBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->equalityset,
01664       lp->nbindex,bflag,lp->objrow,lp->rhscol,&s,&found,&(lp->LPS),&pivots0);
01665   lp->pivots[0]=pivots0;
01666 
01667   if (!found){
01668      lp->se=s;
01669      goto _L99;
01670      /* No LP basis is found, and thus Inconsistent.  
01671      Output the evidence column. */
01672   }
01673 
01674   stop=dd_FALSE;
01675   do {   /* Criss-Cross Method */
01676 #if !defined GMPRATIONAL
01677     if (pivots1>maxpivots) {
01678       *err=dd_LPCycling;
01679       fprintf(stderr,"max number %ld of pivots performed by the criss-cross method. Most likely due to the floating-point arithmetics error.\n", maxpivots);
01680       goto _L99;  /* failure due to max no. of pivots performed */
01681     }
01682 #endif
01683 
01684     dd_SelectCrissCrossPivot(lp->m,lp->d,lp->A,lp->B,bflag,
01685        lp->objrow,lp->rhscol,&r,&s,&chosen,&(lp->LPS));
01686     if (chosen) {
01687       dd_GaussianColumnPivot2(lp->m,lp->d,lp->A,lp->B,lp->nbindex,bflag,r,s);
01688       pivots1++;
01689     } else {
01690       switch (lp->LPS){
01691         case dd_Inconsistent: lp->re=r;
01692         case dd_DualInconsistent: lp->se=s;
01693 
01694         default: break;
01695       }
01696       stop=dd_TRUE;
01697     }
01698   } while(!stop);
01699   
01700 _L99:
01701   lp->pivots[1]=pivots1;
01702   dd_statCCpivots+=pivots1;
01703   dd_SetSolutions(lp->m,lp->d,lp->A,lp->B,
01704    lp->objrow,lp->rhscol,lp->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lp->nbindex,lp->re,lp->se,bflag);
01705   free(nbtemp);
01706 }
01707 
01708 void dd_SetSolutions(dd_rowrange m_size,dd_colrange d_size,
01709    dd_Amatrix A,dd_Bmatrix T,
01710    dd_rowrange objrow,dd_colrange rhscol,dd_LPStatusType LPS,
01711    mytype *optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset, dd_colindex nbindex,
01712    dd_rowrange re,dd_colrange se,dd_rowindex bflag)
01713 /* 
01714 Assign the solution vectors to sol,dsol,*optvalue after solving
01715 the LP.
01716 */
01717 {
01718   dd_rowrange i;
01719   dd_colrange j;
01720   mytype x,sw;
01721   int localdebug=dd_FALSE;
01722   
01723   dd_init(x); dd_init(sw);
01724   switch (LPS){
01725   case dd_Optimal:
01726     for (j=1;j<=d_size; j++) {
01727       dd_set(sol[j-1],T[j-1][rhscol-1]);
01728       dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
01729       dd_neg(dsol[j-1],x);
01730       dd_TableauEntry(optvalue,m_size,d_size,A,T,objrow,rhscol);
01731       if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]); dd_WriteNumber(stderr, dsol[j-1]); }
01732     }
01733     for (i=1; i<=m_size; i++) {
01734       if (bflag[i]==-1) {  /* i is a basic variable */
01735         dd_TableauEntry(&x,m_size,d_size,A,T,i,rhscol);
01736         if (dd_Positive(x)) set_addelem(posset, i);
01737       }
01738     }
01739 
01740     break;
01741   case dd_Inconsistent:
01742     if (localdebug) fprintf(stderr,"DualSimplexSolve: LP is inconsistent.\n");
01743     for (j=1;j<=d_size; j++) {
01744       dd_set(sol[j-1],T[j-1][rhscol-1]);
01745       dd_TableauEntry(&x,m_size,d_size,A,T,re,j);
01746       dd_neg(dsol[j-1],x);
01747       if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]); dd_WriteNumber(stderr,dsol[j-1]);}
01748     }
01749     break;
01750   case dd_DualInconsistent:
01751     for (j=1;j<=d_size; j++) {
01752       dd_set(sol[j-1],T[j-1][se-1]);
01753       dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
01754       dd_neg(dsol[j-1],x);
01755       if (localdebug) {fprintf(stderr,"dsol[%ld]= \n",nbindex[j]);dd_WriteNumber(stderr,dsol[j-1]);}
01756     }
01757     if (localdebug) printf( "DualSimplexSolve: LP is dual inconsistent.\n");
01758     break;
01759 
01760   case dd_StrucDualInconsistent:
01761     dd_TableauEntry(&x,m_size,d_size,A,T,objrow,se);
01762     if (dd_Positive(x)) dd_set(sw,dd_one);
01763     else dd_neg(sw,dd_one);
01764     for (j=1;j<=d_size; j++) {
01765       dd_mul(sol[j-1],sw,T[j-1][se-1]);
01766       dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
01767       dd_neg(dsol[j-1],x);
01768       if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]);dd_WriteNumber(stderr,dsol[j-1]);}
01769     }
01770     if (localdebug) fprintf(stderr,"DualSimplexSolve: LP is dual inconsistent.\n");
01771     break;
01772 
01773   default:break;
01774   }
01775   dd_clear(x); dd_clear(sw);
01776 }
01777 
01778 
01779 void dd_RandomPermutation2(dd_rowindex OV,long t,unsigned int seed)
01780 {
01781   long k,j,ovj;
01782   double u,xk,r,rand_max=(double) RAND_MAX;
01783   int localdebug=dd_FALSE;
01784 
01785   srand(seed);
01786   for (j=t; j>1 ; j--) {
01787     r=rand();
01788     u=r/rand_max;
01789     xk=(double)(j*u +1);
01790     k=(long)xk;
01791     if (localdebug) fprintf(stderr,"u=%g, k=%ld, r=%g, randmax= %g\n",u,k,r,rand_max);
01792     ovj=OV[j];
01793     OV[j]=OV[k];
01794     OV[k]=ovj;
01795     if (localdebug) fprintf(stderr,"row %ld is exchanged with %ld\n",j,k); 
01796   }
01797 }
01798 
01799 void dd_ComputeRowOrderVector2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,
01800     dd_rowindex OV,dd_RowOrderType ho,unsigned int rseed)
01801 {
01802   long i,itemp;
01803   
01804   OV[0]=0;
01805   switch (ho){
01806   case dd_MaxIndex:
01807     for(i=1; i<=m_size; i++) OV[i]=m_size-i+1;
01808     break;
01809 
01810   case dd_LexMin:
01811     for(i=1; i<=m_size; i++) OV[i]=i;
01812     dd_QuickSort(OV,1,m_size,A,d_size);
01813    break;
01814 
01815   case dd_LexMax:
01816     for(i=1; i<=m_size; i++) OV[i]=i;
01817     dd_QuickSort(OV,1,m_size,A,d_size);
01818     for(i=1; i<=m_size/2;i++){   /* just reverse the order */
01819       itemp=OV[i];
01820       OV[i]=OV[m_size-i+1];
01821       OV[m_size-i+1]=itemp;
01822     }
01823     break;
01824 
01825   case dd_RandomRow:
01826     for(i=1; i<=m_size; i++) OV[i]=i;
01827     if (rseed<=0) rseed=1;
01828     dd_RandomPermutation2(OV,m_size,rseed);
01829     break;
01830 
01831   case dd_MinIndex: 
01832     for(i=1; i<=m_size; i++) OV[i]=i;
01833     break;
01834 
01835   default: 
01836     for(i=1; i<=m_size; i++) OV[i]=i;
01837     break;
01838  }
01839 }
01840 
01841 void dd_SelectPreorderedNext2(dd_rowrange m_size,dd_colrange d_size,
01842     rowset excluded,dd_rowindex OV,dd_rowrange *hnext)
01843 {
01844   dd_rowrange i,k;
01845   
01846   *hnext=0;
01847   for (i=1; i<=m_size && *hnext==0; i++){
01848     k=OV[i];
01849     if (!set_member(k,excluded)) *hnext=k ;
01850   }
01851 }
01852 
01853 #ifdef GMPRATIONAL
01854 
01855 ddf_LPObjectiveType Obj2Obj(dd_LPObjectiveType obj)
01856 {
01857    ddf_LPObjectiveType objf=ddf_LPnone;
01858 
01859    switch (obj) {
01860    case dd_LPnone: objf=ddf_LPnone; break;
01861    case dd_LPmax: objf=ddf_LPmax; break;
01862    case dd_LPmin: objf=ddf_LPmin; break;
01863    }
01864    return objf;
01865 }
01866 
01867 ddf_LPPtr dd_LPgmp2LPf(dd_LPPtr lp)
01868 {
01869   dd_rowrange i;
01870   dd_colrange j;
01871   ddf_LPType *lpf;
01872   double val;
01873   dd_boolean localdebug=dd_FALSE;
01874 
01875   if (localdebug) fprintf(stderr,"Converting a GMP-LP to a float-LP.\n");
01876   
01877   lpf=ddf_CreateLPData(Obj2Obj(lp->objective), ddf_Real, lp->m, lp->d);
01878   lpf->Homogeneous = lp->Homogeneous;
01879   lpf->eqnumber=lp->eqnumber;  /* this records the number of equations */
01880 
01881   for (i = 1; i <= lp->m; i++) {
01882     if (set_member(i, lp->equalityset)) set_addelem(lpf->equalityset,i);    
01883           /* it is equality. Its reversed row will not be in this set */
01884       for (j = 1; j <= lp->d; j++) {
01885         val=mpq_get_d(lp->A[i-1][j-1]);
01886         ddf_set_d(lpf->A[i-1][j-1],val);
01887       }  /*of j*/
01888   }  /*of i*/
01889 
01890   return lpf;
01891 }
01892 
01893 
01894 #endif
01895 
01896 
01897 dd_boolean dd_LPSolve(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType *err)
01898 /* 
01899 The current version of dd_LPSolve that solves an LP with floating-arithmetics first
01900 and then with the specified arithimetics if it is GMP.
01901 
01902 When LP is inconsistent then *re returns the evidence row.
01903 When LP is dual-inconsistent then *se returns the evidence column.
01904 */
01905 {
01906   int i;
01907   dd_boolean found=dd_FALSE;
01908 #ifdef GMPRATIONAL
01909   ddf_LPPtr lpf;
01910   ddf_ErrorType errf;
01911   dd_boolean LPScorrect=dd_FALSE;
01912   dd_boolean localdebug=dd_FALSE;
01913 #endif
01914 
01915   *err=dd_NoError;
01916   lp->solver=solver;
01917   
01918    time(&lp->starttime);
01919 
01920 #ifndef GMPRATIONAL
01921   switch (lp->solver) {
01922     case dd_CrissCross:
01923       dd_CrissCrossSolve(lp,err);
01924       break;
01925     case dd_DualSimplex:
01926       dd_DualSimplexSolve(lp,err);
01927       break;
01928   }
01929 #else
01930   lpf=dd_LPgmp2LPf(lp);
01931   switch (lp->solver) {
01932     case dd_CrissCross:
01933       ddf_CrissCrossSolve(lpf,&errf);    /* First, run with double float. */
01934           if (errf==ddf_NoError){   /* 094a:  fix for a bug reported by Dima Pasechnik */
01935         dd_BasisStatus(lpf,lp, &LPScorrect);    /* Check the basis. */
01936           } else {LPScorrect=dd_FALSE;}
01937       if (!LPScorrect) {
01938          if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");
01939          dd_CrissCrossSolve(lp,err);  /* Rerun with GMP if fails. */
01940       } else {
01941          if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");
01942       }
01943       break;
01944     case dd_DualSimplex:
01945       ddf_DualSimplexSolve(lpf,&errf);    /* First, run with double float. */
01946           if (errf==ddf_NoError){   /* 094a:  fix for a bug reported by Dima Pasechnik */
01947         dd_BasisStatus(lpf,lp, &LPScorrect);    /* Check the basis. */
01948           } else {LPScorrect=dd_FALSE;}
01949       if (!LPScorrect){
01950          if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");
01951          dd_DualSimplexSolve(lp,err);  /* Rerun with GMP if fails. */
01952          if (localdebug){
01953             printf("*total number pivots = %ld (ph0 = %ld, ph1 = %ld, ph2 = %ld, ph3 = %ld, ph4 = %ld)\n",
01954                lp->total_pivots,lp->pivots[0],lp->pivots[1],lp->pivots[2],lp->pivots[3],lp->pivots[4]);
01955             ddf_WriteLPResult(stdout, lpf, errf);
01956             dd_WriteLP(stdout, lp);
01957          }
01958       } else {
01959          if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");
01960       }
01961       break;
01962   }
01963   ddf_FreeLPData(lpf);
01964 #endif
01965 
01966   time(&lp->endtime);
01967   lp->total_pivots=0;
01968   for (i=0; i<=4; i++) lp->total_pivots+=lp->pivots[i];
01969   if (*err==dd_NoError) found=dd_TRUE;
01970   return found;
01971 }
01972 
01973 
01974 dd_boolean dd_LPSolve0(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType *err)
01975 /* 
01976 The original version of dd_LPSolve that solves an LP with specified arithimetics.
01977 
01978 When LP is inconsistent then *re returns the evidence row.
01979 When LP is dual-inconsistent then *se returns the evidence column.
01980 */
01981 {
01982   int i;
01983   dd_boolean found=dd_FALSE;
01984 
01985   *err=dd_NoError;
01986   lp->solver=solver;
01987   time(&lp->starttime);
01988 
01989   switch (lp->solver) {
01990     case dd_CrissCross:
01991       dd_CrissCrossSolve(lp,err);
01992       break;
01993     case dd_DualSimplex:
01994       dd_DualSimplexSolve(lp,err);
01995       break;
01996   }
01997 
01998   time(&lp->endtime);
01999   lp->total_pivots=0;
02000   for (i=0; i<=4; i++) lp->total_pivots+=lp->pivots[i];
02001   if (*err==dd_NoError) found=dd_TRUE;
02002   return found;
02003 }
02004 
02005 
02006 dd_LPPtr dd_MakeLPforInteriorFinding(dd_LPPtr lp)
02007 /* Delete the objective row,
02008    add an extra column with -1's to the matrix A,
02009    add an extra row with (bceil, 0,...,0,-1),
02010    add an objective row with (0,...,0,1), and 
02011    rows & columns, and change m_size and d_size accordingly, to output new_A.
02012   This sets up the LP:
02013   maximize      x_{d+1}
02014   s.t.    A x + x_{d+1}  <=  b
02015                 x_{d+1}  <=  bm * bmax,
02016   where bm is set to 2 by default, and bmax=max{1, b[1],...,b[m_size]}.
02017   Note that the equalitions (linearity) in the input lp will be ignored.
02018 */
02019 {
02020   dd_rowrange m;
02021   dd_colrange d;
02022   dd_NumberType numbtype;
02023   dd_LPObjectiveType obj;
02024   dd_LPType *lpnew;
02025   dd_rowrange i; 
02026   dd_colrange j;
02027   mytype bm,bmax,bceil;
02028   int localdebug=dd_FALSE;
02029 
02030   dd_init(bm); dd_init(bmax); dd_init(bceil);
02031   dd_add(bm,dd_one,dd_one); dd_set(bmax,dd_one);
02032   numbtype=lp->numbtype;
02033   m=lp->m+1;
02034   d=lp->d+1;
02035   obj=dd_LPmax;
02036 
02037   lpnew=dd_CreateLPData(obj, numbtype, m, d);
02038 
02039   for (i=1; i<=lp->m; i++) {
02040     if (dd_Larger(lp->A[i-1][lp->rhscol-1],bmax)) 
02041       dd_set(bmax,lp->A[i-1][lp->rhscol-1]);
02042   }
02043   dd_mul(bceil,bm,bmax);
02044   if (localdebug) {fprintf(stderr,"bceil is set to "); dd_WriteNumber(stderr, bceil);}
02045   
02046   for (i=1; i <= lp->m; i++) {
02047     for (j=1; j <= lp->d; j++) {
02048       dd_set(lpnew->A[i-1][j-1],lp->A[i-1][j-1]);
02049     }
02050   }
02051 
02052   for (i=1;i<=lp->m; i++){
02053     dd_neg(lpnew->A[i-1][lp->d],dd_one);  /* new column with all minus one's */
02054   }
02055 
02056   for (j=1;j<=lp->d;j++){
02057     dd_set(lpnew->A[m-2][j-1],dd_purezero);   /* new row (bceil, 0,...,0,-1) */
02058   }
02059   dd_set(lpnew->A[m-2][0],bceil);  /* new row (bceil, 0,...,0,-1) */
02060 
02061   for (j=1;j<= d-1;j++) {
02062     dd_set(lpnew->A[m-1][j-1],dd_purezero);  /* new obj row with (0,...,0,1) */
02063   }
02064   dd_set(lpnew->A[m-1][d-1],dd_one);    /* new obj row with (0,...,0,1) */
02065  
02066   if (localdebug) dd_WriteAmatrix(stderr, lp->A, lp->m, lp->d);
02067   if (localdebug) dd_WriteAmatrix(stderr, lpnew->A, lpnew->m, lpnew->d);
02068   dd_clear(bm); dd_clear(bmax); dd_clear(bceil);
02069 
02070   return lpnew;
02071 }
02072 
02073 
02074 void dd_WriteLPResult(FILE *f,dd_LPPtr lp,dd_ErrorType err)
02075 {
02076   long j;
02077 
02078   fprintf(f,"* cdd LP solver result\n");
02079   
02080   if (err!=dd_NoError) {
02081     dd_WriteErrorMessages(f,err);
02082     goto _L99;
02083   }
02084 
02085   dd_WriteProgramDescription(f);
02086 
02087   fprintf(f,"* #constraints = %ld\n",lp->m-1);
02088   fprintf(f,"* #variables   = %ld\n",lp->d-1);
02089 
02090   switch (lp->solver) {
02091     case dd_DualSimplex:
02092       fprintf(f,"* Algorithm: dual simplex algorithm\n");break; 
02093     case dd_CrissCross:
02094       fprintf(f,"* Algorithm: criss-cross method\n");break;
02095   }
02096 
02097   switch (lp->objective) {
02098     case dd_LPmax:
02099       fprintf(f,"* maximization is chosen\n");break; 
02100     case dd_LPmin:
02101       fprintf(f,"* minimization is chosen\n");break;
02102     case dd_LPnone:
02103       fprintf(f,"* no objective type (max or min) is chosen\n");break;
02104   }
02105   
02106   if (lp->objective==dd_LPmax||lp->objective==dd_LPmin){
02107     fprintf(f,"* Objective function is\n");  
02108     for (j=0; j<lp->d; j++){
02109       if (j>0 && dd_Nonnegative(lp->A[lp->objrow-1][j]) ) fprintf(f," +");
02110       if (j>0 && (j % 5) == 0) fprintf(f,"\n");
02111       dd_WriteNumber(f,lp->A[lp->objrow-1][j]);
02112       if (j>0) fprintf(f," X[%3ld]",j);
02113     }
02114     fprintf(f,"\n");
02115   }
02116 
02117   switch (lp->LPS){
02118   case dd_Optimal:
02119     fprintf(f,"* LP status: a dual pair (x,y) of optimal solutions found.\n");
02120     fprintf(f,"begin\n");
02121     fprintf(f,"  primal_solution\n");
02122     for (j=1; j<lp->d; j++) {
02123       fprintf(f,"  %3ld : ",j);
02124       dd_WriteNumber(f,lp->sol[j]);
02125       fprintf(f,"\n");
02126     }
02127     fprintf(f,"  dual_solution\n");
02128     for (j=1; j<lp->d; j++){
02129       if (lp->nbindex[j+1]>0) {
02130         fprintf(f,"  %3ld : ",lp->nbindex[j+1]);
02131         dd_WriteNumber(f,lp->dsol[j]); fprintf(f,"\n");
02132       }
02133     }
02134     fprintf(f,"  optimal_value : "); dd_WriteNumber(f,lp->optvalue);
02135     fprintf(f,"\nend\n");
02136     break;
02137 
02138   case dd_Inconsistent:
02139     fprintf(f,"* LP status: LP is inconsistent.\n");
02140     fprintf(f,"* The positive combination of original inequalities with\n");
02141     fprintf(f,"* the following coefficients will prove the inconsistency.\n");
02142     fprintf(f,"begin\n");
02143     fprintf(f,"  dual_direction\n");
02144     fprintf(f,"  %3ld : ",lp->re);
02145     dd_WriteNumber(f,dd_one);  fprintf(f,"\n");
02146     for (j=1; j<lp->d; j++){
02147       if (lp->nbindex[j+1]>0) {
02148         fprintf(f,"  %3ld : ",lp->nbindex[j+1]);
02149         dd_WriteNumber(f,lp->dsol[j]); fprintf(f,"\n");
02150       }
02151     }
02152     fprintf(f,"end\n");
02153     break;
02154 
02155   case dd_DualInconsistent: case dd_StrucDualInconsistent:
02156     fprintf(f,"* LP status: LP is dual inconsistent.\n");
02157     fprintf(f,"* The linear combination of columns with\n");
02158     fprintf(f,"* the following coefficients will prove the dual inconsistency.\n");
02159     fprintf(f,"* (It is also an unbounded direction for the primal LP.)\n");
02160     fprintf(f,"begin\n");
02161     fprintf(f,"  primal_direction\n");
02162     for (j=1; j<lp->d; j++) {
02163       fprintf(f,"  %3ld : ",j);
02164       dd_WriteNumber(f,lp->sol[j]);
02165       fprintf(f,"\n");
02166     }
02167     fprintf(f,"end\n");
02168     break;
02169 
02170   default:
02171     break;
02172   }
02173   fprintf(f,"* number of pivot operations = %ld (ph0 = %ld, ph1 = %ld, ph2 = %ld, ph3 = %ld, ph4 = %ld)\n",lp->total_pivots,lp->pivots[0],lp->pivots[1],lp->pivots[2],lp->pivots[3],lp->pivots[4]);
02174   dd_WriteLPTimes(f, lp);
02175 _L99:;
02176 }
02177 
02178 dd_LPPtr dd_CreateLP_H_ImplicitLinearity(dd_MatrixPtr M)
02179 {
02180   dd_rowrange m, i, irev, linc;
02181   dd_colrange d, j;
02182   dd_LPPtr lp;
02183   dd_boolean localdebug=dd_FALSE;
02184 
02185   linc=set_card(M->linset);
02186   m=M->rowsize+1+linc+1; 
02187      /* We represent each equation by two inequalities.
02188         This is not the best way but makes the code simple. */
02189   d=M->colsize+1;
02190   
02191   lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
02192   lp->Homogeneous = dd_TRUE;
02193   lp->objective = dd_LPmax;
02194   lp->eqnumber=linc;  /* this records the number of equations */
02195   lp->redcheck_extensive=dd_FALSE;  /* this is default */
02196 
02197   irev=M->rowsize; /* the first row of the linc reversed inequalities. */
02198   for (i = 1; i <= M->rowsize; i++) {
02199     if (set_member(i, M->linset)) {
02200       irev=irev+1;
02201       set_addelem(lp->equalityset,i);    /* it is equality. */
02202             /* the reversed row irev is not in the equality set. */
02203       for (j = 1; j <= M->colsize; j++) {
02204         dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
02205       }  /*of j*/
02206     } else {
02207       dd_set(lp->A[i-1][d-1],dd_minusone);  /* b_I + A_I x - 1 z >= 0  (z=x_d) */
02208     }
02209     for (j = 1; j <= M->colsize; j++) {
02210       dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
02211       if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
02212     }  /*of j*/
02213   }  /*of i*/
02214   dd_set(lp->A[m-2][0],dd_one);  dd_set(lp->A[m-2][d-1],dd_minusone);
02215       /* make the LP bounded.  */
02216   
02217   dd_set(lp->A[m-1][d-1],dd_one);
02218       /* objective is to maximize z.  */
02219 
02220   if (localdebug) {
02221     fprintf(stderr,"dd_CreateLP_H_ImplicitLinearity: an new lp is\n");
02222     dd_WriteLP(stderr,lp);
02223   }
02224 
02225   return lp;
02226 }
02227 
02228 dd_LPPtr dd_CreateLP_V_ImplicitLinearity(dd_MatrixPtr M)
02229 {
02230   dd_rowrange m, i, irev, linc;
02231   dd_colrange d, j;
02232   dd_LPPtr lp;
02233   dd_boolean localdebug=dd_FALSE;
02234 
02235   linc=set_card(M->linset);
02236   m=M->rowsize+1+linc+1; 
02237      /* We represent each equation by two inequalities.
02238         This is not the best way but makes the code simple. */
02239   d=(M->colsize)+2;  
02240      /* Two more columns.  This is different from the H-reprentation case */
02241   
02242 /* The below must be modified for V-representation!!!  */
02243 
02244   lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
02245   lp->Homogeneous = dd_FALSE;
02246   lp->objective = dd_LPmax;
02247   lp->eqnumber=linc;  /* this records the number of equations */
02248   lp->redcheck_extensive=dd_FALSE;  /* this is default */
02249 
02250   irev=M->rowsize; /* the first row of the linc reversed inequalities. */
02251   for (i = 1; i <= M->rowsize; i++) {
02252     dd_set(lp->A[i-1][0],dd_purezero);  /* It is almost completely degerate LP */
02253     if (set_member(i, M->linset)) {
02254       irev=irev+1;
02255       set_addelem(lp->equalityset,i);    /* it is equality. */
02256             /* the reversed row irev is not in the equality set. */
02257       for (j = 2; j <= (M->colsize)+1; j++) {
02258         dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-2]);
02259       }  /*of j*/
02260       if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
02261     } else {
02262       dd_set(lp->A[i-1][d-1],dd_minusone);  /* b_I x_0 + A_I x - 1 z >= 0 (z=x_d) */
02263     }
02264     for (j = 2; j <= (M->colsize)+1; j++) {
02265       dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-2]);
02266     }  /*of j*/
02267   }  /*of i*/
02268   dd_set(lp->A[m-2][0],dd_one);  dd_set(lp->A[m-2][d-1],dd_minusone);
02269       /* make the LP bounded.  */
02270   dd_set(lp->A[m-1][d-1],dd_one);
02271       /* objective is to maximize z.  */
02272 
02273   if (localdebug) {
02274     fprintf(stderr,"dd_CreateLP_V_ImplicitLinearity: an new lp is\n");
02275     dd_WriteLP(stderr,lp);
02276   }
02277 
02278   return lp;
02279 }
02280 
02281 
02282 dd_LPPtr dd_CreateLP_H_Redundancy(dd_MatrixPtr M, dd_rowrange itest)
02283 {
02284   dd_rowrange m, i, irev, linc;
02285   dd_colrange d, j;
02286   dd_LPPtr lp;
02287   dd_boolean localdebug=dd_FALSE;
02288 
02289   linc=set_card(M->linset);
02290   m=M->rowsize+1+linc; 
02291      /* We represent each equation by two inequalities.
02292         This is not the best way but makes the code simple. */
02293   d=M->colsize;
02294   
02295   lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
02296   lp->Homogeneous = dd_TRUE;
02297   lp->objective = dd_LPmin;
02298   lp->eqnumber=linc;  /* this records the number of equations */
02299   lp->redcheck_extensive=dd_FALSE;  /* this is default */
02300 
02301   irev=M->rowsize; /* the first row of the linc reversed inequalities. */
02302   for (i = 1; i <= M->rowsize; i++) {
02303     if (set_member(i, M->linset)) {
02304       irev=irev+1;
02305       set_addelem(lp->equalityset,i);    /* it is equality. */
02306             /* the reversed row irev is not in the equality set. */
02307       for (j = 1; j <= M->colsize; j++) {
02308         dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
02309       }  /*of j*/
02310       if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
02311     }
02312     for (j = 1; j <= M->colsize; j++) {
02313       dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
02314       if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
02315     }  /*of j*/
02316   }  /*of i*/
02317   for (j = 1; j <= M->colsize; j++) {
02318     dd_set(lp->A[m-1][j-1],M->matrix[itest-1][j-1]);
02319       /* objective is to violate the inequality in question.  */
02320   }  /*of j*/
02321   dd_add(lp->A[itest-1][0],lp->A[itest-1][0],dd_one); /* relax the original inequality by one */
02322 
02323   return lp;
02324 }
02325 
02326 
02327 dd_LPPtr dd_CreateLP_V_Redundancy(dd_MatrixPtr M, dd_rowrange itest)
02328 {
02329   dd_rowrange m, i, irev, linc;
02330   dd_colrange d, j;
02331   dd_LPPtr lp;
02332   dd_boolean localdebug=dd_FALSE;
02333 
02334   linc=set_card(M->linset);
02335   m=M->rowsize+1+linc; 
02336      /* We represent each equation by two inequalities.
02337         This is not the best way but makes the code simple. */
02338   d=(M->colsize)+1;  
02339      /* One more column.  This is different from the H-reprentation case */
02340   
02341 /* The below must be modified for V-representation!!!  */
02342 
02343   lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
02344   lp->Homogeneous = dd_FALSE;
02345   lp->objective = dd_LPmin;
02346   lp->eqnumber=linc;  /* this records the number of equations */
02347   lp->redcheck_extensive=dd_FALSE;  /* this is default */
02348 
02349   irev=M->rowsize; /* the first row of the linc reversed inequalities. */
02350   for (i = 1; i <= M->rowsize; i++) {
02351     if (i==itest){
02352       dd_set(lp->A[i-1][0],dd_one); /* this is to make the LP bounded, ie. the min >= -1 */
02353     } else {
02354       dd_set(lp->A[i-1][0],dd_purezero);  /* It is almost completely degerate LP */
02355     }
02356     if (set_member(i, M->linset)) {
02357       irev=irev+1;
02358       set_addelem(lp->equalityset,i);    /* it is equality. */
02359             /* the reversed row irev is not in the equality set. */
02360       for (j = 2; j <= (M->colsize)+1; j++) {
02361         dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-2]);
02362       }  /*of j*/
02363       if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
02364     }
02365     for (j = 2; j <= (M->colsize)+1; j++) {
02366       dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-2]);
02367     }  /*of j*/
02368   }  /*of i*/
02369   for (j = 2; j <= (M->colsize)+1; j++) {
02370     dd_set(lp->A[m-1][j-1],M->matrix[itest-1][j-2]);
02371       /* objective is to violate the inequality in question.  */
02372   }  /*of j*/
02373   dd_set(lp->A[m-1][0],dd_purezero);   /* the constant term for the objective is zero */
02374 
02375   if (localdebug) dd_WriteLP(stdout, lp);
02376 
02377   return lp;
02378 }
02379 
02380 
02381 dd_LPPtr dd_CreateLP_V_SRedundancy(dd_MatrixPtr M, dd_rowrange itest)
02382 {
02383 /*
02384      V-representation (=boundary problem)
02385        g* = maximize  
02386          1^T b_{I-itest} x_0 + 1^T A_{I-itest}    (the sum of slacks)
02387        subject to
02388          b_itest x_0     + A_itest x      =  0 (the point has to lie on the boundary)
02389          b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators in one side)
02390          1^T b_{I-itest} x_0 + 1^T A_{I-itest} x <=  1 (to make an LP bounded)
02391          b_L x_0         + A_L x = 0.  (linearity generators)
02392          
02393     The redundant row is strongly redundant if and only if g* is zero.
02394 */
02395 
02396   dd_rowrange m, i, irev, linc;
02397   dd_colrange d, j;
02398   dd_LPPtr lp;
02399   dd_boolean localdebug=dd_FALSE;
02400 
02401   linc=set_card(M->linset);
02402   m=M->rowsize+1+linc+2; 
02403      /* We represent each equation by two inequalities.
02404         This is not the best way but makes the code simple.
02405         Two extra constraints are for the first equation and the bouding inequality.
02406         */
02407   d=(M->colsize)+1;  
02408      /* One more column.  This is different from the H-reprentation case */
02409   
02410 /* The below must be modified for V-representation!!!  */
02411 
02412   lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
02413   lp->Homogeneous = dd_FALSE;
02414   lp->objective = dd_LPmax;
02415   lp->eqnumber=linc;  /* this records the number of equations */
02416 
02417   irev=M->rowsize; /* the first row of the linc reversed inequalities. */
02418   for (i = 1; i <= M->rowsize; i++) {
02419     if (i==itest){
02420       dd_set(lp->A[i-1][0],dd_purezero);  /* this is a half of the boundary constraint. */
02421     } else {
02422       dd_set(lp->A[i-1][0],dd_purezero);  /* It is almost completely degerate LP */
02423     }
02424     if (set_member(i, M->linset) || i==itest) {
02425       irev=irev+1;
02426       set_addelem(lp->equalityset,i);    /* it is equality. */
02427             /* the reversed row irev is not in the equality set. */
02428       for (j = 2; j <= (M->colsize)+1; j++) {
02429         dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-2]);
02430       }  /*of j*/
02431       if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
02432     }
02433     for (j = 2; j <= (M->colsize)+1; j++) {
02434       dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-2]);
02435       dd_add(lp->A[m-1][j-1],lp->A[m-1][j-1],lp->A[i-1][j-1]);  /* the objective is the sum of all ineqalities */
02436     }  /*of j*/
02437   }  /*of i*/
02438   for (j = 2; j <= (M->colsize)+1; j++) {
02439     dd_neg(lp->A[m-2][j-1],lp->A[m-1][j-1]);
02440       /* to make an LP bounded.  */
02441   }  /*of j*/
02442   dd_set(lp->A[m-2][0],dd_one);   /* the constant term for the bounding constraint is 1 */
02443 
02444   if (localdebug) dd_WriteLP(stdout, lp);
02445 
02446   return lp;
02447 }
02448 
02449 dd_boolean dd_Redundant(dd_MatrixPtr M, dd_rowrange itest, dd_Arow certificate, dd_ErrorType *error)  
02450   /* 092 */
02451 {
02452   /* Checks whether the row itest is redundant for the representation.
02453      All linearity rows are not checked and considered NONredundant. 
02454      This code works for both H- and V-representations.  A certificate is
02455      given in the case of non-redundancy, showing a solution x violating only the itest
02456      inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
02457      separates the itest from the rest.  More explicitly, the LP to be setup is
02458 
02459      H-representation
02460        f* = minimize  
02461          b_itest     + A_itest x
02462        subject to
02463          b_itest + 1 + A_itest x     >= 0 (relaxed inequality to make an LP bounded)
02464          b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
02465          b_L         + A_L x = 0.  (linearity)
02466 
02467      V-representation (=separation problem)
02468        f* = minimize  
02469          b_itest x_0     + A_itest x
02470        subject to
02471          b_itest x_0     + A_itest x     >= -1 (to make an LP bounded)
02472          b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators except for itest in one side)
02473          b_L x_0         + A_L x = 0.  (linearity generators)
02474     
02475     Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
02476     and the row indices of input is partitioned into I and L where L is the set of linearity.
02477     In both cases, the itest data is nonredundant if and only if the optimal value f* is negative.
02478     The certificate has dimension one more for V-representation case.
02479   */
02480 
02481   dd_colrange j;
02482   dd_LPPtr lp;
02483   dd_LPSolutionPtr lps;
02484   dd_ErrorType err=dd_NoError;
02485   dd_boolean answer=dd_FALSE,localdebug=dd_FALSE;
02486 
02487   *error=dd_NoError;
02488   if (set_member(itest, M->linset)){
02489     if (localdebug) printf("The %ld th row is linearity and redundancy checking is skipped.\n",itest);
02490     goto _L99;
02491   }
02492   
02493   /* Create an LP data for redundancy checking */
02494   if (M->representation==dd_Generator){
02495     lp=dd_CreateLP_V_Redundancy(M, itest);
02496   } else {
02497     lp=dd_CreateLP_H_Redundancy(M, itest);
02498   }
02499 
02500   dd_LPSolve(lp,dd_choiceRedcheckAlgorithm,&err);
02501   if (err!=dd_NoError){
02502     *error=err;
02503     goto _L999;
02504   } else {
02505     lps=dd_CopyLPSolution(lp);
02506 
02507     for (j=0; j<lps->d; j++) {
02508       dd_set(certificate[j], lps->sol[j]);
02509     }
02510 
02511     if (dd_Negative(lps->optvalue)){
02512       answer=dd_FALSE;
02513       if (localdebug) fprintf(stderr,"==> %ld th row is nonredundant.\n",itest);
02514     } else {
02515       answer=dd_TRUE;
02516       if (localdebug) fprintf(stderr,"==> %ld th row is redundant.\n",itest);
02517     }
02518     dd_FreeLPSolution(lps);
02519   }
02520   _L999:
02521   dd_FreeLPData(lp);
02522 _L99:
02523   return answer;
02524 }
02525 
02526 dd_boolean dd_RedundantExtensive(dd_MatrixPtr M, dd_rowrange itest, dd_Arow certificate, 
02527 dd_rowset *redset,dd_ErrorType *error)  
02528   /* 094 */
02529 {
02530   /* This uses the same LP construction as dd_Reduandant.  But, while it is checking
02531      the redundancy of itest, it also tries to find some other variable that are
02532      redundant (i.e. forced to be nonnegative).  This is expensive as it used
02533      the complete tableau information at each DualSimplex pivot.  The redset must
02534      be initialized before this function is called.
02535   */
02536 
02537   dd_colrange j;
02538   dd_LPPtr lp;
02539   dd_LPSolutionPtr lps;
02540   dd_ErrorType err=dd_NoError;
02541   dd_boolean answer=dd_FALSE,localdebug=dd_FALSE;
02542 
02543   *error=dd_NoError;
02544   if (set_member(itest, M->linset)){
02545     if (localdebug) printf("The %ld th row is linearity and redundancy checking is skipped.\n",itest);
02546     goto _L99;
02547   }
02548   
02549   /* Create an LP data for redundancy checking */
02550   if (M->representation==dd_Generator){
02551     lp=dd_CreateLP_V_Redundancy(M, itest);
02552   } else {
02553     lp=dd_CreateLP_H_Redundancy(M, itest);
02554   }
02555   
02556   lp->redcheck_extensive=dd_TRUE;
02557 
02558   dd_LPSolve0(lp,dd_DualSimplex,&err);
02559   if (err!=dd_NoError){
02560     *error=err;
02561     goto _L999;
02562   } else {
02563     set_copy(*redset,lp->redset_extra);
02564     set_delelem(*redset, itest);  
02565     /* itest row might be redundant in the lp but this has nothing to do with its redundancy
02566     in the original system M.   Thus we must delete it.  */
02567     if (localdebug){
02568       fprintf(stderr, "dd_RedundantExtensive: checking for %ld, extra redset with cardinality %ld (%ld)\n",itest,set_card(*redset),set_card(lp->redset_extra)); 
02569       set_fwrite(stderr, *redset); fprintf(stderr, "\n");
02570     }
02571     lps=dd_CopyLPSolution(lp);
02572 
02573     for (j=0; j<lps->d; j++) {
02574       dd_set(certificate[j], lps->sol[j]);
02575     }
02576 
02577     if (dd_Negative(lps->optvalue)){
02578       answer=dd_FALSE;
02579       if (localdebug) fprintf(stderr,"==> %ld th row is nonredundant.\n",itest);
02580     } else {
02581       answer=dd_TRUE;
02582       if (localdebug) fprintf(stderr,"==> %ld th row is redundant.\n",itest);
02583     }
02584     dd_FreeLPSolution(lps);
02585   }
02586   _L999:
02587   dd_FreeLPData(lp);
02588 _L99:
02589   return answer;
02590 }
02591 
02592 dd_rowset dd_RedundantRows(dd_MatrixPtr M, dd_ErrorType *error)  /* 092 */
02593 {
02594   dd_rowrange i,m;
02595   dd_colrange d;
02596   dd_rowset redset;
02597   dd_MatrixPtr Mcopy;
02598   dd_Arow cvec; /* certificate */  
02599   dd_boolean localdebug=dd_FALSE;
02600 
02601   m=M->rowsize;
02602   if (M->representation==dd_Generator){
02603     d=(M->colsize)+1;
02604   } else {
02605     d=M->colsize;
02606   }
02607   Mcopy=dd_MatrixCopy(M);
02608   dd_InitializeArow(d,&cvec); 
02609   set_initialize(&redset, m);
02610   for (i=m; i>=1; i--) {
02611     if (dd_Redundant(Mcopy, i, cvec, error)) {
02612       if (localdebug) printf("dd_RedundantRows: the row %ld is redundant.\n", i);
02613       set_addelem(redset, i);
02614       dd_MatrixRowRemove(&Mcopy, i);
02615     } else {
02616       if (localdebug) printf("dd_RedundantRows: the row %ld is essential.\n", i);
02617     }
02618     if (*error!=dd_NoError) goto _L99;
02619   }
02620 _L99:
02621   dd_FreeMatrix(Mcopy);
02622   dd_FreeArow(d, cvec);
02623   return redset;
02624 }
02625 
02626 
02627 dd_boolean dd_MatrixRedundancyRemove(dd_MatrixPtr *M, dd_rowset *redset,dd_rowindex *newpos, dd_ErrorType *error) /* 094 */
02628 {
02629   /* It returns the set of all redundant rows.  This should be called after all
02630      implicit linearity are recognized with dd_MatrixCanonicalizeLinearity.
02631   */
02632 
02633  
02634   dd_rowrange i,k,m,m1;
02635   dd_colrange d;
02636   dd_rowset redset1;
02637   dd_rowindex newpos1;
02638   dd_MatrixPtr M1=NULL;
02639   dd_Arow cvec; /* certificate */ 
02640   dd_boolean success=dd_FALSE, localdebug=dd_FALSE;
02641 
02642   m=(*M)->rowsize;
02643   set_initialize(redset, m);
02644   M1=dd_MatrixSortedUniqueCopy(*M,newpos);
02645   for (i=1; i<=m; i++){
02646     if ((*newpos)[i]<=0) set_addelem(*redset,i);
02647     if (localdebug) printf(" %ld:%ld",i,(*newpos)[i]);
02648   }
02649   if (localdebug) printf("\n");
02650   
02651   if ((*M)->representation==dd_Generator){
02652     d=((*M)->colsize)+1;
02653   } else {
02654     d=(*M)->colsize;
02655   }
02656   m1=M1->rowsize;
02657   if (localdebug){
02658     fprintf(stderr,"dd_MatrixRedundancyRemove: By sorting, %ld rows have been removed.  The remaining has %ld rows.\n",m-m1,m1);
02659     /* dd_WriteMatrix(stdout,M1);  */
02660   }
02661   dd_InitializeArow(d,&cvec); 
02662   set_initialize(&redset1, M1->rowsize);
02663   k=1;
02664   do {
02665     if (dd_RedundantExtensive(M1, k, cvec, &redset1,error)) {
02666       set_addelem(redset1, k);
02667       dd_MatrixRowsRemove2(&M1,redset1,&newpos1);
02668       for (i=1; i<=m; i++){
02669         if ((*newpos)[i]>0){
02670           if  (set_member((*newpos)[i],redset1)){
02671             set_addelem(*redset,i);
02672             (*newpos)[i]=0;  /* now the original row i is recognized redundant and removed from M1 */
02673           } else {
02674             (*newpos)[i]=newpos1[(*newpos)[i]];  /* update the new pos vector */
02675           }
02676         }
02677       }
02678       set_free(redset1);
02679       set_initialize(&redset1, M1->rowsize); 
02680       if (localdebug) {
02681         printf("dd_MatrixRedundancyRemove: the row %ld is redundant. The new matrix has %ld rows.\n", k, M1->rowsize);
02682         /* dd_WriteMatrix(stderr, M1);  */
02683       }
02684       free(newpos1);
02685     } else {
02686       if (set_card(redset1)>0) {
02687         dd_MatrixRowsRemove2(&M1,redset1,&newpos1);
02688         for (i=1; i<=m; i++){
02689           if ((*newpos)[i]>0){
02690             if  (set_member((*newpos)[i],redset1)){
02691               set_addelem(*redset,i);
02692               (*newpos)[i]=0;  /* now the original row i is recognized redundant and removed from M1 */
02693             } else {
02694               (*newpos)[i]=newpos1[(*newpos)[i]];  /* update the new pos vector */
02695             }
02696           }
02697         }
02698         set_free(redset1);
02699         set_initialize(&redset1, M1->rowsize);
02700         free(newpos1);
02701       }
02702       if (localdebug) {
02703         printf("dd_MatrixRedundancyRemove: the row %ld is essential. The new matrix has %ld rows.\n", k, M1->rowsize);
02704         /* dd_WriteMatrix(stderr, M1);  */
02705       }
02706       k=k+1;
02707     }
02708     if (*error!=dd_NoError) goto _L99;
02709   } while  (k<=M1->rowsize);
02710   if (localdebug) dd_WriteMatrix(stderr, M1);
02711   success=dd_TRUE;
02712   
02713 _L99:
02714   dd_FreeMatrix(*M);
02715   *M=M1;
02716   dd_FreeArow(d, cvec);
02717   set_free(redset1);
02718   return success;
02719 }
02720 
02721 
02722 dd_boolean dd_SRedundant(dd_MatrixPtr M, dd_rowrange itest, dd_Arow certificate, dd_ErrorType *error)  
02723   /* 093a */
02724 {
02725   /* Checks whether the row itest is strongly redundant for the representation.
02726      A row is strongly redundant in H-representation if every point in
02727      the polyhedron satisfies it with strict inequality.
02728      A row is strongly redundant in V-representation if this point is in
02729      the interior of the polyhedron.
02730      
02731      All linearity rows are not checked and considered NOT strongly redundant. 
02732      This code works for both H- and V-representations.  A certificate is
02733      given in the case of non-redundancy, showing a solution x violating only the itest
02734      inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
02735      separates the itest from the rest.  More explicitly, the LP to be setup is
02736 
02737      H-representation
02738        f* = minimize  
02739          b_itest     + A_itest x
02740        subject to
02741          b_itest + 1 + A_itest x     >= 0 (relaxed inequality to make an LP bounded)
02742          b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
02743          b_L         + A_L x = 0.  (linearity)
02744 
02745      V-representation (=separation problem)
02746        f* = minimize  
02747          b_itest x_0     + A_itest x
02748        subject to
02749          b_itest x_0     + A_itest x     >= -1 (to make an LP bounded)
02750          b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators except for itest in one side)
02751          b_L x_0         + A_L x = 0.  (linearity generators)
02752     
02753     Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
02754     and the row indices of input is partitioned into I and L where L is the set of linearity.
02755     In H-representation, the itest data is strongly redundant if and only if the optimal value f* is positive.
02756     In V-representation, the itest data is redundant if and only if the optimal value f* is zero (as the LP
02757     is homogeneous and the optimal value is always non-positive).  To recognize strong redundancy, one
02758     can set up a second LP
02759     
02760      V-representation (=boundary problem)
02761        g* = maximize  
02762          1^T b_{I-itest} x_0 + 1^T A_{I-itest}    (the sum of slacks)
02763        subject to
02764          b_itest x_0     + A_itest x      =  0 (the point has to lie on the boundary)
02765          b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators in one side)
02766          1^T b_{I-itest} x_0 + 1^T A_{I-itest} x <=  1 (to make an LP bounded)
02767          b_L x_0         + A_L x = 0.  (linearity generators)
02768          
02769     The redundant row is strongly redundant if and only if g* is zero.
02770 
02771     The certificate has dimension one more for V-representation case.
02772   */
02773 
02774   dd_colrange j;
02775   dd_LPPtr lp;
02776   dd_LPSolutionPtr lps;
02777   dd_ErrorType err=dd_NoError;
02778   dd_boolean answer=dd_FALSE,localdebug=dd_FALSE;
02779 
02780   *error=dd_NoError;
02781   if (set_member(itest, M->linset)){
02782     if (localdebug) printf("The %ld th row is linearity and strong redundancy checking is skipped.\n",itest);
02783     goto _L99;
02784   }
02785   
02786   /* Create an LP data for redundancy checking */
02787   if (M->representation==dd_Generator){
02788     lp=dd_CreateLP_V_Redundancy(M, itest);
02789   } else {
02790     lp=dd_CreateLP_H_Redundancy(M, itest);
02791   }
02792 
02793   dd_LPSolve(lp,dd_choiceRedcheckAlgorithm,&err);
02794   if (err!=dd_NoError){
02795     *error=err;
02796     goto _L999;
02797   } else {
02798     lps=dd_CopyLPSolution(lp);
02799 
02800     for (j=0; j<lps->d; j++) {
02801       dd_set(certificate[j], lps->sol[j]);
02802     }
02803 
02804     if (localdebug){
02805       printf("Optimum value:");
02806       dd_WriteNumber(stdout, lps->optvalue);
02807       printf("\n");
02808     }
02809 
02810     if (M->representation==dd_Inequality){
02811        if (dd_Positive(lps->optvalue)){
02812           answer=dd_TRUE;
02813           if (localdebug) fprintf(stderr,"==> %ld th inequality is strongly redundant.\n",itest);
02814         } else {
02815           answer=dd_FALSE;
02816           if (localdebug) fprintf(stderr,"==> %ld th inequality is not strongly redundant.\n",itest);
02817         } 
02818     } else {
02819        if (dd_Negative(lps->optvalue)){
02820           answer=dd_FALSE;
02821           if (localdebug) fprintf(stderr,"==> %ld th point is not strongly redundant.\n",itest);
02822         } else {
02823           /* for V-representation, we have to solve another LP */
02824           dd_FreeLPData(lp);
02825           dd_FreeLPSolution(lps);
02826           lp=dd_CreateLP_V_SRedundancy(M, itest);
02827           dd_LPSolve(lp,dd_DualSimplex,&err);
02828           lps=dd_CopyLPSolution(lp);
02829           if (localdebug) dd_WriteLPResult(stdout,lp,err);
02830           if (dd_Positive(lps->optvalue)){
02831             answer=dd_FALSE;
02832             if (localdebug) fprintf(stderr,"==> %ld th point is not strongly redundant.\n",itest);
02833           } else {
02834             answer=dd_TRUE;
02835             if (localdebug) fprintf(stderr,"==> %ld th point is strongly redundant.\n",itest);
02836           }
02837        }
02838     } 
02839     dd_FreeLPSolution(lps);
02840   }
02841   _L999:
02842   dd_FreeLPData(lp);
02843 _L99:
02844   return answer;
02845 }
02846 
02847 dd_rowset dd_SRedundantRows(dd_MatrixPtr M, dd_ErrorType *error)  /* 093a */
02848 {
02849   dd_rowrange i,m;
02850   dd_colrange d;
02851   dd_rowset redset;
02852   dd_MatrixPtr Mcopy;
02853   dd_Arow cvec; /* certificate */  
02854   dd_boolean localdebug=dd_FALSE;
02855 
02856   m=M->rowsize;
02857   if (M->representation==dd_Generator){
02858     d=(M->colsize)+1;
02859   } else {
02860     d=M->colsize;
02861   }
02862   Mcopy=dd_MatrixCopy(M);
02863   dd_InitializeArow(d,&cvec); 
02864   set_initialize(&redset, m);
02865   for (i=m; i>=1; i--) {
02866     if (dd_SRedundant(Mcopy, i, cvec, error)) {
02867       if (localdebug) printf("dd_SRedundantRows: the row %ld is strongly redundant.\n", i);
02868       set_addelem(redset, i);
02869       dd_MatrixRowRemove(&Mcopy, i);
02870     } else {
02871       if (localdebug) printf("dd_SRedundantRows: the row %ld is not strongly redundant.\n", i);
02872     }
02873     if (*error!=dd_NoError) goto _L99;
02874   }
02875 _L99:
02876   dd_FreeMatrix(Mcopy);
02877   dd_FreeArow(d, cvec);
02878   return redset;
02879 }
02880 
02881 dd_rowset dd_RedundantRowsViaShooting(dd_MatrixPtr M, dd_ErrorType *error)  /* 092 */
02882 {
02883   /* 
02884      For H-representation only and not quite reliable,
02885      especially when floating-point arithmetic is used.
02886      Use the ordinary (slower) method dd_RedundantRows.
02887   */
02888 
02889   dd_rowrange i,m, ired, irow=0;
02890   dd_colrange j,k,d;
02891   dd_rowset redset;
02892   dd_rowindex rowflag; 
02893     /* ith comp is negative if the ith inequality (i-1 st row) is redundant.
02894                    zero     if it is not decided.
02895                    k > 0    if it is nonredundant and assigned to the (k-1)th row of M1.
02896     */
02897   dd_MatrixPtr M1;
02898   dd_Arow shootdir, cvec=NULL;
02899   dd_LPPtr lp0, lp;
02900   dd_LPSolutionPtr lps; 
02901   dd_ErrorType err;
02902   dd_LPSolverType solver=dd_DualSimplex; 
02903   dd_boolean localdebug=dd_FALSE;
02904 
02905   m=M->rowsize;
02906   d=M->colsize;
02907   M1=dd_CreateMatrix(m,d);
02908   M1->rowsize=0;  /* cheat the rowsize so that smaller matrix can be stored */
02909   set_initialize(&redset, m);
02910   dd_InitializeArow(d, &shootdir);
02911   dd_InitializeArow(d, &cvec);
02912 
02913   rowflag=(long *)calloc(m+1, sizeof(long)); 
02914 
02915   /* First find some (likely) nonredundant inequalities by Interior Point Find. */
02916   lp0=dd_Matrix2LP(M, &err);
02917   lp=dd_MakeLPforInteriorFinding(lp0);
02918   dd_FreeLPData(lp0); 
02919   dd_LPSolve(lp, solver, &err);  /* Solve the LP */
02920   lps=dd_CopyLPSolution(lp);
02921 
02922   if (dd_Positive(lps->optvalue)){
02923     /* An interior point is found.  Use rayshooting to find some nonredundant
02924        inequalities. */
02925     for (j=1; j<d; j++){
02926       for (k=1; k<=d; k++) dd_set(shootdir[k-1], dd_purezero);
02927       dd_set(shootdir[j], dd_one);  /* j-th unit vector */
02928       ired=dd_RayShooting(M, lps->sol, shootdir);
02929       if (localdebug) printf("nonredundant row %3ld found by shooting.\n", ired);
02930       if (ired>0 && rowflag[ired]<=0) {
02931         irow++;
02932         rowflag[ired]=irow;
02933         for (k=1; k<=d; k++) dd_set(M1->matrix[irow-1][k-1], M->matrix[ired-1][k-1]); 
02934       }
02935         
02936       dd_neg(shootdir[j], dd_one);  /* negative of the j-th unit vector */
02937       ired=dd_RayShooting(M, lps->sol, shootdir);
02938       if (localdebug) printf("nonredundant row %3ld found by shooting.\n", ired);
02939       if (ired>0 && rowflag[ired]<=0) {
02940         irow++;
02941         rowflag[ired]=irow;
02942         for (k=1; k<=d; k++) dd_set(M1->matrix[irow-1][k-1], M->matrix[ired-1][k-1]); 
02943       }
02944     }
02945 
02946     M1->rowsize=irow;
02947     if (localdebug) {
02948       printf("The initial nonredundant set is:");
02949       for (i=1; i<=m; i++) if (rowflag[i]>0) printf(" %ld", i);
02950       printf("\n");
02951     }
02952     
02953     i=1;
02954     while(i<=m){
02955       if (rowflag[i]==0){ /* the ith inequality is not yet checked */
02956         if (localdebug) fprintf(stderr, "Checking redundancy of %ld th inequality\n", i);
02957         irow++;  M1->rowsize=irow;
02958         for (k=1; k<=d; k++) dd_set(M1->matrix[irow-1][k-1], M->matrix[i-1][k-1]);
02959         if (!dd_Redundant(M1, irow, cvec, &err)){
02960           for (k=1; k<=d; k++) dd_sub(shootdir[k-1], cvec[k-1], lps->sol[k-1]); 
02961           ired=dd_RayShooting(M, lps->sol, shootdir);
02962           rowflag[ired]=irow;
02963           for (k=1; k<=d; k++) dd_set(M1->matrix[irow-1][k-1], M->matrix[ired-1][k-1]);
02964           if (localdebug) {
02965             fprintf(stderr, "The %ld th inequality is nonredundant for the subsystem\n", i);
02966             fprintf(stderr, "The nonredundancy of %ld th inequality is found by shooting.\n", ired);
02967           }
02968         } else {
02969           if (localdebug) fprintf(stderr, "The %ld th inequality is redundant for the subsystem and thus for the whole.\n", i);
02970           rowflag[i]=-1;
02971           set_addelem(redset, i);
02972           i++;
02973         }
02974       } else {
02975         i++;
02976       }
02977     } /* endwhile */
02978   } else {
02979     /* No interior point is found.  Apply the standard LP technique.  */
02980     redset=dd_RedundantRows(M, error);
02981   }
02982 
02983   dd_FreeLPData(lp);
02984   dd_FreeLPSolution(lps);
02985 
02986   M1->rowsize=m; M1->colsize=d;  /* recover the original sizes */
02987   dd_FreeMatrix(M1);
02988   dd_FreeArow(d, shootdir);
02989   dd_FreeArow(d, cvec);
02990   free(rowflag);
02991   return redset;
02992 }
02993 
02994 dd_SetFamilyPtr dd_Matrix2Adjacency(dd_MatrixPtr M, dd_ErrorType *error)  /* 093 */
02995 {
02996   /* This is to generate the (facet) graph of a polyheron (H) V-represented by M using LPs.
02997      Since it does not use the representation conversion, it should work for a large
02998      scale problem.
02999   */
03000   dd_rowrange i,m;
03001   dd_colrange d;
03002   dd_rowset redset;
03003   dd_MatrixPtr Mcopy;
03004   dd_SetFamilyPtr F=NULL;
03005 
03006   m=M->rowsize;
03007   d=M->colsize;
03008   if (m<=0 ||d<=0) {
03009     *error=dd_EmptyRepresentation;
03010     goto _L999;
03011   }
03012   Mcopy=dd_MatrixCopy(M);
03013   F=dd_CreateSetFamily(m, m);
03014   for (i=1; i<=m; i++) {
03015     if (!set_member(i, M->linset)){
03016       set_addelem(Mcopy->linset, i);
03017       redset=dd_RedundantRows(Mcopy, error);  /* redset should contain all nonadjacent ones */
03018       set_uni(redset, redset, Mcopy->linset); /* all linearity elements should be nonadjacent */
03019       set_compl(F->set[i-1], redset); /* set the adjacency list of vertex i */
03020       set_delelem(Mcopy->linset, i);
03021       set_free(redset);
03022       if (*error!=dd_NoError) goto _L99;
03023     }
03024   }
03025 _L99:
03026   dd_FreeMatrix(Mcopy);
03027 _L999:
03028   return F;
03029 }
03030 
03031 dd_SetFamilyPtr dd_Matrix2WeakAdjacency(dd_MatrixPtr M, dd_ErrorType *error)  /* 093a */
03032 {
03033   /* This is to generate the weak-adjacency (facet) graph of a polyheron (H) V-represented by M using LPs.
03034      Since it does not use the representation conversion, it should work for a large
03035      scale problem.
03036   */
03037   dd_rowrange i,m;
03038   dd_colrange d;
03039   dd_rowset redset;
03040   dd_MatrixPtr Mcopy;
03041   dd_SetFamilyPtr F=NULL;
03042 
03043   m=M->rowsize;
03044   d=M->colsize;
03045   if (m<=0 ||d<=0) {
03046     *error=dd_EmptyRepresentation;
03047     goto _L999;
03048   }
03049   Mcopy=dd_MatrixCopy(M);
03050   F=dd_CreateSetFamily(m, m);
03051   for (i=1; i<=m; i++) {
03052     if (!set_member(i, M->linset)){
03053       set_addelem(Mcopy->linset, i);
03054       redset=dd_SRedundantRows(Mcopy, error);  /* redset should contain all weakly nonadjacent ones */
03055       set_uni(redset, redset, Mcopy->linset); /* all linearity elements should be nonadjacent */
03056       set_compl(F->set[i-1], redset); /* set the adjacency list of vertex i */
03057       set_delelem(Mcopy->linset, i);
03058       set_free(redset);
03059       if (*error!=dd_NoError) goto _L99;
03060     }
03061   }
03062 _L99:
03063   dd_FreeMatrix(Mcopy);
03064 _L999:
03065   return F;
03066 }
03067 
03068 
03069 dd_boolean dd_ImplicitLinearity(dd_MatrixPtr M, dd_rowrange itest, dd_Arow certificate, dd_ErrorType *error)  
03070   /* 092 */
03071 {
03072   /* Checks whether the row itest is implicit linearity for the representation.
03073      All linearity rows are not checked and considered non implicit linearity (dd_FALSE). 
03074      This code works for both H- and V-representations.  A certificate is
03075      given in the case of dd_FALSE, showing a feasible solution x satisfying the itest
03076      strict inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
03077      separates the itest from the rest.  More explicitly, the LP to be setup is
03078      the same thing as redundancy case but with maximization:
03079 
03080      H-representation
03081        f* = maximize  
03082          b_itest     + A_itest x
03083        subject to
03084          b_itest + 1 + A_itest x     >= 0 (relaxed inequality. This is not necessary but kept for simplicity of the code)
03085          b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
03086          b_L         + A_L x = 0.  (linearity)
03087 
03088      V-representation (=separation problem)
03089        f* = maximize  
03090          b_itest x_0     + A_itest x
03091        subject to
03092          b_itest x_0     + A_itest x     >= -1 (again, this is not necessary but kept for simplicity.)
03093          b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators except for itest in one side)
03094          b_L x_0         + A_L x = 0.  (linearity generators)
03095     
03096     Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
03097     and the row indices of input is partitioned into I and L where L is the set of linearity.
03098     In both cases, the itest data is implicit linearity if and only if the optimal value f* is nonpositive.
03099     The certificate has dimension one more for V-representation case.
03100   */
03101 
03102   dd_colrange j;
03103   dd_LPPtr lp;
03104   dd_LPSolutionPtr lps;
03105   dd_ErrorType err=dd_NoError;
03106   dd_boolean answer=dd_FALSE,localdebug=dd_FALSE;
03107 
03108   *error=dd_NoError;
03109   if (set_member(itest, M->linset)){
03110     if (localdebug) printf("The %ld th row is linearity and redundancy checking is skipped.\n",itest);
03111     goto _L99;
03112   }
03113   
03114   /* Create an LP data for redundancy checking */
03115   if (M->representation==dd_Generator){
03116     lp=dd_CreateLP_V_Redundancy(M, itest);
03117   } else {
03118     lp=dd_CreateLP_H_Redundancy(M, itest);
03119   }
03120 
03121   lp->objective = dd_LPmax;  /* the lp->objective is set by CreateLP* to LPmin */
03122   dd_LPSolve(lp,dd_choiceRedcheckAlgorithm,&err);
03123   if (err!=dd_NoError){
03124     *error=err;
03125     goto _L999;
03126   } else {
03127     lps=dd_CopyLPSolution(lp);
03128 
03129     for (j=0; j<lps->d; j++) {
03130       dd_set(certificate[j], lps->sol[j]);
03131     }
03132 
03133     if (lps->LPS==dd_Optimal && dd_EqualToZero(lps->optvalue)){
03134       answer=dd_TRUE;
03135       if (localdebug) fprintf(stderr,"==> %ld th data is an implicit linearity.\n",itest);
03136     } else {
03137       answer=dd_FALSE;
03138       if (localdebug) fprintf(stderr,"==> %ld th data is not an implicit linearity.\n",itest);
03139     }
03140     dd_FreeLPSolution(lps);
03141   }
03142   _L999:
03143   dd_FreeLPData(lp);
03144 _L99:
03145   return answer;
03146 }
03147 
03148 
03149 int dd_FreeOfImplicitLinearity(dd_MatrixPtr M, dd_Arow certificate, dd_rowset *imp_linrows, dd_ErrorType *error)  
03150   /* 092 */
03151 {
03152   /* Checks whether the matrix M constains any implicit linearity at all.
03153   It returns 1 if it is free of any implicit linearity.  This means that 
03154   the present linearity rows define the linearity correctly.  It returns
03155   nonpositive values otherwise.  
03156 
03157 
03158      H-representation
03159        f* = maximize    z
03160        subject to
03161          b_I  + A_I x - 1 z >= 0 
03162          b_L  + A_L x = 0  (linearity)
03163          z <= 1.
03164 
03165      V-representation (=separation problem)
03166        f* = maximize    z
03167        subject to
03168          b_I x_0 + A_I x - 1 z >= 0 (all nonlinearity generators in one side)
03169          b_L x_0 + A_L x  = 0  (linearity generators)
03170          z <= 1.
03171     
03172     Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
03173     and the row indices of input is partitioned into I and L where L is the set of linearity.
03174     In both cases, any implicit linearity exists if and only if the optimal value f* is nonpositive.
03175     The certificate has dimension one more for V-representation case.
03176   */
03177 
03178   dd_LPPtr lp;
03179   dd_rowrange i,m;
03180   dd_colrange j,d1;
03181   dd_ErrorType err=dd_NoError;
03182   dd_Arow cvec; /* certificate for implicit linearity */
03183 
03184   int answer=0,localdebug=dd_FALSE;
03185 
03186   *error=dd_NoError;
03187   /* Create an LP data for redundancy checking */
03188   if (M->representation==dd_Generator){
03189     lp=dd_CreateLP_V_ImplicitLinearity(M);
03190   } else {
03191     lp=dd_CreateLP_H_ImplicitLinearity(M);
03192   }
03193 
03194   dd_LPSolve(lp,dd_choiceRedcheckAlgorithm,&err);
03195   if (err!=dd_NoError){
03196     *error=err;
03197     goto _L999;
03198   } else {
03199 
03200     for (j=0; j<lp->d; j++) {
03201       dd_set(certificate[j], lp->sol[j]);
03202     }
03203 
03204     if (localdebug) dd_WriteLPResult(stderr,lp,err);
03205     
03206     /* *posset contains a set of row indices that are recognized as nonlinearity.  */
03207     if (localdebug) {
03208       fprintf(stderr,"==> The following variables are not implicit linearity:\n");
03209       set_fwrite(stderr, lp->posset_extra);
03210       fprintf(stderr,"\n");
03211     }
03212     
03213     if (M->representation==dd_Generator){
03214       d1=(M->colsize)+1;
03215     } else {
03216       d1=M->colsize;
03217     }
03218     m=M->rowsize;
03219     dd_InitializeArow(d1,&cvec);
03220     set_initialize(imp_linrows,m);
03221 
03222     if (lp->LPS==dd_Optimal){
03223       if (dd_Positive(lp->optvalue)){
03224         answer=1;
03225         if (localdebug) fprintf(stderr,"==> The matrix has no implicit linearity.\n");
03226       } else if (dd_Negative(lp->optvalue)) {
03227           answer=-1;
03228           if (localdebug) fprintf(stderr,"==> The matrix defines the trivial system.\n");
03229         } else {
03230             answer=0;
03231             if (localdebug) fprintf(stderr,"==> The matrix has some implicit linearity.\n");
03232           }
03233     } else {
03234           answer=-2;
03235           if (localdebug) fprintf(stderr,"==> The LP fails.\n");
03236     }
03237     if (answer==0){
03238       /* List the implicit linearity rows */
03239       for (i=m; i>=1; i--) {
03240         if (!set_member(i,lp->posset_extra)) {
03241           if (dd_ImplicitLinearity(M, i, cvec, error)) {
03242             set_addelem(*imp_linrows, i);
03243             if (localdebug) {
03244               fprintf(stderr," row %ld is implicit linearity\n",i);
03245               fprintf(stderr,"\n");
03246             }
03247           }
03248           if (*error!=dd_NoError) goto _L999;
03249         }
03250       }
03251     }  /* end of if (answer==0) */
03252     if (answer==-1) {      
03253       for (i=m; i>=1; i--) set_addelem(*imp_linrows, i);
03254     } /* all rows are considered implicit linearity */
03255 
03256     dd_FreeArow(d1,cvec);
03257   }
03258   _L999:
03259   dd_FreeLPData(lp);
03260 
03261   return answer;
03262 }
03263 
03264 
03265 dd_rowset dd_ImplicitLinearityRows(dd_MatrixPtr M, dd_ErrorType *error)  /* 092 */
03266 {
03267   dd_colrange d;
03268   dd_rowset imp_linset;
03269   dd_Arow cvec; /* certificate */
03270   int foi;
03271   dd_boolean localdebug=dd_FALSE;
03272 
03273   if (M->representation==dd_Generator){
03274     d=(M->colsize)+2;
03275   } else {
03276     d=M->colsize+1;
03277   }
03278 
03279   dd_InitializeArow(d,&cvec);
03280   if (localdebug) fprintf(stdout, "\ndd_ImplicitLinearityRows: Check whether the system contains any implicit linearity.\n");
03281   foi=dd_FreeOfImplicitLinearity(M, cvec, &imp_linset, error);
03282   if (localdebug){
03283     switch (foi) {
03284       case 1:
03285         fprintf(stdout, "  It is free of implicit linearity.\n");
03286         break;
03287       
03288       case 0:
03289         fprintf(stdout, "  It is not free of implicit linearity.\n");
03290         break;
03291 
03292     case -1:
03293         fprintf(stdout, "  The input system is trivial (i.e. the empty H-polytope or the V-rep of the whole space.\n");
03294         break;
03295     
03296     default:
03297         fprintf(stdout, "  The LP was not solved correctly.\n");
03298         break;
03299     
03300     }
03301   }
03302   
03303   if (localdebug){
03304     fprintf(stderr, "  Implicit linearity rows are:\n");
03305     set_fwrite(stderr,imp_linset);
03306     fprintf(stderr, "\n");  
03307   }
03308   dd_FreeArow(d, cvec);
03309   return imp_linset;
03310 }
03311 
03312 dd_boolean dd_MatrixCanonicalizeLinearity(dd_MatrixPtr *M, dd_rowset *impl_linset,dd_rowindex *newpos, 
03313 dd_ErrorType *error) /* 094 */
03314 {
03315 /* This is to recongnize all implicit linearities, and put all linearities at the top of
03316    the matrix.    All implicit linearities will be returned by *impl_linset.
03317 */
03318   dd_rowrange rank;
03319   dd_rowset linrows,ignoredrows,basisrows;
03320   dd_colset ignoredcols,basiscols;
03321   dd_rowrange i,k,m;
03322   dd_rowindex newpos1;
03323   dd_boolean success=dd_FALSE;
03324   
03325   linrows=dd_ImplicitLinearityRows(*M, error);
03326   if (*error!=dd_NoError) goto _L99;
03327   
03328   m=(*M)->rowsize;
03329       
03330   set_uni((*M)->linset, (*M)->linset, linrows); 
03331       /* add the implicit linrows to the explicit linearity rows */
03332 
03333   /* To remove redundancy of the linearity part, 
03334      we need to compute the rank and a basis of the linearity part. */
03335   set_initialize(&ignoredrows,  (*M)->rowsize);
03336   set_initialize(&ignoredcols,  (*M)->colsize);
03337   set_compl(ignoredrows,  (*M)->linset);
03338   rank=dd_MatrixRank(*M,ignoredrows,ignoredcols,&basisrows,&basiscols);
03339   set_diff(ignoredrows,  (*M)->linset, basisrows);
03340   dd_MatrixRowsRemove2(M,ignoredrows,newpos);
03341   
03342   dd_MatrixShiftupLinearity(M,&newpos1); 
03343  
03344   for (i=1; i<=m; i++){
03345     k=(*newpos)[i];
03346     if (k>0) {
03347       (*newpos)[i]=newpos1[k];
03348     }
03349   }
03350   
03351   *impl_linset=linrows;
03352   success=dd_TRUE;
03353   free(newpos1);
03354   set_free(basisrows);
03355   set_free(basiscols);
03356   set_free(ignoredrows);
03357   set_free(ignoredcols);
03358 _L99:
03359   return success;
03360 }
03361 
03362 dd_boolean dd_MatrixCanonicalize(dd_MatrixPtr *M, dd_rowset *impl_linset, dd_rowset *redset, 
03363 dd_rowindex *newpos, dd_ErrorType *error) /* 094 */
03364 {
03365 /* This is to find a canonical representation of a matrix *M by 
03366    recognizing all implicit linearities and all redundancies.  
03367    All implicit linearities will be returned by *impl_linset and
03368    redundancies will be returned by *redset.
03369 */
03370   dd_rowrange i,k,m;
03371   dd_rowindex newpos1,revpos;
03372   dd_rowset redset1;
03373   dd_boolean success=dd_TRUE;
03374   
03375   m=(*M)->rowsize;
03376   set_initialize(redset, m);
03377   revpos=(long *)calloc(m+1,sizeof(long*));
03378   
03379   success=dd_MatrixCanonicalizeLinearity(M, impl_linset, newpos, error);
03380 
03381   if (!success) goto _L99;  
03382   
03383   for (i=1; i<=m; i++){
03384     k=(*newpos)[i];
03385     if (k>0) revpos[k]=i;  /* inverse of *newpos[] */
03386   }
03387  
03388   success=dd_MatrixRedundancyRemove(M, &redset1, &newpos1, error);  /* 094 */
03389   
03390   if (!success) goto _L99;
03391 
03392   for (i=1; i<=m; i++){
03393     k=(*newpos)[i];
03394     if (k>0) {
03395       (*newpos)[i]=newpos1[k];
03396       if (newpos1[k]<0) (*newpos)[i]=-revpos[-newpos1[k]];  /* update the certificate of its duplicate removal. */
03397       if (set_member(k,redset1)) set_addelem(*redset, i);
03398     }
03399   }
03400 
03401 _L99:
03402   set_free(redset1);
03403   free(newpos1);
03404   free(revpos);
03405   return success;
03406 }
03407 
03408 
03409 dd_boolean dd_ExistsRestrictedFace(dd_MatrixPtr M, dd_rowset R, dd_rowset S, dd_ErrorType *err)
03410 /* 0.94 */
03411 {
03412 /* This function checkes if there is a point that satifies all the constraints of
03413 the matrix M (interpreted as an H-representation) with additional equality contraints
03414 specified by R and additional strict inequality constraints specified by S.
03415 The set S is supposed to be disjoint from both R and M->linset.   When it is not,
03416 the set S will be considered as S\(R U M->linset).
03417 */
03418   dd_boolean answer=dd_FALSE;
03419   dd_LPPtr lp=NULL;
03420 
03421 /*
03422   printf("\n--- ERF ---\n");
03423   printf("R = "); set_write(R); 
03424   printf("S = "); set_write(S);
03425 */
03426   
03427   lp=dd_Matrix2Feasibility2(M, R, S, err);
03428 
03429   if (*err!=dd_NoError) goto _L99;
03430  
03431 /* Solve the LP by cdd LP solver. */
03432   dd_LPSolve(lp, dd_DualSimplex, err);  /* Solve the LP */
03433   if (*err!=dd_NoError) goto _L99;
03434   if (lp->LPS==dd_Optimal && dd_Positive(lp->optvalue)) {
03435     answer=dd_TRUE;
03436   } 
03437 
03438   dd_FreeLPData(lp);
03439 _L99:
03440   return answer;
03441 }
03442 
03443 dd_boolean dd_ExistsRestrictedFace2(dd_MatrixPtr M, dd_rowset R, dd_rowset S, dd_LPSolutionPtr *lps, dd_ErrorType *err)
03444 /* 0.94 */
03445 {
03446 /* This function checkes if there is a point that satifies all the constraints of
03447 the matrix M (interpreted as an H-representation) with additional equality contraints
03448 specified by R and additional strict inequality constraints specified by S.
03449 The set S is supposed to be disjoint from both R and M->linset.   When it is not,
03450 the set S will be considered as S\(R U M->linset).
03451 
03452 This function returns a certificate of the answer in terms of the associated LP solutions.
03453 */
03454   dd_boolean answer=dd_FALSE;
03455   dd_LPPtr lp=NULL;
03456 
03457 /*
03458   printf("\n--- ERF ---\n");
03459   printf("R = "); set_write(R); 
03460   printf("S = "); set_write(S);
03461 */
03462   
03463   lp=dd_Matrix2Feasibility2(M, R, S, err);
03464 
03465   if (*err!=dd_NoError) goto _L99;
03466  
03467 /* Solve the LP by cdd LP solver. */
03468   dd_LPSolve(lp, dd_DualSimplex, err);  /* Solve the LP */
03469   if (*err!=dd_NoError) goto _L99;
03470   if (lp->LPS==dd_Optimal && dd_Positive(lp->optvalue)) {
03471     answer=dd_TRUE;
03472   } 
03473 
03474 
03475   (*lps)=dd_CopyLPSolution(lp);
03476   dd_FreeLPData(lp);
03477 _L99:
03478   return answer;
03479 }
03480 
03481 dd_boolean dd_FindRelativeInterior(dd_MatrixPtr M, dd_rowset *ImL, dd_rowset *Lbasis, dd_LPSolutionPtr *lps, dd_ErrorType *err) 
03482 /* 0.94 */
03483 {
03484 /* This function computes a point in the relative interior of the H-polyhedron given by M.
03485 Even the representation is V-representation, it simply interprete M as H-representation.
03486 lps returns the result of solving an LP whose solution is a relative interior point.
03487 ImL returns all row indices of M that are implicit linearities, i.e. their inqualities
03488 are satisfied by equality by all points in the polyhedron.  Lbasis returns a row basis
03489 of the submatrix of M consisting of all linearities and implicit linearities.  This means
03490 that the dimension of the polyhedron is M->colsize - set_card(Lbasis) -1.
03491 */
03492 
03493   dd_rowset S;
03494   dd_colset T, Lbasiscols;
03495   dd_boolean success=dd_FALSE;
03496   dd_rowrange i;
03497   dd_colrange rank;
03498   
03499 
03500   *ImL=dd_ImplicitLinearityRows(M, err);
03501 
03502   if (*err!=dd_NoError) goto _L99;
03503 
03504   set_initialize(&S, M->rowsize);   /* the empty set */
03505   for (i=1; i <=M->rowsize; i++) {
03506         if (!set_member(i, M->linset) && !set_member(i, *ImL)){
03507           set_addelem(S, i);  /* all nonlinearity rows go to S  */
03508         }
03509   }
03510   if (dd_ExistsRestrictedFace2(M, *ImL, S, lps, err)){
03511     /* printf("a relative interior point found\n"); */
03512     success=dd_TRUE;
03513   }
03514   
03515   set_initialize(&T,  M->colsize); /* empty set */
03516   rank=dd_MatrixRank(M,S,T,Lbasis,&Lbasiscols); /* the rank of the linearity submatrix of M.  */
03517 
03518   set_free(S);
03519   set_free(T);
03520   set_free(Lbasiscols);
03521   
03522 _L99:
03523   return success;
03524 }
03525 
03526 
03527 dd_rowrange dd_RayShooting(dd_MatrixPtr M, dd_Arow p, dd_Arow r)
03528 {
03529 /* 092, find the first inequality "hit" by a ray from an intpt.  */
03530   dd_rowrange imin=-1,i,m;
03531   dd_colrange j, d;
03532   dd_Arow vecmin, vec;
03533   mytype min,t1,t2,alpha, t1min;  
03534   dd_boolean started=dd_FALSE;
03535   dd_boolean localdebug=dd_FALSE;
03536 
03537   m=M->rowsize;
03538   d=M->colsize;
03539   if (!dd_Equal(dd_one, p[0])){
03540     fprintf(stderr, "Warning: RayShooting is called with a point with first coordinate not 1.\n");
03541     dd_set(p[0],dd_one);
03542   }
03543   if (!dd_EqualToZero(r[0])){
03544     fprintf(stderr, "Warning: RayShooting is called with a direction with first coordinate not 0.\n");
03545     dd_set(r[0],dd_purezero);
03546   }
03547 
03548   dd_init(alpha); dd_init(min); dd_init(t1); dd_init(t2); dd_init(t1min);
03549   dd_InitializeArow(d,&vecmin);
03550   dd_InitializeArow(d,&vec);
03551 
03552   for (i=1; i<=m; i++){
03553     dd_InnerProduct(t1, d, M->matrix[i-1], p);
03554     if (dd_Positive(t1)) {
03555       dd_InnerProduct(t2, d, M->matrix[i-1], r);
03556       dd_div(alpha, t2, t1);
03557       if (!started){
03558         imin=i;  dd_set(min, alpha);
03559         dd_set(t1min, t1);  /* store the denominator. */
03560         started=dd_TRUE;
03561         if (localdebug) {
03562           fprintf(stderr," Level 1: imin = %ld and min = ", imin);
03563           dd_WriteNumber(stderr, min);
03564           fprintf(stderr,"\n");
03565         }
03566       } else {
03567         if (dd_Smaller(alpha, min)){
03568           imin=i;  dd_set(min, alpha);
03569           dd_set(t1min, t1);  /* store the denominator. */
03570           if (localdebug) {
03571             fprintf(stderr," Level 2: imin = %ld and min = ", imin);
03572             dd_WriteNumber(stderr, min);
03573             fprintf(stderr,"\n");
03574           }
03575         } else {
03576           if (dd_Equal(alpha, min)) { /* tie break */
03577             for (j=1; j<= d; j++){
03578               dd_div(vecmin[j-1], M->matrix[imin-1][j-1], t1min);
03579               dd_div(vec[j-1], M->matrix[i-1][j-1], t1);
03580             }
03581             if (dd_LexSmaller(vec,vecmin, d)){
03582               imin=i;  dd_set(min, alpha);
03583               dd_set(t1min, t1);  /* store the denominator. */
03584               if (localdebug) {
03585                 fprintf(stderr," Level 3: imin = %ld and min = ", imin);
03586                 dd_WriteNumber(stderr, min);
03587                 fprintf(stderr,"\n");
03588               }
03589             }
03590           }
03591         }
03592       }       
03593     }
03594   }
03595 
03596   dd_clear(alpha); dd_clear(min); dd_clear(t1); dd_clear(t2); dd_clear(t1min);
03597   dd_FreeArow(d, vecmin);
03598   dd_FreeArow(d, vec);
03599   return imin;
03600 }
03601 
03602 #ifdef GMPRATIONAL
03603 void dd_BasisStatusMaximize(dd_rowrange m_size,dd_colrange d_size,
03604     dd_Amatrix A,dd_Bmatrix T,dd_rowset equalityset,
03605     dd_rowrange objrow,dd_colrange rhscol,ddf_LPStatusType LPS,
03606     mytype *optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset, ddf_colindex nbindex,
03607     ddf_rowrange re,ddf_colrange se,long *pivots, int *found, int *LPScorrect)
03608 /*  This is just to check whether the status LPS of the basis given by 
03609 nbindex with extra certificates se or re is correct.  It is done
03610 by recomputing the basis inverse matrix T.  It does not solve the LP
03611 when the status *LPS is undecided.  Thus the input is
03612 m_size, d_size, A, equalityset, LPS, nbindex, re and se.
03613 Other values will be recomputed from scratch.
03614 
03615 The main purpose of the function is to verify the correctness
03616 of the result of floating point computation with the GMP rational
03617 arithmetics.
03618 */
03619 {
03620   long pivots0,pivots1;
03621   dd_rowrange i,is;
03622   dd_colrange s,j;
03623   static dd_rowindex bflag;
03624   static long mlast=0;
03625   static dd_rowindex OrderVector;  /* the permutation vector to store a preordered row indices */
03626   unsigned int rseed=1;
03627   mytype val;
03628   dd_colindex nbtemp;
03629   dd_LPStatusType ddlps;
03630   dd_boolean localdebug=dd_FALSE;
03631 
03632   if (dd_debug) localdebug=dd_debug;
03633   if (localdebug){
03634      printf("\nEvaluating dd_BasisStatusMaximize:\n");
03635   }
03636   dd_init(val);
03637   nbtemp=(long *) calloc(d_size+1,sizeof(long*));
03638   for (i=0; i<= 4; i++) pivots[i]=0;
03639   if (bflag==NULL || mlast!=m_size){
03640      if (mlast!=m_size && mlast>0) {
03641        free(bflag);   /* called previously with different m_size */
03642        free(OrderVector);
03643      }
03644      bflag=(long *) calloc(m_size+1,sizeof(long*));
03645      OrderVector=(long *)calloc(m_size+1,sizeof(long*)); 
03646      /* initialize only for the first time or when a larger space is needed */
03647       mlast=m_size;
03648   }
03649 
03650   /* Initializing control variables. */
03651   dd_ComputeRowOrderVector2(m_size,d_size,A,OrderVector,dd_MinIndex,rseed);
03652 
03653   pivots1=0;
03654 
03655   dd_ResetTableau(m_size,d_size,T,nbtemp,bflag,objrow,rhscol);
03656 
03657   if (localdebug){
03658      printf("\nnbindex:");
03659      for (j=1; j<=d_size; j++) printf(" %ld", nbindex[j]);
03660      printf("\n");
03661      printf("re = %ld,   se=%ld\n", re, se);
03662   }
03663   
03664   is=nbindex[se];
03665   if (localdebug) printf("se=%ld,  is=%ld\n", se, is);
03666   
03667   dd_FindLPBasis2(m_size,d_size,A,T,OrderVector, equalityset,nbindex,bflag,
03668       objrow,rhscol,&s,found,&pivots0);
03669       
03670   pivots[4]=pivots0;  /*GMP postopt pivots */
03671   dd_statBSpivots+=pivots0;
03672 
03673   if (!(*found)){
03674     if (localdebug) {
03675        printf("dd_BasisStatusMaximize: a specified basis DOES NOT exist.\n");
03676     }
03677 
03678        goto _L99;
03679      /* No speficied LP basis is found. */
03680   }
03681 
03682   if (localdebug) {
03683     printf("dd_BasisStatusMaximize: a specified basis exists.\n");
03684     if (m_size <=100 && d_size <=30)
03685     dd_WriteSignTableau(stdout,m_size,d_size,A,T,nbindex,bflag);
03686   }
03687 
03688   /* Check whether a recomputed basis is of the type specified by LPS */
03689   *LPScorrect=dd_TRUE;
03690   switch (LPS){
03691      case dd_Optimal: 
03692        for (i=1; i<=m_size; i++) {
03693          if (i!=objrow && bflag[i]==-1) {  /* i is a basic variable */
03694             dd_TableauEntry(&val,m_size,d_size,A,T,i,rhscol);
03695             if (dd_Negative(val)) {
03696                if (localdebug) printf("RHS entry for %ld is negative\n", i);
03697                *LPScorrect=dd_FALSE;
03698                break;
03699             }
03700           } else if (bflag[i] >0) { /* i is nonbasic variable */
03701             dd_TableauEntry(&val,m_size,d_size,A,T,objrow,bflag[i]);
03702             if (dd_Positive(val)) {
03703                if (localdebug) printf("Reduced cost entry for %ld is positive\n", i);
03704                *LPScorrect=dd_FALSE;
03705                break;
03706             }
03707           }
03708        };
03709        break;
03710      case dd_Inconsistent: 
03711        for (j=1; j<=d_size; j++){
03712           dd_TableauEntry(&val,m_size,d_size,A,T,re,j);
03713           if (j==rhscol){
03714              if (dd_Nonnegative(val)){
03715                if (localdebug) printf("RHS entry for %ld is nonnegative\n", re);
03716                *LPScorrect=dd_FALSE;
03717                break;             
03718              }
03719            } else if (dd_Positive(val)){
03720                if (localdebug) printf("the row entry for(%ld, %ld) is positive\n", re, j);
03721                *LPScorrect=dd_FALSE;
03722                break;             
03723            }
03724        };
03725        break;
03726      case dd_DualInconsistent:
03727         for (i=1; i<=m_size; i++){
03728           dd_TableauEntry(&val,m_size,d_size,A,T,i,bflag[is]);
03729           if (i==objrow){
03730              if (dd_Nonpositive(val)){
03731                if (localdebug) printf("Reduced cost entry for %ld is nonpositive\n", bflag[is]);
03732                *LPScorrect=dd_FALSE;
03733                break;             
03734              }
03735            } else if (dd_Negative(val)){
03736                if (localdebug) printf("the column entry for(%ld, %ld) is positive\n", i, bflag[is]);
03737                *LPScorrect=dd_FALSE;
03738                break;             
03739            }
03740        };
03741        break;
03742 ;
03743      default: break;
03744   }
03745 
03746   ddlps=LPSf2LPS(LPS);
03747 
03748   dd_SetSolutions(m_size,d_size,A,T,
03749    objrow,rhscol,ddlps,optvalue,sol,dsol,posset,nbindex,re,se,bflag);
03750 
03751   
03752 _L99:
03753   dd_clear(val);
03754   free(nbtemp);
03755 }
03756 
03757 void dd_BasisStatusMinimize(dd_rowrange m_size,dd_colrange d_size,
03758     dd_Amatrix A,dd_Bmatrix T,dd_rowset equalityset,
03759     dd_rowrange objrow,dd_colrange rhscol,ddf_LPStatusType LPS,
03760     mytype *optvalue,dd_Arow sol,dd_Arow dsol, dd_rowset posset, ddf_colindex nbindex,
03761     ddf_rowrange re,ddf_colrange se,long *pivots, int *found, int *LPScorrect)
03762 {
03763    dd_colrange j;
03764    
03765    for (j=1; j<=d_size; j++) dd_neg(A[objrow-1][j-1],A[objrow-1][j-1]);
03766    dd_BasisStatusMaximize(m_size,d_size,A,T,equalityset, objrow,rhscol,
03767      LPS,optvalue,sol,dsol,posset,nbindex,re,se,pivots,found,LPScorrect);
03768    dd_neg(*optvalue,*optvalue);
03769    for (j=1; j<=d_size; j++){
03770      dd_neg(dsol[j-1],dsol[j-1]);
03771      dd_neg(A[objrow-1][j-1],A[objrow-1][j-1]);
03772    }
03773 }
03774 #endif
03775 
03776 /* end of cddlp.c */
03777 
03778 #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