Public Member Functions | Data Fields | Private Member Functions | Private Attributes
simplex Class Reference

Linear Programming / Linear Optimization using Simplex - Algorithm. More...

#include <mpr_numeric.h>

Public Member Functions

 simplex (int rows, int cols)
 #rows should be >= m+2, #cols >= n+1 More...
 
 ~simplex ()
 
BOOLEAN mapFromMatrix (matrix m)
 
matrix mapToMatrix (matrix m)
 
intvecposvToIV ()
 
intveczrovToIV ()
 
void compute ()
 

Data Fields

int m
 
int n
 
int m1
 
int m2
 
int m3
 
int icase
 
int * izrov
 
int * iposv
 
mprfloat ** LiPM
 

Private Member Functions

 simplex (const simplex &)
 
void simp1 (mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax)
 
void simp2 (mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1)
 
void simp3 (mprfloat **a, int i1, int k1, int ip, int kp)
 

Private Attributes

int LiPM_cols
 
int LiPM_rows
 

Detailed Description

Linear Programming / Linear Optimization using Simplex - Algorithm.

On output, the tableau LiPM is indexed by two arrays of integers. ipsov[j] contains, for j=1..m, the number i whose original variable is now represented by row j+1 of LiPM. (left-handed vars in solution) (first row is the one with the objective function) izrov[j] contains, for j=1..n, the number i whose original variable x_i is now a right-handed variable, rep. by column j+1 of LiPM. These vars are all zero in the solution. The meaning of n<i<n+m1+m2 is the same as above.

Definition at line 194 of file mpr_numeric.h.

Constructor & Destructor Documentation

simplex::simplex ( int  rows,
int  cols 
)

#rows should be >= m+2, #cols >= n+1

Definition at line 985 of file mpr_numeric.cc.

986  : LiPM_cols(cols), LiPM_rows(rows)
987 {
988  int i;
989 
992 
993  LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) ); // LP matrix
994  for( i= 0; i < LiPM_rows; i++ )
995  {
996  // Mem must be allocated aligned, also for type double!
997  LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) );
998  }
999 
1000  iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
1001  izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
1002 
1003  m=n=m1=m2=m3=icase=0;
1004 
1005 #ifdef mprDEBUG_ALL
1006  Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols);
1007 #endif
1008 }
#define Print
Definition: emacs.cc:83
int LiPM_cols
Definition: mpr_numeric.h:225
#define omAlloc(size)
Definition: omAllocDecl.h:210
double mprfloat
Definition: mpr_global.h:17
int * iposv
Definition: mpr_numeric.h:203
int i
Definition: cfEzgcd.cc:123
mprfloat ** LiPM
Definition: mpr_numeric.h:205
int LiPM_rows
Definition: mpr_numeric.h:225
int icase
Definition: mpr_numeric.h:201
int * izrov
Definition: mpr_numeric.h:203
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define omAlloc0Aligned
Definition: omAllocDecl.h:274
simplex::~simplex ( )

Definition at line 1010 of file mpr_numeric.cc.

1011 {
1012  // clean up
1013  int i;
1014  for( i= 0; i < LiPM_rows; i++ )
1015  {
1016  omFreeSize( (void *) LiPM[i], LiPM_cols * sizeof(mprfloat) );
1017  }
1018  omFreeSize( (void *) LiPM, LiPM_rows * sizeof(mprfloat *) );
1019 
1020  omFreeSize( (void *) iposv, 2*LiPM_rows*sizeof(int) );
1021  omFreeSize( (void *) izrov, 2*LiPM_rows*sizeof(int) );
1022 }
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
int LiPM_cols
Definition: mpr_numeric.h:225
double mprfloat
Definition: mpr_global.h:17
int * iposv
Definition: mpr_numeric.h:203
int i
Definition: cfEzgcd.cc:123
mprfloat ** LiPM
Definition: mpr_numeric.h:205
int LiPM_rows
Definition: mpr_numeric.h:225
int * izrov
Definition: mpr_numeric.h:203
simplex::simplex ( const simplex )
private

