lemon/lp_glpk.cc
author Balazs Dezso <deba@inf.elte.hu>
Tue, 02 Dec 2008 21:40:33 +0100
changeset 458 7afc121e0689
child 459 ed54c0d13df0
permissions -rw-r--r--
Port LP and MIP solvers from SVN -r3509 (#44)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2008
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 ///\file
    20 ///\brief Implementation of the LEMON-GLPK lp solver interface.
    21 
    22 #include <lemon/lp_glpk.h>
    23 //#include <iostream>
    24 
    25 extern "C" {
    26 #include <glpk.h>
    27 }
    28 
    29 #if GLP_MAJOR_VERSION > 4 || (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION > 15)
    30 #define LEMON_glp(func) (glp_##func)
    31 #define LEMON_lpx(func) (lpx_##func)
    32 
    33 #define LEMON_GLP(def) (GLP_##def)
    34 #define LEMON_LPX(def) (LPX_##def)
    35 
    36 #else
    37 
    38 #define LEMON_glp(func) (lpx_##func)
    39 #define LEMON_lpx(func) (lpx_##func)
    40 
    41 #define LEMON_GLP(def) (LPX_##def)
    42 #define LEMON_LPX(def) (LPX_##def)
    43 
    44 #endif
    45 
    46 namespace lemon {
    47 
    48   LpGlpk::LpGlpk() : Parent() {
    49     solved = false;
    50     rows = _lp_bits::LpId(1);
    51     cols = _lp_bits::LpId(1);
    52     lp = LEMON_glp(create_prob)();
    53     LEMON_glp(create_index)(lp);
    54     messageLevel(0);
    55   }
    56 
    57   LpGlpk::LpGlpk(const LpGlpk &glp) : Parent() {
    58     solved = false;
    59     rows = _lp_bits::LpId(1);
    60     cols = _lp_bits::LpId(1);
    61     lp = LEMON_glp(create_prob)();
    62     LEMON_glp(create_index)(lp);
    63     messageLevel(0);
    64     //Coefficient matrix, row bounds
    65     LEMON_glp(add_rows)(lp, LEMON_glp(get_num_rows)(glp.lp));
    66     LEMON_glp(add_cols)(lp, LEMON_glp(get_num_cols)(glp.lp));
    67     int len;
    68     std::vector<int> ind(1+LEMON_glp(get_num_cols)(glp.lp));
    69     std::vector<Value> val(1+LEMON_glp(get_num_cols)(glp.lp));
    70     for (int i=1;i<=LEMON_glp(get_num_rows)(glp.lp);++i)
    71       {
    72         len=LEMON_glp(get_mat_row)(glp.lp,i,&*ind.begin(),&*val.begin());
    73         LEMON_glp(set_mat_row)(lp, i,len,&*ind.begin(),&*val.begin());
    74         LEMON_glp(set_row_bnds)(lp,i,
    75                                 LEMON_glp(get_row_type)(glp.lp,i),
    76                                 LEMON_glp(get_row_lb)(glp.lp,i),
    77                                 LEMON_glp(get_row_ub)(glp.lp,i));
    78       }
    79 
    80     //Objective function, coloumn bounds
    81     LEMON_glp(set_obj_dir)(lp, LEMON_glp(get_obj_dir)(glp.lp));
    82     //Objectif function's constant term treated separately
    83     LEMON_glp(set_obj_coef)(lp,0,LEMON_glp(get_obj_coef)(glp.lp,0));
    84     for (int i=1;i<=LEMON_glp(get_num_cols)(glp.lp);++i)
    85       {
    86         LEMON_glp(set_obj_coef)(lp,i,
    87                                 LEMON_glp(get_obj_coef)(glp.lp,i));
    88         LEMON_glp(set_col_bnds)(lp,i,
    89                                 LEMON_glp(get_col_type)(glp.lp,i),
    90                                 LEMON_glp(get_col_lb)(glp.lp,i),
    91                                 LEMON_glp(get_col_ub)(glp.lp,i));
    92       }
    93     rows = glp.rows;
    94     cols = glp.cols;
    95   }
    96 
    97   LpGlpk::~LpGlpk() {
    98     LEMON_glp(delete_prob)(lp);
    99   }
   100 
   101   int LpGlpk::_addCol() {
   102     int i=LEMON_glp(add_cols)(lp, 1);
   103     LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), 0.0, 0.0);
   104     solved = false;
   105     return i;
   106   }
   107 
   108   ///\e
   109 
   110 
   111   LpSolverBase* LpGlpk::_newLp()
   112   {
   113     LpGlpk* newlp = new LpGlpk;
   114     return newlp;
   115   }
   116 
   117   ///\e
   118 
   119   LpSolverBase* LpGlpk::_copyLp()
   120   {
   121     LpGlpk *newlp = new LpGlpk(*this);
   122     return newlp;
   123   }
   124 
   125   int LpGlpk::_addRow() {
   126     int i=LEMON_glp(add_rows)(lp, 1);
   127     solved = false;
   128     return i;
   129   }
   130 
   131 
   132   void LpGlpk::_eraseCol(int i) {
   133     int ca[2];
   134     ca[1]=i;
   135     LEMON_glp(del_cols)(lp, 1, ca);
   136     solved = false;
   137   }
   138 
   139   void LpGlpk::_eraseRow(int i) {
   140     int ra[2];
   141     ra[1]=i;
   142     LEMON_glp(del_rows)(lp, 1, ra);
   143     solved = false;
   144   }
   145 
   146   void LpGlpk::_getColName(int c, std::string & name) const
   147   {
   148 
   149     const char *n = LEMON_glp(get_col_name)(lp,c);
   150     name = n?n:"";
   151   }
   152 
   153 
   154   void LpGlpk::_setColName(int c, const std::string & name)
   155   {
   156     LEMON_glp(set_col_name)(lp,c,const_cast<char*>(name.c_str()));
   157 
   158   }
   159 
   160   int LpGlpk::_colByName(const std::string& name) const
   161   {
   162     int k = LEMON_glp(find_col)(lp, const_cast<char*>(name.c_str()));
   163     return k > 0 ? k : -1;
   164   }
   165 
   166 
   167   void LpGlpk::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
   168   {
   169     std::vector<int> indices;
   170     std::vector<Value> values;
   171 
   172     indices.push_back(0);
   173     values.push_back(0);
   174 
   175     for(ConstRowIterator it=b; it!=e; ++it) {
   176       indices.push_back(it->first);
   177       values.push_back(it->second);
   178     }
   179 
   180     LEMON_glp(set_mat_row)(lp, i, values.size() - 1,
   181                                 &indices[0], &values[0]);
   182 
   183     solved = false;
   184   }
   185 
   186   void LpGlpk::_getRowCoeffs(int ix, RowIterator b) const
   187   {
   188     int length = LEMON_glp(get_mat_row)(lp, ix, 0, 0);
   189 
   190     std::vector<int> indices(length + 1);
   191     std::vector<Value> values(length + 1);
   192 
   193     LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
   194 
   195     for (int i = 1; i <= length; ++i) {
   196       *b = std::make_pair(indices[i], values[i]);
   197       ++b;
   198     }
   199   }
   200 
   201   void LpGlpk::_setColCoeffs(int ix, ConstColIterator b, ConstColIterator e) {
   202 
   203     std::vector<int> indices;
   204     std::vector<Value> values;
   205 
   206     indices.push_back(0);
   207     values.push_back(0);
   208 
   209     for(ConstColIterator it=b; it!=e; ++it) {
   210       indices.push_back(it->first);
   211       values.push_back(it->second);
   212     }
   213 
   214     LEMON_glp(set_mat_col)(lp, ix, values.size() - 1,
   215                                 &indices[0], &values[0]);
   216 
   217     solved = false;
   218   }
   219 
   220   void LpGlpk::_getColCoeffs(int ix, ColIterator b) const
   221   {
   222     int length = LEMON_glp(get_mat_col)(lp, ix, 0, 0);
   223 
   224     std::vector<int> indices(length + 1);
   225     std::vector<Value> values(length + 1);
   226 
   227     LEMON_glp(get_mat_col)(lp, ix, &indices[0], &values[0]);
   228 
   229     for (int i = 1; i <= length; ++i) {
   230       *b = std::make_pair(indices[i], values[i]);
   231       ++b;
   232     }
   233   }
   234 
   235   void LpGlpk::_setCoeff(int ix, int jx, Value value)
   236   {
   237 
   238     if (LEMON_glp(get_num_cols)(lp) < LEMON_glp(get_num_rows)(lp)) {
   239 
   240       int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0);
   241 
   242       std::vector<int> indices(length + 2);
   243       std::vector<Value> values(length + 2);
   244 
   245       LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
   246 
   247       //The following code does not suppose that the elements of the
   248       //array indices are sorted
   249       bool found=false;
   250       for (int i = 1; i <= length; ++i) {
   251         if (indices[i]==jx){
   252           found=true;
   253           values[i]=value;
   254           break;
   255         }
   256       }
   257       if (!found){
   258         ++length;
   259         indices[length]=jx;
   260         values[length]=value;
   261       }
   262 
   263       LEMON_glp(set_mat_row)(lp, ix, length, &indices[0], &values[0]);
   264 
   265     } else {
   266 
   267       int length=LEMON_glp(get_mat_col)(lp, jx, 0, 0);
   268 
   269       std::vector<int> indices(length + 2);
   270       std::vector<Value> values(length + 2);
   271 
   272       LEMON_glp(get_mat_col)(lp, jx, &indices[0], &values[0]);
   273 
   274       //The following code does not suppose that the elements of the
   275       //array indices are sorted
   276       bool found=false;
   277       for (int i = 1; i <= length; ++i) {
   278         if (indices[i]==ix){
   279           found=true;
   280           values[i]=value;
   281           break;
   282         }
   283       }
   284       if (!found){
   285         ++length;
   286         indices[length]=ix;
   287         values[length]=value;
   288       }
   289 
   290       LEMON_glp(set_mat_col)(lp, jx, length, &indices[0], &values[0]);
   291     }
   292 
   293     solved = false;
   294   }
   295 
   296   LpGlpk::Value LpGlpk::_getCoeff(int ix, int jx) const
   297   {
   298 
   299     int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0);
   300 
   301     std::vector<int> indices(length + 1);
   302     std::vector<Value> values(length + 1);
   303 
   304     LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
   305 
   306     //The following code does not suppose that the elements of the
   307     //array indices are sorted
   308     for (int i = 1; i <= length; ++i) {
   309       if (indices[i]==jx){
   310         return values[i];
   311       }
   312     }
   313     return 0;
   314 
   315   }
   316 
   317 
   318   void LpGlpk::_setColLowerBound(int i, Value lo)
   319   {
   320     if (lo==INF) {
   321       //FIXME error
   322     }
   323     int b=LEMON_glp(get_col_type)(lp, i);
   324     double up=LEMON_glp(get_col_ub)(lp, i);
   325     if (lo==-INF) {
   326       switch (b) {
   327       case LEMON_GLP(FR):
   328       case LEMON_GLP(LO):
   329         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up);
   330         break;
   331       case LEMON_GLP(UP):
   332         break;
   333       case LEMON_GLP(DB):
   334       case LEMON_GLP(FX):
   335         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
   336         break;
   337       default: ;
   338         //FIXME error
   339       }
   340     } else {
   341       switch (b) {
   342       case LEMON_GLP(FR):
   343       case LEMON_GLP(LO):
   344         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up);
   345         break;
   346       case LEMON_GLP(UP):
   347       case LEMON_GLP(DB):
   348       case LEMON_GLP(FX):
   349         if (lo==up)
   350           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up);
   351         else
   352           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up);
   353         break;
   354       default: ;
   355         //FIXME error
   356       }
   357     }
   358 
   359     solved = false;
   360   }
   361 
   362   LpGlpk::Value LpGlpk::_getColLowerBound(int i) const
   363   {
   364     int b=LEMON_glp(get_col_type)(lp, i);
   365       switch (b) {
   366       case LEMON_GLP(LO):
   367       case LEMON_GLP(DB):
   368       case LEMON_GLP(FX):
   369         return LEMON_glp(get_col_lb)(lp, i);
   370       default: ;
   371         return -INF;
   372       }
   373   }
   374 
   375   void LpGlpk::_setColUpperBound(int i, Value up)
   376   {
   377     if (up==-INF) {
   378       //FIXME error
   379     }
   380     int b=LEMON_glp(get_col_type)(lp, i);
   381     double lo=LEMON_glp(get_col_lb)(lp, i);
   382     if (up==INF) {
   383       switch (b) {
   384       case LEMON_GLP(FR):
   385       case LEMON_GLP(LO):
   386         break;
   387       case LEMON_GLP(UP):
   388         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up);
   389         break;
   390       case LEMON_GLP(DB):
   391       case LEMON_GLP(FX):
   392         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up);
   393         break;
   394       default: ;
   395         //FIXME error
   396       }
   397     } else {
   398       switch (b) {
   399       case LEMON_GLP(FR):
   400         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
   401         break;
   402       case LEMON_GLP(UP):
   403         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
   404         break;
   405       case LEMON_GLP(LO):
   406       case LEMON_GLP(DB):
   407       case LEMON_GLP(FX):
   408         if (lo==up)
   409           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up);
   410         else
   411           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up);
   412         break;
   413       default: ;
   414         //FIXME error
   415       }
   416     }
   417 
   418     solved = false;
   419   }
   420 
   421   LpGlpk::Value LpGlpk::_getColUpperBound(int i) const
   422   {
   423     int b=LEMON_glp(get_col_type)(lp, i);
   424       switch (b) {
   425       case LEMON_GLP(UP):
   426       case LEMON_GLP(DB):
   427       case LEMON_GLP(FX):
   428         return LEMON_glp(get_col_ub)(lp, i);
   429       default: ;
   430         return INF;
   431       }
   432   }
   433 
   434   void LpGlpk::_setRowBounds(int i, Value lb, Value ub)
   435   {
   436     //Bad parameter
   437     if (lb==INF || ub==-INF) {
   438       //FIXME error
   439     }
   440 
   441     if (lb == -INF){
   442       if (ub == INF){
   443         LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FR), lb, ub);
   444       }
   445       else{
   446         LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(UP), lb, ub);
   447       }
   448     }
   449     else{
   450       if (ub==INF){
   451         LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(LO), lb, ub);
   452 
   453       }
   454       else{
   455         if (lb == ub){
   456           LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FX), lb, ub);
   457         }
   458         else{
   459           LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(DB), lb, ub);
   460         }
   461       }
   462     }
   463 
   464     solved = false;
   465   }
   466 
   467   void LpGlpk::_getRowBounds(int i, Value &lb, Value &ub) const
   468   {
   469 
   470     int b=LEMON_glp(get_row_type)(lp, i);
   471     switch (b) {
   472     case LEMON_GLP(FR):
   473     case LEMON_GLP(UP):
   474       lb = -INF;
   475         break;
   476     default:
   477       lb=LEMON_glp(get_row_lb)(lp, i);
   478     }
   479 
   480     switch (b) {
   481     case LEMON_GLP(FR):
   482     case LEMON_GLP(LO):
   483       ub = INF;
   484         break;
   485     default:
   486       ub=LEMON_glp(get_row_ub)(lp, i);
   487     }
   488 
   489   }
   490 
   491   void LpGlpk::_setObjCoeff(int i, Value obj_coef)
   492   {
   493     //i=0 means the constant term (shift)
   494     LEMON_glp(set_obj_coef)(lp, i, obj_coef);
   495 
   496     solved = false;
   497   }
   498 
   499   LpGlpk::Value LpGlpk::_getObjCoeff(int i) const {
   500     //i=0 means the constant term (shift)
   501     return LEMON_glp(get_obj_coef)(lp, i);
   502   }
   503 
   504   void LpGlpk::_clearObj()
   505   {
   506     for (int i=0;i<=LEMON_glp(get_num_cols)(lp);++i){
   507       LEMON_glp(set_obj_coef)(lp, i, 0);
   508     }
   509 
   510     solved = false;
   511   }
   512 
   513   LpGlpk::SolveExitStatus LpGlpk::_solve()
   514   {
   515     // A way to check the problem to be solved
   516     //LEMON_glp(write_cpxlp(lp,"naittvan.cpx");
   517 
   518     LEMON_lpx(std_basis)(lp);
   519     int i =  LEMON_lpx(simplex)(lp);
   520 
   521     switch (i) {
   522     case LEMON_LPX(E_OK):
   523       solved = true;
   524       return SOLVED;
   525     default:
   526       return UNSOLVED;
   527     }
   528   }
   529 
   530   LpGlpk::Value LpGlpk::_getPrimal(int i) const
   531   {
   532     return LEMON_glp(get_col_prim)(lp,i);
   533   }
   534 
   535   LpGlpk::Value LpGlpk::_getDual(int i) const
   536   {
   537     return LEMON_glp(get_row_dual)(lp,i);
   538   }
   539 
   540   LpGlpk::Value LpGlpk::_getPrimalValue() const
   541   {
   542     return LEMON_glp(get_obj_val)(lp);
   543   }
   544   bool LpGlpk::_isBasicCol(int i) const
   545   {
   546     return (LEMON_glp(get_col_stat)(lp, i)==LEMON_GLP(BS));
   547   }
   548 
   549 
   550   LpGlpk::SolutionStatus LpGlpk::_getPrimalStatus() const
   551   {
   552     if (!solved) return UNDEFINED;
   553     int stat=  LEMON_lpx(get_status)(lp);
   554     switch (stat) {
   555     case LEMON_LPX(UNDEF)://Undefined (no solve has been run yet)
   556       return UNDEFINED;
   557     case LEMON_LPX(NOFEAS)://There is no feasible solution (primal, I guess)
   558     case LEMON_LPX(INFEAS)://Infeasible
   559       return INFEASIBLE;
   560     case LEMON_LPX(UNBND)://Unbounded
   561       return INFINITE;
   562     case LEMON_LPX(FEAS)://Feasible
   563       return FEASIBLE;
   564     case LEMON_LPX(OPT)://Feasible
   565       return OPTIMAL;
   566     default:
   567       return UNDEFINED; //to avoid gcc warning
   568       //FIXME error
   569     }
   570   }
   571 
   572   LpGlpk::SolutionStatus LpGlpk::_getDualStatus() const
   573   {
   574     if (!solved) return UNDEFINED;
   575     switch (LEMON_lpx(get_dual_stat)(lp)) {
   576     case LEMON_LPX(D_UNDEF)://Undefined (no solve has been run yet)
   577       return UNDEFINED;
   578     case LEMON_LPX(D_NOFEAS)://There is no dual feasible solution
   579 //    case LEMON_LPX(D_INFEAS://Infeasible
   580       return INFEASIBLE;
   581     case LEMON_LPX(D_FEAS)://Feasible
   582       switch (LEMON_lpx(get_status)(lp)) {
   583       case LEMON_LPX(NOFEAS):
   584         return INFINITE;
   585       case LEMON_LPX(OPT):
   586         return OPTIMAL;
   587       default:
   588         return FEASIBLE;
   589       }
   590     default:
   591       return UNDEFINED; //to avoid gcc warning
   592       //FIXME error
   593     }
   594   }
   595 
   596   LpGlpk::ProblemTypes LpGlpk::_getProblemType() const
   597   {
   598     if (!solved) return UNKNOWN;
   599       //int stat=  LEMON_glp(get_status(lp);
   600     int statp=  LEMON_lpx(get_prim_stat)(lp);
   601     int statd=  LEMON_lpx(get_dual_stat)(lp);
   602     if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_FEAS))
   603         return PRIMAL_DUAL_FEASIBLE;
   604     if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_NOFEAS))
   605         return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
   606     if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_FEAS))
   607         return PRIMAL_INFEASIBLE_DUAL_FEASIBLE;
   608     if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_NOFEAS))
   609         return PRIMAL_DUAL_INFEASIBLE;
   610     //In all other cases
   611     return UNKNOWN;
   612   }
   613 
   614   void LpGlpk::_setMax()
   615   {
   616     solved = false;
   617     LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MAX));
   618   }
   619 
   620   void LpGlpk::_setMin()
   621   {
   622     solved = false;
   623     LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MIN));
   624   }
   625 
   626   bool LpGlpk::_isMax() const
   627   {
   628     return (LEMON_glp(get_obj_dir)(lp)==LEMON_GLP(MAX));
   629   }
   630 
   631 
   632 
   633   void LpGlpk::messageLevel(int m)
   634   {
   635     LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_MSGLEV), m);
   636   }
   637 
   638   void LpGlpk::presolver(bool b)
   639   {
   640     LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_PRESOL), b);
   641   }
   642 
   643 
   644 } //END OF NAMESPACE LEMON