nbindex_ref=(long*) calloc(d_size+1,sizeof(long));    for (j=1; j<=d_size; j++){ dd_init(rcost[j-1]);}  }  d_last=d_size;  *err=dd_NoError; *lps=dd_LPSundecided; *s=0;  local_m_size=m_size+1;  /* increase m_size by 1 */  ms=0;  /* ms will be the index of column which has the largest reduced cost */  for (j=1; j<=d_size; j++){    if (j!=rhscol){      dd_TableauEntry(&(rcost[j-1]),local_m_size,d_size,A,T,objrow,j);      if (dd_Larger(rcost[j-1],maxcost)) {dd_set(maxcost,rcost[j-1]); ms = j;}    }  }  if (dd_Positive(maxcost)) dualfeasible=dd_FALSE;  if (!dualfeasible){    for (j=1; j<=d_size; j++){      dd_set(A[local_m_size-1][j-1], dd_purezero);      for (l=1; l<=d_size; l++){        if (nbindex[l]>0) {          dd_set_si(scaling,l+10);          dd_mul(svalue,A[nbindex[l]-1][j-1],scaling);           dd_sub(A[local_m_size-1][j-1],A[local_m_size-1][j-1],svalue);           /* To make the auxiliary row (0,-11,-12,...,-d-10).             It is likely to be better than  (0, -1, -1, ..., -1)             to avoid a degenerate LP.  Version 093c. */        }      }    }        if (localdebug){      fprintf(stderr,"\ndd_FindDualFeasibleBasis: curruent basis is not dual feasible.\n");      fprintf(stderr,"because of the column %ld assoc. with var %ld   dual cost =",       ms,nbindex[ms]);      dd_WriteNumber(stderr, maxcost);      if (localdebug) {        if (m_size <=100 && d_size <=30){          printf("\ndd_FindDualFeasibleBasis: the starting dictionary.\n");          dd_WriteSignTableau(stdout,m_size+1,d_size,A,T,nbindex,bflag);        }      }    }        ms=0;      /* Ratio Test: ms will be now the index of column which has the largest reduced cost         over the auxiliary row entry */    for (j=1; j<=d_size; j++){      if ((j!=rhscol) && dd_Positive(rcost[j-1])){        dd_TableauEntry(&axvalue,local_m_size,d_size,A,T,local_m_size,j);        if (dd_Nonnegative(axvalue)) {          *err=dd_NumericallyInconsistent;            /* This should not happen as they are set negative above.  Quit the phase I.*/          if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected.\n");          goto _L99;        }        dd_neg(axvalue,axvalue);        dd_div(axvalue,rcost[j-1],axvalue);  /* axvalue is the negative of ratio that is to be maximized. */        if (dd_Larger(axvalue,maxratio)) {          dd_set(maxratio,axvalue);           ms = j;        }      }    }    if (ms==0) {      *err=dd_NumericallyInconsistent; /* This should not happen. Quit the phase I.*/      if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected.\n");      goto _L99;    }    /* Pivot on (local_m_size,ms) so that the dual basic solution becomes feasible */    dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,local_m_size,ms);    pivots_p1=pivots_p1+1;    if (localdebug) {      printf("\ndd_FindDualFeasibleBasis: Pivot on %ld %ld.\n",local_m_size,ms);    }  for (j=1; j<=d_size; j++) nbindex_ref[j]=nbindex[j];     /* set the reference basis to be the current feasible basis. */  if (localdebug){    fprintf(stderr, "Store the current feasible basis:");    for (j=1; j<=d_size; j++) fprintf(stderr, " %ld", nbindex_ref[j]);    fprintf(stderr, "\n");       if (m_size <=100 && d_size <=30)      dd_WriteSignTableau2(stdout,m_size+1,d_size,A,T,nbindex_ref,nbindex,bflag);  }    phase1=dd_TRUE; stop=dd_FALSE;    do {   /* Dual Simplex Phase I */      chosen=dd_FALSE; LPSphase1=dd_LPSundecided;      if (pivots_p1>maxpivots) {        *err=dd_LPCycling;        fprintf(stderr,"max number %ld of pivots performed in Phase I. Switch to the anticycling phase.\n", maxpivots);        goto _L99;  /* failure due to max no. of pivots performed */      }      dd_SelectDualSimplexPivot(local_m_size,d_size,phase1,A,T,OV,nbindex_ref,nbindex,bflag,        objrow,rhscol,lexicopivot,&r_val,&s_val,&chosen,&LPSphase1);      if (!chosen) {        /* The current dictionary is terminal.  There are two cases:           dd_TableauEntry(local_m_size,d_size,A,T,objrow,ms) is negative or zero.           The first case implies dual infeasible,           and the latter implies dual feasible but local_m_size is still in nonbasis.           We must pivot in the auxiliary variable local_m_size.         */        dd_TableauEntry(&x,local_m_size,d_size,A,T,objrow,ms);        if (dd_Negative(x)){          *err=dd_NoError; *lps=dd_DualInconsistent;  *s=ms;        }        if (localdebug) {          fprintf(stderr,"\ndd_FindDualFeasibleBasis: the auxiliary variable was forced to enter the basis (# pivots = %ld).\n",pivots_p1);          fprintf(stderr," -- objrow %ld, ms %ld entry: ",objrow,ms);          dd_WriteNumber(stderr, x); fprintf(stderr,"\n");          if (dd_Negative(x)){            fprintf(stderr,"->The basis is dual inconsistent. Terminate.\n");          } else {            fprintf(stderr,"->The basis is feasible. Go to phase II.\n");          }        }        dd_init(minval);        r_val=0;        for (i=1; i<=local_m_size; i++){          if (bflag[i]<0) {              /* i is basic and not the objective variable */            dd_TableauEntry(&val,local_m_size,d_size,A,T,i,ms);  /* auxiliary column*/            if (dd_Smaller(val, minval)) {              r_val=i;              dd_set(minval,val);            }          }        }        dd_clear(minval);                if (r_val==0) {          *err=dd_NumericallyInconsistent; /* This should not happen. Quit the phase I.*/          if (localdebug) fprintf(stderr,"dd_FindDualFeasibleBasis: Numerical Inconsistency detected (r_val is 0).\n");          goto _L99;        }        dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,r_val,ms);        pivots_p1=pivots_p1+1;        if (localdebug) {          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);          if (m_size <=100 && d_size <=30)            dd_WriteSignTableau2(stdout,m_size+1,d_size,A,T,nbindex_ref,nbindex,bflag);        }        stop=dd_TRUE;      } else {        dd_GaussianColumnPivot2(local_m_size,d_size,A,T,nbindex,bflag,r_val,s_val);          pivots_p1=pivots_p1+1;        if (localdebug) {          printf("\ndd_FindDualFeasibleBasis: make a %ld-th pivot on %ld  %ld\n",pivots_p1,r_val,s_val);          if (m_size <=100 && d_size <=30)            dd_WriteSignTableau2(stdout,local_m_size,d_size,A,T,nbindex_ref,nbindex,bflag);        }        if (bflag[local_m_size]<0) {          stop=dd_TRUE;           if (localdebug)             fprintf(stderr,"\nDualSimplex Phase I: the auxiliary variable entered the basis (# pivots = %ld).\nGo to phase II\n",pivots_p1);        }      }    } while(!stop);  }_L99:  *pivot_no=pivots_p1;  dd_statDS1pivots+=pivots_p1;  dd_clear(x); dd_clear(val); dd_clear(maxcost); dd_clear(maxratio);  dd_clear(scaling); dd_clear(svalue); dd_clear(axvalue);}void dd_DualSimplexMinimize(dd_LPPtr lp,dd_ErrorType *err){   dd_colrange j;   *err=dd_NoError;   for (j=1; j<=lp->d; j++)     dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);   dd_DualSimplexMaximize(lp,err);   dd_neg(lp->optvalue,lp->optvalue);   for (j=1; j<=lp->d; j++){     dd_neg(lp->dsol[j-1],lp->dsol[j-1]);     dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);   }}void dd_DualSimplexMaximize(dd_LPPtr lp,dd_ErrorType *err)/* When LP is inconsistent then lp->re returns the evidence row.When LP is dual-inconsistent then lp->se returns the evidence column.*/{  int stop,chosen,phase1,found;  long pivots_ds=0,pivots_p0=0,pivots_p1=0,pivots_pc=0,maxpivots,maxpivfactor=20;  dd_boolean localdebug=dd_FALSE,localdebug1=dd_FALSE;#if !defined GMPRATIONAL  long maxccpivots,maxccpivfactor=100;     /* criss-cross should not cycle, but with floating-point arithmetics, it happens       (very rarely).  Jorg Rambau reported such an LP, in August 2003.  Thanks Jorg!    */#endif  dd_rowrange i,r;  dd_colrange j,s;  static dd_rowindex bflag;  static long mlast=0,nlast=0;  static dd_rowindex OrderVector;  /* the permutation vector to store a preordered row indeces */  static dd_colindex nbindex_ref; /* to be used to store the initial feasible basis for lexico rule */  double redpercent=0,redpercent_prev=0,redgain=0;  unsigned int rseed=1;    /* *err=dd_NoError; */  if (dd_debug) localdebug=dd_debug;  set_emptyset(lp->redset_extra);  for (i=0; i<= 4; i++) lp->pivots[i]=0;  maxpivots=maxpivfactor*lp->d;  /* maximum pivots to be performed before cc pivot is applied. */#if !defined GMPRATIONAL  maxccpivots=maxccpivfactor*lp->d;  /* maximum pivots to be performed with emergency cc pivots. */#endif  if (mlast!=lp->m || nlast!=lp->d){     if (mlast>0) { /* called previously with different lp->m */       free(OrderVector);       free(bflag);       free(nbindex_ref);     }     OrderVector=(long *)calloc(lp->m+1,sizeof(*OrderVector));     bflag=(long *) calloc(lp->m+2,sizeof(*bflag));  /* one more element for an auxiliary variable  */     nbindex_ref=(long*) calloc(lp->d+1,sizeof(long));     mlast=lp->m;nlast=lp->d;  }  /* Initializing control variables. */  dd_ComputeRowOrderVector2(lp->m,lp->d,lp->A,OrderVector,dd_MinIndex,rseed);  lp->re=0; lp->se=0;    dd_ResetTableau(lp->m,lp->d,lp->B,lp->nbindex,bflag,lp->objrow,lp->rhscol);     dd_FindLPBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->equalityset,lp->nbindex,bflag,      lp->objrow,lp->rhscol,&s,&found,&(lp->LPS),&pivots_p0);  lp->pivots[0]=pivots_p0;  if (!found){     lp->se=s;     goto _L99;     /* No LP basis is found, and thus Inconsistent.       Output the evidence column. */  }  dd_FindDualFeasibleBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->nbindex,bflag,      lp->objrow,lp->rhscol,lp->lexicopivot,&s, err,&(lp->LPS),&pivots_p1, maxpivots);  lp->pivots[1]=pivots_p1;  for (j=1; j<=lp->d; j++) nbindex_ref[j]=lp->nbindex[j];     /* set the reference basis to be the current feasible basis. */  if (localdebug){    fprintf(stderr, "dd_DualSimplexMaximize: Store the current feasible basis:");    for (j=1; j<=lp->d; j++) fprintf(stderr, " %ld", nbindex_ref[j]);    fprintf(stderr, "\n");    if (lp->m <=100 && lp->d <=30)      dd_WriteSignTableau2(stdout,lp->m+1,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag);   }    if (*err==dd_LPCycling || *err==dd_NumericallyInconsistent){    if (localdebug) fprintf(stderr, "Phase I failed and thus switch to the Criss-Cross method\n");    dd_CrissCrossMaximize(lp,err);    return;  }  if (lp->LPS==dd_DualInconsistent){     lp->se=s;     goto _L99;     /* No dual feasible basis is found, and thus DualInconsistent.       Output the evidence column. */  }  /* Dual Simplex Method */  stop=dd_FALSE;  do {    chosen=dd_FALSE; lp->LPS=dd_LPSundecided; phase1=dd_FALSE;    if (pivots_ds<maxpivots) {      dd_SelectDualSimplexPivot(lp->m,lp->d,phase1,lp->A,lp->B,OrderVector,nbindex_ref,lp->nbindex,bflag,        lp->objrow,lp->rhscol,lp->lexicopivot,&r,&s,&chosen,&(lp->LPS));    }    if (chosen) {      pivots_ds=pivots_ds+1;      if (lp->redcheck_extensive) {        dd_GetRedundancyInformation(lp->m,lp->d,lp->A,lp->B,lp->nbindex, bflag, lp->redset_extra);        set_uni(lp->redset_accum, lp->redset_accum,lp->redset_extra);        redpercent=100*(double)set_card(lp->redset_extra)/(double)lp->m;        redgain=redpercent-redpercent_prev;        redpercent_prev=redpercent;        if (localdebug1){          fprintf(stderr,"\ndd_DualSimplexMaximize: Phase II pivot %ld on (%ld, %ld).\n",pivots_ds,r,s);          fprintf(stderr,"  redundancy %f percent: redset size = %ld\n",redpercent,set_card(lp->redset_extra));        }      }    }    if (!chosen && lp->LPS==dd_LPSundecided) {        if (localdebug1){         fprintf(stderr,"Warning: an emergency CC pivot in Phase II is performed\n");         /* In principle this should not be executed because we already have dual feasibility            attained and dual simplex pivot should have been chosen.  This might occur            under floating point computation, or the case of cycling.         */      if (localdebug && lp->m <=100 && lp->d <=30){          fprintf(stderr,"\ndd_DualSimplexMaximize: The current dictionary.\n");          dd_WriteSignTableau2(stdout,lp->m,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag);      }    }#if !defined GMPRATIONAL      if (pivots_pc>maxccpivots) {        *err=dd_LPCycling;        stop=dd_TRUE;        goto _L99;      }#endif            dd_SelectCrissCrossPivot(lp->m,lp->d,lp->A,lp->B,bflag,        lp->objrow,lp->rhscol,&r,&s,&chosen,&(lp->LPS));      if (chosen) pivots_pc=pivots_pc+1;    }    if (chosen) {      dd_GaussianColumnPivot2(lp->m,lp->d,lp->A,lp->B,lp->nbindex,bflag,r,s);      if (localdebug && lp->m <=100 && lp->d <=30){          fprintf(stderr,"\ndd_DualSimplexMaximize: The current dictionary.\n");          dd_WriteSignTableau2(stdout,lp->m,lp->d,lp->A,lp->B,nbindex_ref,lp->nbindex,bflag);      }    } else {      switch (lp->LPS){        case dd_Inconsistent: lp->re=r;        case dd_DualInconsistent: lp->se=s;        default: break;      }      stop=dd_TRUE;    }  } while(!stop);_L99:   lp->pivots[2]=pivots_ds;  lp->pivots[3]=pivots_pc;  dd_statDS2pivots+=pivots_ds;  dd_statACpivots+=pivots_pc;  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);}void dd_CrissCrossMinimize(dd_LPPtr lp,dd_ErrorType *err){   dd_colrange j;   *err=dd_NoError;   for (j=1; j<=lp->d; j++)     dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);   dd_CrissCrossMaximize(lp,err);   dd_neg(lp->optvalue,lp->optvalue);   for (j=1; j<=lp->d; j++){     dd_neg(lp->dsol[j-1],lp->dsol[j-1]);     dd_neg(lp->A[lp->objrow-1][j-1],lp->A[lp->objrow-1][j-1]);   }}void dd_CrissCrossMaximize(dd_LPPtr lp,dd_ErrorType *err)/* When LP is inconsistent then lp->re returns the evidence row.When LP is dual-inconsistent then lp->se returns the evidence column.*/{  int stop,chosen,found;  long pivots0,pivots1;#if !defined GMPRATIONAL  long maxpivots,maxpivfactor=1000;     /* criss-cross should not cycle, but with floating-point arithmetics, it happens       (very rarely).  Jorg Rambau reported such an LP, in August 2003.  Thanks Jorg!    */#endif  dd_rowrange i,r;  dd_colrange s;  static dd_rowindex bflag;  static long mlast=0;  static dd_rowindex OrderVector;  /* the permutation vector to store a preordered row indeces */  unsigned int rseed=1;  dd_colindex nbtemp;  *err=dd_NoError;#if !defined GMPRATIONAL  maxpivots=maxpivfactor*lp->d;  /* maximum pivots to be performed when floating-point arithmetics is used. */#endif  nbtemp=(long *) calloc(lp->d+1,sizeof(long*));  for (i=0; i<= 4; i++) lp->pivots[i]=0;