cddlp.c~

CheckMate is a MATLAB-based tool for modeling, simulating and investigating properties of hybrid dyn

  dd_GaussianColumnPivot(m_size,d_size,A,T,r,s);
  entering=nbindex[s];
  bflag[r]=s;     /* the nonbasic variable r corresponds to column s */
  nbindex[s]=r;   /* the nonbasic variable on s column is r */

  if (entering>0) bflag[entering]=-1;
     /* original variables have negative index and should not affect the row index */

  if (localdebug) {
    fprintf(stderr,"dd_GaussianColumnPivot2\n");
    fprintf(stderr," pivot: (leaving, entering) = (%ld, %ld)\n", r,entering);
    fprintf(stderr, " bflag[%ld] is set to %ld\n", r, s);
  }
}


void dd_ResetTableau(dd_rowrange m_size,dd_colrange d_size,dd_Bmatrix T,
    dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol)
{
  dd_rowrange i;
  dd_colrange j;
  
  /* Initialize T and nbindex */
  for (j=1; j<=d_size; j++) nbindex[j]=-j;
  nbindex[rhscol]=0; 
    /* RHS is already in nonbasis and is considered to be associated
       with the zero-th row of input. */
  dd_SetToIdentity(d_size,T);
  
  /* Set the bflag according to nbindex */
  for (i=1; i<=m_size; i++) bflag[i]=-1;  
    /* all basic variables have index -1 */
  bflag[objrow]= 0; 
    /* bflag of the objective variable is 0,
       different from other basic variables which have -1 */
  for (j=1; j<=d_size; j++) if (nbindex[j]>0) bflag[nbindex[j]]=j;
    /* bflag of a nonbasic variable is its column number */

}

void dd_SelectCrissCrossPivot(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
    dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
    dd_rowrange *r,dd_colrange *s,
    int *selected,dd_LPStatusType *lps)
{
  int colselected=dd_FALSE,rowselected=dd_FALSE;
  dd_rowrange i;
  mytype val;
  
  dd_init(val);
  *selected=dd_FALSE;
  *lps=dd_LPSundecided;
  while ((*lps==dd_LPSundecided) && (!rowselected) && (!colselected)) {
    for (i=1; i<=m_size; i++) {
      if (i!=objrow && bflag[i]==-1) {  /* i is a basic variable */
        dd_TableauEntry(&val,m_size,d_size,A,T,i,rhscol);
        if (dd_Negative(val)) {
          rowselected=dd_TRUE;
          *r=i;
          break;
        }
      }
      else if (bflag[i] >0) { /* i is nonbasic variable */
        dd_TableauEntry(&val,m_size,d_size,A,T,objrow,bflag[i]);
        if (dd_Positive(val)) {
          colselected=dd_TRUE;
          *s=bflag[i];
          break;
        }
      }
    }
    if  ((!rowselected) && (!colselected)) {
      *lps=dd_Optimal;
      return;
    }
    else if (rowselected) {
     for (i=1; i<=m_size; i++) {
       if (bflag[i] >0) { /* i is nonbasic variable */
          dd_TableauEntry(&val,m_size,d_size,A,T,*r,bflag[i]);
          if (dd_Positive(val)) {
            colselected=dd_TRUE;
            *s=bflag[i];
            *selected=dd_TRUE;
            break;
          }
        }
      }
    }
    else if (colselected) {
      for (i=1; i<=m_size; i++) {
        if (i!=objrow && bflag[i]==-1) {  /* i is a basic variable */
          dd_TableauEntry(&val,m_size,d_size,A,T,i,*s);
          if (dd_Negative(val)) {
            rowselected=dd_TRUE;
            *r=i;
            *selected=dd_TRUE;
            break;
          }
        }
      }
    }
    if (!rowselected) {
      *lps=dd_DualInconsistent;
    }
    else if (!colselected) {
      *lps=dd_Inconsistent;
    }
  }
  dd_clear(val);
}

