/* cddlp.c:  dual simplex method c-code   written by Komei Fukuda, fukuda@ifor.math.ethz.ch   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 1typedef 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 GMPRATIONALvoid 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 *);#endifvoid 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;  }}