Member Function Documentation

void simplex::compute ( )

Definition at line 1108 of file mpr_numeric.cc.

1109 {
1110  int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
1111  int *l1,*l2,*l3;
1112  mprfloat q1, bmax;
1113 
1114  if ( m != (m1+m2+m3) )
1115  {
1116  // error: bad input
1117  error(WarnS("simplex::compute: Bad input constraint counts!");)
1118  icase=-2;
1119  return;
1120  }
1121 
1122  l1= (int *) omAlloc0( (n+1) * sizeof(int) );
1123  l2= (int *) omAlloc0( (m+1) * sizeof(int) );
1124  l3= (int *) omAlloc0( (m+1) * sizeof(int) );
1125 
1126  nl1= n;
1127  for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k;
1128  nl2=m;
1129  for ( i=1; i<=m; i++ )
1130  {
1131  if ( LiPM[i+1][1] < 0.0 )
1132  {
1133  // error: bad input
1134  error(WarnS("simplex::compute: Bad input tableau!");)
1135  error(Warn("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);)
1136  icase=-2;
1137  // free mem l1,l2,l3;
1138  omFreeSize( (void *) l3, (m+1) * sizeof(int) );
1139  omFreeSize( (void *) l2, (m+1) * sizeof(int) );
1140  omFreeSize( (void *) l1, (n+1) * sizeof(int) );
1141  return;
1142  }
1143  l2[i]= i;
1144  iposv[i]= n+i;
1145  }
1146  for ( i=1; i<=m2; i++) l3[i]= 1;
1147  ir= 0;
1148  if (m2+m3)
1149  {
1150  ir=1;
1151  for ( k=1; k <= (n+1); k++ )
1152  {
1153  q1=0.0;
1154  for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k];
1155  LiPM[m+2][k]= -q1;
1156  }
1157 
1158  do
1159  {
1160  simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax);
1161  if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS )
1162  {
1163  icase= -1; // no solution found
1164  // free mem l1,l2,l3;
1165  omFreeSize( (void *) l3, (m+1) * sizeof(int) );
1166  omFreeSize( (void *) l2, (m+1) * sizeof(int) );
1167  omFreeSize( (void *) l1, (n+1) * sizeof(int) );
1168  return;
1169  }
1170  else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS )
1171  {
1172  m12= m1+m2+1;
1173  if ( m12 <= m )
1174  {
1175  for ( ip= m12; ip <= m; ip++ )
1176  {
1177  if ( iposv[ip] == (ip+n) )
1178  {
1179  simp1(LiPM,ip,l1,nl1,1,&kp,&bmax);
1180  if ( fabs(bmax) >= SIMPLEX_EPS)
1181  goto one;
1182  }
1183  }
1184  }
1185  ir= 0;
1186  --m12;
1187  if ( m1+1 <= m12 )
1188  for ( i=m1+1; i <= m12; i++ )
1189  if ( l3[i-m1] == 1 )
1190  for ( k=1; k <= n+1; k++ )
1191  LiPM[i+1][k] = -(LiPM[i+1][k]);
1192  break;
1193  }
1194  //#if DEBUG
1195  //print_bmat( a, m+2, n+3);
1196  //#endif
1197  simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1198  if ( ip == 0 )
1199  {
1200  icase = -1; // no solution found
1201  // free mem l1,l2,l3;
1202  omFreeSize( (void *) l3, (m+1) * sizeof(int) );
1203  omFreeSize( (void *) l2, (m+1) * sizeof(int) );
1204  omFreeSize( (void *) l1, (n+1) * sizeof(int) );
1205  return;
1206  }
1207  one: simp3(LiPM,m+1,n,ip,kp);
1208  // #if DEBUG
1209  // print_bmat(a,m+2,n+3);
1210  // #endif
1211  if ( iposv[ip] >= (n+m1+m2+1))
1212  {
1213  for ( k= 1; k <= nl1; k++ )
1214  if ( l1[k] == kp ) break;
1215  --nl1;
1216  for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1];
1217  ++(LiPM[m+2][kp+1]);
1218  for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1219  }
1220  else
1221  {
1222  if ( iposv[ip] >= (n+m1+1) )
1223  {
1224  kh= iposv[ip]-m1-n;
1225  if ( l3[kh] )
1226  {
1227  l3[kh]= 0;
1228  ++(LiPM[m+2][kp+1]);
1229  for ( i=1; i<= m+2; i++ )
1230  LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1231  }
1232  }
1233  }
1234  is= izrov[kp];
1235  izrov[kp]= iposv[ip];
1236  iposv[ip]= is;
1237  } while (ir);
1238  }
1239  /* end of phase 1, have feasible sol, now optimize it */
1240  loop
1241  {
1242  // #if DEBUG
1243  // print_bmat( a, m+1, n+5);
1244  // #endif
1245  simp1(LiPM,0,l1,nl1,0,&kp,&bmax);
1246  if (bmax <= /*SIMPLEX_EPS*/0.0)
1247  {
1248  icase=0; // finite solution found
1249  // free mem l1,l2,l3
1250  omFreeSize( (void *) l3, (m+1) * sizeof(int) );
1251  omFreeSize( (void *) l2, (m+1) * sizeof(int) );
1252  omFreeSize( (void *) l1, (n+1) * sizeof(int) );
1253  return;
1254  }
1255  simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1256  if (ip == 0)
1257  {
1258  //printf("Unbounded:");
1259  // #if DEBUG
1260  // print_bmat( a, m+1, n+1);
1261  // #endif
1262  icase=1; /* unbounded */
1263  // free mem
1264  omFreeSize( (void *) l3, (m+1) * sizeof(int) );
1265  omFreeSize( (void *) l2, (m+1) * sizeof(int) );
1266  omFreeSize( (void *) l1, (n+1) * sizeof(int) );
1267  return;
1268  }
1269  simp3(LiPM,m,n,ip,kp);
1270  is= izrov[kp];
1271  izrov[kp]= iposv[ip];
1272  iposv[ip]= is;
1273  }/*for ;;*/
1274 }
loop
Definition: myNF.cc:98
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
int k
Definition: cfEzgcd.cc:93
#define WarnS
Definition: emacs.cc:81
double mprfloat
Definition: mpr_global.h:17
int * iposv
Definition: mpr_numeric.h:203
int i
Definition: cfEzgcd.cc:123
void simp2(mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1)
#define SIMPLEX_EPS
Definition: mpr_numeric.h:181
#define error(a)
Definition: mpr_numeric.cc:979
mprfloat ** LiPM
Definition: mpr_numeric.h:205
int icase
Definition: mpr_numeric.h:201
void simp1(mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax)
int * izrov
Definition: mpr_numeric.h:203
#define omAlloc0(size)
Definition: omAllocDecl.h:211
void simp3(mprfloat **a, int i1, int k1, int ip, int kp)
#define Warn
Definition: emacs.cc:80
BOOLEAN simplex::mapFromMatrix ( matrix  m)