void dd_CrissCrossSolve(dd_LPPtr lp, dd_ErrorType *err)
{
  switch (lp->objective) {
    case dd_LPmax:
         dd_CrissCrossMaximize(lp,err);
      break;
      
    case dd_LPmin:
         dd_CrissCrossMinimize(lp,err);
      break;

    case dd_LPnone: *err=dd_NoLPObjective; break;
  }

}

void dd_DualSimplexSolve(dd_LPPtr lp, dd_ErrorType *err)
{
  switch (lp->objective) {
    case dd_LPmax:
         dd_DualSimplexMaximize(lp,err);
      break;
      
    case dd_LPmin:
         dd_DualSimplexMinimize(lp,err);
      break;

    case dd_LPnone: *err=dd_NoLPObjective; break;
  }
}

#ifdef GMPRATIONAL

dd_LPStatusType LPSf2LPS(ddf_LPStatusType lpsf)
{
   dd_LPStatusType lps=dd_LPSundecided;

   switch (lpsf) {
   case ddf_LPSundecided: lps=dd_LPSundecided; break;
   case ddf_Optimal: lps=dd_Optimal; break;
   case ddf_Inconsistent: lps=dd_Inconsistent; break; 
   case ddf_DualInconsistent: lps=dd_DualInconsistent; break;
   case ddf_StrucInconsistent: lps=dd_StrucInconsistent; break; 
   case ddf_StrucDualInconsistent: lps=dd_StrucDualInconsistent; break;
   case ddf_Unbounded: lps=dd_Unbounded; break;
   case ddf_DualUnbounded: lps=dd_DualUnbounded; break;
   }
   return lps;
}


void dd_BasisStatus(ddf_LPPtr lpf, dd_LPPtr lp, dd_boolean *LPScorrect)
{
  int i;
  dd_colrange j;
  dd_boolean basisfound; 
 
  switch (lp->objective) {
    case dd_LPmax:
      dd_BasisStatusMaximize(lp->m,lp->d,lp->A,lp->B,lp->equalityset,lp->objrow,lp->rhscol,
           lpf->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lpf->nbindex,lpf->re,lpf->se,lp->pivots, 
           &basisfound, LPScorrect);
      if (*LPScorrect) {
         /* printf("BasisStatus Check: the current basis is verified with GMP\n"); */
         lp->LPS=LPSf2LPS(lpf->LPS);
         lp->re=lpf->re;
         lp->se=lpf->se;
         for (j=1; j<=lp->d; j++) lp->nbindex[j]=lpf->nbindex[j]; 
      }
      for (i=1; i<=5; i++) lp->pivots[i-1]+=lpf->pivots[i-1]; 
      break;
    case dd_LPmin:
      dd_BasisStatusMinimize(lp->m,lp->d,lp->A,lp->B,lp->equalityset,lp->objrow,lp->rhscol,
           lpf->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lpf->nbindex,lpf->re,lpf->se,lp->pivots, 
           &basisfound, LPScorrect);
      if (*LPScorrect) {
         /* printf("BasisStatus Check: the current basis is verified with GMP\n"); */
         lp->LPS=LPSf2LPS(lpf->LPS);
         lp->re=lpf->re;
         lp->se=lpf->se;
         for (j=1; j<=lp->d; j++) lp->nbindex[j]=lpf->nbindex[j]; 
      }
      for (i=1; i<=5; i++) lp->pivots[i-1]+=lpf->pivots[i-1]; 
      break;
    case dd_LPnone:  break;
   }      
}
#endif

