if you have a nuclear power plant and wanna compute small magic squares, then let's do it
5 #include <lemon/time_measure.h>
6 #include <lp_solver_base.h>
10 using namespace lemon;
13 for n=3,4 , the program is very fast
14 for n=5, with glpk, the run takes hours
17 int main(int, char **) {
19 const double row_sum=(1.0+n*n)*n/2;
22 typedef LPGLPK LPSolver;
23 typedef LPSolver::Col Col;
25 typedef std::map<std::pair<int, int>, Col> Coords;
27 // we create a new variable for each entry
28 // of the magic square
29 for (int i=1; i<=n; ++i) {
30 for (int j=1; j<=n; ++j) {
32 x[std::make_pair(i,j)]=col;
33 lp.setColLowerBound(col, 1.0);
34 lp.setColUpperBound(col, double(n*n));
37 LPSolver::Expression expr3, expr4;
38 for (int i=1; i<=n; ++i) {
39 LPSolver::Expression expr1, expr2;
40 for (int j=1; j<=n; ++j) {
41 expr1+=x[std::make_pair(i, j)];
42 expr2+=x[std::make_pair(j, i)];
44 // sum of rows and columns
45 lp.addRow(expr1==row_sum);
46 lp.addRow(expr2==row_sum);
47 expr3+=x[std::make_pair(i, i)];
48 expr4+=x[std::make_pair(i, (n+1)-i)];
50 // sum of the diagonal entries
51 lp.addRow(expr3==row_sum);
52 lp.addRow(expr4==row_sum);
54 cout << "elapsed time: " << ts << endl;
55 for (int i=1; i<=n; ++i) {
56 for (int j=1; j<=n; ++j) {
57 cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)])
61 // we make new binary variables for each pair of
62 // entries of the square to achieve that
63 // the values of different entries are different
65 for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) {
66 Coords::const_iterator jt=it; ++jt;
67 for(; jt!=x.end(); ++jt) {
68 Col col1=(*it).second;
69 Col col2=(*jt).second;
71 lp.setColLowerBound(col, 0.0);
72 lp.setColUpperBound(col, 1.0);
73 lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0);
77 cout << "elapsed time: " << ts << endl;
79 // let's solve the integer problem
81 cout << "elapsed time: " << ts << endl;
82 for (int i=1; i<=n; ++i) {
83 for (int j=1; j<=n; ++j) {
84 cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)])