1 | // -*- c++ -*- |
---|
2 | #include <iostream> |
---|
3 | #include <fstream> |
---|
4 | |
---|
5 | #include <lemon/time_measure.h> |
---|
6 | #include <lp_solver_glpk.h> |
---|
7 | |
---|
8 | using std::cout; |
---|
9 | using std::endl; |
---|
10 | using namespace lemon; |
---|
11 | |
---|
12 | /* |
---|
13 | On an 1537Mhz PC, the run times with |
---|
14 | glpk are the following. |
---|
15 | for n=3,4, some secondes |
---|
16 | for n=5, 25 hours |
---|
17 | */ |
---|
18 | |
---|
19 | int main(int, char **) { |
---|
20 | const int n=4; |
---|
21 | const double row_sum=(1.0+n*n)*n/2; |
---|
22 | Timer ts; |
---|
23 | ts.reset(); |
---|
24 | typedef LpGlpk LPSolver; |
---|
25 | typedef LPSolver::Col Col; |
---|
26 | LPSolver lp; |
---|
27 | typedef std::map<std::pair<int, int>, Col> Coords; |
---|
28 | Coords x; |
---|
29 | // we create a new variable for each entry |
---|
30 | // of the magic square |
---|
31 | for (int i=1; i<=n; ++i) { |
---|
32 | for (int j=1; j<=n; ++j) { |
---|
33 | Col col=lp.addCol(); |
---|
34 | x[std::make_pair(i,j)]=col; |
---|
35 | lp.setColLowerBound(col, 1.0); |
---|
36 | lp.setColUpperBound(col, double(n*n)); |
---|
37 | } |
---|
38 | } |
---|
39 | LPSolver::Expression expr3, expr4; |
---|
40 | for (int i=1; i<=n; ++i) { |
---|
41 | LPSolver::Expression expr1, expr2; |
---|
42 | for (int j=1; j<=n; ++j) { |
---|
43 | expr1+=x[std::make_pair(i, j)]; |
---|
44 | expr2+=x[std::make_pair(j, i)]; |
---|
45 | } |
---|
46 | // sum of rows and columns |
---|
47 | lp.addRow(expr1==row_sum); |
---|
48 | lp.addRow(expr2==row_sum); |
---|
49 | expr3+=x[std::make_pair(i, i)]; |
---|
50 | expr4+=x[std::make_pair(i, (n+1)-i)]; |
---|
51 | } |
---|
52 | // sum of the diagonal entries |
---|
53 | lp.addRow(expr3==row_sum); |
---|
54 | lp.addRow(expr4==row_sum); |
---|
55 | lp.solveSimplex(); |
---|
56 | cout << "elapsed time: " << ts << endl; |
---|
57 | for (int i=1; i<=n; ++i) { |
---|
58 | for (int j=1; j<=n; ++j) { |
---|
59 | cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)]) |
---|
60 | << endl; |
---|
61 | } |
---|
62 | } |
---|
63 | // we make new binary variables for each pair of |
---|
64 | // entries of the square to achieve that |
---|
65 | // the values of different entries are different |
---|
66 | lp.setMIP(); |
---|
67 | for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) { |
---|
68 | Coords::const_iterator jt=it; ++jt; |
---|
69 | for(; jt!=x.end(); ++jt) { |
---|
70 | Col col1=(*it).second; |
---|
71 | Col col2=(*jt).second; |
---|
72 | Col col=lp.addCol(); |
---|
73 | lp.setColLowerBound(col, 0.0); |
---|
74 | lp.setColUpperBound(col, 1.0); |
---|
75 | lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0); |
---|
76 | lp.setColInt(col); |
---|
77 | } |
---|
78 | } |
---|
79 | cout << "elapsed time: " << ts << endl; |
---|
80 | lp.solveSimplex(); |
---|
81 | // let's solve the integer problem |
---|
82 | lp.solveBandB(); |
---|
83 | cout << "elapsed time: " << ts << endl; |
---|
84 | for (int i=1; i<=n; ++i) { |
---|
85 | for (int j=1; j<=n; ++j) { |
---|
86 | cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)]) |
---|
87 | << endl; |
---|
88 | } |
---|
89 | } |
---|
90 | } |
---|