cddlp.c~

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

  if (bflag==NULL || mlast!=lp->m){
     if (mlast!=lp->m && mlast>0) {
       free(bflag);   /* called previously with different lp->m */
       free(OrderVector);
     }
     bflag=(long *) calloc(lp->m+1,sizeof(long*));
     OrderVector=(long *)calloc(lp->m+1,sizeof(long*)); 
     /* initialize only for the first time or when a larger space is needed */
     
     mlast=lp->m;
  }
  /* Initializing control variables. */
  dd_ComputeRowOrderVector2(lp->m,lp->d,lp->A,OrderVector,dd_MinIndex,rseed);

  lp->re=0; lp->se=0; pivots1=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),&pivots0);
  lp->pivots[0]=pivots0;

  if (!found){
     lp->se=s;
     goto _L99;
     /* No LP basis is found, and thus Inconsistent.  
     Output the evidence column. */
  }

  stop=dd_FALSE;
  do {   /* Criss-Cross Method */
#if !defined GMPRATIONAL
    if (pivots1>maxpivots) {
      *err=dd_LPCycling;
      fprintf(stderr,"max number %ld of pivots performed by the criss-cross method. Most likely due to the floating-point arithmetics error.\n", maxpivots);
      goto _L99;  /* failure due to max no. of pivots performed */
    }
#endif

    dd_SelectCrissCrossPivot(lp->m,lp->d,lp->A,lp->B,bflag,
       lp->objrow,lp->rhscol,&r,&s,&chosen,&(lp->LPS));
    if (chosen) {
      dd_GaussianColumnPivot2(lp->m,lp->d,lp->A,lp->B,lp->nbindex,bflag,r,s);
      pivots1++;
    } 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[1]=pivots1;
  dd_statCCpivots+=pivots1;
  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);
  free(nbtemp);
}

void dd_SetSolutions(dd_rowrange m_size,dd_colrange d_size,
   dd_Amatrix A,dd_Bmatrix T,
   dd_rowrange objrow,dd_colrange rhscol,dd_LPStatusType LPS,
   mytype *optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset, dd_colindex nbindex,
   dd_rowrange re,dd_colrange se,dd_rowindex bflag)
/* 
Assign the solution vectors to sol,dsol,*optvalue after solving
the LP.
*/
{
  dd_rowrange i;
  dd_colrange j;
  mytype x,sw;
  int localdebug=dd_FALSE;
  
  dd_init(x); dd_init(sw);
  switch (LPS){
  case dd_Optimal:
    for (j=1;j<=d_size; j++) {
      dd_set(sol[j-1],T[j-1][rhscol-1]);
      dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
      dd_neg(dsol[j-1],x);
      dd_TableauEntry(optvalue,m_size,d_size,A,T,objrow,rhscol);
      if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]); dd_WriteNumber(stderr, dsol[j-1]); }
    }
    for (i=1; i<=m_size; i++) {
      if (bflag[i]==-1) {  /* i is a basic variable */
        dd_TableauEntry(&x,m_size,d_size,A,T,i,rhscol);
        if (dd_Positive(x)) set_addelem(posset, i);
      }
    }

    break;
  case dd_Inconsistent:
    if (localdebug) fprintf(stderr,"DualSimplexSolve: LP is inconsistent.\n");
    for (j=1;j<=d_size; j++) {
      dd_set(sol[j-1],T[j-1][rhscol-1]);
      dd_TableauEntry(&x,m_size,d_size,A,T,re,j);
      dd_neg(dsol[j-1],x);
      if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]); dd_WriteNumber(stderr,dsol[j-1]);}
    }
    break;
  case dd_DualInconsistent:
    for (j=1;j<=d_size; j++) {
      dd_set(sol[j-1],T[j-1][se-1]);
      dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
      dd_neg(dsol[j-1],x);
      if (localdebug) {fprintf(stderr,"dsol[%ld]= \n",nbindex[j]);dd_WriteNumber(stderr,dsol[j-1]);}
    }
    if (localdebug) printf( "DualSimplexSolve: LP is dual inconsistent.\n");
    break;

  case dd_StrucDualInconsistent:
    dd_TableauEntry(&x,m_size,d_size,A,T,objrow,se);
    if (dd_Positive(x)) dd_set(sw,dd_one);
    else dd_neg(sw,dd_one);
    for (j=1;j<=d_size; j++) {
      dd_mul(sol[j-1],sw,T[j-1][se-1]);
      dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);
      dd_neg(dsol[j-1],x);
      if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]);dd_WriteNumber(stderr,dsol[j-1]);}
    }
    if (localdebug) fprintf(stderr,"DualSimplexSolve: LP is dual inconsistent.\n");
    break;

  default:break;
  }
  dd_clear(x); dd_clear(sw);
}


