My Project
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() [1/2]

simplex::simplex ( int  rows,
int  cols 
)

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

Definition at line 971 of file mpr_numeric.cc.

972 : LiPM_cols(cols), LiPM_rows(rows)
973{
974 int i;
975
978
979 LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) ); // LP matrix
980 for( i= 0; i < LiPM_rows; i++ )
981 {
982 // Mem must be allocated aligned, also for type double!
983 LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) );
984 }
985
986 iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
987 izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
988
989 m=n=m1=m2=m3=icase=0;
990
991#ifdef mprDEBUG_ALL
992 Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols);
993#endif
994}
int i
Definition: cfEzgcd.cc:132
mprfloat ** LiPM
Definition: mpr_numeric.h:205
int * iposv
Definition: mpr_numeric.h:203
int LiPM_rows
Definition: mpr_numeric.h:225
int * izrov
Definition: mpr_numeric.h:203
int icase
Definition: mpr_numeric.h:201
int LiPM_cols
Definition: mpr_numeric.h:225
#define Print
Definition: emacs.cc:80
double mprfloat
Definition: mpr_global.h:17
#define omAlloc0Aligned
Definition: omAllocDecl.h:274
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omAlloc0(size)
Definition: omAllocDecl.h:211

◆ ~simplex()

simplex::~simplex ( )

Definition at line 996 of file mpr_numeric.cc.

997{
998 // clean up
999 int i;
1000 for( i= 0; i < LiPM_rows; i++ )
1001 {
1002 omFreeSize( (void *) LiPM[i], LiPM_cols * sizeof(mprfloat) );
1003 }
1004 omFreeSize( (void *) LiPM, LiPM_rows * sizeof(mprfloat *) );
1005
1006 omFreeSize( (void *) iposv, 2*LiPM_rows*sizeof(int) );
1007 omFreeSize( (void *) izrov, 2*LiPM_rows*sizeof(int) );
1008}
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260

◆ simplex() [2/2]

simplex::simplex ( const simplex )
private

Member Function Documentation

◆ compute()

void simplex::compute ( )

Definition at line 1094 of file mpr_numeric.cc.

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

◆ mapFromMatrix()

BOOLEAN simplex::mapFromMatrix ( matrix  m)

Definition at line 1010 of file mpr_numeric.cc.

1011{
1012 int i,j;
1013// if ( MATROWS( m ) > LiPM_rows || MATCOLS( m ) > LiPM_cols ) {
1014// WarnS("");
1015// return FALSE;
1016// }
1017
1018 number coef;
1019 for ( i= 1; i <= MATROWS( mm ); i++ )
1020 {
1021 for ( j= 1; j <= MATCOLS( mm ); j++ )
1022 {
1023 if ( MATELEM(mm,i,j) != NULL )
1024 {
1025 coef= pGetCoeff( MATELEM(mm,i,j) );
1026 if ( coef != NULL && !nIsZero(coef) )
1027 LiPM[i][j]= (double)(*(gmp_float*)coef);
1028 //#ifdef mpr_DEBUG_PROT
1029 //Print("%f ",LiPM[i][j]);
1030 //#endif
1031 }
1032 }
1033 // PrintLn();
1034 }
1035
1036 return TRUE;
1037}
#define NULL
Definition: auxiliary.h:104
#define TRUE
Definition: auxiliary.h:100
int j
Definition: facHensel.cc:110
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define nIsZero(n)
Definition: numbers.h:19

◆ mapToMatrix()

matrix simplex::mapToMatrix ( matrix  m)

Definition at line 1039 of file mpr_numeric.cc.

1040{
1041 int i,j;
1042// if ( MATROWS( mm ) < LiPM_rows-3 || MATCOLS( m ) < LiPM_cols-2 ) {
1043// WarnS("");
1044// return NULL;
1045// }
1046
1047//Print(" %d x %d\n",MATROWS( mm ),MATCOLS( mm ));
1048
1049 number coef;
1050 gmp_float * bla;
1051 for ( i= 1; i <= MATROWS( mm ); i++ )
1052 {
1053 for ( j= 1; j <= MATCOLS( mm ); j++ )
1054 {
1055 pDelete( &(MATELEM(mm,i,j)) );
1056 MATELEM(mm,i,j)= NULL;
1057//Print(" %3.0f ",LiPM[i][j]);
1058 if ( LiPM[i][j] != 0.0 )
1059 {
1060 bla= new gmp_float(LiPM[i][j]);
1061 coef= (number)bla;
1062 MATELEM(mm,i,j)= pOne();
1063 pSetCoeff( MATELEM(mm,i,j), coef );
1064 }
1065 }
1066//PrintLn();
1067 }
1068
1069 return mm;
1070}
#define pDelete(p_ptr)
Definition: polys.h:186
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define pOne()
Definition: polys.h:315

◆ posvToIV()

intvec * simplex::posvToIV ( )

Definition at line 1072 of file mpr_numeric.cc.

