lemon/lp_glpk.cc
author kpeter
Mon, 18 Feb 2008 03:32:06 +0000
changeset 2575 e866e288cba6
parent 2441 d8d6ab871608
child 2591 3b4d5bc3b4fb
permissions -rw-r--r--
Major improvements in NetworkSimplex.

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