Definition at line 1024 of file mpr_numeric.cc.

1025 {
1026  int i,j;
1027 // if ( MATROWS( m ) > LiPM_rows || MATCOLS( m ) > LiPM_cols ) {
1028 // WarnS("");
1029 // return FALSE;
1030 // }
1031 
1032  number coef;
1033  for ( i= 1; i <= MATROWS( mm ); i++ )
1034  {
1035  for ( j= 1; j <= MATCOLS( mm ); j++ )
1036  {
1037  if ( MATELEM(mm,i,j) != NULL )
1038  {
1039  coef= pGetCoeff( MATELEM(mm,i,j) );
1040  if ( coef != NULL && !nIsZero(coef) )
1041  LiPM[i][j]= (double)(*(gmp_float*)coef);
1042  //#ifdef mpr_DEBUG_PROT
1043  //Print("%f ",LiPM[i][j]);
1044  //#endif
1045  }
1046  }
1047  // PrintLn();
1048  }
1049 
1050  return TRUE;
1051 }
#define TRUE
Definition: auxiliary.h:144
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
int j
Definition: myNF.cc:70
int i
Definition: cfEzgcd.cc:123
#define MATCOLS(i)
Definition: matpol.h:28
#define nIsZero(n)
Definition: numbers.h:19
mprfloat ** LiPM
Definition: mpr_numeric.h:205
#define NULL
Definition: omList.c:10
#define MATROWS(i)
Definition: matpol.h:27
#define MATELEM(mat, i, j)
Definition: matpol.h:29
matrix simplex::mapToMatrix ( matrix  m)