1073{
1074 int i;
1075 intvec * iv = new intvec( m );
1076 for ( i= 1; i <= m; i++ )
1077 {
1078 IMATELEM(*iv,i,1)= iposv[i];
1079 }
1080 return iv;
1081}
Definition: intvec.h:23
#define IMATELEM(M, I, J)
Definition: intvec.h:85

◆ simp1()

void simplex::simp1 ( mprfloat **  a,
int  mm,
int  ll[],
int  nll,
int  iabf,
int *  kp,
mprfloat bmax 
)
private

Definition at line 1262 of file mpr_numeric.cc.

1263{
1264 int k;
1265 mprfloat test;
1266
1267 if( nll <= 0)
1268 { /* init'tion: fixed */
1269 *bmax = 0.0;
1270 return;
1271 }
1272 *kp=ll[1];
1273 *bmax=a[mm+1][*kp+1];
1274 for (k=2;k<=nll;k++)
1275 {
1276 if (iabf == 0)
1277 {
1278 test=a[mm+1][ll[k]+1]-(*bmax);
1279 if (test > 0.0)
1280 {
1281 *bmax=a[mm+1][ll[k]+1];
1282 *kp=ll[k];
1283 }
1284 }
1285 else
1286 { /* abs values: have fixed it */
1287 test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
1288 if (test > 0.0)
1289 {
1290 *bmax=a[mm+1][ll[k]+1];
1291 *kp=ll[k];
1292 }
1293 }
1294 }
1295}
CanonicalForm test
Definition: cfModGcd.cc:4098

◆ simp2()

void simplex::simp2 ( mprfloat **  a,
int  n,
int  l2[],
int  nl2,
int *  ip,
int  kp,
mprfloat q1 
)
private

Definition at line 1297 of file mpr_numeric.cc.

1298{
1299 int k,ii,i;
1300 mprfloat qp,q0,q;
1301
1302 *ip= 0;
1303 for ( i=1; i <= nl2; i++ )
1304 {
1305 if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS )
1306 {
1307 *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1];
1308 *ip= l2[i];
1309 for ( i= i+1; i <= nl2; i++ )
1310 {
1311 ii= l2[i];
1312 if (a[ii+1][kp+1] < -SIMPLEX_EPS)
1313 {
1314 q= -a[ii+1][1] / a[ii+1][kp+1];
1315 if (q - *q1 < -SIMPLEX_EPS)
1316 {
1317 *ip=ii;
1318 *q1=q;
1319 }
1320 else if (q - *q1 < SIMPLEX_EPS)
1321 {
1322 for ( k=1; k<= nn; k++ )
1323 {
1324 qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
1325 q0= -a[ii+1][k+1]/a[ii+1][kp+1];
1326 if ( q0 != qp ) break;
1327 }
1328 if ( q0 < qp ) *ip= ii;
1329 }
1330 }
1331 }
1332 }
1333 }
1334}

◆ simp3()

void simplex::simp3 ( mprfloat **  a,
int  i1,
int  k1,
int  ip,
int  kp 
)
private

Definition at line 1336 of file mpr_numeric.cc.

1337{
1338 int kk,ii;
1339 mprfloat piv;
1340
1341 piv= 1.0 / a[ip+1][kp+1];
1342 for ( ii=1; ii <= i1+1; ii++ )
1343 {
1344 if ( ii -1 != ip )
1345 {
1346 a[ii][kp+1] *= piv;
1347 for ( kk=1; kk <= k1+1; kk++ )
1348 if ( kk-1 != kp )
1349 a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
1350 }
1351 }
1352 for ( kk=1; kk<= k1+1; kk++ )
1353 if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
1354 a[ip+1][kp+1]= piv;
1355}

◆ zrovToIV()

intvec * simplex::zrovToIV ( )

Definition at line 1083 of file mpr_numeric.cc.

1084{
1085 int i;
1086 intvec * iv = new intvec( n );
1087 for ( i= 1; i <= n; i++ )
1088 {
1089 IMATELEM(*iv,i,1)= izrov[i];
1090 }
1091 return iv;
1092}

Field Documentation

◆ icase

int simplex::icase

Definition at line 201 of file mpr_numeric.h.

◆ iposv

int * simplex::iposv

Definition at line 203 of file mpr_numeric.h.

◆ izrov

int* simplex::izrov

Definition at line 203 of file mpr_numeric.h.

◆ LiPM

mprfloat** simplex::LiPM

Definition at line 205 of file mpr_numeric.h.

◆ LiPM_cols

int simplex::LiPM_cols
private

Definition at line 225 of file mpr_numeric.h.

◆ LiPM_rows

int simplex::LiPM_rows
private

Definition at line 225 of file mpr_numeric.h.

◆ m

int simplex::m

Definition at line 198 of file mpr_numeric.h.

◆ m1

int simplex::m1

Definition at line 200 of file mpr_numeric.h.

◆ m2

int simplex::m2

Definition at line 200 of file mpr_numeric.h.

◆ m3

int simplex::m3

Definition at line 200 of file mpr_numeric.h.

◆ n

int simplex::n

Definition at line 199 of file mpr_numeric.h.


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