fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <string.h>
  5. #include <assert.h>
  6.  
  7. #define M 20
  8. #define N 20
  9.  
  10. static const double epsilon   = 1.0e-8;
  11. int equal(double a, double b) { return fabs(a-b) < epsilon; }
  12.  
  13. typedef struct {
  14.   int m, n; // m=rows, n=columns, mat[m x n]
  15.   double mat[M][N];
  16. } Tableau;
  17.  
  18. void nl(int k){ int j; for(j=0;j<k;j++) putchar('-'); putchar('\n'); }
  19.  
  20. void print_tableau(Tableau *tab, const char* mes) {
  21.   static int counter=0;
  22.   int i, j;
  23.   printf("\n%d. Tableau %s:\n", ++counter, mes);
  24.   nl(70);
  25.  
  26.   printf("%-6s%5s", "col:", "b[i]");
  27.   for(j=1;j<tab->n; j++) { printf("    x%d,", j); } printf("\n");
  28.  
  29.   for(i=0;i<tab->m; i++) {
  30.     if (i==0) printf("max:"); else
  31.     printf("b%d: ", i);
  32.     for(j=0;j<tab->n; j++) {
  33.       if (equal((int)tab->mat[i][j], tab->mat[i][j]))
  34.         printf(" %6d", (int)tab->mat[i][j]);
  35.       else
  36.         printf(" %6.2lf", tab->mat[i][j]);
  37.       }
  38.     printf("\n");
  39.   }
  40.   nl(70);
  41. }
  42.  
  43. /* Example input file for read_tableau:
  44.      4 5
  45.       0   -0.5  -3 -1  -4 
  46.      40    1     1  1   1 
  47.      10   -2    -1  1   1 
  48.      10    0     1  0  -1  
  49. */
  50. void read_tableau(Tableau *tab, const char * filename) {
  51.   int err, i, j;
  52.   FILE * fp;
  53.  
  54.   fp  = fopen(filename, "r" );
  55.   if( !fp ) {
  56.     printf("Cannot read %s\n", filename); exit(1);
  57.   }
  58.   memset(tab, 0, sizeof(*tab));
  59.   err = fscanf(fp, "%d %d", &tab->m, &tab->n);
  60.   if (err == 0 || err == EOF) {
  61.     printf("Cannot read m or n\n"); exit(1);
  62.   }
  63.   for(i=0;i<tab->m; i++) {
  64.     for(j=0;j<tab->n; j++) {
  65.       err = fscanf(fp, "%lf", &tab->mat[i][j]);
  66.       if (err == 0 || err == EOF) {
  67.         printf("Cannot read A[%d][%d]\n", i, j); exit(1);
  68.       }
  69.     }
  70.   }
  71.   printf("Read tableau [%d rows x %d columns] from file '%s'.\n",
  72.     tab->m, tab->n, filename);
  73.   fclose(fp);
  74. }
  75.  
  76. void pivot_on(Tableau *tab, int row, int col) {
  77.   int i, j;
  78.   double pivot;
  79.  
  80.   pivot = tab->mat[row][col];
  81.   assert(pivot>0);
  82.   for(j=0;j<tab->n;j++)
  83.     tab->mat[row][j] /= pivot;
  84.   assert( equal(tab->mat[row][col], 1. ));
  85.  
  86.   for(i=0; i<tab->m; i++) { // foreach remaining row i do
  87.     double multiplier = tab->mat[i][col];
  88.     if(i==row) continue;
  89.     for(j=0; j<tab->n; j++) { // r[i] = r[i] - z * r[row];
  90.       tab->mat[i][j] -= multiplier * tab->mat[row][j];
  91.     }
  92.   }
  93. }
  94.  
  95. // Find pivot_col = most negative column in mat[0][1..n]
  96. int find_pivot_column(Tableau *tab) {
  97.   int j, pivot_col = 1;
  98.   double lowest = tab->mat[0][pivot_col];
  99.   for(j=1; j<tab->n; j++) {
  100.     if (tab->mat[0][j] < lowest) {
  101.       lowest = tab->mat[0][j];
  102.       pivot_col = j;
  103.     }
  104.   }
  105.   printf("Most negative column in row[0] is col %d = %g.\n", pivot_col, lowest);
  106.   if( lowest >= 0 ) {
  107.     return -1; // All positive columns in row[0], this is optimal.
  108.   }
  109.   return pivot_col;
  110. }
  111.  
  112. // Find the pivot_row, with smallest positive ratio = col[0] / col[pivot]
  113. int find_pivot_row(Tableau *tab, int pivot_col) {
  114.   int i, pivot_row = 0;
  115.   double min_ratio = -1;
  116.   printf("Ratios A[row_i,0]/A[row_i,%d] = [",pivot_col);
  117.   for(i=1;i<tab->m;i++){
  118.     double ratio = tab->mat[i][0] / tab->mat[i][pivot_col];
  119.     printf("%3.2lf, ", ratio);
  120.     if ( (ratio > 0  && ratio < min_ratio ) || min_ratio < 0 ) {
  121.       min_ratio = ratio;
  122.       pivot_row = i;
  123.     }
  124.   }
  125.   printf("].\n");
  126.   if (min_ratio == -1)
  127.     return -1; // Unbounded.
  128.   printf("Found pivot A[%d,%d], min positive ratio=%g in row=%d.\n",
  129.       pivot_row, pivot_col, min_ratio, pivot_row);
  130.   return pivot_row;
  131. }
  132.  
  133. void add_slack_variables(Tableau *tab) {
  134.   int i, j;
  135.   for(i=1; i<tab->m; i++) {
  136.     for(j=1; j<tab->m; j++)
  137.       tab->mat[i][j + tab->n -1] = (i==j);
  138.   }
  139.   tab->n += tab->m -1;
  140. }
  141.  
  142. void check_b_positive(Tableau *tab) {
  143.   int i;
  144.   for(i=1; i<tab->m; i++)
  145.     assert(tab->mat[i][0] >= 0);
  146. }
  147.  
  148. // Given a column of identity matrix, find the row containing 1.
  149. // return -1, if the column as not from an identity matrix.
  150. int find_basis_variable(Tableau *tab, int col) {
  151.   int i, xi=-1;
  152.   for(i=1; i < tab->m; i++) {
  153.     if (equal( tab->mat[i][col],1) ) {
  154.       if (xi == -1)
  155.         xi=i;   // found first '1', save this row number.
  156.       else
  157.         return -1; // found second '1', not an identity matrix.
  158.  
  159.     } else if (!equal( tab->mat[i][col],0) ) {
  160.       return -1; // not an identity matrix column.
  161.     }
  162.   }
  163.   return xi;
  164. }
  165.  
  166. void print_optimal_vector(Tableau *tab, char *message) {
  167.   int j, xi;
  168.   printf("%s at ", message);
  169.   for(j=1;j<tab->n;j++) { // for each column.
  170.     xi = find_basis_variable(tab, j);
  171.     if (xi != -1)
  172.       printf("x%d=%3.2lf, ", j, tab->mat[xi][0] );
  173.     else
  174.       printf("x%d=0, ", j);
  175.   }
  176.   printf("\n");
  177. }
  178.  
  179. void simplex(Tableau *tab) {
  180.   int loop=0;
  181.   add_slack_variables(tab);
  182.   check_b_positive(tab);
  183.   print_tableau(tab,"Padded with slack variables");
  184.   while( ++loop ) {
  185.     int pivot_col, pivot_row;
  186.  
  187.     pivot_col = find_pivot_column(tab);
  188.     if( pivot_col < 0 ) {
  189.       printf("Found optimal value=A[0,0]=%3.2lf (no negatives in row 0).\n",
  190.         tab->mat[0][0]);
  191.       print_optimal_vector(tab, "Optimal vector");
  192.       break;
  193.     }
  194.     printf("Entering variable x%d to be made basic, so pivot_col=%d.\n",
  195.       pivot_col, pivot_col);
  196.  
  197.     pivot_row = find_pivot_row(tab, pivot_col);
  198.     if (pivot_row < 0) {
  199.       printf("unbounded (no pivot_row).\n");
  200.       break;
  201.     }
  202.     printf("Leaving variable x%d, so pivot_row=%d\n", pivot_row, pivot_row);
  203.  
  204.     pivot_on(tab, pivot_row, pivot_col);
  205.     print_tableau(tab,"After pivoting");
  206.     print_optimal_vector(tab, "Basic feasible solution");
  207.  
  208.     if(loop > 20) {
  209.       printf("Too many iterations > %d.\n", loop);
  210.       break;
  211.     }
  212.   }
  213. }
  214.  
  215. Tableau tab  = { 4, 5, {                     // Size of tableau [4 rows x 5 columns ]
  216.     {  0.0 , -0.5 , -3.0 ,-1.0 , -4.0,   },  // Max: z = 0.5*x + 3*y + z + 4*w,
  217.     { 40.0 ,  1.0 ,  1.0 , 1.0 ,  1.0,   },  //    x + y + z + w <= 40 .. b1
  218.     { 10.0 , -2.0 , -1.0 , 1.0 ,  1.0,   },  //  -2x - y + z + w <= 10 .. b2
  219.     { 10.0 ,  0.0 ,  1.0 , 0.0 , -1.0,   },  //        y     - w <= 10 .. b3
  220.   }
  221. };
  222.  
  223. int main(int argc, char *argv[]){
  224.   if (argc > 1) { // usage: cmd datafile
  225.     read_tableau(&tab, argv[1]);
  226.   }
  227.   print_tableau(&tab,"Initial");
  228.   simplex(&tab);
  229.   return 0;
  230. }