void dd_RandomPermutation2(dd_rowindex OV,long t,unsigned int seed)
{
  long k,j,ovj;
  double u,xk,r,rand_max=(double) RAND_MAX;
  int localdebug=dd_FALSE;

  srand(seed);
  for (j=t; j>1 ; j--) {
    r=rand();
    u=r/rand_max;
    xk=(double)(j*u +1);
    k=(long)xk;
    if (localdebug) fprintf(stderr,"u=%g, k=%ld, r=%g, randmax= %g\n",u,k,r,rand_max);
    ovj=OV[j];
    OV[j]=OV[k];
    OV[k]=ovj;
    if (localdebug) fprintf(stderr,"row %ld is exchanged with %ld\n",j,k); 
  }
}

void dd_ComputeRowOrderVector2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,
    dd_rowindex OV,dd_RowOrderType ho,unsigned int rseed)
{
  long i,itemp;
  
  OV[0]=0;
  switch (ho){
  case dd_MaxIndex:
    for(i=1; i<=m_size; i++) OV[i]=m_size-i+1;
    break;

  case dd_LexMin:
    for(i=1; i<=m_size; i++) OV[i]=i;
    dd_QuickSort(OV,1,m_size,A,d_size);
   break;

  case dd_LexMax:
    for(i=1; i<=m_size; i++) OV[i]=i;
    dd_QuickSort(OV,1,m_size,A,d_size);
    for(i=1; i<=m_size/2;i++){   /* just reverse the order */
      itemp=OV[i];
      OV[i]=OV[m_size-i+1];
      OV[m_size-i+1]=itemp;
    }
    break;

  case dd_RandomRow:
    for(i=1; i<=m_size; i++) OV[i]=i;
    if (rseed<=0) rseed=1;
    dd_RandomPermutation2(OV,m_size,rseed);
    break;

  case dd_MinIndex: 
    for(i=1; i<=m_size; i++) OV[i]=i;
    break;

  default: 
    for(i=1; i<=m_size; i++) OV[i]=i;
    break;
 }
}

void dd_SelectPreorderedNext2(dd_rowrange m_size,dd_colrange d_size,
    rowset excluded,dd_rowindex OV,dd_rowrange *hnext)
{
  dd_rowrange i,k;
  
  *hnext=0;
  for (i=1; i<=m_size && *hnext==0; i++){
    k=OV[i];
    if (!set_member(k,excluded)) *hnext=k ;
  }
}

#ifdef GMPRATIONAL

ddf_LPObjectiveType Obj2Obj(dd_LPObjectiveType obj)
{
   ddf_LPObjectiveType objf=ddf_LPnone;

   switch (obj) {
   case dd_LPnone: objf=ddf_LPnone; break;
   case dd_LPmax: objf=ddf_LPmax; break;
   case dd_LPmin: objf=ddf_LPmin; break;
   }
   return objf;
}

ddf_LPPtr dd_LPgmp2LPf(dd_LPPtr lp)
{
  dd_rowrange i;
  dd_colrange j;
  ddf_LPType *lpf;
  double val;
  dd_boolean localdebug=dd_FALSE;

  if (localdebug) fprintf(stderr,"Converting a GMP-LP to a float-LP.\n");
  
  lpf=ddf_CreateLPData(Obj2Obj(lp->objective), ddf_Real, lp->m, lp->d);
  lpf->Homogeneous = lp->Homogeneous;
  lpf->eqnumber=lp->eqnumber;  /* this records the number of equations */

  for (i = 1; i <= lp->m; i++) {
    if (set_member(i, lp->equalityset)) set_addelem(lpf->equalityset,i);    
          /* it is equality. Its reversed row will not be in this set */
      for (j = 1; j <= lp->d; j++) {
        val=mpq_get_d(lp->A[i-1][j-1]);
        ddf_set_d(lpf->A[i-1][j-1],val);
      }  /*of j*/
  }  /*of i*/

  return lpf;
}


#endif


