cddlp.c

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

/* cddlp.c:  dual simplex method c-code
   written by Komei Fukuda, [e-mail]
   Version 0.94a, Aug. 24, 2005
*/

/* cddlp.c : C-Implementation of the dual simplex method for
   solving an LP: max/min  A_(m-1).x subject to  x in P, where
   P= {x :  A_i.x >= 0, i=0,...,m-2, and  x_0=1}, and
   A_i is the i-th row of an m x n matrix A.
   Please read COPYING (GNU General Public Licence) and
   the manual cddlibman.tex for detail.
*/

#include "setoper.h"  /* set operation library header (Ver. May 18, 2000 or later) */
#include "cdd.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>

#if defined GMPRATIONAL
#include "cdd_f.h"
#endif

#define dd_CDDLPVERSION  "Version 0.94b (August 25, 2005)"

#define dd_FALSE 0
#define dd_TRUE 1

typedef set_type rowset;  /* set_type defined in setoper.h */
typedef set_type colset;

void dd_CrissCrossSolve(dd_LPPtr lp,dd_ErrorType *);
void dd_DualSimplexSolve(dd_LPPtr lp,dd_ErrorType *);
void dd_CrissCrossMinimize(dd_LPPtr,dd_ErrorType *);
void dd_CrissCrossMaximize(dd_LPPtr,dd_ErrorType *);
void dd_DualSimplexMinimize(dd_LPPtr,dd_ErrorType *);
void dd_DualSimplexMaximize(dd_LPPtr,dd_ErrorType *);
void dd_FindLPBasis(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowindex,dd_rowset,
    dd_colindex,dd_rowindex,dd_rowrange,dd_colrange,
    dd_colrange *,int *,dd_LPStatusType *,long *);
void dd_FindDualFeasibleBasis(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowindex,
    dd_colindex,long *,dd_rowrange,dd_colrange,dd_boolean,
    dd_colrange *,dd_ErrorType *,dd_LPStatusType *,long *, long maxpivots);


#ifdef GMPRATIONAL
void dd_BasisStatus(ddf_LPPtr lpf, dd_LPPtr lp, dd_boolean*);
void dd_BasisStatusMinimize(dd_rowrange,dd_colrange, dd_Amatrix,dd_Bmatrix,dd_rowset,
    dd_rowrange,dd_colrange,ddf_LPStatusType,mytype *,dd_Arow,dd_Arow,dd_rowset,ddf_colindex,
    ddf_rowrange,ddf_colrange,long *, int *, int *);
void dd_BasisStatusMaximize(dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,dd_rowset,
    dd_rowrange,dd_colrange,ddf_LPStatusType,mytype *,dd_Arow,dd_Arow,dd_rowset,ddf_colindex,
    ddf_rowrange,ddf_colrange,long *, int *, int *);
#endif

void dd_WriteBmatrix(FILE *f,dd_colrange d_size,dd_Bmatrix T);
void dd_SetNumberType(char *line,dd_NumberType *number,dd_ErrorType *Error);
void dd_ComputeRowOrderVector2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,
    dd_rowindex OV,dd_RowOrderType ho,unsigned int rseed);
void dd_SelectPreorderedNext2(dd_rowrange m_size,dd_colrange d_size,
    rowset excluded,dd_rowindex OV,dd_rowrange *hnext);
void dd_SetSolutions(dd_rowrange,dd_colrange,
   dd_Amatrix,dd_Bmatrix,dd_rowrange,dd_colrange,dd_LPStatusType,
   mytype *,dd_Arow,dd_Arow,dd_rowset,dd_colindex,dd_rowrange,dd_colrange,dd_rowindex);
   
void dd_WriteTableau(FILE *,dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,
  dd_colindex,dd_rowindex);

void dd_WriteSignTableau(FILE *,dd_rowrange,dd_colrange,dd_Amatrix,dd_Bmatrix,
  dd_colindex,dd_rowindex);


