SoPlex Doxygen Documentation
spxsolver.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2012 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file spxsolver.h
17  * @brief main LP solver class
18  */
19 #ifndef _SPXSOLVER_H_
20 #define _SPXSOLVER_H_
21 
22 #include <assert.h>
23 #include <iostream>
24 #include <iomanip>
25 #include <sstream>
26 
27 #include "spxdefines.h"
28 #include "timer.h"
29 #include "spxlp.h"
30 #include "spxbasis.h"
31 #include "array.h"
32 #include "random.h"
33 #include "unitvector.h"
34 #include "updatevector.h"
35 
36 #define SPARSITYTHRESHOLD 0.5 /**< percentage of basic infeasibilities that is considered sparse */
37 #define DENSEROUNDS 5 /**< number of refactorization until sparsity is tested again */
38 #define SPARSITY_TRADEOFF 0.8 /**< threshold to decide whether slack or structural variables enter the basis
39  * slack variables are preferred if their pricing value is not worse than
40  * SPARSITY_TRADEOFF * the best structural value
41  */
42 
43 namespace soplex
44 {
45 class SPxPricer;
46 class SPxRatioTester;
47 class SPxStarter;
48 class DVector_exact;
49 
50 /**@brief Sequential object-oriented SimPlex.
51  @ingroup Algo
52 
53  SPxSolver is an LP solver class using the revised Simplex algorithm. It
54  provids two basis representations, namely a column basis and a row basis
55  (see #Representation). For both representations, a primal and
56  dual algorithm is available (see \ref Type).
57 
58  In addition, SPxSolver can be custumized with various respects:
59  - pricing algorithms using SPxPricer
60  - ratio test using class SPxRatioTester
61  - computation of a start basis using class SPxStarter
62  - preprocessing of the LP using class SPxSimplifier
63  - termination criteria by overriding
64 
65  SPxSolver is derived from SPxLP that is used to store the LP to be solved.
66  Hence, the LPs solved with SPxSolver have the general format
67 
68  \f[
69  \begin{array}{rl}
70  \hbox{max} & \mbox{maxObj}^T x \\
71  \hbox{s.t.} & \mbox{lhs} \le Ax \le \mbox{rhs} \\
72  & \mbox{low} \le x \le \mbox{up}
73  \end{array}
74  \f]
75 
76  Also, SPxLP provide all manipulation methods for the LP. They allow
77  SPxSolver to be used within cutting plane algorithms.
78 */
79 class SPxSolver : public SPxLP, protected SPxBasis
80 {
81  friend class SPxFastRT;
82  friend class SPxBoundFlippingRT;
83 
84 public:
85 
86  //-----------------------------
87  /**@name Data Types */
88  //@{
89  /// LP basis representation.
90  /** Solving LPs with the Simplex algorithm requires the definition of a
91  * \em basis. A basis can be defined as a set of column vectors or a
92  * set of row vectors building a nonsingular matrix. We will refer to
93  * the first case as the \em columnwise representation and the latter
94  * case will be called the \em rowwise representation.
95  *
96  * Type Representation determines the representation of SPxSolver, i.e.
97  * a columnwise (#COLUMN == 1) or rowwise (#ROW == -1) one.
98  */
99  enum Representation
100  {
101  ROW = -1, ///< rowwise representation.
102  COLUMN = 1 ///< columnwise representation.
103  };
104 
105  /// Algorithmic type.
106  /** SPxSolver uses the reviesed Simplex algorithm to solve LPs.
107  * Mathematically, one distinguishes the \em primal from the
108  * \em dual algorihm. Algorithmically, these relate to the two
109  * types #ENTER or #LEAVE. How they relate, depends on the chosen
110  * basis representation. This is desribed by the following table:
111  *
112  * <TABLE>
113  * <TR><TD>&nbsp;</TD><TD>ENTER </TD><TD>LEAVE </TD></TR>
114  * <TR><TD>ROW </TD><TD>DUAL </TD><TD>PRIMAL</TD></TR>
115  * <TR><TD>COLUMN</TD><TD>PRIMAL</TD><TD>DUAL </TD></TR>
116  * </TABLE>
117  */
118  enum Type
119  {
120  /// Entering Simplex.
121  /** The Simplex loop for the entering Simplex can be sketched
122  * as follows:
123  * - \em Pricing : Select a variable to #ENTER the basis.
124  * - \em Ratio-Test : Select variable to #LEAVE the
125  * basis such that the basis remains feasible.
126  * - Perform the basis update.
127  */
128  ENTER = -1,
129  /// Leaving Simplex.
130  /** The Simplex loop for the leaving Simplex can be sketched
131  * as follows:
132  * - \em Pricing: Select a variable to #LEAVE the basis.
133  * - \em Ratio-Test: Select variable to #ENTER the
134  * basis such that the basis remains priced.
135  * - Perform the basis update.
136  */
137  LEAVE = 1
138  };
139 
140  /// Pricing type.
141  /** In case of the #ENTER%ing Simplex algorithm, for performance
142  * reasons it may be advisable not to compute and maintain up to
143  * date vectors #pVec() and #test() and instead compute only some
144  * of its elements explicitely. This is controled by the #Pricing type.
145  */
146  enum Pricing
147  {
148  /// Full pricing.
149  /** If #FULL pricing in selected for the #ENTER%ing Simplex,
150  * vectors #pVec() and #test() are kept up to date by
151  * SPxSolver. An SPxPricer only needs to select an Id such
152  * that the #test() or #coTest() value is < 0.
153  */
154  FULL,
155  /// Partial pricing.
156  /** When #PARTIAL pricing in selected for the #ENTER%ing
157  * Simplex, vectors #pVec() and #test() are not set up and
158  * updated by SPxSolver. However, vectors #coPvec() and
159  * #coTest() are still kept up to date by SPxSolver.
160  * An SPxPricer object needs to compute the values for
161  * #pVec() and #test() itself in order to select an
162  * appropriate pivot with #test() < 0. Methods \ref computePvec(int)
163  * "computePvec(i)" and \ref computeTest(int) "computeTest(i)"
164  * will assist the used to do so. Note
165  * that it may be feasible for a pricer to return an Id with
166  * #test() > 0; such will be rejected by SPxSolver.
167  */
168  PARTIAL
169  };
170 
171  enum VarStatus
172  {
173  ON_UPPER, ///< variable set to its upper bound.
174  ON_LOWER, ///< variable set to its lower bound.
175  FIXED, ///< variable fixed to identical bounds.
176  ZERO, ///< free variable fixed to zero.
177  BASIC, ///< variable is basic.
178  UNDEFINED ///< nothing known about basis status (possibly due to a singular basis in transformed problem)
179  };
180 
181  /**@todo In spxchange, change the status to
182  if (m_status > 0) m_status = REGULAR;
183  */
184  enum Status
185  {
186 
187  ERROR = -13, ///< an error occured.
188  NO_RATIOTESTER = -12, ///< No ratiotester loaded
189  NO_PRICER = -11, ///< No pricer loaded
190  NO_SOLVER = -10, ///< No linear solver loaded
191  NOT_INIT = -9, ///< not initialised error
192  ABORT_CYCLING = -8, ///< solve() aborted due to detection of cycling.
193  ABORT_TIME = -7, ///< solve() aborted due to time limit.
194  ABORT_ITER = -6, ///< solve() aborted due to iteration limit.
195  ABORT_VALUE = -5, ///< solve() aborted due to objective limit.
196  SINGULAR = -4, ///< Basis is singular, numerical troubles?
197  NO_PROBLEM = -3, ///< No Problem has been loaded.
198  REGULAR = -2, ///< LP has a usable Basis (maybe LP is changed).
199  RUNNING = -1, ///< algorithm is running
200  UNKNOWN = 0, ///< nothing known on loaded problem.
201  OPTIMAL = 1, ///< LP has been solved to optimality.
202  UNBOUNDED = 2, ///< LP has been proven to be primal unbounded.
203  INFEASIBLE = 3 ///< LP has been proven to be primal infeasible.
204  };
205 
206  //@}
207 
208 private:
209 
210  //-----------------------------
211  /**@name Private data */
212  //@{
213  Type theType; ///< entering or leaving algortihm.
214  Pricing thePricing; ///< full or partial pricing.
215  Representation theRep; ///< row or column representation.
216  Timer theTime; ///< time spent in last call to method solve()
217  Real theCumulativeTime; ///< cumulative time spent in all calls to method solve()
218  int maxIters; ///< maximum allowed iterations.
219  int maxRefines; ///< maximum allowed refinement rounds.
220  Real maxTime; ///< maximum allowed time.
221  Real objLimit; ///< objective value limit.
222  Status m_status; ///< status of algorithm.
223 
224  Real m_entertol; ///< feasibility tolerance maintained during entering algorithm
225  Real m_leavetol; ///< feasibility tolerance maintained during leaving algorithm
226  Real m_irthreshold; ///< iterative refinement threshold
227  Real theShift; ///< sum of all shifts applied to any bound.
228  Real lastShift; ///< for forcing feasibility.
229  int m_maxCycle; ///< maximum steps before cycling is detected.
230  int m_numCycle; ///< actual number of degenerate steps so far.
231  bool initialized; ///< true, if all vectors are setup.
232 
233  Vector* solveVector2; ///< when 2 systems are to solve at a time
234  SSVector* solveVector2rhs; ///< when 2 systems are to solve at a time
235  Vector* solveVector3; ///< when 3 systems are to be solved at a time; typically reserved for bound flipping ratio test (basic solution will be modified!)
236  SSVector* solveVector3rhs; ///< when 3 systems are to be solved at a time; typically reserved for bound flipping ratio test (basic solution will be modified!)
237  Vector* coSolveVector2; ///< when 2 systems are to solve at a time
238  SSVector* coSolveVector2rhs; ///< when 2 systems are to solve at a time
239 
240  bool freePricer; ///< true iff thepricer should be freed inside of object
241  bool freeRatioTester; ///< true iff theratiotester should be freed inside of object
242  bool freeStarter; ///< true iff thestarter should be freed inside of object
243 
244  /* Store the index of a leaving variable if only an instable entering variable has been found.
245  instableLeave == true iff this instable basis change should be performed.
246  (see spxsolve.cpp and leave.cpp) */
247  int instableLeaveNum;
248  bool instableLeave;
249  Real instableLeaveVal;
250  //@}
251 
252 protected:
253 
254  //-----------------------------
255  /**@name Protected data */
256  //@{
257  Array < UnitVector > unitVecs; ///< array of unit vectors
258  const SVSet* thevectors; ///< the LP vectors according to representation
259  const SVSet* thecovectors; ///< the LP coVectors according to representation
260 
261  DVector primRhs; ///< rhs vector for computing the primal vector
262  UpdateVector primVec; ///< primal vector
263  DVector dualRhs; ///< rhs vector for computing the dual vector
264  UpdateVector dualVec; ///< dual vector
265  UpdateVector addVec; ///< storage for thePvec = &addVec
266 
267  DVector theURbound; ///< Upper Row Feasibility bound
268  DVector theLRbound; ///< Lower Row Feasibility bound
269  DVector theUCbound; ///< Upper Column Feasibility bound
270  DVector theLCbound; ///< Lower Column Feasibility bound
271 
272  /** In entering Simplex algorithm, the ratio test must ensure that all
273  * \em basic variables remain within their feasibility bounds. To give fast
274  * acces to them, the bounds of basic variables are copied into the
275  * following two vectors.
276  */
277  DVector theUBbound; ///< Upper Basic Feasibility bound
278  DVector theLBbound; ///< Lower Basic Feasibility bound
279 
280  DVector* theFrhs;
281  UpdateVector* theFvec;
282 
283  DVector* theCoPrhs;
284  UpdateVector* theCoPvec;
285  UpdateVector* thePvec;
286 
287  UpdateVector* theRPvec; ///< row pricing vector
288  UpdateVector* theCPvec; ///< column pricing vector
289 
290  // The following vectors serve for the virtualization of shift bounds
291  //@todo In prinziple this schould be references.
292  DVector* theUbound; ///< Upper bound for vars
293  DVector* theLbound; ///< Lower bound for vars
294  DVector* theCoUbound; ///< Upper bound for covars
295  DVector* theCoLbound; ///< Lower bound for covars
296 
297  // The following vectors serve for the virtualization of testing vectors
298  DVector theCoTest;
299  DVector theTest;
300 
301  DSVector primalRay; ///< stores primal ray in case of unboundedness
302  DSVector dualFarkas; ///< stores dual farkas proof in case of infeasibility
303 
304  int leaveCount; ///< number of LEAVE iterations
305  int enterCount; ///< number of ENTER iterations
306 
307  int boundflips; ///< number of performed bound flips
308  int totalboundflips; ///< total number of bound flips
309 
310  SPxPricer* thepricer;
311  SPxRatioTester* theratiotester;
312  SPxStarter* thestarter;
313  //@}
314 
315  //-----------------------------
316  /**@name Precision */
317  //@{
318  /// is the solution precise enough, or should we increase delta() ?
319  virtual bool precisionReached(Real& newpricertol) const;
320  //@}
321 
322 public:
323 
324  /** For the leaving Simplex algorithm this vector contains the indices of infeasible basic variables;
325  * for the entering Simplex algorithm this vector contains the indices of infeasible slack variables.
326  */
327  DIdxSet infeasibilities;
328  /**For the entering Simplex algorithm these vectors contains the indices of infeasible basic variables.
329  */
330  DIdxSet infeasibilitiesCo;
331 
332  /** Binary vectors to store whether basic indices are infeasible
333  * the i-th entry equals false, if the i-th basic variable is not infeasible
334  * the i-th entry equals true, if the i-th basic variable is infeasible
335  */
336  Array<bool> isInfeasible; ///< belongs to \ref soplex::SPxSolver::infeasibilities "infeasibilities" in the leaving and entering Simplex
337  Array<bool> isInfeasibleCo; ///< belongs to \ref soplex::SPxSolver::infeasibilitiesCo "infeasibilitiesCo" in the entering Simplex
338 
339  /// These values enable or disable sparse pricing
340  bool sparsePricingLeave; ///< true if sparsePricing is turned on in the leaving Simplex
341  bool sparsePricingEnter; ///< true if sparsePricing is turned on in the entering Simplex for slack variables
342  bool sparsePricingEnterCo; ///< true if sparsePricing is turned on in the entering Simplex
343 
344  int remainingRoundsLeave; ///< number of dense rounds/refactorizations until sparsePricing is enabled again
345  int remainingRoundsEnter;
346  int remainingRoundsEnterCo;
347 
348  int sparsityThresholdLeave; ///< maximum number of infeasibilities that is considered sparse for leaving Simplex
349  int sparsityThresholdEnter; ///< maximum number of infeasibilities that is considered sparse for entering Simplex (dim)
350  int sparsityThresholdEnterCo; ///< maximum number of infeasibilities that is considered sparse for entering Simplex (coDim)
351 
352  //-----------------------------
353  /**@name Access */
354  //@{
355  /// return the version of SPxSolver as number like 123 for 1.2.3
356  int version() const
357  {
358  return SOPLEX_VERSION;
359  }
360  /// return the internal subversion of SPxSolver as number
361  int subversion() const
362  {
363  return SOPLEX_SUBVERSION;
364  }
365  /// return the current basis representation.
366  Representation rep() const
367  {
368  return theRep;
369  }
370 
371  /// return current Type.
372  Type type() const
373  {
374  return theType;
375  }
376 
377  /// return current Pricing.
378  Pricing pricing() const
379  {
380  return thePricing;
381  }
382 
383  /// return current starter.
384  SPxStarter* starter() const
385  {
386  return thestarter;
387  }
388  //@}
389 
390  //-----------------------------
391  /**@name Setup
392  * Before solving an LP with an instance of SPxSolver,
393  * the following steps must be performed:
394  *
395  * -# Load the LP by copying an external LP or reading it from an
396  * input stream.
397  * -# Setup the pricer to use by loading an \ref soplex::SPxPricer
398  * "SPxPricer" object (if not already done in a previous call).
399  * -# Setup the ratio test method to use by loading an
400  * \ref soplex::SPxRatioTester "SPxRatioTester" object
401  * (if not already done in a previous call).
402  * -# Setup the linear system solver to use by loading an
403  * \ref soplex::SLinSolver "SLinSolver" object
404  * (if not already done in a previous call).
405  * -# Optionally setup an start basis generation method by loading an
406  * \ref soplex::SPxStarter "SPxStarter" object.
407  * -# Optionally setup a start basis by loading a
408  * \ref soplex::SPxBasis::Desc "SPxBasis::Desc" object.
409  * -# Optionally switch to another basis
410  * \ref soplex::SPxSolver::Representation "Representation"
411  * by calling method \ref soplex::SPxSolver::setRep() "setRep()".
412  * -# Optionally switch to another algorithm
413  * \ref soplex::SPxSolver::Type "Type"
414  * by calling method \ref soplex::SPxSolver::setType() "setType()".
415  *
416  * Now the solver is ready for execution. If the loaded LP is to be solved
417  * again from scratch, this can be done with method
418  * \ref soplex::SPxSolver::reLoad() "reLoad()". Finally,
419  * \ref soplex::SPxSolver::clear() "clear()" removes the LP from the solver.
420  */
421  //@{
422  /// read LP from input stream.
423  virtual bool read(std::istream& in, NameSet* rowNames = 0,
424  NameSet* colNames = 0, DIdxSet* intVars = 0);
425 
426  /// copy LP.
427  virtual void loadLP(const SPxLP& LP);
428  /// setup linear solver to use. If \p destroy is true, \p slusolver will be freed in destructor.
429  virtual void setSolver(SLinSolver* slu, const bool destroy = false);
430  /// setup pricer to use. If \p destroy is true, \p pricer will be freed in destructor.
431  virtual void setPricer(SPxPricer* pricer, const bool destroy = false);
432  /// setup ratio-tester to use. If \p destroy is true, \p tester will be freed in destructor.
433  virtual void setTester(SPxRatioTester* tester, const bool destroy = false);
434  /// setup starting basis generator to use. If \p destroy is true, \p starter will be freed in destructor.
435  virtual void setStarter(SPxStarter* starter, const bool destroy = false);
436  /// set a start basis.
437  virtual void loadBasis(const SPxBasis::Desc&);
438 
439  /// initialize #ROW or #COLUMN representation.
440  void initRep (Representation p_rep);
441  /// switch to #ROW or #COLUMN representation if not already used.
442  void setRep (Representation p_rep);
443  /// set \ref soplex::SPxSolver::LEAVE "LEAVE" or \ref soplex::SPxSolver::ENTER "ENTER" algorithm.
444  void setType(Type tp);
445  /// set \ref soplex::SPxSolver::FULL "FULL" or \ref soplex::SPxSolver::PARTIAL "PARTIAL" pricing.
446  void setPricing(Pricing pr);
447 
448  /// reload LP.
449  virtual void reLoad();
450 
451  /// clear all data in solver.
452  virtual void clear();
453 
454  /** Load basis from \p filename in MPS format. If \p rowNames and \p
455  * colNames are \c NULL, default names are used for the constraints and
456  * variables.
457  */
458  virtual bool readBasisFile(const char* filename,
459  const NameSet* rowNames, const NameSet* colNames);
460 
461  /** Write basis to \p filename in MPS format. If \p rowNames and \p
462  * colNames are \c NULL, default names are used for the constraints and
463  * variables.
464  */
465  virtual bool writeBasisFile(const char* filename,
466  const NameSet* rowNames, const NameSet* colNames) const;
467 
468  /** Write current LP, basis, and parameter settings.
469  * LP is written in MPS format to "\p filename".mps, basis is written in "\p filename".bas, and parameters
470  * are written to "\p filename".set. If \p rowNames and \p colNames are \c NULL, default names are used for
471  * the constraints and variables.
472  */
473  virtual bool writeState(const char* filename,
474  const NameSet* rowNames = NULL, const NameSet* colNames = NULL) const;
475 
476  //@}
477 
478  /**@name Solving LPs */
479  //@{
480  /// solve loaded LP.
481  /** Solves the loaded LP by calling the floating point routine fpsolve().
482  * If feastol() or opttol() are below irtol(), iterative refinement is applied.
483  */
484  virtual Status solve();
485 
486  /// Status of solution process.
487  Status status() const;
488 
489  /// current objective value.
490  /**@return Objective value of the current solution vector
491  * (see #getPrimal()).
492  */
493  virtual Real value() const;
494 
495 #if 0
496  /// returns dualsol^T b + min{(objvec^T - dualsol^T A) x} calculated in interval arithmetics
497  Real provedBound(Vector& dualsol, const Vector& objvec) const;
498 
499  /// proved dual bound for objective value.
500  virtual Real provedDualbound() const;
501 
502  /// returns whether an infeasible LP is proven to be infeasible.
503  virtual bool isProvenInfeasible() const;
504 #endif
505 
506  /// get solution vector for primal variables.
507  /** This method returns the Status of the basis.
508  * If it is #REGULAR or better,
509  * the primal solution vector of the current basis will be copied
510  * to the argument \p vector. Hence, \p vector must be of dimension
511  * #nCols().
512  *
513  * @throw SPxStatusException if not initialized
514  */
515  virtual Status getPrimal(Vector& vector) const;
516 
517  /// get vector of slack variables.
518  /** This method returns the Status of the basis.
519  * If it is #REGULAR or better,
520  * the slack variables of the current basis will be copied
521  * to the argument \p vector. Hence, \p vector must be of dimension
522  * #nRows().
523  *
524  * @warning Because SPxSolver supports range constraints as its
525  * default, slack variables are defined in a nonstandard way:
526  * Let \em x be the current solution vector and \em A the constraint
527  * matrix. Then the vector of slack variables is defined as
528  * \f$s = Ax\f$.
529  *
530  * @throw SPxStatusException if no problem loaded
531  */
532  virtual Status getSlacks (Vector& vector) const;
533 
534  /// get current solution vector for dual variables.
535  /** This method returns the Status of the basis.
536  * If it is #REGULAR or better,
537  * the vector of dual variables of the current basis will be copied
538  * to the argument \p vector. Hence, \p vector must be of dimension
539  * #nRows().
540  *
541  * @warning Even though mathematically, each range constraint would
542  * account for two dual variables (one for each inequaility), only
543  * #nRows() dual variables are setup via the following
544  * construction: Given a range constraint, there are three possible
545  * situations:
546  * - None of its inequalities is tight: The dual variables
547  * for both are 0. However, when shifting (see below)
548  * occurs, it may be set to a value other than 0, which
549  * models a perturbed objective vector.
550  * - Both of its inequalities are tight: In this case the
551  * range constraint models an equality and we adopt the
552  * standard definition.
553  * - One of its inequalities is tight while the other is not:
554  * In this case only the dual variable for the tight
555  * constraint is given with the standard definition, while
556  * the other constraint is implicitely set to 0.
557  *
558  * @throw SPxStatusException if no problem loaded
559  */
560  virtual Status getDual (Vector& vector) const;
561 
562  /// get vector of reduced costs.
563  /** This method returns the Status of the basis.
564  * If it is #REGULAR or better,
565  * the vector of reduced costs of the current basis will be copied
566  * to the argument \p vector. Hence, \p vector must be of dimension
567  * #nCols().
568  *
569  * Let \em d denote the vector of dual variables, as defined above,
570  * and \em A the LPs constraint matrix. Then the reduced cost vector
571  * \em r is defined as \f$r^T = c^T - d^TA\f$.
572  *
573  * @throw SPxStatusException if no problem loaded
574  */
575  virtual Status getRedCost (Vector& vector) const;
576 
577  /// get primal ray in case of unboundedness.
578  /// @throw SPxStatusException if no problem loaded
579  virtual Status getPrimalray (Vector& vector) const;
580 
581  /// get dual farkas proof of infeasibility.
582  /// @throw SPxStatusException if no problem loaded
583  virtual Status getDualfarkas (Vector& vector) const;
584 
585  /// Termination criterion.
586  /** This method is called in each Simplex iteration to determine, if
587  * the algorithm is to terminate. In this case a nonzero value is
588  * returned.
589  *
590  * This method is declared virtual to allow for implementation of
591  * other stopping criteria or using it as callback method within the
592  * Simplex loop, by overriding the method in a derived class.
593  * However, all implementations must terminate with the
594  * statement \c return SPxSolver::#terminate(), if no own termination
595  * criteria is encountered.
596  *
597  * Note, that the Simplex loop stopped even when #terminate()
598  * returns 0, if the LP has been solved to optimality (i.e. no
599  * further pricing succeeds and no shift is present).
600  */
601  virtual bool terminate ();
602  //@}
603 
604  //-----------------------------
605  /**@name Control Parameters */
606  //@{
607  /// values \f$|x| < \epsilon\f$ are considered to be 0.
608  /** if you want another value for epsilon, use
609  * \ref soplex::Param::setEpsilon() "Param::setEpsilon()".
610  */
611  Real epsilon() const
612  {
613  return primVec.delta().getEpsilon();
614  }
615  /// feasibility tolerance maintained by ratio test during ENTER algorithm.
616  Real entertol() const
617  {
618  assert(m_entertol > 0.0);
619 
620  return m_entertol;
621  }
622  /// feasibility tolerance maintained by ratio test during LEAVE algorithm.
623  Real leavetol() const
624  {
625  assert(m_leavetol > 0.0);
626 
627  return m_leavetol;
628  }
629  /// allowed primal feasibility tolerance.
630  Real feastol() const
631  {
632  assert(m_entertol >= 0.0);
633  assert(m_leavetol >= 0.0);
634 
635  return theRep == COLUMN ? m_entertol : m_leavetol;
636  }
637  /// allowed optimality, i.e., dual feasibility tolerance.
638  Real opttol() const
639  {
640  assert(m_entertol >= 0.0);
641  assert(m_leavetol >= 0.0);
642 
643  return theRep == COLUMN ? m_leavetol : m_entertol;
644  }
645  /// guaranteed primal and dual bound violation for optimal solution, returning the maximum of feastol() and opttol(), i.e., the less tight tolerance.
646  Real delta() const
647  {
648  assert(m_entertol >= 0.0);
649  assert(m_leavetol >= 0.0);
650 
651  return m_entertol > m_leavetol ? m_entertol : m_leavetol;
652  }
653  /// iterative refinement threshold: if feastol() or opttol() are below this value, iterative refinement is applied.
654  Real irthreshold() const
655  {
656  assert(m_irthreshold > 0.0);
657 
658  return m_irthreshold;
659  }
660  /// set parameter \p feastol.
661  void setFeastol(Real d);
662  /// set parameter \p opttol.
663  void setOpttol(Real d);
664  /// set parameter \p delta, i.e., set \p feastol and \p opttol to same value.
665  void setDelta(Real d);
666  /// set parameter \p irthreshold.
667  void setIrthreshold(Real d);
668 
669  /** SPxSolver considers a Simplex step as degenerate if the
670  * steplength does not exceed #epsilon(). Cycling occurs if only
671  * degenerate steps are taken. To prevent this situation, SPxSolver
672  * perturbs the problem such that nondegenerate steps are ensured.
673  *
674  * maxCycle() controls how agressive such perturbation is
675  * performed, since no more than maxCycle() degenerate steps are
676  * accepted before perturbing the LP. The current number of consecutive
677  * degenerate steps is counted by numCycle().
678  */
679  /// maximum number of degenerate simplex steps before we detect cycling.
680  int maxCycle() const
681  {
682  return m_maxCycle;
683  }
684  /// actual number of degenerate simplex steps encountered so far.
685  int numCycle() const
686  {
687  return m_numCycle;
688  }
689  //@}
690 
691 private:
692 
693  //-----------------------------
694  /**@name Private helpers */
695  //@{
696  /// solve loaded LP using floating point arithmetic.
697  /** Solves the loaded LP by processing the Simplex iteration until
698  * the termination criteria is fullfilled (see #terminate()).
699  * The SPxStatus of the solver will indicate the reason for termination.
700  *
701  * @throw SPxStatusException if either no problem, solver, pricer
702  * or ratiotester loaded or if solve is still running when it shouldn't be
703  */
704  Status fpsolve();
705  ///
706  void localAddRows(int start);
707  ///
708  void localAddCols(int start);
709  /// apply iterative refinement until irfeastol and iropttol are reached or modified problem is not solved to
710  /// optimality; returns true if and only if precision has been reached
711  bool refine(
712  Real irfeastol, /**< primal feasibility tolerance */
713  Real iropttol, /**< dual feasibility tolerance */
714  Vector_exact& primal_ex, /**< buffer to return refined primal solution values */
715  Vector_exact& slack_ex, /**< buffer to return refined slack values */
716  Vector_exact& dual_ex, /**< buffer to return refined dual solution values */
717  Vector_exact& redcost_ex, /**< buffer to return refined reduced cost values */
718  int maxitersround /**< iteration limit per refinement round */
719  );
720  ///
721  void setPrimal(Vector& p_vector);
722  ///
723  void setSlacks(Vector& p_vector);
724  ///
725  void setDual(Vector& p_vector);
726  ///
727  void setRedCost(Vector& p_vector);
728  //@}
729 
730 protected:
731 
732  //-----------------------------
733  /**@name Protected helpers */
734  //@{
735  ///
736  virtual void addedRows(int n);
737  ///
738  virtual void addedCols(int n);
739  ///
740  virtual void doRemoveRow(int i);
741  ///
742  virtual void doRemoveRows(int perm[]);
743  ///
744  virtual void doRemoveCol(int i);
745  ///
746  virtual void doRemoveCols(int perm[]);
747  //@}
748 
749 public:
750 
751  //-----------------------------
752  /**@name Modification */
753  //@{
754  ///
755  virtual void changeObj(const Vector& newObj);
756  ///
757  virtual void changeObj(int i, Real newVal);
758  ///
759  virtual void changeObj(SPxColId p_id, Real p_newVal)
760  {
761  changeObj(number(p_id), p_newVal);
762  }
763  ///
764  virtual void changeLower(const Vector& newLower);
765  ///
766  virtual void changeLower(int i, Real newLower);
767  ///
768  virtual void changeLower(SPxColId p_id, Real p_newLower)
769  {
770  changeLower(number(p_id), p_newLower);
771  }
772  ///
773  virtual void changeUpper(const Vector& newUpper);
774  ///
775  virtual void changeUpper(int i, Real newUpper);
776  ///
777  virtual void changeUpper(SPxColId p_id, Real p_newUpper)
778  {
779  changeUpper(number(p_id), p_newUpper);
780  }
781  ///
782  virtual void changeBounds(const Vector& newLower, const Vector& newUpper);
783  ///
784  virtual void changeBounds(int i, Real newLower, Real newUpper);
785  ///
786  virtual void changeBounds(
787  SPxColId p_id, Real p_newLower, Real p_newUpper)
788  {
789  changeBounds(number(p_id), p_newLower, p_newUpper);
790  }
791  ///
792  virtual void changeLhs(const Vector& newLhs);
793  ///
794  virtual void changeLhs(int i, Real newLhs);
795  ///
796  virtual void changeLhs(SPxRowId p_id, Real p_newLhs)
797  {
798  changeLhs(number(p_id), p_newLhs);
799  }
800  ///
801  virtual void changeRhs(const Vector& newRhs);
802  ///
803  virtual void changeRhs(int i, Real newRhs);
804  ///
805  virtual void changeRhs(SPxRowId p_id, Real p_newRhs)
806  {
807  changeRhs(number(p_id), p_newRhs);
808  }
809  ///
810  virtual void changeRange(const Vector& newLhs, const Vector& newRhs);
811  ///
812  virtual void changeRange(int i, Real newLhs, Real newRhs);
813  ///
814  virtual void changeRange(
815  SPxRowId p_id, Real p_newLhs, Real p_newRhs)
816  {
817  changeRange(number(p_id), p_newLhs, p_newRhs);
818  }
819  ///
820  virtual void changeRow(int i, const LPRow& newRow);
821  ///
822  virtual void changeRow(SPxRowId p_id, const LPRow& p_newRow)
823  {
824  changeRow(number(p_id), p_newRow);
825  }
826  ///
827  virtual void changeCol(int i, const LPCol& newCol);
828  ///
829  virtual void changeCol(SPxColId p_id, const LPCol& p_newCol)
830  {
831  changeCol(number(p_id), p_newCol);
832  }
833  ///
834  virtual void changeElement(int i, int j, Real val);
835  ///
836  virtual void changeElement(
837  SPxRowId rid, SPxColId cid, Real val)
838  {
839  changeElement(number(rid), number(cid), val);
840  }
841  ///
842  virtual void changeSense(SPxSense sns);
843  //@}
844 
845  //------------------------------------
846  /**@name Dimension and codimension */
847  //@{
848  /// dimension of basis matrix.
849  int dim() const
850  {
851  return thecovectors->num();
852  }
853  /// codimension.
854  int coDim() const
855  {
856  return thevectors->num();
857  }
858  //@}
859 
860  //------------------------------------
861  /**@name Variables and Covariables
862  * Class SPxLP introduces \ref soplex::SPxId "SPxIds" to identify
863  * row or column data of an LP. SPxSolver uses this concept to
864  * access data with respect to the chosen representation.
865  */
866  //@{
867  /// id of \p i 'th vector.
868  /** The \p i 'th Id is the \p i 'th SPxRowId for a rowwise and the
869  * \p i 'th SPxColId for a columnwise basis represenation. Hence,
870  * 0 <= i < #coDim().
871  */
872  SPxId id(int i) const
873  {
874  if (rep() == ROW)
875  {
876  SPxRowId rid = SPxLP::rId(i);
877  return SPxId(rid);
878  }
879  else
880  {
881  SPxColId cid = SPxLP::cId(i);
882  return SPxId(cid);
883  }
884  }
885 
886  /// id of \p i 'th covector.
887  /** The \p i 'th #coId() is the \p i 'th SPxColId for a rowwise and the
888  * \p i 'th SPxRowId for a columnwise basis represenation. Hence,
889  * 0 <= i < #dim().
890  */
891  SPxId coId(int i) const
892  {
893  if (rep() == ROW)
894  {
895  SPxColId cid = SPxLP::cId(i);
896  return SPxId(cid);
897  }
898  else
899  {
900  SPxRowId rid = SPxLP::rId(i);
901  return SPxId(rid);
902  }
903  }
904 
905  /// Is \p p_id an SPxId ?
906  /** This method returns wheather or not \p p_id identifies a vector
907  * with respect to the chosen representation.
908  */
909  int isId(const SPxId& p_id) const
910  {
911  return p_id.info * theRep > 0;
912  }
913 
914  /// Is \p p_id a CoId.
915  /** This method returns wheather or not \p p_id identifies a coVector
916  * with respect to the chosen representation.
917  */
918  int isCoId(const SPxId& p_id) const
919  {
920  return p_id.info * theRep < 0;
921  }
922  //@}
923 
924  //------------------------------------
925  /**@name Vectors and Covectors */
926  //@{
927  /// \p i 'th vector.
928  /**@return a reference to the \p i 'th, 0 <= i < #coDim(), vector of
929  * the loaded LP (with respect to the chosen representation).
930  */
931  const SVector& vector(int i) const
932  {
933  return (*thevectors)[i];
934  }
935 
936  ///
937  const SVector& vector(const SPxRowId& rid) const
938  {
939  assert(rid.isValid());
940  return (rep() == ROW)
941  ? (*thevectors)[number(rid)]
942  : static_cast<const SVector&>(unitVecs[number(rid)]);
943  }
944  ///
945  const SVector& vector(const SPxColId& cid) const
946  {
947  assert(cid.isValid());
948  return (rep() == COLUMN)
949  ? (*thevectors)[number(cid)]
950  : static_cast<const SVector&>(unitVecs[number(cid)]);
951  }
952 
953  /// vector associated to \p p_id.
954  /**@return Returns a reference to the vector of the loaded LP corresponding
955  * to \p id (with respect to the chosen representation). If \p p_id is
956  * an id, a vector of the constraint matrix is returned, otherwise
957  * the corresponding unit vector (of the slack variable or bound
958  * inequality) is returned.
959  * @todo The implementation does not exactly look like it will do
960  * what is promised in the describtion.
961  */
962  const SVector& vector(const SPxId& p_id) const
963  {
964  assert(p_id.isValid());
965 
966  return p_id.isSPxRowId()
967  ? vector(SPxRowId(p_id))
968  : vector(SPxColId(p_id));
969  }
970 
971  /// \p i 'th covector of LP.
972  /**@return a reference to the \p i 'th, 0 <= i < #dim(), covector of
973  * the loaded LP (with respect to the chosen representation).
974  */
975  const SVector& coVector(int i) const
976  {
977  return (*thecovectors)[i];
978  }
979  ///
980  const SVector& coVector(const SPxRowId& rid) const
981  {
982  assert(rid.isValid());
983  return (rep() == COLUMN)
984  ? (*thecovectors)[number(rid)]
985  : static_cast<const SVector&>(unitVecs[number(rid)]);
986  }
987  ///
988  const SVector& coVector(const SPxColId& cid) const
989  {
990  assert(cid.isValid());
991  return (rep() == ROW)
992  ? (*thecovectors)[number(cid)]
993  : static_cast<const SVector&>(unitVecs[number(cid)]);
994  }
995  /// coVector associated to \p p_id.
996  /**@return a reference to the covector of the loaded LP
997  * corresponding to \p p_id (with respect to the chosen
998  * representation). If \p p_id is a coid, a covector of the constraint
999  * matrix is returned, otherwise the corresponding unit vector is
1000  * returned.
1001  */
1002  const SVector& coVector(const SPxId& p_id) const
1003  {
1004  assert(p_id.isValid());
1005  return p_id.isSPxRowId()
1006  ? coVector(SPxRowId(p_id))
1007  : coVector(SPxColId(p_id));
1008  }
1009  /// return \p i 'th unit vector.
1010  const SVector& unitVector(int i) const
1011  {
1012  return unitVecs[i];
1013  }
1014  //@}
1015 
1016  //------------------------------------
1017  /**@name Variable status
1018  * The Simplex basis assigns a \ref soplex::SPxBasis::Desc::Status
1019  * "Status" to each variable and covariable. Depending on the
1020  * representation, the status indicates that the corresponding
1021  * vector is in the basis matrix or not.
1022  */
1023  //@{
1024  /// Status of \p i 'th variable.
1025  SPxBasis::Desc::Status varStatus(int i) const
1026  {
1027  return desc().status(i);
1028  }
1029 
1030  /// Status of \p i 'th covariable.
1031  SPxBasis::Desc::Status covarStatus(int i) const
1032  {
1033  return desc().coStatus(i);
1034  }
1035 
1036  /// does \p stat describe a basic index ?
1037  int isBasic(SPxBasis::Desc::Status stat) const
1038  {
1039  return (stat * rep() > 0);
1040  }
1041 
1042  /// is the \p p_id 'th vector basic ?
1043  int isBasic(const SPxId& p_id) const
1044  {
1045  assert(p_id.isValid());
1046  return p_id.isSPxRowId()
1047  ? isBasic(SPxRowId(p_id))
1048  : isBasic(SPxColId(p_id));
1049  }
1050 
1051  /// is the \p rid 'th vector basic ?
1052  int isBasic(const SPxRowId& rid) const
1053  {
1054  return isBasic(desc().rowStatus(number(rid)));
1055  }
1056 
1057  /// is the \p cid 'th vector basic ?
1058  int isBasic(const SPxColId& cid) const
1059  {
1060  return isBasic(desc().colStatus(number(cid)));
1061  }
1062 
1063  /// is the \p i 'th row vector basic ?
1064  int isRowBasic(int i) const
1065  {
1066  return isBasic(desc().rowStatus(i));
1067  }
1068 
1069  /// is the \p i 'th column vector basic ?
1070  int isColBasic(int i) const
1071  {
1072  return isBasic(desc().colStatus(i));
1073  }
1074 
1075  /// is the \p i 'th vector basic ?
1076  int isBasic(int i) const
1077  {
1078  return isBasic(desc().status(i));
1079  }
1080 
1081  /// is the \p i 'th covector basic ?
1082  int isCoBasic(int i) const
1083  {
1084  return isBasic(desc().coStatus(i));
1085  }
1086  //@}
1087 
1088  /// feasibility vector.
1089  /** This method return the \em feasibility vector. If it satisfies its
1090  * bound, the basis is called feasible (independently of the chosen
1091  * representation). The feasibility vector has dimension #dim().
1092  *
1093  * For the entering Simplex, #fVec is kept within its bounds. In
1094  * contrast to this, the pricing of the leaving Simplex selects an
1095  * element of #fVec, that violates its bounds.
1096  */
1097  UpdateVector& fVec() const
1098  {
1099  return *theFvec;
1100  }
1101  /// right-hand side vector for \ref soplex::SPxSolver::fVec "fVec"
1102  /** The feasibility vector is computed by solving a linear system with the
1103  * basis matrix. The right-hand side vector of this system is referred
1104  * to as \em feasibility, \em right-hand \em side \em vector #fRhs().
1105  *
1106  * For a row basis, #fRhs() is the objective vector (ignoring shifts).
1107  * For a column basis, it is the sum of all nonbasic vectors scaled by
1108  * the factor of their bound.
1109  */
1110  const Vector& fRhs() const
1111  {
1112  return *theFrhs;
1113  }
1114  /// upper bound for \ref soplex::SPxSolver::fVec "fVec".
1115  const Vector& ubBound() const
1116  {
1117  return theUBbound;
1118  }
1119  /// upper bound for #fVec, writable.
1120  /** This method returns the upper bound for the feasibility vector.
1121  * It may only be called for the #ENTER%ing Simplex.
1122  *
1123  * For the #ENTER%ing Simplex algorithms, the feasibility vector is
1124  * maintained to fullfill its bounds. As #fVec itself, also its
1125  * bounds depend on the chosen representation. Further, they may
1126  * need to be shifted (see below).
1127  */
1128  Vector& ubBound()
1129  {
1130  return theUBbound;
1131  }
1132  /// lower bound for \ref soplex::SPxSolver::fVec "fVec".
1133  const Vector& lbBound() const
1134  {
1135  return theLBbound;
1136  }
1137  /// lower bound for #fVec, writable.
1138  /** This method returns the lower bound for the feasibility vector.
1139  * It may only be called for the #ENTER%ing Simplex.
1140  *
1141  * For the #ENTER%ing Simplex algorithms, the feasibility vector is
1142  * maintained to fullfill its bounds. As #fVec itself, also its
1143  * bound depend on the chosen representation. Further, they may
1144  * need to be shifted (see below).
1145  */
1146  Vector& lbBound()
1147  {
1148  return theLBbound;
1149  }
1150 
1151  /// Violations of \ref soplex::SPxSolver::fVec "fVec"
1152  /** For the leaving Simplex algorithm, pricing involves selecting a
1153  * variable from #fVec that violates its bounds that is to leave
1154  * the basis. When a SPxPricer is called to select such a
1155  * leaving variable, #fTest() contains the vector of violations:
1156  * For #fTest()[i] < 0, the \c i 'th basic variable violates one of
1157  * its bounds by the given value. Otherwise no bound is violated.
1158  */
1159  const Vector& fTest() const
1160  {
1161  assert(type() == LEAVE);
1162  return theCoTest;
1163  }
1164 
1165  /// copricing vector.
1166  /** The copricing vector #coPvec along with the pricing vector
1167  * #pVec are used for pricing in the #ENTER%ing Simplex algorithm,
1168  * i.e. one variable is selected, that violates its bounds. In
1169  * contrast to this, the #LEAVE%ing Simplex algorithm keeps both
1170  * vectors within their bounds.
1171  */
1172  UpdateVector& coPvec() const
1173  {
1174  return *theCoPvec;
1175  }
1176 
1177  /// Right-hand side vector for \ref soplex::SPxSolver::coPvec "coPvec".
1178  /** The vector #coPvec is computed by solving a linear system with the
1179  * basis matrix and #coPrhs as the right-hand side vector. For
1180  * column basis representation, #coPrhs is build up of the
1181  * objective vector elements of all basic variables. For a row
1182  * basis, it consists of the thight bounds of all basic
1183  * constraints.
1184  */
1185  const Vector& coPrhs() const
1186  {
1187  return *theCoPrhs;
1188  }
1189 
1190  ///
1191  const Vector& ucBound() const
1192  {
1193  assert(theType == LEAVE);
1194  return *theCoUbound;
1195  }
1196  /// upper bound for #coPvec.
1197  /** This method returns the upper bound for #coPvec. It may only be
1198  * called for the leaving Simplex algorithm.
1199  *
1200  * For the leaving Simplex algorithms #coPvec is maintained to
1201  * fullfill its bounds. As #coPvec itself, also its bound depend
1202  * on the chosen representation. Further, they may need to be
1203  * shifted (see below).
1204  */
1205  Vector& ucBound()
1206  {
1207  assert(theType == LEAVE);
1208  return *theCoUbound;
1209  }
1210 
1211  ///
1212  const Vector& lcBound() const
1213  {
1214  assert(theType == LEAVE);
1215  return *theCoLbound;
1216  }
1217  /// lower bound for #coPvec.
1218  /** This method returns the lower bound for #coPvec. It may only be
1219  * called for the leaving Simplex algorithm.
1220  *
1221  * For the leaving Simplex algorithms #coPvec is maintained to
1222  * fullfill its bounds. As #coPvec itself, also its bound depend
1223  * on the chosen representation. Further, they may need to be
1224  * shifted (see below).
1225  */
1226  Vector& lcBound()
1227  {
1228  assert(theType == LEAVE);
1229  return *theCoLbound;
1230  }
1231 
1232  /// violations of \ref soplex::SPxSolver::coPvec "coPvec".
1233  /** In entering Simplex pricing selects checks vectors #coPvec()
1234  * and #pVec() for violation of its bounds. #coTest() contains
1235  * the violations for #coPvec() which are indicated by a negative
1236  * value. That is, if #coTest()[i] < 0, the \p i 'th element of #coPvec()
1237  * is violated by -#coTest()[i].
1238  */
1239  const Vector& coTest() const
1240  {
1241  assert(type() == ENTER);
1242  return theCoTest;
1243  }
1244  /// pricing vector.
1245  /** The pricing vector #pVec is the product of #coPvec with the
1246  * constraint matrix. As #coPvec, also #pVec is maintained within
1247  * its bound for the leaving Simplex algorithm, while the bounds
1248  * are tested for the entering Simplex. #pVec is of dimension
1249  * #coDim(). Vector #pVec() is only up to date for #LEAVE%ing
1250  * Simplex or #FULL pricing in #ENTER%ing Simplex.
1251  */
1252  UpdateVector& pVec() const
1253  {
1254  return *thePvec;
1255  }
1256  ///
1257  const Vector& upBound() const
1258  {
1259  assert(theType == LEAVE);
1260  return *theUbound;
1261  }
1262  /// upper bound for #pVec.
1263  /** This method returns the upper bound for #pVec. It may only be
1264  * called for the leaving Simplex algorithm.
1265  *
1266  * For the leaving Simplex algorithms #pVec is maintained to
1267  * fullfill its bounds. As #pVec itself, also its bound depend
1268  * on the chosen representation. Further, they may need to be
1269  * shifted (see below).
1270  */
1271  Vector& upBound()
1272  {
1273  assert(theType == LEAVE);
1274  return *theUbound;
1275  }
1276 
1277  ///
1278  const Vector& lpBound() const
1279  {
1280  assert(theType == LEAVE);
1281  return *theLbound;
1282  }
1283  /// lower bound for #pVec.
1284  /** This method returns the lower bound for #pVec. It may only be
1285  * called for the leaving Simplex algorithm.
1286  *
1287  * For the leaving Simplex algorithms #pVec is maintained to
1288  * fullfill its bounds. As #pVec itself, also its bound depend
1289  * on the chosen representation. Further, they may need to be
1290  * shifted (see below).
1291  */
1292  Vector& lpBound()
1293  {
1294  assert(theType == LEAVE);
1295  return *theLbound;
1296  }
1297 
1298  /// Violations of \ref soplex::SPxSolver::pVec "pVec".
1299  /** In entering Simplex pricing selects checks vectors #coPvec()
1300  * and #pVec() for violation of its bounds. Vector #test()
1301  * contains the violations for #pVec(), i.e., if #test()[i] < 0,
1302  * the i'th element of #pVec() is violated by #test()[i].
1303  * Vector #test() is only up to date for #FULL pricing.
1304  */
1305  const Vector& test() const
1306  {
1307  assert(type() == ENTER);
1308  return theTest;
1309  }
1310 
1311  /// compute and return \ref soplex::SPxSolver::pVec() "pVec()"[i].
1312  Real computePvec(int i);
1313  /// compute entire \ref soplex::SPxSolver::pVec() "pVec()".
1314  void computePvec();
1315  /// compute and return \ref soplex::SPxSolver::test() "test()"[i] in \ref soplex::SPxSolver::ENTER "ENTER"ing Simplex.
1316  Real computeTest(int i);
1317  /// compute test vector in \ref soplex::SPxSolver::ENTER "ENTER"ing Simplex.
1318  void computeTest();
1319 
1320  //------------------------------------
1321  /**@name Shifting
1322  * The task of the ratio test (implemented in SPxRatioTester classes)
1323  * is to select a variable for the basis update, such that the basis
1324  * remains priced (i.e. both, the pricing and copricing vectors satisfy
1325  * their bounds) or feasible (i.e. the feasibility vector satisfies its
1326  * bounds). However, this can lead to numerically instable basis matrices
1327  * or -- after accumulation of various errors -- even to a singular basis
1328  * matrix.
1329  *
1330  * The key to overcome this problem is to allow the basis to become "a
1331  * bit" infeasible or unpriced, in order provide a better choice for the
1332  * ratio test to select a stable variable. This is equivalent to enlarging
1333  * the bounds by a small amount. This is referred to as \em shifting.
1334  *
1335  * These methods serve for shifting feasibility bounds, either in order
1336  * to maintain numerical stability or initially for computation of
1337  * phase 1. The sum of all shifts applied to any bound is stored in
1338  * \ref soplex::SPxSolver::theShift "theShift".
1339  *
1340  * The following methods are used to shift individual bounds. They are
1341  * mainly intended for stable implenentations of SPxRatioTester.
1342  */
1343  //@{
1344  /// Perform initial shifting to optain an feasible or pricable basis.
1345  void shiftFvec();
1346  /// Perform initial shifting to optain an feasible or pricable basis.
1347  void shiftPvec();
1348 
1349  /// shift \p i 'th \ref soplex::SPxSolver::ubBound "ubBound" to \p to.
1350  void shiftUBbound(int i, Real to)
1351  {
1352  assert(theType == ENTER);
1353  theShift += to - theUBbound[i];
1354  theUBbound[i] = to;
1355  }
1356  /// shift \p i 'th \ref soplex::SPxSolver::lbBound "lbBound" to \p to.
1357  void shiftLBbound(int i, Real to)
1358  {
1359  assert(theType == ENTER);
1360  theShift += theLBbound[i] - to;
1361  theLBbound[i] = to;
1362  }
1363  /// shift \p i 'th \ref soplex::SPxSolver::upBound "upBound" to \p to.
1364  void shiftUPbound(int i, Real to)
1365  {
1366  assert(theType == LEAVE);
1367  theShift += to - (*theUbound)[i];
1368  (*theUbound)[i] = to;
1369  }
1370  /// shift \p i 'th \ref soplex::SPxSolver::lpBound "lpBound" to \p to.
1371  void shiftLPbound(int i, Real to)
1372  {
1373  assert(theType == LEAVE);
1374  theShift += (*theLbound)[i] - to;
1375  (*theLbound)[i] = to;
1376  }
1377  /// shift \p i 'th \ref soplex::SPxSolver::ucBound "ucBound" to \p to.
1378  void shiftUCbound(int i, Real to)
1379  {
1380  assert(theType == LEAVE);
1381  theShift += to - (*theCoUbound)[i];
1382  (*theCoUbound)[i] = to;
1383  }
1384  /// shift \p i 'th \ref soplex::SPxSolver::lcBound "lcBound" to \p to.
1385  void shiftLCbound(int i, Real to)
1386  {
1387  assert(theType == LEAVE);
1388  theShift += (*theCoLbound)[i] - to;
1389  (*theCoLbound)[i] = to;
1390  }
1391  ///
1392  void testBounds() const;
1393 
1394  /// total current shift amount.
1395  virtual Real shift() const
1396  {
1397  return theShift;
1398  }
1399  /// remove shift as much as possible.
1400  virtual void unShift(void);
1401 
1402  /// get violation of constraints.
1403  virtual void qualConstraintViolation(Real& maxviol, Real& sumviol) const;
1404  /// get violations of bounds.
1405  virtual void qualBoundViolation(Real& maxviol, Real& sumviol) const;
1406  /// get the residuum |Ax-b|.
1407  virtual void qualSlackViolation(Real& maxviol, Real& sumviol) const;
1408  /// get violation of optimality criterion.
1409  virtual void qualRedCostViolation(Real& maxviol, Real& sumviol) const;
1410  //@}
1411 
1412 private:
1413 
1414  //------------------------------------
1415  /**@name Perturbation */
1416  //@{
1417  ///
1418  void perturbMin(
1419  const UpdateVector& vec, Vector& low, Vector& up, Real eps, Real delta,
1420  int start = 0, int incr = 1);
1421  ///
1422  void perturbMax(
1423  const UpdateVector& vec, Vector& low, Vector& up, Real eps, Real delta,
1424  int start = 0, int incr = 1);
1425  ///
1426  Real perturbMin(const UpdateVector& uvec,
1427  Vector& low, Vector& up, Real eps, Real delta,
1428  const SPxBasis::Desc::Status* stat, int start, int incr) const;
1429  ///
1430  Real perturbMax(const UpdateVector& uvec,
1431  Vector& low, Vector& up, Real eps, Real delta,
1432  const SPxBasis::Desc::Status* stat, int start, int incr) const;
1433  //@}
1434 
1435  //------------------------------------
1436  /**@name The Simplex Loop
1437  * We now present a set of methods that may be usefull when implementing
1438  * own SPxPricer or SPxRatioTester classes. Here is, how
1439  * SPxSolver will call methods from its loaded SPxPricer and
1440  * SPxRatioTester.
1441  *
1442  * For the entering Simplex:
1443  * -# \ref soplex::SPxPricer::selectEnter() "SPxPricer::selectEnter()"
1444  * -# \ref soplex::SPxRatioTester::selectLeave() "SPxRatioTester::selectLeave()"
1445  * -# \ref soplex::SPxPricer::entered4() "SPxPricer::entered4()"
1446  *
1447  * For the leaving Simplex:
1448  * -# \ref soplex::SPxPricer::selectLeave() "SPxPricer::selectLeave()"
1449  * -# \ref soplex::SPxRatioTester::selectEnter() "SPxRatioTester::selectEnter()"
1450  * -# \ref soplex::SPxPricer::left4() "SPxPricer::left4()"
1451  */
1452  //@{
1453 public:
1454  /// Setup vectors to be solved within Simplex loop.
1455  /** Load vector \p y to be #solve%d with the basis matrix during the
1456  * #LEAVE Simplex. The system will be solved after #SPxSolver%'s call
1457  * to SPxRatioTester. The system will be solved along with
1458  * another system. Solving two linear system at a time has
1459  * performance advantages over solving the two linear systems
1460  * seperately.
1461  */
1462  void setup4solve(Vector* p_y, SSVector* p_rhs)
1463  {
1464  assert(type() == LEAVE);
1465  solveVector2 = p_y;
1466  solveVector2rhs = p_rhs;
1467  }
1468  /// Setup vectors to be solved within Simplex loop.
1469  /** Load a second additional vector \p y2 to be #solve%d with the
1470  * basis matrix during the #LEAVE Simplex. The system will be
1471  * solved after #SPxSolver%'s call to SPxRatioTester.
1472  * The system will be solved along with at least one
1473  * other system. Solving several linear system at a time has
1474  * performance advantages over solving them seperately.
1475  */
1476  void setup4solve2(Vector* p_y2, SSVector* p_rhs2)
1477  {
1478  assert(type() == LEAVE);
1479  solveVector3 = p_y2;
1480  solveVector3rhs = p_rhs2;
1481  }
1482  /// Setup vectors to be cosolved within Simplex loop.
1483  /** Load vector \p y to be #coSolve%d with the basis matrix during
1484  * the #ENTER Simplex. The system will be solved after #SPxSolver%'s
1485  * call to SPxRatioTester. The system will be solved along
1486  * with another system. Solving two linear system at a time has
1487  * performance advantages over solving the two linear systems
1488  * seperately.
1489  */
1490  void setup4coSolve(Vector* p_y, SSVector* p_rhs)
1491  {
1492  assert(type() == ENTER);
1493  coSolveVector2 = p_y;
1494  coSolveVector2rhs = p_rhs;
1495  }
1496  /// maximal infeasibility of basis
1497  /** This method is called before concluding optimality. Since it is
1498  * possible that some stable implementation of class
1499  * SPxRatioTester yielded a slightly infeasible (or unpriced)
1500  * basis, this must be checked before terminating with an optimal
1501  * solution.
1502  */
1503  virtual Real maxInfeas() const;
1504 
1505  /// Return current basis.
1506  /**@note The basis can be used to solve linear systems or use
1507  * any other of its (const) methods. It is, however, encuraged
1508  * to use methods #setup4solve() and #setup4coSolve() for solving
1509  * systems, since this is likely to have perfomance advantages.
1510  */
1511  const SPxBasis& basis() const
1512  {
1513  return *this;
1514  }
1515  ///
1516  SPxBasis& basis()
1517  {
1518  return *this;
1519  }
1520  /// return loaded SPxPricer.
1521  const SPxPricer* pricer() const
1522  {
1523  return thepricer;
1524  }
1525  /// return loaded SLinSolver.
1526  const SLinSolver* slinSolver() const
1527  {
1528  return SPxBasis::factor;
1529  }
1530  /// return loaded SPxRatioTester.
1531  const SPxRatioTester* ratiotester() const
1532  {
1533  return theratiotester;
1534  }
1535 
1536 protected:
1537  /// Factorize basis matrix.
1538  /// @throw SPxStatusException if loaded matrix is singular
1539  virtual void factorize();
1540 
1541 private:
1542 
1543  /** let index \p i leave the basis and manage entering of another one.
1544  @returns \c false if LP is unbounded/infeasible. */
1545  bool leave(int i);
1546  /** let id enter the basis and manage leaving of another one.
1547  @returns \c false if LP is unbounded/infeasible. */
1548  bool enter(SPxId& id);
1549 
1550  /// test coVector \p i with status \p stat.
1551  Real coTest(int i, SPxBasis::Desc::Status stat) const;
1552  /// compute coTest vector.
1553  void computeCoTest();
1554  /// recompute coTest vector.
1555  void updateCoTest();
1556 
1557  /// test vector \p i with status \p stat.
1558  Real test(int i, SPxBasis::Desc::Status stat) const;
1559  /// recompute test vector.
1560  void updateTest();
1561 
1562  /// compute basis feasibility test vector.
1563  void computeFtest();
1564  /// update basis feasibility test vector.
1565  void updateFtest();
1566  //@}
1567 
1568  //------------------------------------
1569  /**@name Parallelization
1570  * In this section we present the methods, that are provided in order to
1571  * allow a parallel version to be implemented as a derived class, thereby
1572  * inheriting most of the code of SPxSolver.
1573  *
1574  * @par Initialization
1575  * These methods are used to setup all the vectors used in the Simplex
1576  * loop, that where described in the previous sectios.
1577  */
1578  //@{
1579 public:
1580  /// intialize data structures.
1581  /** If SPxSolver is not \ref isInitialized() "initialized", the method
1582  * #solve() calls #init() to setup all vectors and internal data structures.
1583  * Most of the other methods within this section are called by #init().
1584  *
1585  * Derived classes should add the initialization of additional
1586  * data structures by overriding this method. Don't forget,
1587  * however, to call SPxSolver::init().
1588  */
1589  virtual void init();
1590 
1591 protected:
1592 
1593  /// has the internal data been initialized?
1594  /** As long as an instance of SPxSolver is not initialized, no member
1595  * contains setup data. Initialization is performed via method
1596  * #init(). Afterwards all data structures are kept up to date (even
1597  * for all manipulation methods), until #unInit() is called. However,
1598  * some manipulation methods call #unInit() themselfs.
1599  */
1600  bool isInitialized() const
1601  {
1602  return initialized;
1603  }
1604  /// unintialize data structures.
1605  virtual void unInit()
1606  {
1607  initialized = false;
1608  }
1609  /// setup all vecs fresh
1610  virtual void reinitializeVecs();
1611  /// reset dimensions of vectors according to loaded LP.
1612  virtual void reDim();
1613  /// compute feasibility vector from scratch.
1614  void computeFrhs();
1615  ///
1616  virtual void computeFrhsXtra();
1617  ///
1618  virtual void computeFrhs1(const Vector&, const Vector&);
1619  ///
1620  void computeFrhs2(const Vector&, const Vector&);
1621  /// compute \ref soplex::SPxSolver::theCoPrhs "theCoPrhs" for entering Simplex.
1622  virtual void computeEnterCoPrhs();
1623  ///
1624  void computeEnterCoPrhs4Row(int i, int n);
1625  ///
1626  void computeEnterCoPrhs4Col(int i, int n);
1627  /// compute \ref soplex::SPxSolver::theCoPrhs "theCoPrhs" for leaving Simplex.
1628  virtual void computeLeaveCoPrhs();
1629  ///
1630  void computeLeaveCoPrhs4Row(int i, int n);
1631  ///
1632  void computeLeaveCoPrhs4Col(int i, int n);
1633 
1634  /// Compute part of objective value.
1635  /** This method is called from #value() in order to compute the part of
1636  * the objective value resulting form nonbasic variables for #COLUMN
1637  * Representation.
1638  */
1639  Real nonbasicValue() const;
1640 
1641  /// Get pointer to the \p id 'th vector
1642  virtual const SVector* enterVector(const SPxId& p_id)
1643  {
1644  assert(p_id.isValid());
1645  return p_id.isSPxRowId()
1646  ? &vector(SPxRowId(p_id)) : &vector(SPxColId(p_id));
1647  }
1648  ///
1649  virtual void getLeaveVals(int i,
1650  SPxBasis::Desc::Status& leaveStat, SPxId& leaveId,
1651  Real& leaveMax, Real& leavebound, int& leaveNum);
1652  ///
1653  virtual void getLeaveVals2(Real leaveMax, SPxId enterId,
1654  Real& enterBound, Real& newUBbound,
1655  Real& newLBbound, Real& newCoPrhs);
1656  ///
1657  virtual void getEnterVals(SPxId id, Real& enterTest,
1658  Real& enterUB, Real& enterLB, Real& enterVal, Real& enterMax,
1659  Real& enterPric, SPxBasis::Desc::Status& enterStat, Real& enterRO);
1660  ///
1661  virtual void getEnterVals2(int leaveIdx,
1662  Real enterMax, Real& leaveBound);
1663  ///
1664  virtual void ungetEnterVal(SPxId enterId, SPxBasis::Desc::Status enterStat,
1665  Real leaveVal, const SVector& vec);
1666  ///
1667  virtual void rejectEnter(SPxId enterId,
1668  Real enterTest, SPxBasis::Desc::Status enterStat);
1669  ///
1670  virtual void rejectLeave(int leaveNum, SPxId leaveId,
1671  SPxBasis::Desc::Status leaveStat, const SVector* newVec = 0);
1672  ///
1673  virtual void setupPupdate(void);
1674  ///
1675  virtual void doPupdate(void);
1676  ///
1677  virtual void clearUpdateVecs(void);
1678  ///
1679  virtual void perturbMinEnter(void);
1680  /// perturb basis bounds.
1681  virtual void perturbMaxEnter(void);
1682  ///
1683  virtual void perturbMinLeave(void);
1684  /// perturb nonbasic bounds.
1685  virtual void perturbMaxLeave(void);
1686  //@}
1687 
1688  //------------------------------------
1689  /** The following methods serve for initializing the bounds for dual or
1690  * primal Simplex algorithm of entering or leaving type.
1691  */
1692  //@{
1693  ///
1694  void clearDualBounds(SPxBasis::Desc::Status, Real&, Real&) const;
1695  ///
1696  void setDualColBounds();
1697  ///
1698  void setDualRowBounds();
1699  /// setup feasibility bounds for entering algorithm
1700  void setPrimalBounds();
1701  ///
1702  void setEnterBound4Col(int, int);
1703  ///
1704  void setEnterBound4Row(int, int);
1705  ///
1706  virtual void setEnterBounds();
1707  ///
1708  void setLeaveBound4Row(int i, int n);
1709  ///
1710  void setLeaveBound4Col(int i, int n);
1711  ///
1712  virtual void setLeaveBounds();
1713  //@}
1714 
1715 public:
1716 
1717  //------------------------------------
1718  /** Limits and status inquiry */
1719  //@{
1720  /// set time limit.
1721  virtual void setTerminationTime(Real time = infinity);
1722  /// return time limit.
1723  virtual Real terminationTime() const;
1724  /// set iteration limit.
1725  virtual void setTerminationIter(int iteration = -1);
1726  /// return iteration limit.
1727  virtual int terminationIter() const;
1728  /// set refinement limit.
1729  virtual void setMaxRefinements(int p_maxrefinements);
1730  /// return refinement limit.
1731  virtual int maxRefinements() const;
1732  /// set objective limit.
1733  virtual void setTerminationValue(Real value = infinity);
1734  /// return objective limit.
1735  virtual Real terminationValue() const;
1736  /// get objective value of current solution.
1737  virtual Real objValue() const
1738  {
1739  return value();
1740  }
1741  /// get all results of last solve.
1742  Status
1743  getResult( Real* value = 0, Vector* primal = 0,
1744  Vector* slacks = 0, Vector* dual = 0,
1745  Vector* reduCost = 0) const;
1746 
1747 protected:
1748 
1749  /**@todo put the following basis methods near the variable status methods!*/
1750  /// converts basis status to VarStatus
1751  VarStatus basisStatusToVarStatus( SPxBasis::Desc::Status stat ) const;
1752 
1753  /// converts VarStatus to basis status for rows
1754  SPxBasis::Desc::Status varStatusToBasisStatusRow( int row, VarStatus stat )
1755  const;
1756 
1757  /// converts VarStatus to basis status for columns
1758  SPxBasis::Desc::Status varStatusToBasisStatusCol( int col, VarStatus stat )
1759  const;
1760 
1761 public:
1762 
1763  /// gets basis status for a single row
1764  VarStatus getBasisRowStatus( int row ) const;
1765 
1766  /// gets basis status for a single column
1767  VarStatus getBasisColStatus( int col ) const;
1768 
1769  /// get current basis, and return solver status.
1770  Status getBasis(VarStatus rows[], VarStatus cols[]) const;
1771 
1772  /// gets basis status
1773  SPxBasis::SPxStatus getBasisStatus() const
1774  {
1775  return SPxBasis::status();
1776  }
1777 
1778  /// check a given basis for validity.
1779  bool isBasisValid(DataArray<VarStatus> rows, DataArray<VarStatus> cols);
1780 
1781  /// set the lp solver's basis.
1782  void setBasis(const VarStatus rows[], const VarStatus cols[]);
1783 
1784  /// set the lp solver's basis status.
1785  void setBasisStatus( SPxBasis::SPxStatus stat )
1786  {
1787  if( m_status == OPTIMAL )
1788  m_status = UNKNOWN;
1789  SPxBasis::setStatus( stat );
1790  }
1791  /// reset cumulative time counter to zero.
1792  void resetCumulativeTime()
1793  {
1794  theCumulativeTime = 0.0;
1795  }
1796 
1797  /// get number of iterations of current solution.
1798  int iterations() const
1799  {
1800  return basis().iteration();
1801  }
1802 
1803  /// time spent in last call to method solve().
1804  Real time() const
1805  {
1806  return theTime.userTime();
1807  }
1808  /// cumulative time spent in all calls to method solve().
1809  Real cumulativeTime() const
1810  {
1811  return theCumulativeTime;
1812  }
1813  /// return const lp's rows if available.
1814  const LPRowSet& rows() const
1815  {
1816  return *lprowset();
1817  }
1818 
1819  /// return const lp's cols if available.
1820  const LPColSet& cols() const
1821  {
1822  return *lpcolset();
1823  }
1824 
1825  /// copy lower bound vector to \p p_low.
1826  void getLower(Vector& p_low) const
1827  {
1828  p_low = SPxLP::lower();
1829  }
1830  /// copy upper bound vector to \p p_up.
1831  void getUpper(Vector& p_up) const
1832  {
1833  p_up = SPxLP::upper();
1834  }
1835 
1836  /// copy lhs value vector to \p p_lhs.
1837  void getLhs(Vector& p_lhs) const
1838  {
1839  p_lhs = SPxLP::lhs();
1840  }
1841 
1842  /// copy rhs value vector to \p p_rhs.
1843  void getRhs(Vector& p_rhs) const
1844  {
1845  p_rhs = SPxLP::rhs();
1846  }
1847 
1848  /// optimization sense.
1849  SPxSense sense() const
1850  {
1851  return spxSense();
1852  }
1853 
1854  /// returns statistical information in form of a string.
1855  std::string statistics() const
1856  {
1857  std::stringstream s;
1858  s << basis().statistics()
1859  << "Solution time : " << std::setw(10) << std::fixed << std::setprecision(2) << time() << std::endl
1860  << "Iterations : " << std::setw(10) << iterations() << std::endl;
1861 
1862  return s.str();
1863  }
1864  //@}
1865 
1866  //------------------------------------
1867  /** Mapping between numbers and Ids */
1868  //@{
1869  /// RowId of \p i 'th inequality.
1870  SPxRowId rowId(int i) const
1871  {
1872  return rId(i);
1873  }
1874  /// ColId of \p i 'th column.
1875  SPxColId colId(int i) const
1876  {
1877  return cId(i);
1878  }
1879  //@}
1880 
1881  //------------------------------------
1882  /** Constructors / destructors */
1883  //@{
1884  /// default constructor.
1885  explicit
1886  SPxSolver( Type type = LEAVE,
1887  Representation rep = ROW );
1888  // virtual destructor
1889  virtual ~SPxSolver();
1890  //@}
1891 
1892  //------------------------------------
1893  /** Miscellaneous */
1894  //@{
1895  /// check consistency.
1896  bool isConsistent() const;
1897  //@}
1898 
1899  //------------------------------------
1900  /** assignment operator and copy constructor */
1901  //@{
1902  /// assignment operator
1903  SPxSolver& operator=(const SPxSolver& base);
1904  /// copy constructor
1905  SPxSolver(const SPxSolver& base);
1906  //@}
1907 
1908  void testVecs();
1909 };
1910 
1911 //
1912 // Auxiliary functions.
1913 //
1914 
1915 /// Pretty-printing of variable status.
1916 std::ostream& operator<<( std::ostream& os,
1917  const SPxSolver::VarStatus& status );
1918 
1919 /// Pretty-printing of solver status.
1920 std::ostream& operator<<( std::ostream& os,
1921  const SPxSolver::Status& status );
1922 
1923 /// Pretty-printing of algorithm.
1924 std::ostream& operator<<( std::ostream& os,
1925  const SPxSolver::Type& status );
1926 
1927 /// Pretty-printing of representation.
1928 std::ostream& operator<<( std::ostream& os,
1929  const SPxSolver::Representation& status );
1930 
1931 
1932 } // namespace soplex
1933 #endif // _SPXSOLVER_H_
1934 
1935 //-----------------------------------------------------------------------------
1936 //Emacs Local Variables:
1937 //Emacs mode:c++
1938 //Emacs c-basic-offset:3
1939 //Emacs tab-width:8
1940 //Emacs indent-tabs-mode:nil
1941 //Emacs End:
1942 //-----------------------------------------------------------------------------
© 2003-2013 by Zuse Institute Berlin (ZIB), Imprint Generated on Wed Jan 9 2013 for SoPlex by doxygen