Success #stdin #stdout 0.01s 5292KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1. Tableau Initial:
----------------------------------------------------------------------
col:   b[i]    x1,    x2,    x3,    x4,
max:      0  -0.50     -3     -1     -4
b1:      40      1      1      1      1
b2:      10     -2     -1      1      1
b3:      10      0      1      0     -1
----------------------------------------------------------------------

2. Tableau Padded with slack variables:
----------------------------------------------------------------------
col:   b[i]    x1,    x2,    x3,    x4,    x5,    x6,    x7,
max:      0  -0.50     -3     -1     -4      0      0      0
b1:      40      1      1      1      1      1      0      0
b2:      10     -2     -1      1      1      0      1      0
b3:      10      0      1      0     -1      0      0      1
----------------------------------------------------------------------
Most negative column in row[0] is col 4 = -4.
Entering variable x4 to be made basic, so pivot_col=4.
Ratios A[row_i,0]/A[row_i,4] = [40.00, 10.00, -10.00, ].
Found pivot A[2,4], min positive ratio=10 in row=2.
Leaving variable x2, so pivot_row=2

3. Tableau After pivoting:
----------------------------------------------------------------------
col:   b[i]    x1,    x2,    x3,    x4,    x5,    x6,    x7,
max:     40  -8.50     -7      3      0      0      4      0
b1:      30      3      2      0      0      1     -1      0
b2:      10     -2     -1      1      1      0      1      0
b3:      20     -2      0      1      0      0      1      1
----------------------------------------------------------------------
Basic feasible solution at x1=0, x2=0, x3=0, x4=10.00, x5=30.00, x6=0, x7=20.00, 
Most negative column in row[0] is col 1 = -8.5.
Entering variable x1 to be made basic, so pivot_col=1.
Ratios A[row_i,0]/A[row_i,1] = [10.00, -5.00, -10.00, ].
Found pivot A[1,1], min positive ratio=10 in row=1.
Leaving variable x1, so pivot_row=1

