§1 LP INTRODUCTION 1
October 3, 2008 at 23:08
1. Introduction. I’m writing this program in order to gain personal experience implementing the
simplex algorithm—even though I know that commercial codes do the job much, much better. My aim
is to learn, to make the logic “crystal clear” if not lightning fast, and perhaps also to watch and interact
with this magical process.
The computational task to be solved has been traditionally called Linear Programming, and it takes many
forms. The particular case considered here is to maximize b
1
v
1
+ · · · + b
n
v
n
subject to the constraints v
j
≥ 0
for 1 ≤ j ≤ n and a
i1
v
1
+ · · · + a
in
v
n
≤ c
i
for 1 ≤ i ≤ m, where a
ij
is a given m × n matrix of integers and
the other parameters b
j
, c
i
are integers with c
i
≥ 0. For example, what is the maximum of v
1
+ 3v
2
− v
3
, if
we require that
3v
1
− v
2
+ 5v
3
≤ 8, −v
1
+ 2v
2
− v
3
≤ 1, 2v
1
+ 4v
2
+ v
3
≤ 5,
and v
1
, v
2
, v
3
≥ 0? [I took this example from Muroga’s book on threshold logic (1971).] The fact that the
constants c
i
on the right-hand side are nonnegative means that the trivial values v
1
= · · · = v
n
= 0 always
satisfy the constraints; hence the maximum is always nonnegative.
The algorithm below also solves a “dual” problem as a special bonus. Namely, it tells us how to minimize
the quantity c
1
u
1
+ · · · + c
m
u
m
subject to u
i
≥ 0 and a
1j
u
1
+ · · · + a
mj
u
m
≥ b
j
for 1 ≤ j ≤ n.
There may no values (u
1
, . . . , u
m
) that meet those dual constraints. In such a case, the algorithm will
effectively prove the impossibility, and it will also demonstrate that the maximum in the original problem
is +∞. (For example, suppose m = n = 1, b
1
= 1, c
1
= 0, and a
11
= −1. The problem of maximizing
v
1
subject to −v
1
≤ 0 and v
1
≥ 0 obviously has +∞ as its answer; and +∞ is also the “minimum” of the
quantity 0 taken over all u
1
such that u
1
≥ 0 and −u
1
≥ 1, because no such numbers u
1
exist.)
The first line of the standard input file should contain the integers b
1
b
2
. . . b
n
in decimal notation,
separated by spaces. That first line should be followed by m further lines that each contain n + 1 integer
values c
i
a
i1
a
i2
. . . a
in
, for 1 ≤ i ≤ m.
To enhance this learning experience, I’m solving the problem both with floating-point arithmetic and with
an all-integer method that produces rational numbers as output. Hopefully the two answers will agree. But
the all-integer method might overflow, and the floating-point method may suffer from rounding errors.
#define maxm 10 /∗ of course these limits can be raised if desired ∗/
#define maxn 100 /∗ (up to a point) ∗/
#define buf
size BUFSIZ
#include <stdio.h>
#include <ctype.h>
h Include tricky code for zapping 2 i
typedef long intword; /∗ will be long long on my other computer ∗/
char buf [buf
size ];
intword a[maxm + 1][maxm + maxn + 1]; /∗ integer work area ∗/
intword denom [maxm + maxn + 1]; /∗ scale factors ∗/
double aa [maxm + 1][maxm + maxn + 1]; /∗ floating-point work area ∗/
double trial [maxm + maxn + 1]; /∗ pivot testing area ∗/
int verbose = 1; /∗ can be set positive for extra output ∗/
int count ; /∗ the number of steps taken so far ∗/
int p[maxm + 1], q[maxm + maxn + 1]; /∗ current basis and inverse ∗/
main ( )
{
register intword h, i, j, k, l, m, n, s;
register double z;
h Check the zap trick 3 i;
h Read the input matrix 4 i;
h Solve the problem 13 i;
}