dd_LPSolutionPtr dd_CopyLPSolution(dd_LPPtr lp)
{
  dd_LPSolutionPtr lps;
  dd_colrange j;
  long i;

  lps=(dd_LPSolutionPtr) calloc(1,sizeof(dd_LPSolutionType));
  for (i=1; i<=dd_filenamelen; i++) lps->filename[i-1]=lp->filename[i-1];
  lps->objective=lp->objective;
  lps->solver=lp->solver; 
  lps->m=lp->m;
  lps->d=lp->d;
  lps->numbtype=lp->numbtype;

  lps->LPS=lp->LPS;  /* the current solution status */
  dd_init(lps->optvalue);
  dd_set(lps->optvalue,lp->optvalue);  /* optimal value */
  dd_InitializeArow(lp->d+1,&(lps->sol));
  dd_InitializeArow(lp->d+1,&(lps->dsol));
  lps->nbindex=(long*) calloc((lp->d)+1,sizeof(long));  /* dual solution */
  for (j=0; j<=lp->d; j++){
    dd_set(lps->sol[j],lp->sol[j]);
    dd_set(lps->dsol[j],lp->dsol[j]);
    lps->nbindex[j]=lp->nbindex[j];
  }
  lps->pivots[0]=lp->pivots[0];
  lps->pivots[1]=lp->pivots[1];
  lps->pivots[2]=lp->pivots[2];
  lps->pivots[3]=lp->pivots[3];
  lps->pivots[4]=lp->pivots[4];
  lps->total_pivots=lp->total_pivots;

  return lps;
}


dd_LPPtr dd_CreateLPData(dd_LPObjectiveType obj,
   dd_NumberType nt,dd_rowrange m,dd_colrange d)
{
  dd_LPType *lp;

  lp=(dd_LPPtr) calloc(1,sizeof(dd_LPType));
  lp->solver=dd_choiceLPSolverDefault;  /* set the default lp solver */
  lp->d=d;
  lp->m=m;
  lp->numbtype=nt;
  lp->objrow=m;
  lp->rhscol=1L;
  lp->objective=dd_LPnone;
  lp->LPS=dd_LPSundecided;
  lp->eqnumber=0;  /* the number of equalities */

  lp->nbindex=(long*) calloc(d+1,sizeof(long));
  lp->given_nbindex=(long*) calloc(d+1,sizeof(long));
  set_initialize(&(lp->equalityset),m);  
    /* i must be in the set iff i-th row is equality . */

  lp->redcheck_extensive=dd_FALSE; /* this is on only for RedundantExtensive */
  lp->ired=0;
  set_initialize(&(lp->redset_extra),m);  
    /* i is in the set if i-th row is newly recognized redundant (during the checking the row ired). */
  set_initialize(&(lp->redset_accum),m);  
    /* i is in the set if i-th row is recognized redundant (during the checking the row ired). */
   set_initialize(&(lp->posset_extra),m);  
    /* i is in the set if i-th row is recognized non-linearity (during the course of computation). */
 lp->lexicopivot=dd_choiceLexicoPivotQ;  /* dd_choice... is set in dd_set_global_constants() */

  lp->m_alloc=lp->m+2;
  lp->d_alloc=lp->d+2;
  lp->objective=obj;
  dd_InitializeBmatrix(lp->d_alloc,&(lp->B));
  dd_InitializeAmatrix(lp->m_alloc,lp->d_alloc,&(lp->A));
  dd_InitializeArow(lp->d_alloc,&(lp->sol));
  dd_InitializeArow(lp->d_alloc,&(lp->dsol));
  dd_init(lp->optvalue);
  return lp;
}