void dd_FindLPBasis(dd_rowrange m_size,dd_colrange d_size,
    dd_Amatrix A, dd_Bmatrix T,dd_rowindex OV,dd_rowset equalityset, dd_colindex nbindex,
    dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
    dd_colrange *cs,int *found,dd_LPStatusType *lps,long *pivot_no)
{ 
  /* Find a LP basis using Gaussian pivots.
     If the problem has an LP basis,
     the procedure returns *found=dd_TRUE,*lps=LPSundecided and an LP basis.
     If the constraint matrix A (excluding the rhs and objective) is not
     column indepent, there are two cases.  If the dependency gives a dual
     inconsistency, this returns *found=dd_FALSE, *lps=dd_StrucDualInconsistent and 
     the evidence column *s.  Otherwise, this returns *found=dd_TRUE, 
     *lps=LPSundecided and an LP basis of size less than d_size.  Columns j
     that do not belong to the basis (i.e. cannot be chosen as pivot because
     they are all zero) will be indicated in nbindex vector: nbindex[j] will
     be negative and set to -j.
  */
  int chosen,stop;
  long pivots_p0=0,rank;
  colset ColSelected;
  rowset RowSelected;
  mytype val;

  dd_rowrange r;
  dd_colrange j,s;

  dd_init(val);
  *found=dd_FALSE; *cs=0; rank=0;
  *lps=dd_LPSundecided;

  set_initialize(&RowSelected,m_size);
  set_initialize(&ColSelected,d_size);
  set_addelem(RowSelected,objrow);
  set_addelem(ColSelected,rhscol);

  stop=dd_FALSE;
  do {   /* Find a LP basis */
    dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset,
      m_size,RowSelected,ColSelected,&r,&s,&chosen);
    if (chosen) {
      set_addelem(RowSelected,r);
      set_addelem(ColSelected,s);
      dd_GaussianColumnPivot2(m_size,d_size,A,T,nbindex,bflag,r,s);
      pivots_p0++;
      rank++;
    } else {
      for (j=1;j<=d_size  && *lps==dd_LPSundecided; j++) {
        if (j!=rhscol && nbindex[j]<0){
          dd_TableauEntry(&val,m_size,d_size,A,T,objrow,j);
          if (dd_Nonzero(val)){  /* dual inconsistent */
            *lps=dd_StrucDualInconsistent;
            *cs=j;
            /* dual inconsistent because the nonzero reduced cost */
          }
        }
      }
      if (*lps==dd_LPSundecided) *found=dd_TRUE;  
         /* dependent columns but not dual inconsistent. */
      stop=dd_TRUE;
    }
    if (rank==d_size-1) {
      stop = dd_TRUE;
      *found=dd_TRUE;
    }
  } while (!stop);

  *pivot_no=pivots_p0;
  dd_statBApivots+=pivots_p0;
  set_free(RowSelected);
  set_free(ColSelected);
  dd_clear(val);
}


