lemon/lp_cplex.cc
author deba
Mon, 18 Dec 2006 10:12:07 +0000
changeset 2330 9dccb1abc721
parent 2312 07e46cbb7d85
child 2361 f2ef1aa8189a
permissions -rw-r--r--
Better handling of inexact computation.
We do not use tolerance for excess, just for edges
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2006
     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 #include <iostream>
    20 #include<lemon/lp_cplex.h>
    21 
    22 ///\file
    23 ///\brief Implementation of the LEMON-CPLEX lp solver interface.
    24 namespace lemon {
    25   
    26   LpCplex::LpCplex() : LpSolverBase() {
    27 
    28     //    env = CPXopenCPLEXdevelop(&status);     
    29     env = CPXopenCPLEX(&status);     
    30     lp = CPXcreateprob(env, &status, "LP problem");
    31   }
    32   
    33   LpCplex::~LpCplex() {
    34     CPXfreeprob(env,&lp); 
    35     CPXcloseCPLEX(&env); 
    36   }
    37   
    38   LpSolverBase &LpCplex::_newLp() 
    39   {
    40     //The first approach opens a new environment
    41     LpCplex* newlp=new LpCplex();
    42     return *newlp;
    43   }
    44 
    45   LpSolverBase &LpCplex::_copyLp() {
    46     ///\bug FixID data is not copied!
    47     //The first approach opens a new environment
    48     LpCplex* newlp=new LpCplex();
    49     //The routine CPXcloneprob can be used to create a new CPLEX problem 
    50     //object and copy all the problem data from an existing problem 
    51     //object to it. Solution and starting information is not copied.
    52     newlp->lp = CPXcloneprob(env, lp, &status);
    53     return *newlp;
    54   }
    55 
    56   int LpCplex::_addCol()
    57   {
    58     int i = CPXgetnumcols(env, lp);
    59     Value lb[1],ub[1];
    60     lb[0]=-INF;//-CPX_INFBOUND;
    61     ub[0]=INF;//CPX_INFBOUND;
    62     status = CPXnewcols(env, lp, 1, NULL, lb, ub, NULL, NULL);
    63     return i;
    64   }
    65 
    66   
    67   int LpCplex::_addRow() 
    68   {
    69     //We want a row that is not constrained
    70     char sense[1];
    71     sense[0]='L';//<= constraint
    72     Value rhs[1];
    73     rhs[0]=INF;
    74     int i = CPXgetnumrows(env, lp);
    75     status = CPXnewrows(env, lp, 1, rhs, sense, NULL, NULL);
    76     return i;
    77   }
    78 
    79 
    80   void LpCplex::_eraseCol(int i) {
    81     CPXdelcols(env, lp, i, i);
    82   }
    83   
    84   void LpCplex::_eraseRow(int i) {
    85     CPXdelrows(env, lp, i, i);
    86   }
    87   
    88   void LpCplex::_getColName(int col, std::string &name)
    89   {
    90     ///\bug Untested
    91     int storespace;
    92     CPXgetcolname(env, lp, 0, 0, 0, &storespace, col, col);
    93 
    94     storespace *= -1;
    95     char buf[storespace];
    96     char *names[1];
    97     int dontcare;
    98     ///\bug return code unchecked for error
    99     CPXgetcolname(env, lp, names, buf, storespace, &dontcare, col, col);
   100     name = names[0];
   101   }
   102   
   103   void LpCplex::_setColName(int col, const std::string &name)
   104   {
   105     ///\bug Untested
   106     char *names[1];
   107     names[0] = const_cast<char*>(name.c_str());
   108     ///\bug return code unchecked for error
   109     CPXchgcolname(env, lp, 1, &col, names);    
   110   }
   111   
   112   ///\warning Data at index 0 is ignored in the arrays.
   113   void LpCplex::_setRowCoeffs(int i, LpRowIterator b, LpRowIterator e)
   114   {
   115     std::vector<int> indices;
   116     std::vector<int> rowlist;
   117     std::vector<Value> values;
   118 
   119     for(LpRowIterator it=b; it!=e; ++it) {
   120       indices.push_back(it->first);
   121       values.push_back(it->second);
   122       rowlist.push_back(i);
   123     }
   124 
   125     status = CPXchgcoeflist(env, lp, values.size(), 
   126 			    &rowlist[0], &indices[0], &values[0]); 
   127   }
   128   
   129   void LpCplex::_setColCoeffs(int i, LpColIterator b, LpColIterator e)
   130   {
   131     std::vector<int> indices;
   132     std::vector<int> collist;
   133     std::vector<Value> values;
   134 
   135     for(LpColIterator it=b; it!=e; ++it) {
   136       indices.push_back(it->first);
   137       values.push_back(it->second);
   138       collist.push_back(i);
   139     }
   140 
   141     status = CPXchgcoeflist(env, lp, values.size(), 
   142 			    &indices[0], &collist[0], &values[0]); 
   143   }
   144   
   145   void LpCplex::_setCoeff(int row, int col, Value value) 
   146   {
   147     CPXchgcoef(env, lp, row, col, value);
   148   }
   149 
   150   void LpCplex::_setColLowerBound(int i, Value value)
   151   {
   152     int indices[1];
   153     indices[0]=i;
   154     char lu[1];
   155     lu[0]='L';
   156     Value bd[1];
   157     bd[0]=value;
   158     status = CPXchgbds(env, lp, 1, indices, lu, bd);
   159  
   160   }
   161   
   162   void LpCplex::_setColUpperBound(int i, Value value)
   163   {
   164     int indices[1];
   165     indices[0]=i;
   166     char lu[1];
   167     lu[0]='U';
   168     Value bd[1];
   169     bd[0]=value;
   170     status = CPXchgbds(env, lp, 1, indices, lu, bd);
   171   }
   172 
   173   //This will be easier to implement
   174   void LpCplex::_setRowBounds(int i, Value lb, Value ub)
   175   {
   176     //Bad parameter
   177     if (lb==INF || ub==-INF) {
   178       //FIXME error
   179     }
   180     
   181     int cnt=1;
   182     int indices[1];
   183     indices[0]=i;
   184     char sense[1];
   185 
   186     if (lb==-INF){
   187       sense[0]='L';
   188       CPXchgsense(env, lp, cnt, indices, sense);
   189       CPXchgcoef(env, lp, i, -1, ub);
   190       
   191     }
   192     else{
   193       if (ub==INF){
   194 	sense[0]='G';
   195 	CPXchgsense(env, lp, cnt, indices, sense);
   196 	CPXchgcoef(env, lp, i, -1, lb);
   197       }
   198       else{
   199 	if (lb == ub){
   200 	  sense[0]='E';
   201 	  CPXchgsense(env, lp, cnt, indices, sense);
   202 	  CPXchgcoef(env, lp, i, -1, lb);
   203 	}
   204 	else{
   205 	  sense[0]='R';
   206 	  CPXchgsense(env, lp, cnt, indices, sense);
   207 	  CPXchgcoef(env, lp, i, -1, lb);
   208 	  CPXchgcoef(env, lp, i, -2, ub-lb);	  
   209 	}
   210       }
   211     }
   212   }
   213 
   214 //   void LpCplex::_setRowLowerBound(int i, Value value)
   215 //   {
   216 //     //Not implemented, obsolete
   217 //   }
   218   
   219 //   void LpCplex::_setRowUpperBound(int i, Value value)
   220 //   {
   221 //     //Not implemented, obsolete
   222 // //     //TODO Ezt kell meg megirni
   223 // //     //type of the problem
   224 // //     char sense[1];
   225 // //     status = CPXgetsense(env, lp, sense, i, i);
   226 // //     Value rhs[1];
   227 // //     status = CPXgetrhs(env, lp, rhs, i, i);
   228 
   229 // //     switch (sense[0]) {
   230 // //     case 'L'://<= constraint
   231 // //       break;
   232 // //     case 'E'://= constraint
   233 // //       break;
   234 // //     case 'G'://>= constraint
   235 // //       break;
   236 // //     case 'R'://ranged constraint
   237 // //       break;
   238 // //     default: ;
   239 // //       //FIXME error
   240 // //     }
   241 
   242 // //     status = CPXchgcoef(env, lp, i, -2, value_rng);
   243 //   }
   244   
   245   void LpCplex::_setObjCoeff(int i, Value obj_coef)
   246   {
   247     CPXchgcoef(env, lp, -1, i, obj_coef);
   248   }
   249 
   250   void LpCplex::_clearObj()
   251   {
   252     for (int i=0;i< CPXgetnumcols(env, lp);++i){
   253       CPXchgcoef(env, lp, -1, i, 0);
   254     }
   255     
   256   }
   257   // The routine returns zero unless an error occurred during the
   258   // optimization. Examples of errors include exhausting available
   259   // memory (CPXERR_NO_MEMORY) or encountering invalid data in the
   260   // CPLEX problem object (CPXERR_NO_PROBLEM). Exceeding a
   261   // user-specified CPLEX limit, or proving the model infeasible or
   262   // unbounded, are not considered errors. Note that a zero return
   263   // value does not necessarily mean that a solution exists. Use query
   264   // routines CPXsolninfo, CPXgetstat, and CPXsolution to obtain
   265   // further information about the status of the optimization.
   266   LpCplex::SolveExitStatus LpCplex::_solve()
   267   {
   268     //CPX_PARAM_LPMETHOD 
   269     status = CPXlpopt(env, lp);
   270     //status = CPXprimopt(env, lp);
   271 #if CPX_VERSION >= 800
   272     if (status)
   273     {
   274       return UNSOLVED;
   275     }
   276     else
   277     {
   278       switch (CPXgetstat(env, lp))
   279       {
   280         case CPX_STAT_OPTIMAL:
   281         case CPX_STAT_INFEASIBLE:
   282         case CPX_STAT_UNBOUNDED:
   283           return SOLVED;
   284         default:
   285           return UNSOLVED;
   286       }
   287     }
   288 #else
   289     if (status == 0){
   290       //We want to exclude some cases
   291       switch (CPXgetstat(env, lp)){
   292       case CPX_OBJ_LIM:
   293       case CPX_IT_LIM_FEAS:
   294       case CPX_IT_LIM_INFEAS:               
   295       case CPX_TIME_LIM_FEAS:
   296       case CPX_TIME_LIM_INFEAS:
   297 	return UNSOLVED;
   298       default:
   299 	return SOLVED; 
   300       }
   301     }
   302     else{
   303       return UNSOLVED;
   304     }
   305 #endif
   306   }
   307 
   308   LpCplex::Value LpCplex::_getPrimal(int i)
   309   {
   310     Value x;
   311     CPXgetx(env, lp, &x, i, i);
   312     return x;
   313   }
   314 
   315   LpCplex::Value LpCplex::_getDual(int i)
   316   {
   317     Value y;
   318     CPXgetpi(env, lp, &y, i, i);
   319     return y;
   320   }
   321   
   322   LpCplex::Value LpCplex::_getPrimalValue()
   323   {
   324     Value objval;
   325     //method = CPXgetmethod (env, lp);
   326     //printf("CPXgetprobtype %d \n",CPXgetprobtype(env,lp));
   327     status = CPXgetobjval(env, lp, &objval);
   328     //printf("Objective value: %g \n",objval);
   329     return objval;
   330   }
   331   bool LpCplex::_isBasicCol(int i) {
   332     int cstat[CPXgetnumcols(env, lp)];
   333     CPXgetbase(env, lp, cstat, NULL);
   334     return (cstat[i]==CPX_BASIC);
   335   }  
   336 
   337 //7.5-os cplex statusai (Vigyazat: a 9.0-asei masok!)
   338 // This table lists the statuses, returned by the CPXgetstat() routine, for solutions to LP problems or mixed integer problems. If no solution exists, the return value is zero.
   339 
   340 // For Simplex, Barrier  
   341 // 1  	CPX_OPTIMAL  
   342 // 	 Optimal solution found  
   343 // 2  	CPX_INFEASIBLE  
   344 // 	 Problem infeasible  
   345 // 3    CPX_UNBOUNDED  
   346 // 	 Problem unbounded  
   347 // 4  	CPX_OBJ_LIM  
   348 // 	 Objective limit exceeded in Phase II  
   349 // 5  	CPX_IT_LIM_FEAS  
   350 // 	 Iteration limit exceeded in Phase II  
   351 // 6  	CPX_IT_LIM_INFEAS  
   352 // 	 Iteration limit exceeded in Phase I  
   353 // 7  	CPX_TIME_LIM_FEAS  
   354 // 	 Time limit exceeded in Phase II  
   355 // 8  	CPX_TIME_LIM_INFEAS  
   356 // 	 Time limit exceeded in Phase I  
   357 // 9  	CPX_NUM_BEST_FEAS  
   358 // 	 Problem non-optimal, singularities in Phase II  
   359 // 10 	CPX_NUM_BEST_INFEAS  
   360 // 	 Problem non-optimal, singularities in Phase I  
   361 // 11 	CPX_OPTIMAL_INFEAS  
   362 // 	 Optimal solution found, unscaled infeasibilities  
   363 // 12 	CPX_ABORT_FEAS  
   364 // 	 Aborted in Phase II  
   365 // 13 	CPX_ABORT_INFEAS  
   366 // 	 Aborted in Phase I  
   367 // 14  	CPX_ABORT_DUAL_INFEAS  
   368 // 	 Aborted in barrier, dual infeasible  
   369 // 15  	CPX_ABORT_PRIM_INFEAS  
   370 // 	 Aborted in barrier, primal infeasible  
   371 // 16  	CPX_ABORT_PRIM_DUAL_INFEAS  
   372 // 	 Aborted in barrier, primal and dual infeasible  
   373 // 17  	CPX_ABORT_PRIM_DUAL_FEAS  
   374 // 	 Aborted in barrier, primal and dual feasible  
   375 // 18  	CPX_ABORT_CROSSOVER  
   376 // 	 Aborted in crossover  
   377 // 19  	CPX_INForUNBD  
   378 // 	 Infeasible or unbounded  
   379 // 20   CPX_PIVOT
   380 //       User pivot used
   381 //
   382 //     Ezeket hova tegyem:
   383 // ??case CPX_ABORT_DUAL_INFEAS           
   384 // ??case CPX_ABORT_CROSSOVER             
   385 // ??case CPX_INForUNBD                   
   386 // ??case CPX_PIVOT              
   387          
   388 //Some more interesting stuff:
   389 
   390 // CPX_PARAM_LPMETHOD  1062  int  LPMETHOD
   391 // 0 Automatic 
   392 // 1 Primal Simplex 
   393 // 2 Dual Simplex 
   394 // 3 Network Simplex 
   395 // 4 Standard Barrier 
   396 // Default: 0 
   397 // Description: Method for linear optimization. 
   398 // Determines which algorithm is used when CPXlpopt() (or "optimize" in the Interactive Optimizer) is called. Currently the behavior of the "Automatic" setting is that CPLEX simply invokes the dual simplex method, but this capability may be expanded in the future so that CPLEX chooses the method based on problem characteristics 
   399   //Hulye cplex
   400   void statusSwitch(CPXENVptr env,int& stat){
   401 #if CPX_VERSION < 900
   402     int lpmethod;
   403     CPXgetintparam (env,CPX_PARAM_LPMETHOD,&lpmethod);
   404     if (lpmethod==2){
   405       if (stat==CPX_UNBOUNDED){
   406 	stat=CPX_INFEASIBLE;
   407       }
   408       else{
   409 	if (stat==CPX_INFEASIBLE)
   410 	  stat=CPX_UNBOUNDED;
   411       }
   412     }
   413 #endif
   414   }
   415 
   416   LpCplex::SolutionStatus LpCplex::_getPrimalStatus()
   417   {
   418     //Unboundedness not treated well: the following is from cplex 9.0 doc
   419     // About Unboundedness
   420 
   421     // The treatment of models that are unbounded involves a few
   422     // subtleties. Specifically, a declaration of unboundedness means that
   423     // ILOG CPLEX has determined that the model has an unbounded
   424     // ray. Given any feasible solution x with objective z, a multiple of
   425     // the unbounded ray can be added to x to give a feasible solution
   426     // with objective z-1 (or z+1 for maximization models). Thus, if a
   427     // feasible solution exists, then the optimal objective is
   428     // unbounded. Note that ILOG CPLEX has not necessarily concluded that
   429     // a feasible solution exists. Users can call the routine CPXsolninfo
   430     // to determine whether ILOG CPLEX has also concluded that the model
   431     // has a feasible solution.
   432 
   433     int stat = CPXgetstat(env, lp);
   434 #if CPX_VERSION >= 800
   435     switch (stat)
   436     {
   437       case CPX_STAT_OPTIMAL:
   438         return OPTIMAL;
   439       case CPX_STAT_UNBOUNDED:
   440         return INFINITE;
   441       case CPX_STAT_INFEASIBLE:
   442         return INFEASIBLE;
   443       default:
   444         return UNDEFINED;
   445     }
   446 #else
   447     statusSwitch(env,stat);
   448     //CPXgetstat(env, lp);
   449     //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL);
   450     switch (stat) {
   451     case 0:
   452       return UNDEFINED; //Undefined
   453     case CPX_OPTIMAL://Optimal
   454       return OPTIMAL;
   455     case CPX_UNBOUNDED://Unbounded
   456       return INFEASIBLE;//In case of dual simplex
   457       //return INFINITE;
   458     case CPX_INFEASIBLE://Infeasible 
   459  //    case CPX_IT_LIM_INFEAS:
   460 //     case CPX_TIME_LIM_INFEAS:
   461 //     case CPX_NUM_BEST_INFEAS:             
   462 //     case CPX_OPTIMAL_INFEAS:              
   463 //     case CPX_ABORT_INFEAS:                
   464 //     case CPX_ABORT_PRIM_INFEAS:           
   465 //     case CPX_ABORT_PRIM_DUAL_INFEAS:      
   466       return INFINITE;//In case of dual simplex
   467       //return INFEASIBLE;
   468 //     case CPX_OBJ_LIM:                    
   469 //     case CPX_IT_LIM_FEAS:             
   470 //     case CPX_TIME_LIM_FEAS:                
   471 //     case CPX_NUM_BEST_FEAS:                
   472 //     case CPX_ABORT_FEAS:                  
   473 //     case CPX_ABORT_PRIM_DUAL_FEAS:        
   474 //       return FEASIBLE;
   475     default:
   476       return UNDEFINED; //Everything else comes here
   477       //FIXME error
   478     }
   479 #endif
   480   }
   481 
   482 //9.0-as cplex verzio statusai
   483 // CPX_STAT_ABORT_DUAL_OBJ_LIM
   484 // CPX_STAT_ABORT_IT_LIM
   485 // CPX_STAT_ABORT_OBJ_LIM
   486 // CPX_STAT_ABORT_PRIM_OBJ_LIM
   487 // CPX_STAT_ABORT_TIME_LIM
   488 // CPX_STAT_ABORT_USER
   489 // CPX_STAT_FEASIBLE_RELAXED
   490 // CPX_STAT_INFEASIBLE
   491 // CPX_STAT_INForUNBD
   492 // CPX_STAT_NUM_BEST
   493 // CPX_STAT_OPTIMAL
   494 // CPX_STAT_OPTIMAL_FACE_UNBOUNDED
   495 // CPX_STAT_OPTIMAL_INFEAS
   496 // CPX_STAT_OPTIMAL_RELAXED
   497 // CPX_STAT_UNBOUNDED
   498 
   499   LpCplex::SolutionStatus LpCplex::_getDualStatus()
   500   {
   501     int stat = CPXgetstat(env, lp);
   502 #if CPX_VERSION >= 800
   503     switch (stat)
   504     {
   505       case CPX_STAT_OPTIMAL:
   506         return OPTIMAL;
   507       case CPX_STAT_UNBOUNDED:
   508         return INFEASIBLE;
   509       default:
   510         return UNDEFINED;
   511     }
   512 #else
   513     statusSwitch(env,stat);
   514     switch (stat) {
   515     case 0:
   516       return UNDEFINED; //Undefined
   517     case CPX_OPTIMAL://Optimal
   518       return OPTIMAL;
   519     case CPX_UNBOUNDED:
   520      return INFEASIBLE;
   521     default:
   522       return UNDEFINED; //Everything else comes here
   523       //FIXME error
   524     }
   525 #endif
   526   }
   527 
   528   LpCplex::ProblemTypes LpCplex::_getProblemType()
   529   {
   530     int stat = CPXgetstat(env, lp);
   531 #if CPX_VERSION >= 800
   532     switch (stat)
   533     {
   534       case CPX_STAT_OPTIMAL:
   535 	return PRIMAL_DUAL_FEASIBLE;
   536       case CPX_STAT_UNBOUNDED:
   537  	return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
   538       default:
   539         return UNKNOWN;
   540     }
   541 #else
   542     switch (stat) {
   543     case CPX_OPTIMAL://Optimal
   544 	return PRIMAL_DUAL_FEASIBLE;
   545     case CPX_UNBOUNDED:
   546  	return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
   547 // 	return PRIMAL_INFEASIBLE_DUAL_FEASIBLE;
   548 // 	return PRIMAL_DUAL_INFEASIBLE;
   549 
   550 //Seems to be that this is all we can say for sure
   551     default:
   552 	//In all other cases
   553 	return UNKNOWN;
   554       //FIXME error
   555     }
   556 #endif
   557   }
   558 
   559   void LpCplex::_setMax()
   560   {
   561     CPXchgobjsen(env, lp, CPX_MAX);
   562    }
   563   void LpCplex::_setMin()
   564   {
   565     CPXchgobjsen(env, lp, CPX_MIN);
   566    }
   567   
   568 } //namespace lemon
   569