Scippy

SoPlex

Sequential object-oriented simPlex

spxvecs.cpp
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-2019 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 #include <assert.h>
17 #include <iostream>
18 
19 #include "soplex/spxdefines.h"
20 #include "soplex/spxsolver.h"
21 #include "soplex/exceptions.h"
22 
23 namespace soplex
24 {
25 /** Initialize Vectors
26 
27  Computing the right hand side vector for the feasibility vector depends on
28  the chosen representation and type of the basis.
29 
30  In columnwise case, |theFvec| = \f$x_B = A_B^{-1} (- A_N x_N)\f$, where \f$x_N\f$
31  are either upper or lower bounds for the nonbasic variables (depending on
32  the variables |Status|). If these values remain unchanged throughout the
33  simplex algorithm, they may be taken directly from LP. However, in the
34  entering type algorith they are changed and, hence, retreived from the
35  column or row upper or lower bound vectors.
36 
37  In rowwise case, |theFvec| = \f$\pi^T_B = (c^T - 0^T A_N) A_B^{-1}\f$. However,
38  this applies only to leaving type algorithm, where no bounds on dual
39  variables are altered. In entering type algorithm they are changed and,
40  hence, retreived from the column or row upper or lower bound vectors.
41  */
43 {
44 
45  if(rep() == COLUMN)
46  {
47  theFrhs->clear();
48 
49  if(type() == LEAVE)
50  {
52 
53  for(int i = 0; i < nRows(); i++)
54  {
55  Real x;
56 
58 
59  if(!isBasic(stat))
60  {
61  switch(stat)
62  {
63  // columnwise cases:
65  continue;
66 
68  assert(EQ(lhs(i), rhs(i)));
69 
70  //lint -fallthrough
72  x = rhs(i);
73  break;
74 
76  x = lhs(i);
77  break;
78 
79  default:
80  MSG_ERROR(std::cerr << "ESVECS01 ERROR: "
81  << "inconsistent basis must not happen!"
82  << std::endl;)
83  throw SPxInternalCodeException("XSVECS01 This should never happen.");
84  }
85 
86  assert(x < infinity);
87  assert(x > -infinity);
88  (*theFrhs)[i] += x; // slack !
89  }
90  }
91  }
92  else
93  {
96  }
97  }
98  else
99  {
100  assert(rep() == ROW);
101 
102  if(type() == ENTER)
103  {
104  theFrhs->clear();
107  *theFrhs += maxObj();
108  }
109  else
110  {
111  ///@todo put this into a separate method
112  *theFrhs = maxObj();
113  const SPxBasis::Desc& ds = desc();
114 
115  for(int i = 0; i < nRows(); ++i)
116  {
117  SPxBasis::Desc::Status stat = ds.rowStatus(i);
118 
119  if(!isBasic(stat))
120  {
121  Real x;
122 
123  switch(stat)
124  {
126  continue;
127 
131  x = maxRowObj(i);
132  break;
133 
134  default:
135  assert(lhs(i) <= -infinity && rhs(i) >= infinity);
136  x = 0.0;
137  break;
138  }
139 
140  assert(x < infinity);
141  assert(x > -infinity);
142  // assert(x == 0.0);
143 
144  if(x != 0.0)
145  theFrhs->multAdd(x, vector(i));
146  }
147  }
148  }
149  }
150 }
151 
153 {
154 
155  assert(rep() == COLUMN);
156  assert(type() == LEAVE);
157 
158  for(int i = 0; i < nCols(); ++i)
159  {
161 
162  if(!isBasic(stat))
163  {
164  Real x;
165 
166  switch(stat)
167  {
168  // columnwise cases:
170  continue;
171 
173  assert(EQ(SPxLP::lower(i), SPxLP::upper(i)));
174 
175  //lint -fallthrough
177  x = SPxLP::upper(i);
178  break;
179 
181  x = SPxLP::lower(i);
182  break;
183 
184  default:
185  MSG_ERROR(std::cerr << "ESVECS02 ERROR: "
186  << "inconsistent basis must not happen!"
187  << std::endl;)
188  throw SPxInternalCodeException("XSVECS02 This should never happen.");
189  }
190 
191  assert(x < infinity);
192  assert(x > -infinity);
193 
194  if(x != 0.0)
195  theFrhs->multAdd(-x, vector(i));
196  }
197  }
198 }
199 
200 
201 /** This methods subtracts \f$A_N x_N\f$ or \f$\pi_N^T A_N\f$ from |theFrhs| as
202  specified by the |Status| of all nonbasic variables. The values of \f$x_N\f$ or
203  \f$\pi_N\f$ are taken from the passed arrays.
204  */
206  const Vector& ufb, ///< upper feasibility bound for variables
207  const Vector& lfb) ///< lower feasibility bound for variables
208 {
209 
210  const SPxBasis::Desc& ds = desc();
211 
212  for(int i = 0; i < coDim(); ++i)
213  {
214  SPxBasis::Desc::Status stat = ds.status(i);
215 
216  if(!isBasic(stat))
217  {
218  Real x;
219 
220  switch(stat)
221  {
225  continue;
226 
229  x = ufb[i];
230  break;
231 
234  x = lfb[i];
235  break;
236 
238  assert(EQ(lfb[i], ufb[i]));
239 
240  //lint -fallthrough
242  x = lfb[i];
243  break;
244 
245  default:
246  MSG_ERROR(std::cerr << "ESVECS03 ERROR: "
247  << "inconsistent basis must not happen!"
248  << std::endl;)
249  throw SPxInternalCodeException("XSVECS04 This should never happen.");
250  }
251 
252  assert(x < infinity);
253  assert(x > -infinity);
254 
255  if(x != 0.0)
256  theFrhs->multAdd(-x, vector(i));
257  }
258  }
259 }
260 
261 /** This methods subtracts \f$A_N x_N\f$ or \f$\pi_N^T A_N\f$ from |theFrhs| as
262  specified by the |Status| of all nonbasic variables. The values of
263  \f$x_N\f$ or \f$\pi_N\f$ are taken from the passed arrays.
264  */
266  Vector& coufb, ///< upper feasibility bound for covariables
267  Vector& colfb) ///< lower feasibility bound for covariables
268 {
269  const SPxBasis::Desc& ds = desc();
270 
271  for(int i = 0; i < dim(); ++i)
272  {
273  SPxBasis::Desc::Status stat = ds.coStatus(i);
274 
275  if(!isBasic(stat))
276  {
277  Real x;
278 
279  switch(stat)
280  {
284  continue;
285 
286  case SPxBasis::Desc::P_ON_LOWER : // negative slack bounds!
288  x = coufb[i];
289  break;
290 
291  case SPxBasis::Desc::P_ON_UPPER : // negative slack bounds!
293  x = colfb[i];
294  break;
295 
298 
299  if(colfb[i] != coufb[i])
300  {
301  MSG_WARNING((*spxout), (*spxout) << "WSVECS04 Frhs2[" << i << "]: " << stat << " "
302  << colfb[i] << " " << coufb[i]
303  << " shouldn't be" << std::endl;)
304 
305  if(isZero(colfb[i]) || isZero(coufb[i]))
306  colfb[i] = coufb[i] = 0.0;
307  else
308  {
309  Real mid = (colfb[i] + coufb[i]) / 2.0;
310  colfb[i] = coufb[i] = mid;
311  }
312  }
313 
314  assert(EQ(colfb[i], coufb[i]));
315  x = colfb[i];
316  break;
317 
318  default:
319  MSG_ERROR(std::cerr << "ESVECS05 ERROR: "
320  << "inconsistent basis must not happen!"
321  << std::endl;)
322  throw SPxInternalCodeException("XSVECS05 This should never happen.");
323  }
324 
325  assert(x < infinity);
326  assert(x > -infinity);
327 
328  (*theFrhs)[i] -= x; // This is a slack, so no need to multiply a vector.
329  }
330  }
331 }
332 
333 /** Computing the right hand side vector for |theCoPvec| depends on
334  the type of the simplex algorithm. In entering algorithms, the
335  values are taken from the inequality's right handside or the
336  column's objective value.
337 
338  In contrast to this leaving algorithms take the values from vectors
339  |theURbound| and so on.
340 
341  We reflect this difference by providing two pairs of methods
342  |computeEnterCoPrhs(n, stat)| and |computeLeaveCoPrhs(n, stat)|. The first
343  pair operates for entering algorithms, while the second one is intended for
344  leaving algorithms. The return value of these methods is the right hand
345  side value for the \f$n\f$-th row or column id, respectively, if it had the
346  passed |Status| for both.
347 
348  Both methods are again split up into two methods named |...4Row(i,n)| and
349  |...4Col(i,n)|, respectively. They do their job for the |i|-th basis
350  variable, being the |n|-th row or column.
351 */
353 {
354  assert(baseId(i).isSPxRowId());
355  assert(number(SPxRowId(baseId(i))) == n);
356 
357  switch(desc().rowStatus(n))
358  {
359  // rowwise representation:
361  assert(lhs(n) > -infinity);
362  assert(EQ(rhs(n), lhs(n)));
363 
364  //lint -fallthrough
366  assert(rep() == ROW);
367  assert(rhs(n) < infinity);
368  (*theCoPrhs)[i] = rhs(n);
369  break;
370 
372  assert(rep() == ROW);
373  assert(lhs(n) > -infinity);
374  (*theCoPrhs)[i] = lhs(n);
375  break;
376 
377  // columnwise representation:
378  // slacks must be left 0!
379  default:
380  (*theCoPrhs)[i] = maxRowObj(n);
381  break;
382  }
383 }
384 
386 {
387  assert(baseId(i).isSPxColId());
388  assert(number(SPxColId(baseId(i))) == n);
389 
390  switch(desc().colStatus(n))
391  {
392  // rowwise representation:
394  assert(EQ(SPxLP::upper(n), SPxLP::lower(n)));
395  assert(SPxLP::lower(n) > -infinity);
396 
397  //lint -fallthrough
399  assert(rep() == ROW);
400  assert(SPxLP::upper(n) < infinity);
401  (*theCoPrhs)[i] = SPxLP::upper(n);
402  break;
403 
405  assert(rep() == ROW);
406  assert(SPxLP::lower(n) > -infinity);
407  (*theCoPrhs)[i] = SPxLP::lower(n);
408  break;
409 
410  // columnwise representation:
416  assert(rep() == COLUMN);
417  (*theCoPrhs)[i] = maxObj(n);
418  break;
419 
420  default: // variable left 0
421  (*theCoPrhs)[i] = 0;
422  break;
423  }
424 }
425 
427 {
428  assert(type() == ENTER);
429 
430  for(int i = dim() - 1; i >= 0; --i)
431  {
432  SPxId l_id = baseId(i);
433 
434  if(l_id.isSPxRowId())
436  else
438  }
439 }
440 
442 {
443  assert(baseId(i).isSPxRowId());
444  assert(number(SPxRowId(baseId(i))) == n);
445 
446  switch(desc().rowStatus(n))
447  {
450  assert(theLRbound[n] > -infinity);
451  assert(EQ(theURbound[n], theLRbound[n]));
452 
453  //lint -fallthrough
456  assert(theURbound[n] < infinity);
457  (*theCoPrhs)[i] = theURbound[n];
458  break;
459 
462  assert(theLRbound[n] > -infinity);
463  (*theCoPrhs)[i] = theLRbound[n];
464  break;
465 
466  default:
467  (*theCoPrhs)[i] = maxRowObj(n);
468  break;
469  }
470 }
471 
473 {
474  assert(baseId(i).isSPxColId());
475  assert(number(SPxColId(baseId(i))) == n);
476 
477  switch(desc().colStatus(n))
478  {
482  assert(theLCbound[n] > -infinity);
483  assert(EQ(theUCbound[n], theLCbound[n]));
484 
485  //lint -fallthrough
488  assert(theUCbound[n] < infinity);
489  (*theCoPrhs)[i] = theUCbound[n];
490  break;
491 
494  assert(theLCbound[n] > -infinity);
495  (*theCoPrhs)[i] = theLCbound[n];
496  break;
497 
498  default:
499  (*theCoPrhs)[i] = maxObj(n);
500  // (*theCoPrhs)[i] = 0;
501  break;
502  }
503 }
504 
506 {
507  assert(type() == LEAVE);
508 
509  for(int i = dim() - 1; i >= 0; --i)
510  {
511  SPxId l_id = baseId(i);
512 
513  if(l_id.isSPxRowId())
515  else
517  }
518 }
519 
520 
521 /** When computing the copricing vector, we expect the pricing vector to be
522  setup correctly. Then computing the copricing vector is nothing but
523  computing all scalar products of the pricing vector with the vectors of the
524  LPs constraint matrix.
525  */
527 {
528  int i;
529 
530  for(i = coDim() - 1; i >= 0; --i)
531  (*thePvec)[i] = vector(i) * (*theCoPvec);
532 }
533 
535 {
536  SSVector& p = thePvec->delta();
537  SSVector& c = theCoPvec->delta();
538 
539  if(c.isSetup())
540  {
541  if(c.size() < 0.95 * theCoPvec->dim())
543  else
544  p.assign2product(c, *thevectors);
545  }
546  else
547  {
549  }
550 
551  p.setup();
552 }
553 
555 {
556  theCoPvec->update();
557 
558  if(pricing() == FULL)
559  {
560  // thePvec->delta().setup();
561  thePvec->update();
562  }
563 }
564 } // namespace soplex
const VectorBase< R > & rhs() const
Returns right hand side vector.
Definition: spxlpbase.h:221
Exception class for things that should NEVER happen.This class is derived from the SoPlex exception b...
Definition: exceptions.h:109
Status & coStatus(int i)
Definition: spxbasis.h:284
bool isSetup() const
Returns setup status.
Definition: ssvectorbase.h:121
primal variable is fixed to both bounds
Definition: spxbasis.h:190
primal or dual variable is undefined
Definition: spxbasis.h:195
const VectorBase< Real > & upper() const
Returns upper bound vector.
Definition: spxlpbase.h:461
SPxOut * spxout
message handler
Definition: spxsolver.h:463
THREADLOCAL const Real infinity
Definition: spxdefines.cpp:26
DVector * theUbound
Upper bound for vars.
Definition: spxsolver.h:375
Pricing pricing() const
return current Pricing.
Definition: spxsolver.h:518
DVector * theCoUbound
Upper bound for covars.
Definition: spxsolver.h:377
Exception classes for SoPlex.
Status & rowStatus(int i)
Definition: spxbasis.h:239
int number(const SPxRowId &id) const
Returns the row number of the row with identifier id.
Definition: spxlpbase.h:527
virtual void computeFrhsXtra()
Definition: spxvecs.cpp:152
int dim() const
dimension of basis matrix.
Definition: spxsolver.h:1075
Ids for LP columns.Class SPxColId provides DataKeys for the column indices of an SPxLP.
Definition: spxid.h:36
virtual void doPupdate(void)
Definition: spxvecs.cpp:554
rowwise representation.
Definition: spxsolver.h:107
virtual void computeLeaveCoPrhs()
compute theCoPrhs for leaving Simplex.
Definition: spxvecs.cpp:505
dual variable is left free, but unset
Definition: spxbasis.h:191
int nRows() const
Returns number of rows in LP.
Definition: spxlpbase.h:152
Entering Simplex.
Definition: spxsolver.h:134
SSVectorBase< R > & assign2productAndSetup(const SVSetBase< S > &A, SSVectorBase< T > &x)
Assigns SSVectorBase to thereby setting up x.
Definition: basevectors.h:798
primal variable is set to its upper bound
Definition: spxbasis.h:188
Generic Ids for LP rows or columns.Both SPxColIds and SPxRowIds may be treated uniformly as SPxIds: ...
Definition: spxid.h:85
Leaving Simplex.
Definition: spxsolver.h:143
double Real
Definition: spxdefines.h:218
DVector * theFrhs
Definition: spxsolver.h:360
SSVectorBase< R > & assign2product(const SSVectorBase< S > &x, const SVSetBase< T > &A)
Assigns to SSVectorBase.
Definition: basevectors.h:539
DVector * theLbound
Lower bound for vars.
Definition: spxsolver.h:376
UpdateVector * thePvec
Definition: spxsolver.h:368
void update()
Perform the update.
Definition: updatevector.h:147
#define MSG_ERROR(x)
Prints out message x if the verbosity level is at least SPxOut::ERROR.
Definition: spxdefines.h:114
const SVSet * thecovectors
the LP coVectors according to representation
Definition: spxsolver.h:338
void computeEnterCoPrhs4Col(int i, int n)
Definition: spxvecs.cpp:385
dual variable is set to its upper bound
Definition: spxbasis.h:192
main LP solver class
virtual void setupPupdate(void)
Definition: spxvecs.cpp:534
bool isBasic(SPxBasis::Desc::Status stat) const
does stat describe a basic index ?
Definition: spxsolver.h:1263
primal variable is left free, but unset
Definition: spxbasis.h:189
const VectorBase< R > & lhs() const
Returns left hand side vector.
Definition: spxlpbase.h:255
Basis descriptor.
Definition: spxbasis.h:104
DVector theURbound
Upper Row Feasibility bound.
Definition: spxsolver.h:346
const VectorBase< R > & maxRowObj() const
Definition: spxlpbase.h:300
Status & colStatus(int i)
Definition: spxbasis.h:254
SPxId & baseId(int i)
Definition: spxbasis.h:504
DVector * theCoLbound
Lower bound for covars.
Definition: spxsolver.h:378
Debugging, floating point type and parameter definitions.
bool EQ(Real a, Real b, Real eps=Param::epsilon())
returns true iff |a-b| <= eps
Definition: spxdefines.h:376
Full pricing.
Definition: spxsolver.h:160
DVector theLRbound
Lower Row Feasibility bound.
Definition: spxsolver.h:347
DVector theUCbound
Upper Column Feasibility bound.
Definition: spxsolver.h:348
void computeLeaveCoPrhs4Col(int i, int n)
Definition: spxvecs.cpp:472
VectorBase< R > & multAdd(const S &x, const VectorBase< T > &vec)
Addition of scaled vector.
Definition: vectorbase.h:412
void computeEnterCoPrhs4Row(int i, int n)
Definition: spxvecs.cpp:352
int dim() const
Dimension of vector.
Definition: vectorbase.h:217
int size() const
Returns the number of nonzeros.
Definition: ssvectorbase.h:200
Everything should be within this namespace.
void computeLeaveCoPrhs4Row(int i, int n)
Definition: spxvecs.cpp:441
#define MSG_WARNING(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::WARNING.
Definition: spxdefines.h:116
primal variable is set to its lower bound
Definition: spxbasis.h:187
const VectorBase< R > & maxObj() const
Returns objective vector for maximization problem.
Definition: spxlpbase.h:434
void setup()
Initializes nonzero indices for elements with absolute values above epsilon and sets all other elemen...
Definition: ssvectorbase.h:133
DVector theLCbound
Lower Column Feasibility bound.
Definition: spxsolver.h:349
void computeFrhs()
compute feasibility vector from scratch.
Definition: spxvecs.cpp:42
void clear()
Set vector to 0.
Definition: vectorbase.h:262
bool isSPxRowId() const
is id a row id?
Definition: spxid.h:161
dual variable is set to its lower bound
Definition: spxbasis.h:193
Type type() const
return current Type.
Definition: spxsolver.h:512
int coDim() const
codimension.
Definition: spxsolver.h:1080
int nCols() const
Returns number of columns in LP.
Definition: spxlpbase.h:158
const SVSet * thevectors
the LP vectors according to representation
Definition: spxsolver.h:337
dual variable has two bounds
Definition: spxbasis.h:194
Ids for LP rows.Class SPxRowId provides DataKeys for the row indices of an SPxLP. ...
Definition: spxid.h:55
virtual void computeEnterCoPrhs()
compute theCoPrhs for entering Simplex.
Definition: spxvecs.cpp:426
bool isZero(Real a, Real eps=Param::epsilon())
returns true iff |a| <= eps
Definition: spxdefines.h:412
const SVector & vector(int i) const
i &#39;th vector.
Definition: spxsolver.h:1157
SSVector & delta()
update vector , writeable
Definition: updatevector.h:122
Status & status(int i)
Definition: spxbasis.h:269
void computePvec()
compute entire pVec().
Definition: spxvecs.cpp:526
UpdateVector * theCoPvec
Definition: spxsolver.h:366
void computeFrhs2(Vector &, Vector &)
Definition: spxvecs.cpp:265
Status
Status of a variable.
Definition: spxbasis.h:185
const VectorBase< Real > & lower() const
Returns (internal and possibly scaled) lower bound vector.
Definition: spxlpbase.h:488
virtual void computeFrhs1(const Vector &, const Vector &)
Definition: spxvecs.cpp:205
const Desc & desc() const
Definition: spxbasis.h:464
Representation rep() const
return the current basis representation.
Definition: spxsolver.h:506
columnwise representation.
Definition: spxsolver.h:108
SSVectorBase< R > & assign2product4setup(const SVSetBase< S > &A, const SSVectorBase< T > &x)
Assigns SSVectorBase to for a setup x.
Definition: basevectors.h:570