void dd_FindLPBasis2(dd_rowrange m_size,dd_colrange d_size,
    dd_Amatrix A, dd_Bmatrix T,dd_rowindex OV,dd_rowset equalityset, dd_colindex nbindex,
    dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol,
    dd_colrange *cs,int *found,long *pivot_no)
{ 
  /* Similar to dd_FindLPBasis but it is much simpler.  This tries to recompute T for
  the specified basis given by nbindex.  It will return *found=dd_FALSE if the specified
  basis is not a basis.
  */
  int chosen,stop;
  long pivots_p0=0,rank;
  dd_colset ColSelected,DependentCols;
  dd_rowset RowSelected, NopivotRow;
  mytype val;
  dd_boolean localdebug=dd_FALSE;

  dd_rowrange r,negcount=0;
  dd_colrange j,s;

  dd_init(val);
  *found=dd_FALSE; *cs=0; rank=0;

  set_initialize(&RowSelected,m_size);
  set_initialize(&DependentCols,d_size);
  set_initialize(&ColSelected,d_size);
  set_initialize(&NopivotRow,m_size);
  set_addelem(RowSelected,objrow);
  set_addelem(ColSelected,rhscol);
  set_compl(NopivotRow, NopivotRow);  /* set NopivotRow to be the groundset */
  
  for (j=2; j<=d_size; j++) 
    if (nbindex[j]>0) 
       set_delelem(NopivotRow, nbindex[j]);
    else if (nbindex[j]<0){ 
       negcount++;       
       set_addelem(DependentCols, -nbindex[j]); 
       set_addelem(ColSelected, -nbindex[j]); 
    }
 
  set_uni(RowSelected, RowSelected, NopivotRow);  /* RowSelected is the set of rows not allowed to poviot on */

  stop=dd_FALSE;
  do {   /* Find a LP basis */
    dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset, m_size,RowSelected,ColSelected,&r,&s,&chosen);
    if (chosen) {
      set_addelem(RowSelected,r);
      set_addelem(ColSelected,s);

      dd_GaussianColumnPivot2(m_size,d_size,A,T,nbindex,bflag,r,s);
      if (localdebug && m_size <=10){
        dd_WriteBmatrix(stderr,d_size,T);
        dd_WriteTableau(stderr,m_size,d_size,A,T,nbindex,bflag);
      }
      pivots_p0++;
      rank++;
    } else{
      *found=dd_FALSE;   /* cannot pivot on any of the spacified positions. */
      stop=dd_TRUE;
    }
    if (rank==d_size-1-negcount) {
      if (negcount){
        /* Now it tries to pivot on rows that are supposed to be dependent. */ 
        set_diff(ColSelected, ColSelected, DependentCols); 
        dd_SelectPivot2(m_size,d_size,A,T,dd_MinIndex,OV,equalityset, m_size,RowSelected,ColSelected,&r,&s,&chosen);
        if (chosen) *found=dd_FALSE;  /* not supposed to be independent */
        else *found=dd_TRUE;
        if (localdebug){
          printf("Try to check the dependent cols:");
          set_write(DependentCols);
          if (chosen) printf("They are not dependent.  Can still pivot on (%ld, %ld)\n",r, s);
          else printf("They are indeed dependent.\n");
        }
      } else {
        *found=dd_TRUE;
     }   
     stop = dd_TRUE;
    }
  } while (!stop);

  for (j=1; j<=d_size; j++) if (nbindex[j]>0) bflag[nbindex[j]]=j;
  *pivot_no=pivots_p0;
  set_free(RowSelected);
  set_free(ColSelected);
  set_free(NopivotRow);
  set_free(DependentCols);
  dd_clear(val);
}

void dd_FindDualFeasibleBasis(dd_rowrange m_size,dd_colrange d_size,
    dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,
    dd_colindex nbindex,dd_rowindex bflag,dd_rowrange objrow,dd_colrange rhscol, dd_boolean lexicopivot,
    dd_colrange *s,dd_ErrorType *err,dd_LPStatusType *lps,long *pivot_no, long maxpivots)
{ 
  /* Find a dual feasible basis using Phase I of Dual Simplex method.
     If the problem is dual feasible,
     the procedure returns *err=NoError, *lps=LPSundecided and a dual feasible
     basis.   If the problem is dual infeasible, this returns
     *err=NoError, *lps=DualInconsistent and the evidence column *s.
     Caution: matrix A must have at least one extra row:  the row space A[m_size] must
     have been allocated.
  */
  dd_boolean phase1,dualfeasible=dd_TRUE;
  dd_boolean localdebug=dd_FALSE,chosen,stop;
  dd_LPStatusType LPSphase1;
  long pivots_p1=0;
  dd_rowrange i,r_val;
  dd_colrange j,l,ms=0,s_val,local_m_size;
  mytype x,val,maxcost,axvalue,maxratio;
  static dd_colrange d_last=0;
  static dd_Arow rcost;
  static dd_colindex nbindex_ref; /* to be used to store the initial feasible basis for lexico rule */

  mytype scaling,svalue;  /* random scaling mytype value */
  mytype minval;

  if (dd_debug) localdebug=dd_debug;
  dd_init(x); dd_init(val); dd_init(scaling); dd_init(svalue);  dd_init(axvalue);
  dd_init(maxcost);  dd_set(maxcost,dd_minuszero);
  dd_init(maxratio);  dd_set(maxratio,dd_minuszero);
  if (d_last<d_size) {
    if (d_last>0) {
      for (j=1; j<=d_last; j++){ dd_clear(rcost[j-1]);}
      free(rcost);
      free(nbindex_ref);
    }
    rcost=(mytype*) calloc(d_size,sizeof(mytype));