4. Tableau After pivoting:
----------------------------------------------------------------------
col:   b[i]    x1,    x2,    x3,    x4,    x5,    x6,    x7,
max:    125      0  -1.33      3      0   2.83   1.17      0
b1:      10      1   0.67      0      0   0.33  -0.33      0
b2:      30      0   0.33      1      1   0.67   0.33      0
b3:      40      0   1.33      1      0   0.67   0.33      1
----------------------------------------------------------------------
Basic feasible solution at x1=10.00, x2=0, x3=0, x4=30.00, x5=0, x6=0, x7=40.00, 
Most negative column in row[0] is col 2 = -1.33333.
Entering variable x2 to be made basic, so pivot_col=2.
Ratios A[row_i,0]/A[row_i,2] = [15.00, 90.00, 30.00, ].
Found pivot A[1,2], min positive ratio=15 in row=1.
Leaving variable x1, so pivot_row=1

5. Tableau After pivoting:
----------------------------------------------------------------------
col:   b[i]    x1,    x2,    x3,    x4,    x5,    x6,    x7,
max:    145      2      0      3      0   3.50   0.50      0
b1:      15   1.50      1      0      0   0.50  -0.50      0
b2:      25  -0.50      0      1      1   0.50   0.50      0
b3:      20     -2      0      1      0      0      1      1
----------------------------------------------------------------------
Basic feasible solution at x1=0, x2=15.00, x3=0, x4=25.00, x5=0, x6=0, x7=20.00, 
Most negative column in row[0] is col 2 = 0.
Found optimal value=A[0,0]=145.00 (no negatives in row 0).
Optimal vector at x1=0, x2=15.00, x3=0, x4=25.00, x5=0, x6=0, x7=20.00,
https://ideone.com/tbmLHE
language:
C++ (gcc 8.3)
created:
12 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!