dd_LPPtr dd_Matrix2LP(dd_MatrixPtr M, dd_ErrorType *err)
{
  dd_rowrange m, i, irev, linc;
  dd_colrange d, j;
  dd_LPType *lp;
  dd_boolean localdebug=dd_FALSE;

  *err=dd_NoError;
  linc=set_card(M->linset);
  m=M->rowsize+1+linc; 
     /* We represent each equation by two inequalities.
        This is not the best way but makes the code simple. */
  d=M->colsize;
  if (localdebug) fprintf(stderr,"number of equalities = %ld\n", linc);
  
  lp=dd_CreateLPData(M->objective, M->numbtype, m, d);
  lp->Homogeneous = dd_TRUE;
  lp->eqnumber=linc;  /* this records the number of equations */

  irev=M->rowsize; /* the first row of the linc reversed inequalities. */
  for (i = 1; i <= M->rowsize; i++) {
    if (set_member(i, M->linset)) {
      irev=irev+1;
      set_addelem(lp->equalityset,i);    /* it is equality. */
                                         /* the reversed row irev is not in the equality set. */
      for (j = 1; j <= M->colsize; j++) {
        dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
      }  /*of j*/
      if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
    }
    for (j = 1; j <= M->colsize; j++) {
      dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
      if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
    }  /*of j*/
  }  /*of i*/
  for (j = 1; j <= M->colsize; j++) {
    dd_set(lp->A[m-1][j-1],M->rowvec[j-1]);  /* objective row */
  }  /*of j*/

  return lp;
}

dd_LPPtr dd_Matrix2Feasibility(dd_MatrixPtr M, dd_ErrorType *err)
/* Load a matrix to create an LP object for feasibility.   It is 
   essentially the dd_Matrix2LP except that the objject function
   is set to identically ZERO (maximization).
   
*/  
	 /*  094 */
{
  dd_rowrange m, linc;
  dd_colrange j;
  dd_LPType *lp;

  *err=dd_NoError;
  linc=set_card(M->linset);
  m=M->rowsize+1+linc; 
     /* We represent each equation by two inequalities.
        This is not the best way but makes the code simple. */
  
  lp=dd_Matrix2LP(M, err);
  lp->objective = dd_LPmax;   /* since the objective is zero, this is not important */
  for (j = 1; j <= M->colsize; j++) {
    dd_set(lp->A[m-1][j-1],dd_purezero);  /* set the objective to zero. */
  }  /*of j*/

  return lp;
}

dd_LPPtr dd_Matrix2Feasibility2(dd_MatrixPtr M, dd_rowset R, dd_rowset S, dd_ErrorType *err)
/* Load a matrix to create an LP object for feasibility with additional equality and
   strict inequality constraints given by R and S.  There are three types of inequalities:
   
   b_r + A_r x =  0     Linearity (Equations) specified by M
   b_s + A_s x >  0     Strict Inequalities specified by row index set S
   b_t + A_t x >= 0     The rest inequalities in M
   
   Where the linearity is considered here as the union of linearity specified by
   M and the additional set R.  When S contains any linearity rows, those
   rows are considered linearity (equation).  Thus S does not overlide linearity.
   To find a feasible solution, we set an LP
   
   maximize  z
   subject to
   b_r + A_r x     =  0      all r in Linearity
   b_s + A_s x - z >= 0      for all s in S
   b_t + A_t x     >= 0      for all the rest rows t
   1           - z >= 0      to make the LP bounded.
   
   Clearly, the feasibility problem has a solution iff the LP has a positive optimal value. 
   The variable z will be the last variable x_{d+1}.
   
*/  
/*  094 */
{
  dd_rowrange m, i, irev, linc;
  dd_colrange d, j;
  dd_LPType *lp;
  dd_rowset L;
  dd_boolean localdebug=dd_FALSE;

  *err=dd_NoError;
  set_initialize(&L, M->rowsize);
  set_uni(L,M->linset,R);
  linc=set_card(L);
  m=M->rowsize+1+linc+1; 
     /* We represent each equation by two inequalities.
        This is not the best way but makes the code simple. */
  d=M->colsize+1;
  if (localdebug) fprintf(stderr,"number of equalities = %ld\n", linc);
  
  lp=dd_CreateLPData(dd_LPmax, M->numbtype, m, d);
  lp->Homogeneous = dd_TRUE;
  lp->eqnumber=linc;  /* this records the number of equations */

  irev=M->rowsize; /* the first row of the linc reversed inequalities. */
  for (i = 1; i <= M->rowsize; i++) {
    if (set_member(i, L)) {
      irev=irev+1;
      set_addelem(lp->equalityset,i);    /* it is equality. */
                                         /* the reversed row irev is not in the equality set. */
      for (j = 1; j <= M->colsize; j++) {
        dd_neg(lp->A[irev-1][j-1],M->matrix[i-1][j-1]);
      }  /*of j*/
      if (localdebug) fprintf(stderr,"equality row %ld generates the reverse row %ld.\n",i,irev);
    } else if (set_member(i, S)) {
	  dd_set(lp->A[i-1][M->colsize],dd_minusone);
    }
    for (j = 1; j <= M->colsize; j++) {
      dd_set(lp->A[i-1][j-1],M->matrix[i-1][j-1]);
      if (j==1 && i<M->rowsize && dd_Nonzero(M->matrix[i-1][j-1])) lp->Homogeneous = dd_FALSE;
    }  /*of j*/
  }  /*of i*/
  for (j = 1; j <= d; j++) {
    dd_set(lp->A[m-2][j-1],dd_purezero);  /* initialize */
  }  /*of j*/
  dd_set(lp->A[m-2][0],dd_one);  /* the bounding constraint. */
  dd_set(lp->A[m-2][M->colsize],dd_minusone);  /* the bounding constraint. */
  for (j = 1; j <= d; j++) {
    dd_set(lp->A[m-1][j-1],dd_purezero);  /* initialize */
  }  /*of j*/
  dd_set(lp->A[m-1][M->colsize],dd_one);  /* maximize  z */
  
  set_free(L);
  return lp;
}