Definition at line 1053 of file mpr_numeric.cc.

1054 {
1055  int i,j;
1056 // if ( MATROWS( mm ) < LiPM_rows-3 || MATCOLS( m ) < LiPM_cols-2 ) {
1057 // WarnS("");
1058 // return NULL;
1059 // }
1060 
1061 //Print(" %d x %d\n",MATROWS( mm ),MATCOLS( mm ));
1062 
1063  number coef;
1064  gmp_float * bla;
1065  for ( i= 1; i <= MATROWS( mm ); i++ )
1066  {
1067  for ( j= 1; j <= MATCOLS( mm ); j++ )
1068  {
1069  pDelete( &(MATELEM(mm,i,j)) );
1070  MATELEM(mm,i,j)= NULL;
1071 //Print(" %3.0f ",LiPM[i][j]);
1072  if ( LiPM[i][j] != 0.0 )
1073  {
1074  bla= new gmp_float(LiPM[i][j]);
1075  coef= (number)bla;
1076  MATELEM(mm,i,j)= pOne();
1077  pSetCoeff( MATELEM(mm,i,j), coef );
1078  }
1079  }
1080 //PrintLn();
1081  }
1082 
1083  return mm;
1084 }
int j
Definition: myNF.cc:70
int i
Definition: cfEzgcd.cc:123
#define pOne()
Definition: polys.h:286
#define MATCOLS(i)
Definition: matpol.h:28
mprfloat ** LiPM
Definition: mpr_numeric.h:205
#define NULL
Definition: omList.c:10
#define pDelete(p_ptr)
Definition: polys.h:157
#define MATROWS(i)
Definition: matpol.h:27
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define MATELEM(mat, i, j)
Definition: matpol.h:29
intvec * simplex::posvToIV ( )

Definition at line 1086 of file mpr_numeric.cc.

1087 {
1088  int i;
1089  intvec * iv = new intvec( m );
1090  for ( i= 1; i <= m; i++ )
1091  {
1092  IMATELEM(*iv,i,1)= iposv[i];
1093  }
1094  return iv;
1095 }
Definition: intvec.h:16
int * iposv
Definition: mpr_numeric.h:203
int i
Definition: cfEzgcd.cc:123
#define IMATELEM(M, I, J)
Definition: intvec.h:77
void simplex::simp1 ( mprfloat **  a,
int  mm,
int  ll[],
int  nll,
int  iabf,
int *  kp,
mprfloat bmax 
)
private

Definition at line 1276 of file mpr_numeric.cc.

1277 {
1278  int k;
1279  mprfloat test;
1280 
1281  if( nll <= 0)
1282  { /* init'tion: fixed */
1283  *bmax = 0.0;
1284  return;
1285  }
1286  *kp=ll[1];
1287  *bmax=a[mm+1][*kp+1];
1288  for (k=2;k<=nll;k++)
1289  {
1290  if (iabf == 0)
1291  {
1292  test=a[mm+1][ll[k]+1]-(*bmax);
1293  if (test > 0.0)
1294  {
1295  *bmax=a[mm+1][ll[k]+1];
1296  *kp=ll[k];
1297  }
1298  }
1299  else
1300  { /* abs values: have fixed it */
1301  test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
1302  if (test > 0.0)
1303  {
1304  *bmax=a[mm+1][ll[k]+1];
1305  *kp=ll[k];
1306  }
1307  }
1308  }
1309 }
const poly a
Definition: syzextra.cc:212
int k
Definition: cfEzgcd.cc:93
double mprfloat
Definition: mpr_global.h:17
CanonicalForm test
Definition: cfModGcd.cc:4037
void simplex::simp2 ( mprfloat **  a,
int  n,
int  l2[],
int  nl2,
int *  ip,
int  kp,
mprfloat q1 
)
private

