marci@1153: // -*- c++ -*- marci@1153: #include <iostream> marci@1153: #include <fstream> marci@1153: marci@1153: #include <lemon/time_measure.h> athos@1241: #include <lp_solver_glpk.h> marci@1153: marci@1153: using std::cout; marci@1153: using std::endl; marci@1153: using namespace lemon; marci@1153: marci@1153: /* marci@1176: On an 1537Mhz PC, the run times with marci@1176: glpk are the following. marci@1176: for n=3,4, some secondes marci@1176: for n=5, 25 hours marci@1153: */ marci@1153: marci@1153: int main(int, char **) { athos@1243: const int n=3; marci@1153: const double row_sum=(1.0+n*n)*n/2; marci@1153: Timer ts; marci@1153: ts.reset(); athos@1241: typedef LpGlpk LPSolver; marci@1153: typedef LPSolver::Col Col; marci@1153: LPSolver lp; marci@1153: typedef std::map<std::pair<int, int>, Col> Coords; marci@1153: Coords x; marci@1153: // we create a new variable for each entry marci@1153: // of the magic square marci@1153: for (int i=1; i<=n; ++i) { marci@1153: for (int j=1; j<=n; ++j) { marci@1153: Col col=lp.addCol(); marci@1153: x[std::make_pair(i,j)]=col; marci@1153: lp.setColLowerBound(col, 1.0); marci@1153: lp.setColUpperBound(col, double(n*n)); marci@1153: } marci@1153: } marci@1153: LPSolver::Expression expr3, expr4; marci@1153: for (int i=1; i<=n; ++i) { marci@1153: LPSolver::Expression expr1, expr2; marci@1153: for (int j=1; j<=n; ++j) { marci@1153: expr1+=x[std::make_pair(i, j)]; marci@1153: expr2+=x[std::make_pair(j, i)]; marci@1153: } athos@1243: marci@1153: // sum of rows and columns marci@1153: lp.addRow(expr1==row_sum); marci@1153: lp.addRow(expr2==row_sum); athos@1243: cout <<"Expr1: "<<expr1<<endl; athos@1243: cout <<"Expr2: "<<expr2<<endl; athos@1243: marci@1153: expr3+=x[std::make_pair(i, i)]; marci@1153: expr4+=x[std::make_pair(i, (n+1)-i)]; marci@1153: } athos@1243: cout <<"Expr3: "<<expr3<<endl; athos@1243: cout <<"Expr4: "<<expr4<<endl; marci@1153: // sum of the diagonal entries marci@1153: lp.addRow(expr3==row_sum); marci@1153: lp.addRow(expr4==row_sum); marci@1153: lp.solveSimplex(); marci@1153: cout << "elapsed time: " << ts << endl; marci@1153: for (int i=1; i<=n; ++i) { marci@1153: for (int j=1; j<=n; ++j) { marci@1153: cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)]) marci@1153: << endl; marci@1153: } marci@1153: } athos@1243: athos@1243: athos@1243: athos@1243: // // we make new binary variables for each pair of athos@1243: // // entries of the square to achieve that athos@1243: // // the values of different entries are different athos@1243: // lp.setMIP(); athos@1243: // for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) { athos@1243: // Coords::const_iterator jt=it; ++jt; athos@1243: // for(; jt!=x.end(); ++jt) { athos@1243: // Col col1=(*it).second; athos@1243: // Col col2=(*jt).second; athos@1243: // Col col=lp.addCol(); athos@1243: // lp.setColLowerBound(col, 0.0); athos@1243: // lp.setColUpperBound(col, 1.0); athos@1243: // lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0); athos@1243: // lp.setColInt(col); athos@1243: // } athos@1243: // } athos@1243: // cout << "elapsed time: " << ts << endl; athos@1243: // lp.solveSimplex(); athos@1243: // // let's solve the integer problem athos@1243: // lp.solveBandB(); athos@1243: // cout << "elapsed time: " << ts << endl; athos@1243: // for (int i=1; i<=n; ++i) { athos@1243: // for (int j=1; j<=n; ++j) { athos@1243: // cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)]) athos@1243: // << endl; athos@1243: // } athos@1243: // } marci@1153: }