void dd_FreeLPData(dd_LPPtr lp)
{
  if ((lp)!=NULL){
    dd_clear(lp->optvalue);
    dd_FreeArow(lp->d_alloc,lp->dsol);
    dd_FreeArow(lp->d_alloc,lp->sol);
    dd_FreeBmatrix(lp->d_alloc,lp->B);
    dd_FreeAmatrix(lp->m_alloc,lp->d_alloc,lp->A);
    set_free(lp->equalityset);
    set_free(lp->redset_extra);
    set_free(lp->redset_accum);
    set_free(lp->posset_extra);
    free(lp->nbindex);
    free(lp->given_nbindex);
    free(lp);
  }
}

void dd_FreeLPSolution(dd_LPSolutionPtr lps)
{
  if (lps!=NULL){
    free(lps->nbindex);
    dd_FreeArow(lps->d+1,lps->dsol);
    dd_FreeArow(lps->d+1,lps->sol);
    dd_clear(lps->optvalue);
    
    free(lps);
  }
}

int dd_LPReverseRow(dd_LPPtr lp, dd_rowrange i)
{
  dd_colrange j;
  int success=0;

  if (i>=1 && i<=lp->m){
    lp->LPS=dd_LPSundecided;
    for (j=1; j<=lp->d; j++) {
      dd_neg(lp->A[i-1][j-1],lp->A[i-1][j-1]);
      /* negating the i-th constraint of A */
    }
    success=1;
  }
  return success;
}

int dd_LPReplaceRow(dd_LPPtr lp, dd_rowrange i, dd_Arow a)
{
  dd_colrange j;
  int success=0;

  if (i>=1 && i<=lp->m){
    lp->LPS=dd_LPSundecided;
    for (j=1; j<=lp->d; j++) {
      dd_set(lp->A[i-1][j-1],a[j-1]);
      /* replacing the i-th constraint by a */
    }
    success=1;
  }
  return success;
}

dd_Arow dd_LPCopyRow(dd_LPPtr lp, dd_rowrange i)
{
  dd_colrange j;
  dd_Arow a;

  if (i>=1 && i<=lp->m){
    dd_InitializeArow(lp->d, &a);
    for (j=1; j<=lp->d; j++) {
      dd_set(a[j-1],lp->A[i-1][j-1]);
      /* copying the i-th row to a */
    }
  }
  return a;
}


void dd_SetNumberType(char *line,dd_NumberType *number,dd_ErrorType *Error)
{
  if (strncmp(line,"integer",7)==0) {
    *number = dd_Integer;
    return;
  }
  else if (strncmp(line,"rational",8)==0) {
    *number = dd_Rational;
    return;
  }
  else if (strncmp(line,"real",4)==0) {
    *number = dd_Real;
    return;
  }
  else { 
    *number=dd_Unknown;
    *Error=dd_ImproperInputFormat;
  }
}