dd_boolean dd_LPSolve(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType *err)
/* 
The current version of dd_LPSolve that solves an LP with floating-arithmetics first
and then with the specified arithimetics if it is GMP.

When LP is inconsistent then *re returns the evidence row.
When LP is dual-inconsistent then *se returns the evidence column.
*/
{
  int i;
  dd_boolean found=dd_FALSE;
#ifdef GMPRATIONAL
  ddf_LPPtr lpf;
  ddf_ErrorType errf;
  dd_boolean LPScorrect=dd_FALSE;
  dd_boolean localdebug=dd_FALSE;
#endif

  *err=dd_NoError;
  lp->solver=solver;
  
   time(&lp->starttime);

#ifndef GMPRATIONAL
  switch (lp->solver) {
    case dd_CrissCross:
      dd_CrissCrossSolve(lp,err);
      break;
    case dd_DualSimplex:
      dd_DualSimplexSolve(lp,err);
      break;
  }
#else
  lpf=dd_LPgmp2LPf(lp);
  switch (lp->solver) {
    case dd_CrissCross:
      ddf_CrissCrossSolve(lpf,&errf);    /* First, run with double float. */
	  if (errf==ddf_NoError){   /* 094a:  fix for a bug reported by Dima Pasechnik */
        dd_BasisStatus(lpf,lp, &LPScorrect);    /* Check the basis. */
	  } else {LPScorrect=dd_FALSE;}
      if (!LPScorrect) {
         if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");
         dd_CrissCrossSolve(lp,err);  /* Rerun with GMP if fails. */
      } else {
         if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");
      }
      break;
    case dd_DualSimplex:
      ddf_DualSimplexSolve(lpf,&errf);    /* First, run with double float. */
	  if (errf==ddf_NoError){   /* 094a:  fix for a bug reported by Dima Pasechnik */
        dd_BasisStatus(lpf,lp, &LPScorrect);    /* Check the basis. */
	  } else {LPScorrect=dd_FALSE;}
      if (!LPScorrect){
         if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");
         dd_DualSimplexSolve(lp,err);  /* Rerun with GMP if fails. */
         if (localdebug){
            printf("*total number pivots = %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]);
            ddf_WriteLPResult(stdout, lpf, errf);
            dd_WriteLP(stdout, lp);
         }
      } else {
         if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");
      }
      break;
  }
  ddf_FreeLPData(lpf);
#endif

  time(&lp->endtime);
  lp->total_pivots=0;
  for (i=0; i<=4; i++) lp->total_pivots+=lp->pivots[i];
  if (*err==dd_NoError) found=dd_TRUE;
  return found;
}


dd_boolean dd_LPSolve0(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType *err)
/* 
The original version of dd_LPSolve that solves an LP with specified arithimetics.

When LP is inconsistent then *re returns the evidence row.
When LP is dual-inconsistent then *se returns the evidence column.
*/
{
  int i;
  dd_boolean found=dd_FALSE;

  *err=dd_NoError;
  lp->solver=solver;
  time(&lp->starttime);

  switch (lp->solver) {
    case dd_CrissCross:
      dd_CrissCrossSolve(lp,err);
      break;
    case dd_DualSimplex:
      dd_DualSimplexSolve(lp,err);
      break;
  }

  time(&lp->endtime);
  lp->total_pivots=0;
  for (i=0; i<=4; i++) lp->total_pivots+=lp->pivots[i];
  if (*err==dd_NoError) found=dd_TRUE;
  return found;
}


dd_LPPtr dd_MakeLPforInteriorFinding(dd_LPPtr lp)
/* Delete the objective row,
   add an extra column with -1's to the matrix A,
   add an extra row with (bceil, 0,...,0,-1),
   add an objective row with (0,...,0,1), and 
   rows & columns, and change m_size and d_size accordingly, to output new_A.
  This sets up the LP:
  maximize      x_{d+1}
  s.t.    A x + x_{d+1}  <=  b
                x_{d+1}  <=  bm * bmax,
  where bm is set to 2 by default, and bmax=max{1, b[1],...,b[m_size]}.
  Note that the equalitions (linearity) in the input lp will be ignored.
*/
{
  dd_rowrange m;
  dd_colrange d;
  dd_NumberType numbtype;
  dd_LPObjectiveType obj;
  dd_LPType *lpnew;
  dd_rowrange i; 
  dd_colrange j;
  mytype bm,bmax,bceil;
  int localdebug=dd_FALSE;

  dd_init(bm); dd_init(bmax); dd_init(bceil);
  dd_add(bm,dd_one,dd_one); dd_set(bmax,dd_one);
  numbtype=lp->numbtype;
  m=lp->m+1;
  d=lp->d+1;
  obj=dd_LPmax;

  lpnew=dd_CreateLPData(obj, numbtype, m, d);

  for (i=1; i<=lp->m; i++) {
    if (dd_Larger(lp->A[i-1][lp->rhscol-1],bmax)) 
      dd_set(bmax,lp->A[i-1][lp->rhscol-1]);
  }
  dd_mul(bceil,bm,bmax);
  if (lo