Definition at line 1311 of file mpr_numeric.cc.

1312 {
1313  int k,ii,i;
1314  mprfloat qp,q0,q;
1315 
1316  *ip= 0;
1317  for ( i=1; i <= nl2; i++ )
1318  {
1319  if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS )
1320  {
1321  *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1];
1322  *ip= l2[i];
1323  for ( i= i+1; i <= nl2; i++ )
1324  {
1325  ii= l2[i];
1326  if (a[ii+1][kp+1] < -SIMPLEX_EPS)
1327  {
1328  q= -a[ii+1][1] / a[ii+1][kp+1];
1329  if (q - *q1 < -SIMPLEX_EPS)
1330  {
1331  *ip=ii;
1332  *q1=q;
1333  }
1334  else if (q - *q1 < SIMPLEX_EPS)
1335  {
1336  for ( k=1; k<= nn; k++ )
1337  {
1338  qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
1339  q0= -a[ii+1][k+1]/a[ii+1][kp+1];
1340  if ( q0 != qp ) break;
1341  }
1342  if ( q0 < qp ) *ip= ii;
1343  }
1344  }
1345  }
1346  }
1347  }
1348 }
const poly a
Definition: syzextra.cc:212
int k
Definition: cfEzgcd.cc:93
double mprfloat
Definition: mpr_global.h:17
int i
Definition: cfEzgcd.cc:123
#define SIMPLEX_EPS
Definition: mpr_numeric.h:181
void simplex::simp3 ( mprfloat **  a,
int  i1,
int  k1,
int  ip,
int  kp 
)
private

Definition at line 1350 of file mpr_numeric.cc.

1351 {
1352  int kk,ii;
1353  mprfloat piv;
1354 
1355  piv= 1.0 / a[ip+1][kp+1];
1356  for ( ii=1; ii <= i1+1; ii++ )
1357  {
1358  if ( ii -1 != ip )
1359  {
1360  a[ii][kp+1] *= piv;
1361  for ( kk=1; kk <= k1+1; kk++ )
1362  if ( kk-1 != kp )
1363  a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
1364  }
1365  }
1366  for ( kk=1; kk<= k1+1; kk++ )
1367  if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
1368  a[ip+1][kp+1]= piv;
1369 }
const poly a
Definition: syzextra.cc:212
double mprfloat
Definition: mpr_global.h:17
intvec * simplex::zrovToIV ( )

Definition at line 1097 of file mpr_numeric.cc.

1098 {
1099  int i;
1100  intvec * iv = new intvec( n );
1101  for ( i= 1; i <= n; i++ )
1102  {
1103  IMATELEM(*iv,i,1)= izrov[i];
1104  }
1105  return iv;
1106 }
Definition: intvec.h:16
int i
Definition: cfEzgcd.cc:123
int * izrov
Definition: mpr_numeric.h:203
#define IMATELEM(M, I, J)
Definition: intvec.h:77

Field Documentation

int simplex::icase

Definition at line 201 of file mpr_numeric.h.

int * simplex::iposv

Definition at line 203 of file mpr_numeric.h.

int* simplex::izrov

Definition at line 203 of file mpr_numeric.h.

mprfloat** simplex::LiPM

Definition at line 205 of file mpr_numeric.h.

int simplex::LiPM_cols
private

Definition at line 225 of file mpr_numeric.h.

int simplex::LiPM_rows
private

Definition at line 225 of file mpr_numeric.h.

int simplex::m

Definition at line 198 of file mpr_numeric.h.

int simplex::m1

Definition at line 200 of file mpr_numeric.h.

int simplex::m2

Definition at line 200 of file mpr_numeric.h.

int simplex::m3

Definition at line 200 of file mpr_numeric.h.

int simplex::n

Definition at line 199 of file mpr_numeric.h.


The documentation for this class was generated from the following files: