SoPlex Doxygen Documentation
enter.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 // #define DEBUGGING 1
17 
18 /* \SubSection{Updating the Basis for Entering Variables}
19  */
20 #include <assert.h>
21 
22 #include "spxdefines.h"
23 #include "soplex.h"
24 #include "spxratiotester.h"
25 #include "spxout.h"
26 #include "exceptions.h"
27 
28 namespace soplex
29 {
30 
31 /*
32 In the entering simplex algorithms (i.e. iteratively a vector is selected to
33 \em enter the simplex basis as in the dual rowwise and primal columnwise case)
34 let $A$ denote the current basis, $x$ and entering vector and $f$ the
35 feasibility vector. For a feasible basis $l \le f \le u$ holds.
36 For the rowwise case $f$ is obtained by solving $f^T = c^T A^{-1}$,
37 wherease in columnwisecase $f = A^{-1} b$.
38 
39 Let us further consider the rowwise case. Exchanging $x$ with the $i$-th
40 vector of $A$ yields
41 
42 \begin{equation}\label{update.eq}
43  A^{(i)} = E_i A \hbox{, with } E_i = I + e_i (x^T A^{-1} - e_i^T).
44 \end{equation}
45 
46 With $E_i^{-1} = I + e_i \frac{e_i^T - \delta^T}{\delta}$,
47 $\delta^T = x^T A^{-1}$ one gets the new feasibility vector
48 
49 \begin{eqnarray*}
50  (f^{(i)})^T
51  &=& c^T (A^{(i)})^{-1} \\
52  &=& c^T A^{-1} + c^T A^{-1} e_i \frac{e_i^T - \delta^T}{\delta_i} \\
53  &=& f^T + \frac{f_i}{\delta_i} e_i^T - \frac{f_i}{\delta_i} \delta^T. \\
54 \end{eqnarray*}
55 
56 The selection of the leaving vector $i^*$ for the basis must ensure, that for
57 all $j \ne i^*$ $f^{(i^*)}_j$ remains within its bounds $l_j$ and $u_j$.
58  */
59 
60 
61 /*
62  Testing all values of |pVec| against its bounds. If $i$, say, is violated
63  the violation is saved as negative value in |theTest[i]|.
64  */
65 Real SPxSolver::test(int i, SPxBasis::Desc::Status stat) const
66 {
67  METHOD( "SPxSolver::test()" );
68  assert(type() == ENTER);
69  assert(!isBasic(stat));
70 
71  Real x;
72 
73  switch (stat)
74  {
75  case SPxBasis::Desc::D_FREE:
76  case SPxBasis::Desc::D_ON_BOTH:
77  assert(rep() == ROW);
78  x = (*thePvec)[i] - lhs(i);
79  if (x < 0)
80  return x;
81  // no break: next is else case
82  //lint -fallthrough
83  case SPxBasis::Desc::D_ON_LOWER:
84  assert(rep() == ROW);
85  return rhs(i) - (*thePvec)[i];
86  case SPxBasis::Desc::D_ON_UPPER:
87  assert(rep() == ROW);
88  return (*thePvec)[i] - lhs(i);
89 
90  case SPxBasis::Desc::P_ON_UPPER:
91  assert(rep() == COLUMN);
92  return maxObj(i) - (*thePvec)[i];
93  case SPxBasis::Desc::P_ON_LOWER:
94  assert(rep() == COLUMN);
95  return (*thePvec)[i] - maxObj(i);
96  case SPxBasis::Desc::P_FREE :
97  x = maxObj(i) - (*thePvec)[i];
98  return (x < 0) ? x : -x;
99 
100  default:
101  return 0;
102  }
103 }
104 
105 void SPxSolver::computeTest()
106 {
107  METHOD( "SPxSolver::computeTest()" );
108 
109  const SPxBasis::Desc& ds = desc();
110  Real pricingTol = leavetol();
111  infeasibilitiesCo.clear();
112  int ninfeasibilities = 0;
113 
114  for(int i = 0; i < coDim(); ++i)
115  {
116  SPxBasis::Desc::Status stat = ds.status(i);
117 
118  if(isBasic(stat))
119  theTest[i] = 0.0;
120  else
121  {
122  theTest[i] = test(i, stat);
123 
124  if( remainingRoundsEnterCo == 0 )
125  {
126  if( theTest[i] < -pricingTol )
127  {
128  assert(infeasibilitiesCo.size() < infeasibilitiesCo.max());
129  infeasibilitiesCo.addIdx(i);
130  isInfeasibleCo[i] = true;
131  ++ninfeasibilities;
132  }
133  else
134  isInfeasibleCo[i] = false;
135  if( ninfeasibilities > sparsityThresholdEnterCo )
136  {
137  MSG_INFO2( spxout << "IENTER04 too many infeasibilities for sparse pricing"
138  << std::endl; )
139  remainingRoundsEnterCo = DENSEROUNDS;
140  sparsePricingEnterCo = false;
141  ninfeasibilities = 0;
142  }
143  }
144  }
145  }
146  if( ninfeasibilities == 0 && !sparsePricingEnterCo )
147  --remainingRoundsEnterCo;
148  else if( ninfeasibilities <= sparsityThresholdEnterCo && !sparsePricingEnterCo )
149  {
150  std::streamsize prec = spxout.precision();
151  MSG_INFO2( spxout << "IENTER03 sparse pricing active, "
152  << "sparsity: "
153  << std::setw(6) << std::fixed << std::setprecision(4)
154  << (Real) ninfeasibilities/coDim()
155  << std::scientific << std::setprecision(int(prec))
156  << std::endl; )
157  sparsePricingEnterCo = true;
158  }
159 }
160 
161 Real SPxSolver::computePvec(int i)
162 {
163  METHOD( "SPxSolver::computePvec()" );
164 
165  return (*thePvec)[i] = vector(i) * (*theCoPvec);
166 }
167 
168 Real SPxSolver::computeTest(int i)
169 {
170  METHOD( "SPxSolver::computeTest()" );
171  SPxBasis::Desc::Status stat = desc().status(i);
172  if (isBasic(stat))
173  return theTest[i] = 0;
174  else
175  return theTest[i] = test(i, stat);
176 }
177 
178 /*
179  Testing all values of #coPvec# against its bounds. If $i$, say, is violated
180  the violation is saved as negative value in |theCoTest[i]|.
181  */
182 Real SPxSolver::coTest(int i, SPxBasis::Desc::Status stat) const
183 {
184  METHOD( "SPxSolver::coTest()" );
185  assert(type() == ENTER);
186  assert(!isBasic(stat));
187 
188  Real x;
189 
190  switch (stat)
191  {
192  case SPxBasis::Desc::D_FREE:
193  case SPxBasis::Desc::D_ON_BOTH :
194  assert(rep() == ROW);
195  x = (*theCoPvec)[i] - SPxLP::lower(i);
196  if (x < 0)
197  return x;
198  // no break: next is else case
199  //lint -fallthrough
200  case SPxBasis::Desc::D_ON_LOWER:
201  assert(rep() == ROW);
202  return SPxLP::upper(i) - (*theCoPvec)[i];
203  case SPxBasis::Desc::D_ON_UPPER:
204  assert(rep() == ROW);
205  return (*theCoPvec)[i] - SPxLP::lower(i);
206 
207  case SPxBasis::Desc::P_ON_UPPER:
208  assert(rep() == COLUMN);
209  return (*theCoPvec)[i] - 0; // slacks !
210  case SPxBasis::Desc::P_ON_LOWER:
211  assert(rep() == COLUMN);
212  return 0 - (*theCoPvec)[i]; // slacks !
213 
214  default:
215  return 0;
216  }
217 }
218 
219 void SPxSolver::computeCoTest()
220 {
221  METHOD( "SPxSolver::computeCoTest()" );
222  int i;
223  Real pricingTol = leavetol();
224  infeasibilities.clear();
225  int ninfeasibilities = 0;
226  const SPxBasis::Desc& ds = desc();
227 
228  for (i = dim() - 1; i >= 0; --i)
229  {
230  SPxBasis::Desc::Status stat = ds.coStatus(i);
231  if (isBasic(stat))
232  theCoTest[i] = 0;
233  else
234  {
235  theCoTest[i] = coTest(i, stat);
236  if( remainingRoundsEnter == 0 )
237  {
238  if( theCoTest[i] < -pricingTol )
239  {
240  assert(infeasibilities.size() < infeasibilities.max());
241  infeasibilities.addIdx(i);
242  isInfeasible[i] = true;
243  ++ninfeasibilities;
244  }
245  else
246  isInfeasible[i] = false;
247  if( ninfeasibilities > sparsityThresholdEnter )
248  {
249  MSG_INFO2( spxout << "IENTER06 too many infeasibilities for sparse pricing"
250  << std::endl; )
251  remainingRoundsEnter = DENSEROUNDS;
252  sparsePricingEnter = false;
253  ninfeasibilities = 0;
254  }
255  }
256  }
257  }
258  if( ninfeasibilities == 0 && !sparsePricingEnter )
259  --remainingRoundsEnter;
260  else if( ninfeasibilities <= sparsityThresholdEnter && !sparsePricingEnter )
261  {
262  MSG_INFO2( spxout << "IENTER05 sparse pricing active, "
263  << "sparsity: "
264  << std::setw(6) << std::fixed << std::setprecision(4)
265  << (Real) ninfeasibilities/dim()
266  << std::endl; )
267  sparsePricingEnter = true;
268  }
269 }
270 
271 
272 /*
273  The following methods require propersy initialized vectors |fVec| and
274  #coPvec#.
275  */
276 void SPxSolver::updateTest()
277 {
278  METHOD( "SPxSolver::updateTest()" );
279  thePvec->delta().setup();
280 
281  const IdxSet& idx = thePvec->idx();
282  const SPxBasis::Desc& ds = desc();
283  Real pricingTol = epsilon();
284 
285  int i;
286  for (i = idx.size() - 1; i >= 0; --i)
287  {
288  int j = idx.index(i);
289  SPxBasis::Desc::Status stat = ds.status(j);
290  if (!isBasic(stat))
291  {
292  theTest[j] = test(j, stat);
293 
294  if( sparsePricingEnterCo && theTest[j] < -pricingTol )
295  {
296  assert(remainingRoundsEnterCo == 0);
297  if( !isInfeasibleCo[j] )
298  {
299  infeasibilitiesCo.addIdx(j);
300  isInfeasibleCo[j] = true;
301  }
302  }
303  }
304  else
305  theTest[j] = 0;
306  }
307 }
308 
309 void SPxSolver::updateCoTest()
310 {
311  METHOD( "SPxSolver::updateCoTest()" );
312  theCoPvec->delta().setup();
313 
314  const IdxSet& idx = theCoPvec->idx();
315  const SPxBasis::Desc& ds = desc();
316  Real pricingTol = epsilon();
317 
318  int i;
319  for (i = idx.size() - 1; i >= 0; --i)
320  {
321  int j = idx.index(i);
322  SPxBasis::Desc::Status stat = ds.coStatus(j);
323  if (!isBasic(stat))
324  {
325  theCoTest[j] = coTest(j, stat);
326 
327  if( sparsePricingEnter && theCoTest[j] < -pricingTol )
328  {
329  assert(remainingRoundsEnter == 0);
330  if( !isInfeasible[j] )
331  {
332  infeasibilities.addIdx(j);
333  isInfeasible[j] = true;
334  }
335  }
336  }
337  else
338  theCoTest[j] = 0;
339  }
340 }
341 
342 
343 
344 /* \Section{Compute statistics on entering variable}
345  Here is a list of variables relevant when including |Id| to the basis.
346  They are computed by |computeEnterStats()|.
347  */
348 void SPxSolver::getEnterVals
349 (
350  SPxId enterId,
351  Real& enterTest,
352  Real& enterUB,
353  Real& enterLB,
354  Real& enterVal,
355  Real& enterMax,
356  Real& enterPric,
357  SPxBasis::Desc::Status& enterStat,
358  Real& enterRO
359 )
360 {
361  METHOD( "SPxSolver::getEnterVals()" );
362  int enterIdx;
363  SPxBasis::Desc& ds = desc();
364 
365  if (enterId.isSPxColId())
366  {
367  enterIdx = number(SPxColId(enterId));
368  enterStat = ds.colStatus(enterIdx);
369  assert(!isBasic(enterStat));
370 
371  /* For an #Id# to enter the basis we better recompute the Test value.
372  */
373  if (rep() == COLUMN)
374  {
375  computePvec(enterIdx);
376  enterTest = computeTest(enterIdx);
377  theTest[enterIdx] = 0;
378  }
379  else
380  {
381  enterTest = coTest()[enterIdx];
382  theCoTest[enterIdx] = 0;
383  }
384 
385  switch (enterStat)
386  {
387  // primal/columnwise cases:
388  case SPxBasis::Desc::P_ON_UPPER :
389  assert( rep() == COLUMN );
390  enterUB = theUCbound[enterIdx];
391  enterLB = theLCbound[enterIdx];
392  enterVal = enterUB;
393  enterMax = enterLB - enterUB;
394  enterPric = (*thePvec)[enterIdx];
395  enterRO = maxObj(enterIdx);
396  if( enterLB <= -infinity )
397  ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_LOWER;
398  else if( EQ( enterLB, enterUB ) )
399  ds.colStatus(enterIdx) = SPxBasis::Desc::D_FREE;
400  else
401  ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
402  break;
403  case SPxBasis::Desc::P_ON_LOWER :
404  assert( rep() == COLUMN );
405  enterUB = theUCbound[enterIdx];
406  enterLB = theLCbound[enterIdx];
407  enterVal = enterLB;
408  enterMax = enterUB - enterLB;
409  enterPric = (*thePvec)[enterIdx];
410  enterRO = maxObj(enterIdx);
411  if( enterUB >= infinity )
412  ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_UPPER;
413  else if( EQ( enterLB, enterUB ) )
414  ds.colStatus(enterIdx) = SPxBasis::Desc::D_FREE;
415  else
416  ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
417  break;
418  case SPxBasis::Desc::P_FREE :
419  assert( rep() == COLUMN );
420  enterUB = theUCbound[enterIdx];
421  enterLB = theLCbound[enterIdx];
422  enterVal = 0;
423  enterPric = (*thePvec)[enterIdx];
424  enterRO = maxObj(enterIdx);
425  ds.colStatus(enterIdx) = SPxBasis::Desc::D_UNDEFINED;
426  enterMax = (enterRO - enterPric > 0) ? infinity : -infinity;
427  break;
428 
429  // dual/rowwise cases:
430  case SPxBasis::Desc::D_ON_UPPER :
431  assert( rep() == ROW );
432  assert(theUCbound[enterIdx] < infinity);
433  enterUB = theUCbound[enterIdx];
434  enterLB = -infinity;
435  enterMax = -infinity;
436  enterVal = enterUB;
437  enterPric = (*theCoPvec)[enterIdx];
438  enterRO = SPxLP::lower(enterIdx);
439  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
440  break;
441  case SPxBasis::Desc::D_ON_LOWER :
442  assert( rep() == ROW );
443  assert(theLCbound[enterIdx] > -infinity);
444  enterLB = theLCbound[enterIdx];
445  enterUB = infinity;
446  enterMax = infinity;
447  enterVal = enterLB;
448  enterPric = (*theCoPvec)[enterIdx];
449  enterRO = SPxLP::upper(enterIdx);
450  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
451  break;
452  case SPxBasis::Desc::D_FREE:
453  assert( rep() == ROW );
454  assert(SPxLP::lower(enterIdx) == SPxLP::upper(enterIdx));
455  enterUB = infinity;
456  enterLB = -infinity;
457  enterVal = 0;
458  enterRO = SPxLP::upper(enterIdx);
459  enterPric = (*theCoPvec)[enterIdx];
460  if (enterPric > enterRO)
461  enterMax = infinity;
462  else
463  enterMax = -infinity;
464  ds.colStatus(enterIdx) = SPxBasis::Desc::P_FIXED;
465  break;
466  case SPxBasis::Desc::D_ON_BOTH :
467  assert( rep() == ROW );
468  enterPric = (*theCoPvec)[enterIdx];
469  if (enterPric > SPxLP::upper(enterIdx))
470  {
471  enterLB = theLCbound[enterIdx];
472  enterUB = infinity;
473  enterMax = infinity;
474  enterVal = enterLB;
475  enterRO = SPxLP::upper(enterIdx);
476  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
477  }
478  else
479  {
480  enterUB = theUCbound[enterIdx];
481  enterVal = enterUB;
482  enterRO = SPxLP::lower(enterIdx);
483  enterLB = -infinity;
484  enterMax = -infinity;
485  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
486  }
487  break;
488  default:
489  throw SPxInternalCodeException("XENTER01 This should never happen.");
490  }
491  MSG_DEBUG( spxout << "DENTER03 SPxSolver::getEnterVals() : col " << enterIdx
492  << ": " << enterStat
493  << " -> " << ds.colStatus(enterIdx)
494  << std::endl; )
495  }
496 
497  else
498  {
499  assert(enterId.isSPxRowId());
500  enterIdx = number(SPxRowId(enterId));
501  enterStat = ds.rowStatus(enterIdx);
502  assert(!isBasic(enterStat));
503 
504  /* For an #Id# to enter the basis we better recompute the Test value.
505  */
506  if (rep() == ROW)
507  {
508  computePvec(enterIdx);
509  enterTest = computeTest(enterIdx);
510  theTest[enterIdx] = 0;
511  }
512  else
513  {
514  enterTest = coTest()[enterIdx];
515  theCoTest[enterIdx] = 0;
516  }
517 
518  switch (enterStat)
519  {
520  // primal/columnwise cases:
521  case SPxBasis::Desc::P_ON_UPPER :
522  assert( rep() == COLUMN );
523  enterUB = theURbound[enterIdx];
524  enterLB = theLRbound[enterIdx];
525  enterVal = enterLB;
526  enterMax = enterUB - enterLB;
527  enterPric = (*theCoPvec)[enterIdx];
528  enterRO = 0;
529  if( enterUB >= infinity )
530  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_LOWER;
531  else if( EQ( enterLB, enterUB ) )
532  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_FREE;
533  else
534  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
535  break;
536  case SPxBasis::Desc::P_ON_LOWER :
537  assert( rep() == COLUMN );
538  enterUB = theURbound[enterIdx];
539  enterLB = theLRbound[enterIdx];
540  enterVal = enterUB;
541  enterMax = enterLB - enterUB;
542  enterPric = (*theCoPvec)[enterIdx];
543  enterRO = 0;
544  if( enterLB <= -infinity )
545  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_UPPER;
546  else if( EQ( enterLB, enterUB ) )
547  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_FREE;
548  else
549  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
550  break;
551  case SPxBasis::Desc::P_FREE :
552  assert( rep() == COLUMN );
553 #if 1
554  throw SPxInternalCodeException("XENTER02 This should never happen.");
555 #else
556  MSG_ERROR( spxout << "EENTER99 ERROR: not yet debugged!" << std::endl; )
557  enterPric = (*theCoPvec)[enterIdx];
558  enterRO = 0;
559  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_UNDEFINED;
560 #endif
561  break;
562 
563  // dual/rowwise cases:
564  case SPxBasis::Desc::D_ON_UPPER :
565  assert( rep() == ROW );
566  assert(theURbound[enterIdx] < infinity);
567  enterUB = theURbound[enterIdx];
568  enterLB = -infinity;
569  enterVal = enterUB;
570  enterMax = -infinity;
571  enterPric = (*thePvec)[enterIdx];
572  enterRO = lhs(enterIdx);
573  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
574  break;
575  case SPxBasis::Desc::D_ON_LOWER :
576  assert( rep() == ROW );
577  assert(theLRbound[enterIdx] > -infinity);
578  enterLB = theLRbound[enterIdx];
579  enterUB = infinity;
580  enterVal = enterLB;
581  enterMax = infinity;
582  enterPric = (*thePvec)[enterIdx];
583  enterRO = rhs(enterIdx);
584  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
585  break;
586  case SPxBasis::Desc::D_FREE:
587  assert( rep() == ROW );
588  assert(rhs(enterIdx) == lhs(enterIdx));
589  enterUB = infinity;
590  enterLB = -infinity;
591  enterVal = 0;
592  enterPric = (*thePvec)[enterIdx];
593  enterRO = rhs(enterIdx);
594  enterMax = (enterPric > enterRO) ? infinity : -infinity;
595  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_FIXED;
596  break;
597  case SPxBasis::Desc::D_ON_BOTH :
598  assert( rep() == ROW );
599  enterPric = (*thePvec)[enterIdx];
600  if (enterPric > rhs(enterIdx))
601  {
602  enterLB = theLRbound[enterIdx];
603  enterVal = enterLB;
604  enterUB = infinity;
605  enterMax = infinity;
606  enterRO = rhs(enterIdx);
607  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
608  }
609  else
610  {
611  enterUB = theURbound[enterIdx];
612  enterVal = enterUB;
613  enterLB = -infinity;
614  enterMax = -infinity;
615  enterRO = lhs(enterIdx);
616  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
617  }
618  break;
619 
620  default:
621  throw SPxInternalCodeException("XENTER03 This should never happen.");
622  }
623  MSG_DEBUG( spxout << "DENTER05 SPxSolver::getEnterVals() : row "
624  << enterIdx << ": " << enterStat
625  << " -> " << ds.rowStatus(enterIdx)
626  << std::endl; )
627  }
628 }
629 
630 /* process leaving variable
631  */
632 void SPxSolver::getEnterVals2
633 (
634  int leaveIdx,
635  Real enterMax,
636  Real& leavebound
637 )
638 {
639  METHOD( "SPxSolver::getEnterVals2()" );
640  int idx;
641  SPxBasis::Desc& ds = desc();
642  SPxId leftId = baseId(leaveIdx);
643 
644  if (leftId.isSPxRowId())
645  {
646  idx = number(SPxRowId(leftId));
647  SPxBasis::Desc::Status leaveStat = ds.rowStatus(idx);
648 
649  switch (leaveStat)
650  {
651  case SPxBasis::Desc::P_FIXED :
652  assert(rep() == ROW);
653  throw SPxInternalCodeException("XENTER04 This should never happen.");
654  break;
655  case SPxBasis::Desc::P_ON_UPPER :
656  assert(rep() == ROW);
657  leavebound = theLBbound[leaveIdx];
658  theLRbound[idx] = leavebound;
659  ds.rowStatus(idx) = dualRowStatus(idx);
660  break;
661  case SPxBasis::Desc::P_ON_LOWER :
662  assert(rep() == ROW);
663  leavebound = theUBbound[leaveIdx];
664  theURbound[idx] = leavebound;
665  ds.rowStatus(idx) = dualRowStatus(idx);
666  break;
667  case SPxBasis::Desc::P_FREE :
668  assert(rep() == ROW);
669 #if 1
670  throw SPxInternalCodeException("XENTER05 This should never happen.");
671 #else
672  MSG_ERROR( spxout << "EENTER98 ERROR: not yet debugged!" << std::endl; )
673  if ((*theCoPvec)[leaveIdx] - theLBbound[leaveIdx] <
674  theUBbound[leaveIdx] - (*theCoPvec)[leaveIdx])
675  {
676  leavebound = theLBbound[leaveIdx];
677  theLRbound[idx] = leavebound;
678  }
679  else
680  {
681  leavebound = theUBbound[leaveIdx];
682  theURbound[idx] = leavebound;
683  }
684  ds.rowStatus(idx) = SPxBasis::Desc::D_UNDEFINED;
685 #endif
686  break;
687  // primal/columnwise cases:
688  case SPxBasis::Desc::D_UNDEFINED :
689  assert(rep() == COLUMN);
690  throw SPxInternalCodeException("XENTER06 This should never happen.");
691  break;
692  case SPxBasis::Desc::D_FREE :
693  assert(rep() == COLUMN);
694  if (theFvec->delta()[leaveIdx] * enterMax < 0)
695  leavebound = theUBbound[leaveIdx];
696  else
697  leavebound = theLBbound[leaveIdx];
698  theLRbound[idx] = leavebound;
699  theURbound[idx] = leavebound;
700  ds.rowStatus(idx) = SPxBasis::Desc::P_FIXED;
701  break;
702  case SPxBasis::Desc::D_ON_UPPER :
703  assert(rep() == COLUMN);
704  leavebound = theUBbound[leaveIdx];
705  theURbound[idx] = leavebound;
706  ds.rowStatus(idx) = SPxBasis::Desc::P_ON_LOWER;
707  break;
708  case SPxBasis::Desc::D_ON_LOWER :
709  assert(rep() == COLUMN);
710  leavebound = theLBbound[leaveIdx];
711  theLRbound[idx] = leavebound;
712  ds.rowStatus(idx) = SPxBasis::Desc::P_ON_UPPER;
713  break;
714  case SPxBasis::Desc::D_ON_BOTH :
715  assert(rep() == COLUMN);
716  if (enterMax * theFvec->delta()[leaveIdx] < 0)
717  {
718  leavebound = theUBbound[leaveIdx];
719  theURbound[idx] = leavebound;
720  ds.rowStatus(idx) = SPxBasis::Desc::P_ON_LOWER;
721  }
722  else
723  {
724  leavebound = theLBbound[leaveIdx];
725  theLRbound[idx] = leavebound;
726  ds.rowStatus(idx) = SPxBasis::Desc::P_ON_UPPER;
727  }
728  break;
729 
730  default:
731  throw SPxInternalCodeException("XENTER07 This should never happen.");
732  }
733  MSG_DEBUG( spxout << "DENTER06 SPxSolver::getEnterVals2(): row "
734  << idx << ": " << leaveStat
735  << " -> " << ds.rowStatus(idx)
736  << std::endl; )
737  }
738 
739  else
740  {
741  assert(leftId.isSPxColId());
742  idx = number(SPxColId(leftId));
743  SPxBasis::Desc::Status leaveStat = ds.colStatus(idx);
744 
745  switch (leaveStat)
746  {
747  case SPxBasis::Desc::P_ON_UPPER :
748  assert(rep() == ROW);
749  leavebound = theLBbound[leaveIdx];
750  theLCbound[idx] = leavebound;
751  ds.colStatus(idx) = dualColStatus(idx);
752  break;
753  case SPxBasis::Desc::P_ON_LOWER :
754  assert(rep() == ROW);
755  leavebound = theUBbound[leaveIdx];
756  theUCbound[idx] = leavebound;
757  ds.colStatus(idx) = dualColStatus(idx);
758  break;
759  case SPxBasis::Desc::P_FREE :
760  assert(rep() == ROW);
761  if (theFvec->delta()[leaveIdx] * enterMax > 0)
762  {
763  leavebound = theLBbound[leaveIdx];
764  theLCbound[idx] = leavebound;
765  }
766  else
767  {
768  leavebound = theUBbound[leaveIdx];
769  theUCbound[idx] = leavebound;
770  }
771  ds.colStatus(idx) = SPxBasis::Desc::D_UNDEFINED;
772  break;
773  case SPxBasis::Desc::P_FIXED:
774  assert(rep() == ROW);
775  throw SPxInternalCodeException("XENTER08 This should never happen.");
776  break;
777  // primal/columnwise cases:
778  case SPxBasis::Desc::D_FREE :
779  assert(rep() == COLUMN);
780  if (theFvec->delta()[leaveIdx] * enterMax > 0)
781  leavebound = theLBbound[leaveIdx];
782  else
783  leavebound = theUBbound[leaveIdx];
784  theUCbound[idx] =
785  theLCbound[idx] = leavebound;
786  ds.colStatus(idx) = SPxBasis::Desc::P_FIXED;
787  break;
788  case SPxBasis::Desc::D_ON_UPPER :
789  assert(rep() == COLUMN);
790  leavebound = theLBbound[leaveIdx];
791  theLCbound[idx] = leavebound;
792  ds.colStatus(idx) = SPxBasis::Desc::P_ON_LOWER;
793  break;
794  case SPxBasis::Desc::D_ON_LOWER :
795  assert(rep() == COLUMN);
796  leavebound = theUBbound[leaveIdx];
797  theUCbound[idx] = leavebound;
798  ds.colStatus(idx) = SPxBasis::Desc::P_ON_UPPER;
799  break;
800  case SPxBasis::Desc::D_ON_BOTH :
801  case SPxBasis::Desc::D_UNDEFINED :
802  assert(rep() == COLUMN);
803  if (enterMax * theFvec->delta()[leaveIdx] < 0)
804  {
805  leavebound = theUBbound[leaveIdx];
806  theUCbound[idx] = leavebound;
807  ds.colStatus(idx) = SPxBasis::Desc::P_ON_UPPER;
808  }
809  else
810  {
811  leavebound = theLBbound[leaveIdx];
812  theLCbound[idx] = leavebound;
813  ds.colStatus(idx) = SPxBasis::Desc::P_ON_LOWER;
814  }
815  break;
816 
817  default:
818  throw SPxInternalCodeException("XENTER09 This should never happen.");
819  }
820  MSG_DEBUG( spxout << "DENTER07 SPxSolver::getEnterVals2(): col "
821  << idx << ": " << leaveStat
822  << " -> " << ds.colStatus(idx)
823  << std::endl; )
824  }
825 }
826 
827 
828 void
829 SPxSolver::ungetEnterVal(
830  SPxId enterId,
831  SPxBasis::Desc::Status enterStat,
832  Real leaveVal,
833  const SVector& vec
834 )
835 {
836  METHOD( "SPxSolver::ungetEnterVal()" );
837  int enterIdx;
838  SPxBasis::Desc& ds = desc();
839 
840  if (enterId.isSPxColId())
841  {
842  enterIdx = number(SPxColId(enterId));
843  if (enterStat == SPxBasis::Desc::P_ON_UPPER)
844  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
845  else
846  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
847  theFrhs->multAdd(leaveVal, vec);
848  }
849  else
850  {
851  enterIdx = number(SPxRowId(enterId));
852  assert(enterId.isSPxRowId());
853  if (enterStat == SPxBasis::Desc::P_ON_UPPER)
854  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
855  else
856  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
857  (*theFrhs)[enterIdx] += leaveVal;
858  }
859  if (isId(enterId))
860  theTest[enterIdx] = 0;
861  else
862  theCoTest[enterIdx] = 0;
863 }
864 
865 void SPxSolver::rejectEnter(
866  SPxId enterId,
867  Real enterTest,
868  SPxBasis::Desc::Status enterStat
869 )
870 {
871  METHOD( "SPxSolver::rejectEnter()" );
872  int enterIdx = number(enterId);
873  if (isId(enterId))
874  {
875  theTest[enterIdx] = enterTest;
876  desc().status(enterIdx) = enterStat;
877  }
878  else
879  {
880  theCoTest[enterIdx] = enterTest;
881  desc().coStatus(enterIdx) = enterStat;
882  }
883 }
884 
885 bool SPxSolver::enter(SPxId& enterId)
886 {
887  METHOD( "SPxSolver::enter()" );
888  assert(enterId.isValid());
889  assert(type() == ENTER);
890  assert(initialized);
891 
892  SPxId none; // invalid id used when enter fails
893  Real enterTest; // correct test value of entering var
894  Real enterUB; // upper bound of entering variable
895  Real enterLB; // lower bound of entering variable
896  Real enterVal; // current value of entering variable
897  Real enterMax; // maximum value for entering shift
898  Real enterPric; // priced value of entering variable
899  SPxBasis::Desc::Status enterStat; // status of entering variable
900  Real enterRO; // rhs/obj of entering variable
901  const SVector* enterVec = enterVector(enterId);
902 
903  getEnterVals(enterId, enterTest, enterUB, enterLB,
904  enterVal, enterMax, enterPric, enterStat, enterRO);
905 
906  if (enterTest > -epsilon())
907  {
908  rejectEnter(enterId, enterTest, enterStat);
909  change(-1, none, 0);
910 
911  MSG_DEBUG( spxout << "DENTER08 rejecting false enter pivot" << std::endl; )
912 
913  return true;
914  }
915 
916  /* Before performing the actual basis update, we must determine, how this
917  is to be accomplished.
918  */
919  // BH 2005-11-15: Obviously solve4update() is only called if theFvec.delta()
920  // is setup (i.e. the indices of the NZEs are stored within it) and there are
921  // 0 NZEs (???).
922  // In that case theFvec->delta() is set such that
923  // Base * theFvec->delta() = enterVec
924  if (theFvec->delta().isSetup() && theFvec->delta().size() == 0)
925  SPxBasis::solve4update(theFvec->delta(), *enterVec);
926 #ifdef ENABLE_ADDITIONAL_CHECKS
927  else
928  {
929  // BH 2005-11-29: This code block seems to check the assertion
930  // || Base * theFvec->delta() - enterVec ||_2 <= entertol()
931  DVector tmp(dim());
932  // BH 2005-11-15: This cast is necessary since SSVector inherits protected from DVector.
933  tmp = reinterpret_cast<DVector&>(theFvec->delta());
934  multBaseWith(tmp);
935  tmp -= *enterVec;
936  if (tmp.length() > entertol()) {
937  // This happens frequently and does usually not hurt, so print these
938  // warnings only with verbose level INFO2 and higher.
939  MSG_INFO2( spxout << "WENTER09 fVec updated error = "
940  << tmp.length() << std::endl; )
941  }
942  }
943 #endif // ENABLE_ADDITIONAL_CHECKS
944 
945  if (m_numCycle > m_maxCycle)
946  {
947  if (-enterMax > 0)
948  perturbMaxEnter();
949  else
950  perturbMinEnter();
951  }
952 
953  Real leaveVal = -enterMax;
954 
955  int leaveIdx = theratiotester->selectLeave(leaveVal, enterId);
956 
957  /* in row representation, fixed columns and rows should not leave the basis */
958  assert(leaveIdx < 0 || !baseId(leaveIdx).isSPxColId() || desc().colStatus(number(SPxColId(baseId(leaveIdx)))) != SPxBasis::Desc::P_FIXED);
959  assert(leaveIdx < 0 || !baseId(leaveIdx).isSPxRowId() || desc().rowStatus(number(SPxRowId(baseId(leaveIdx)))) != SPxBasis::Desc::P_FIXED);
960 
961  /*
962  We now tried to find a variable to leave the basis. If one has been
963  found, a regular basis update is to be performed.
964  */
965  if (leaveIdx >= 0)
966  {
967  if (fabs(leaveVal) < entertol())
968  {
969  if (theUBbound[leaveIdx] != theLBbound[leaveIdx]
970  && enterStat != Desc::P_FREE && enterStat != Desc::D_FREE)
971  m_numCycle++;
972  }
973  else
974  m_numCycle /= 2;
975 
976  // setup for updating the copricing vector
977  if (coSolveVector2)
978  SPxBasis::coSolve(theCoPvec->delta(), *coSolveVector2, unitVecs[leaveIdx], *coSolveVector2rhs);
979  else
980  SPxBasis::coSolve(theCoPvec->delta(), unitVecs[leaveIdx]);
981 
982  (*theCoPrhs)[leaveIdx] = enterRO;
983  theCoPvec->value() = (enterRO - enterPric) / theFvec->delta()[leaveIdx];
984 
985  if (theCoPvec->value() > epsilon() || theCoPvec->value() < -epsilon())
986  {
987  if (pricing() == FULL)
988  {
989  thePvec->value() = theCoPvec->value();
990  setupPupdate();
991  }
992  doPupdate();
993  }
994  assert(thePvec->isConsistent());
995  assert(theCoPvec->isConsistent());
996 
997  assert(!baseId(leaveIdx).isSPxRowId() || desc().rowStatus(number(SPxRowId(baseId(leaveIdx)))) != SPxBasis::Desc::P_FIXED);
998  assert(!baseId(leaveIdx).isSPxColId() || desc().colStatus(number(SPxColId(baseId(leaveIdx)))) != SPxBasis::Desc::P_FIXED);
999 
1000  Real leavebound; // bound on which leaving variable moves
1001  try
1002  {
1003  getEnterVals2(leaveIdx, enterMax, leavebound);
1004  }
1005  catch( SPxException F )
1006  {
1007  rejectEnter(enterId, enterTest, enterStat);
1008  change(-1, none, 0);
1009  throw F;
1010  }
1011 
1012  // process entering variable
1013  theUBbound[leaveIdx] = enterUB;
1014  theLBbound[leaveIdx] = enterLB;
1015 
1016  // compute tests:
1017  updateCoTest();
1018  if (pricing() == FULL)
1019  updateTest();
1020 
1021 
1022  // update feasibility vectors
1023  theFvec->value() = leaveVal;
1024  theFvec->update();
1025  (*theFvec)[leaveIdx] = enterVal - leaveVal;
1026 
1027  if (leavebound > epsilon() || leavebound < -epsilon())
1028  theFrhs->multAdd(-leavebound, baseVec(leaveIdx));
1029 
1030  if (enterVal > epsilon() || enterVal < -epsilon())
1031  theFrhs->multAdd(enterVal, *enterVec);
1032 
1033 
1034  // change basis matrix
1035  change(leaveIdx, enterId, enterVec, &(theFvec->delta()));
1036  }
1037  /* No leaving vector could be found that would yield a stable pivot step.
1038  */
1039  else if (leaveVal != -enterMax)
1040  {
1041  rejectEnter(enterId, REAL(0.01) * enterTest - REAL(2.0) * leavetol(), enterStat);
1042  change(-1, none, 0);
1043  }
1044  /* No leaving vector has been selected from the basis. However, if the
1045  shift amount for |fVec| is bounded, we are in the case, that the
1046  entering variable is moved from one bound to its other, before any of
1047  the basis feasibility variables reaches their bound. This may only
1048  happen in primal/columnwise case with upper and lower bounds on
1049  variables.
1050  */
1051  else if (leaveVal < infinity && leaveVal > -infinity)
1052  {
1053  assert(rep() == COLUMN);
1054  assert(leaveVal == -enterMax);
1055 
1056  change(-1, enterId, enterVec);
1057 
1058  theFvec->value() = leaveVal;
1059  theFvec->update();
1060 
1061  ungetEnterVal(enterId, enterStat, leaveVal, *enterVec);
1062  }
1063  /* No variable could be selected to leave the basis and even the entering
1064  variable is unbounded --- this is a failure.
1065  */
1066  else
1067  {
1068  /* The following line originally was in the "lastUpdate() > 1" case;
1069  we need it in the INFEASIBLE/UNBOUNDED case, too, to have the
1070  basis descriptor at the correct size.
1071  */
1072  rejectEnter(enterId, enterTest, enterStat);
1073  change(-1, none, 0);
1074 
1075  if (lastUpdate() > 1)
1076  {
1077  MSG_INFO3( spxout << "IENTER01 factorization triggered in "
1078  << "enter() for feasibility test" << std::endl; )
1079  factorize();
1080 
1081  /* after a factorization, the entering column/row might not be infeasible or suboptimal anymore, hence we do
1082  * not try to call leave(leaveIdx), but rather return to the main solving loop and call the pricer again
1083  */
1084  return true;
1085  }
1086 
1087  MSG_INFO3( spxout << "IENTER02 unboundness/infeasiblity found in "
1088  << "enter()" << std::endl; )
1089 
1090  if (rep() == ROW)
1091  {
1092  Real sign;
1093 
1094  dualFarkas.clear();
1095  dualFarkas.setMax(fVec().delta().size() + 1);
1096  sign = (leaveVal > 0 ? -1.0 : 1.0);
1097 
1098  for( int j = 0; j < fVec().delta().size(); ++j )
1099  {
1100  SPxId id = baseId(fVec().idx().index(j));
1101 
1102  if( id.isSPxRowId() )
1103  dualFarkas.add(number(SPxRowId(id)), sign * fVec().delta().value(j));
1104  }
1105 
1106  if( enterId.isSPxRowId() )
1107  dualFarkas.add(number(SPxRowId(enterId)), -sign);
1108 
1109  setBasisStatus(SPxBasis::INFEASIBLE);
1110  }
1111  /**@todo if shift() is not zero, we must not conclude primal unboundedness */
1112  else
1113  {
1114  Real sign;
1115 
1116  primalRay.clear();
1117  primalRay.setMax(fVec().delta().size() + 1);
1118  sign = leaveVal > 0 ? 1.0 : -1.0;
1119 
1120  for( int j = 0; j < fVec().delta().size(); ++j )
1121  {
1122  SPxId i = baseId(fVec().idx().index(j));
1123 
1124  if( i.isSPxColId() )
1125  primalRay.add(number(SPxColId(i)), sign*fVec().delta().value(j));
1126  }
1127 
1128  if( enterId.isSPxColId() )
1129  primalRay.add(number(SPxColId(enterId)), -sign);
1130 
1131  setBasisStatus(SPxBasis::UNBOUNDED);
1132  }
1133  return false;
1134  }
1135  return true;
1136 }
1137 } // namespace soplex
1138 
1139 //-----------------------------------------------------------------------------
1140 //Emacs Local Variables:
1141 //Emacs mode:c++
1142 //Emacs c-basic-offset:3
1143 //Emacs tab-width:8
1144 //Emacs indent-tabs-mode:nil
1145 //Emacs End:
1146 //-----------------------------------------------------------------------------
© 2003-2013 by Zuse Institute Berlin (ZIB), Imprint Generated on Wed Jan 9 2013 for SoPlex by doxygen