SoPlex Doxygen Documentation
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-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 #include <assert.h>
17 #include <iostream>
18 
19 #include "spxdefines.h"
20 #include "spxsolver.h"
21 #include "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  */
42 void SPxSolver::computeFrhs()
43 {
44  METHOD( "SPxSolver::computeFrhs()" );
45 
46  if (rep() == COLUMN)
47  {
48  theFrhs->clear();
49 
50  if (type() == LEAVE)
51  {
52  computeFrhsXtra();
53 
54  for(int i = 0; i < nRows(); i++)
55  {
56  Real x;
57 
58  SPxBasis::Desc::Status stat = desc().rowStatus(i);
59 
60  if (!isBasic(stat))
61  {
62  switch (stat)
63  {
64  // columnwise cases:
65  case SPxBasis::Desc::P_FREE :
66  continue;
67 
68  case (SPxBasis::Desc::P_ON_UPPER + SPxBasis::Desc::P_ON_LOWER) :
69  assert(lhs(i) == rhs(i));
70  //lint -fallthrough
71  case SPxBasis::Desc::P_ON_UPPER :
72  x = rhs(i);
73  break;
74  case SPxBasis::Desc::P_ON_LOWER :
75  x = lhs(i);
76  break;
77 
78  default:
79  MSG_ERROR( spxout << "ESVECS01 ERROR: "
80  << "inconsistent basis must not happen!"
81  << std::endl; )
82  throw SPxInternalCodeException("XSVECS01 This should never happen.");
83  }
84  assert(x < infinity);
85  assert(x > -infinity);
86  (*theFrhs)[i] += x; // slack !
87  }
88  }
89  }
90  else
91  {
92  computeFrhs1(*theUbound, *theLbound);
93  computeFrhs2(*theCoUbound, *theCoLbound);
94  }
95  }
96  else
97  {
98  assert(rep() == ROW);
99 
100  if (type() == ENTER)
101  {
102  theFrhs->clear();
103  computeFrhs1(*theUbound, *theLbound);
104  computeFrhs2(*theCoUbound, *theCoLbound);
105  *theFrhs += maxObj();
106  }
107  else
108  *theFrhs = maxObj();
109  }
110 }
111 
112 void SPxSolver::computeFrhsXtra()
113 {
114  METHOD( "SPxSolver::computeFrhsXtra()" );
115 
116  assert(rep() == COLUMN);
117  assert(type() == LEAVE);
118 
119  for (int i = 0; i < nCols(); ++i)
120  {
121  SPxBasis::Desc::Status stat = desc().colStatus(i);
122 
123  if (!isBasic(stat))
124  {
125  Real x;
126 
127  switch (stat)
128  {
129  // columnwise cases:
130  case SPxBasis::Desc::P_FREE :
131  continue;
132 
133  case (SPxBasis::Desc::P_ON_UPPER + SPxBasis::Desc::P_ON_LOWER) :
134  assert(SPxLP::lower(i) == SPxLP::upper(i));
135  //lint -fallthrough
136  case SPxBasis::Desc::P_ON_UPPER :
137  x = SPxLP::upper(i);
138  break;
139  case SPxBasis::Desc::P_ON_LOWER :
140  x = SPxLP::lower(i);
141  break;
142 
143  default:
144  MSG_ERROR( spxout << "ESVECS02 ERROR: "
145  << "inconsistent basis must not happen!"
146  << std::endl; )
147  throw SPxInternalCodeException("XSVECS02 This should never happen.");
148  }
149  assert(x < infinity);
150  assert(x > -infinity);
151 
152  if (x != 0.0)
153  theFrhs->multAdd(-x, vector(i));
154  }
155  }
156 }
157 
158 
159 /** This methods subtracts \f$A_N x_N\f$ or \f$\pi_N^T A_N\f$ from |theFrhs| as
160  specified by the |Status| of all nonbasic variables. The values of \f$x_N\f$ or
161  \f$\pi_N\f$ are taken from the passed arrays.
162  */
163 void SPxSolver::computeFrhs1(
164  const Vector& ufb, ///< upper feasibility bound for variables
165  const Vector& lfb) ///< lower feasibility bound for variables
166 {
167  METHOD( "SPxSolver::computeFrhs1()" );
168 
169  const SPxBasis::Desc& ds = desc();
170 
171  for (int i = 0; i < coDim(); ++i)
172  {
173  SPxBasis::Desc::Status stat = ds.status(i);
174 
175  if (!isBasic(stat))
176  {
177  Real x;
178 
179  switch (stat)
180  {
181  case SPxBasis::Desc::D_FREE :
182  case SPxBasis::Desc::D_UNDEFINED :
183  case SPxBasis::Desc::P_FREE :
184  continue;
185 
186  case SPxBasis::Desc::P_ON_UPPER :
187  case SPxBasis::Desc::D_ON_UPPER :
188  x = ufb[i];
189  break;
190  case SPxBasis::Desc::P_ON_LOWER :
191  case SPxBasis::Desc::D_ON_LOWER :
192  x = lfb[i];
193  break;
194 
195  case (SPxBasis::Desc::P_ON_UPPER + SPxBasis::Desc::P_ON_LOWER) :
196  assert(lfb[i] == ufb[i]);
197  case (SPxBasis::Desc::D_ON_UPPER + SPxBasis::Desc::D_ON_LOWER) :
198  x = lfb[i];
199  break;
200 
201  default:
202  MSG_ERROR( spxout << "ESVECS03 ERROR: "
203  << "inconsistent basis must not happen!"
204  << std::endl; )
205  throw SPxInternalCodeException("XSVECS04 This should never happen.");
206  }
207  assert(x < infinity);
208  assert(x > -infinity);
209 
210  if (x != 0.0)
211  theFrhs->multAdd(-x, vector(i));
212  }
213  }
214 }
215 
216 /** This methods subtracts \f$A_N x_N\f$ or \f$\pi_N^T A_N\f$ from |theFrhs| as
217  specified by the |Status| of all nonbasic variables. The values of
218  \f$x_N\f$ or \f$\pi_N\f$ are taken from the passed arrays.
219  */
220 void SPxSolver::computeFrhs2(
221  const Vector& coufb, ///< upper feasibility bound for covariables
222  const Vector& colfb) ///< lower feasibility bound for covariables
223 {
224  METHOD( "SPxSolver::computeFrhs2()" );
225  const SPxBasis::Desc& ds = desc();
226 
227  for(int i = 0; i < dim(); ++i)
228  {
229  SPxBasis::Desc::Status stat = ds.coStatus(i);
230 
231  if (!isBasic(stat))
232  {
233  Real x;
234 
235  switch (stat)
236  {
237  case SPxBasis::Desc::D_FREE :
238  case SPxBasis::Desc::D_UNDEFINED :
239  case SPxBasis::Desc::P_FREE :
240  continue;
241 
242  case SPxBasis::Desc::P_ON_LOWER : // negative slack bounds!
243  case SPxBasis::Desc::D_ON_UPPER :
244  x = coufb[i];
245  break;
246  case SPxBasis::Desc::P_ON_UPPER : // negative slack bounds!
247  case SPxBasis::Desc::D_ON_LOWER :
248  x = colfb[i];
249  break;
250  case SPxBasis::Desc::P_FIXED :
251  case SPxBasis::Desc::D_ON_BOTH :
252 
253  if (colfb[i] != coufb[i])
254  {
255  MSG_WARNING( spxout << "WSVECS04 Frhs2: " << stat << " "
256  << colfb[i] << " " << coufb[i]
257  << " shouldn't be" << std::endl; )
258  }
259  //assert(colfb[i] == coufb[i]);
260  x = colfb[i];
261  break;
262 
263  default:
264  MSG_ERROR( spxout << "ESVECS05 ERROR: "
265  << "inconsistent basis must not happen!"
266  << std::endl; )
267  throw SPxInternalCodeException("XSVECS05 This should never happen.");
268  }
269  assert(x < infinity);
270  assert(x > -infinity);
271 
272  (*theFrhs)[i] -= x; // This is a slack, so no need to multiply a vector.
273  }
274  }
275 }
276 
277 /** Computing the right hand side vector for |theCoPvec| depends on
278  the type of the simplex algorithm. In entering algorithms, the
279  values are taken from the inequality's right handside or the
280  column's objective value.
281 
282  In contrast to this leaving algorithms take the values from vectors
283  |theURbound| and so on.
284 
285  We reflect this difference by providing two pairs of methods
286  |computeEnterCoPrhs(n, stat)| and |computeLeaveCoPrhs(n, stat)|. The first
287  pair operates for entering algorithms, while the second one is intended for
288  leaving algorithms. The return value of these methods is the right hand
289  side value for the \f$n\f$-th row or column id, respectively, if it had the
290  passed |Status| for both.
291 
292  Both methods are again split up into two methods named |...4Row(i,n)| and
293  |...4Col(i,n)|, respectively. They do their job for the |i|-th basis
294  variable, being the |n|-th row or column.
295 */
296 void SPxSolver::computeEnterCoPrhs4Row(int i, int n)
297 {
298  METHOD( "SPxSolver::computeEnterCoPrhs4Row()" );
299  assert(baseId(i).isSPxRowId());
300  assert(number(SPxRowId(baseId(i))) == n);
301 
302  switch (desc().rowStatus(n))
303  {
304  // rowwise representation:
305  case SPxBasis::Desc::P_ON_LOWER + SPxBasis::Desc::P_ON_UPPER :
306  assert(lhs(n) > -infinity);
307  assert(rhs(n) == lhs(n));
308  //lint -fallthrough
309  case SPxBasis::Desc::P_ON_UPPER :
310  assert(rep() == ROW);
311  assert(rhs(n) < infinity);
312  (*theCoPrhs)[i] = rhs(n);
313  break;
314  case SPxBasis::Desc::P_ON_LOWER :
315  assert(rep() == ROW);
316  assert(lhs(n) > -infinity);
317  (*theCoPrhs)[i] = lhs(n);
318  break;
319 
320  // columnwise representation:
321  // slacks must be left 0!
322  default:
323  (*theCoPrhs)[i] = 0;
324  break;
325  }
326 }
327 
328 void SPxSolver::computeEnterCoPrhs4Col(int i, int n)
329 {
330  METHOD( "SPxSolver::computeEnterCoPrhs4Col()" );
331  assert(baseId(i).isSPxColId());
332  assert(number(SPxColId(baseId(i))) == n);
333  switch (desc().colStatus(n))
334  {
335  // rowwise representation:
336  case SPxBasis::Desc::P_ON_LOWER + SPxBasis::Desc::P_ON_UPPER :
337  assert(SPxLP::upper(n) == SPxLP::lower(n));
338  assert(SPxLP::lower(n) > -infinity);
339  //lint -fallthrough
340  case SPxBasis::Desc::P_ON_UPPER :
341  assert(rep() == ROW);
342  assert(SPxLP::upper(n) < infinity);
343  (*theCoPrhs)[i] = SPxLP::upper(n);
344  break;
345  case SPxBasis::Desc::P_ON_LOWER :
346  assert(rep() == ROW);
347  assert(SPxLP::lower(n) > -infinity);
348  (*theCoPrhs)[i] = SPxLP::lower(n);
349  break;
350 
351  // columnwise representation:
352  case SPxBasis::Desc::D_UNDEFINED :
353  case SPxBasis::Desc::D_ON_UPPER + SPxBasis::Desc::D_ON_LOWER :
354  case SPxBasis::Desc::D_ON_UPPER :
355  case SPxBasis::Desc::D_ON_LOWER :
356  case SPxBasis::Desc::D_FREE :
357  assert(rep() == COLUMN);
358  (*theCoPrhs)[i] = maxObj(n);
359  break;
360 
361  default: // variable left 0
362  (*theCoPrhs)[i] = 0;
363  break;
364  }
365 }
366 
367 void SPxSolver::computeEnterCoPrhs()
368 {
369  METHOD( "SPxSolver::computeEnterCoPrhs()" );
370  assert(type() == ENTER);
371 
372  for (int i = dim() - 1; i >= 0; --i)
373  {
374  SPxId l_id = baseId(i);
375  if (l_id.isSPxRowId())
376  computeEnterCoPrhs4Row(i, number(SPxRowId(l_id)));
377  else
378  computeEnterCoPrhs4Col(i, number(SPxColId(l_id)));
379  }
380 }
381 
382 void SPxSolver::computeLeaveCoPrhs4Row(int i, int n)
383 {
384  METHOD( "SPxSolver::computeLeaveCoPrhs4Row()" );
385  assert(baseId(i).isSPxRowId());
386  assert(number(SPxRowId(baseId(i))) == n);
387  switch (desc().rowStatus(n))
388  {
389  case SPxBasis::Desc::D_ON_LOWER + SPxBasis::Desc::D_ON_UPPER :
390  case SPxBasis::Desc::P_ON_LOWER + SPxBasis::Desc::P_ON_UPPER :
391  assert(theLRbound[n] > -infinity);
392  assert(theURbound[n] == theLRbound[n]);
393  //lint -fallthrough
394  case SPxBasis::Desc::D_ON_UPPER :
395  case SPxBasis::Desc::P_ON_UPPER :
396  assert(theURbound[n] < infinity);
397  (*theCoPrhs)[i] = theURbound[n];
398  break;
399 
400  case SPxBasis::Desc::D_ON_LOWER :
401  case SPxBasis::Desc::P_ON_LOWER :
402  assert(theLRbound[n] > -infinity);
403  (*theCoPrhs)[i] = theLRbound[n];
404  break;
405 
406  default:
407  (*theCoPrhs)[i] = 0;
408  break;
409  }
410 }
411 
412 void SPxSolver::computeLeaveCoPrhs4Col(int i, int n)
413 {
414  METHOD( "SPxSolver::computeLeaveCoPrhs4Col()" );
415  assert(baseId(i).isSPxColId());
416  assert(number(SPxColId(baseId(i))) == n);
417  switch (desc().colStatus(n))
418  {
419  case SPxBasis::Desc::D_UNDEFINED :
420  case SPxBasis::Desc::D_ON_LOWER + SPxBasis::Desc::D_ON_UPPER :
421  case SPxBasis::Desc::P_ON_LOWER + SPxBasis::Desc::P_ON_UPPER :
422  assert(theLCbound[n] > -infinity);
423  assert(theUCbound[n] == theLCbound[n]);
424  //lint -fallthrough
425  case SPxBasis::Desc::D_ON_LOWER :
426  case SPxBasis::Desc::P_ON_UPPER :
427  assert(theUCbound[n] < infinity);
428  (*theCoPrhs)[i] = theUCbound[n];
429  break;
430  case SPxBasis::Desc::D_ON_UPPER :
431  case SPxBasis::Desc::P_ON_LOWER :
432  assert(theLCbound[n] > -infinity);
433  (*theCoPrhs)[i] = theLCbound[n];
434  break;
435 
436  default:
437  (*theCoPrhs)[i] = maxObj(n);
438  // (*theCoPrhs)[i] = 0;
439  break;
440  }
441 }
442 
443 void SPxSolver::computeLeaveCoPrhs()
444 {
445  METHOD( "SPxSolver::computeLeaveCoPrhs()" );
446  assert(type() == LEAVE);
447 
448  for (int i = dim() - 1; i >= 0; --i)
449  {
450  SPxId l_id = baseId(i);
451  if (l_id.isSPxRowId())
452  computeLeaveCoPrhs4Row(i, number(SPxRowId(l_id)));
453  else
454  computeLeaveCoPrhs4Col(i, number(SPxColId(l_id)));
455  }
456 }
457 
458 
459 /** When computing the copricing vector, we expect the pricing vector to be
460  setup correctly. Then computing the copricing vector is nothing but
461  computing all scalar products of the pricing vector with the vectors of the
462  LPs constraint matrix.
463  */
464 void SPxSolver::computePvec()
465 {
466  METHOD( "SPxSolver::computePvec()" );
467  int i;
468 
469  for (i = coDim() - 1; i >= 0; --i)
470  (*thePvec)[i] = vector(i) * (*theCoPvec);
471 }
472 
473 void SPxSolver::setupPupdate(void)
474 {
475  METHOD( "SPxSolver::setupPupdate()" );
476  SSVector& p = thePvec->delta();
477  SSVector& c = theCoPvec->delta();
478 
479  if (c.isSetup())
480  {
481  if (c.size() < 0.95 * theCoPvec->dim())
482  p.assign2product4setup(*thecovectors, c);
483  else
484  p.assign2product(c, *thevectors);
485  }
486  else
487  {
488  p.assign2productAndSetup(*thecovectors, c);
489  }
490 
491  p.setup();
492 }
493 
494 void SPxSolver::doPupdate(void)
495 {
496  METHOD( "SPxSolver::doPupdate()" );
497  theCoPvec->update();
498  if (pricing() == FULL)
499  {
500  // thePvec->delta().setup();
501  thePvec->update();
502  }
503 }
504 } // namespace soplex
505 
506 //-----------------------------------------------------------------------------
507 //Emacs Local Variables:
508 //Emacs mode:c++
509 //Emacs c-basic-offset:3
510 //Emacs tab-width:8
511 //Emacs indent-tabs-mode:nil
512 //Emacs End:
513 //-----------------------------------------------------------------------------
© 2003-2013 by Zuse Institute Berlin (ZIB), Imprint Generated on Wed Jan 9 2013 for SoPlex by doxygen