Scippy

SoPlex

Sequential object-oriented simPlex

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-2015 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 /* \SubSection{Updating the Basis for Entering Variables}
17  */
18 #include <assert.h>
19 
20 #include "spxdefines.h"
21 #include "spxratiotester.h"
22 #include "spxpricer.h"
23 #include "spxout.h"
24 #include "exceptions.h"
25 
26 namespace soplex
27 {
28 
29 /*
30 In the entering simplex algorithms (i.e. iteratively a vector is selected to
31 \em enter the simplex basis as in the dual rowwise and primal columnwise case)
32 let $A$ denote the current basis, $x$ and entering vector and $f$ the
33 feasibility vector. For a feasible basis $l \le f \le u$ holds.
34 For the rowwise case $f$ is obtained by solving $f^T = c^T A^{-1}$,
35 wherease in columnwisecase $f = A^{-1} b$.
36 
37 Let us further consider the rowwise case. Exchanging $x$ with the $i$-th
38 vector of $A$ yields
39 
40 \begin{equation}\label{update.eq}
41  A^{(i)} = E_i A \hbox{, with } E_i = I + e_i (x^T A^{-1} - e_i^T).
42 \end{equation}
43 
44 With $E_i^{-1} = I + e_i \frac{e_i^T - \delta^T}{\delta}$,
45 $\delta^T = x^T A^{-1}$ one gets the new feasibility vector
46 
47 \begin{eqnarray*}
48  (f^{(i)})^T
49  &=& c^T (A^{(i)})^{-1} \\
50  &=& c^T A^{-1} + c^T A^{-1} e_i \frac{e_i^T - \delta^T}{\delta_i} \\
51  &=& f^T + \frac{f_i}{\delta_i} e_i^T - \frac{f_i}{\delta_i} \delta^T. \\
52 \end{eqnarray*}
53 
54 The selection of the leaving vector $i^*$ for the basis must ensure, that for
55 all $j \ne i^*$ $f^{(i^*)}_j$ remains within its bounds $l_j$ and $u_j$.
56  */
57 
58 
59 /*
60  Testing all values of |pVec| against its bounds. If $i$, say, is violated
61  the violation is saved as negative value in |theTest[i]|.
62  */
64 {
65  assert(type() == ENTER);
66  assert(!isBasic(stat));
67 
68  Real x;
69 
70  switch (stat)
71  {
74  assert(rep() == ROW);
75  x = (*thePvec)[i] - lhs(i);
76  if (x < 0)
77  return x;
78  // no break: next is else case
79  //lint -fallthrough
81  assert(rep() == ROW);
82  return rhs(i) - (*thePvec)[i];
84  assert(rep() == ROW);
85  return (*thePvec)[i] - lhs(i);
86 
88  assert(rep() == COLUMN);
89  return maxObj(i) - (*thePvec)[i];
91  assert(rep() == COLUMN);
92  return (*thePvec)[i] - maxObj(i);
94  x = maxObj(i) - (*thePvec)[i];
95  return (x < 0) ? x : -x;
96 
97  default:
98  return 0;
99  }
100 }
101 
103 {
104 
105  const SPxBasis::Desc& ds = desc();
106  Real pricingTol = leavetol();
108  int ninfeasibilities = 0;
109  int sparsitythreshold = (int) (sparsePricingFactor * coDim());
110 
111  for(int i = 0; i < coDim(); ++i)
112  {
113  SPxBasis::Desc::Status stat = ds.status(i);
114 
115  if(isBasic(stat))
116  {
117  theTest[i] = 0.0;
118  if( remainingRoundsEnterCo == 0 )
120  }
121  else
122  {
123  assert(!isBasic(stat));
124  theTest[i] = test(i, stat);
125 
126  if( remainingRoundsEnterCo == 0 )
127  {
128  if( theTest[i] < -pricingTol )
129  {
133  ++ninfeasibilities;
134  }
135  else
137  if( ninfeasibilities > sparsitythreshold)
138  {
139  MSG_INFO2( (*spxout), (*spxout) << " --- using dense pricing"
140  << std::endl; )
142  sparsePricingEnterCo = false;
143  ninfeasibilities = 0;
144  }
145  }
146  }
147  }
148  if( ninfeasibilities == 0 && !sparsePricingEnterCo )
150  else if( ninfeasibilities <= sparsitythreshold && !sparsePricingEnterCo )
151  {
152  MSG_INFO2( (*spxout),
153  std::streamsize prec = spxout->precision();
154  if( hyperPricingEnter )
155  (*spxout) << " --- using hypersparse pricing, ";
156  else
157  (*spxout) << " --- using sparse pricing, ";
158  (*spxout) << "sparsity: "
159  << std::setw(6) << std::fixed << std::setprecision(4)
160  << (Real) ninfeasibilities/coDim()
161  << std::scientific << std::setprecision(int(prec))
162  << std::endl;
163  )
164  sparsePricingEnterCo = true;
165  }
166 }
167 
169 {
170 
171  return (*thePvec)[i] = vector(i) * (*theCoPvec);
172 }
173 
175 {
176  SPxBasis::Desc::Status stat = desc().status(i);
177  if (isBasic(stat))
178  return theTest[i] = 0;
179  else
180  return theTest[i] = test(i, stat);
181 }
182 
183 /*
184  Testing all values of #coPvec# against its bounds. If $i$, say, is violated
185  the violation is saved as negative value in |theCoTest[i]|.
186  */
188 {
189  assert(type() == ENTER);
190  assert(!isBasic(stat));
191 
192  Real x;
193 
194  switch (stat)
195  {
198  assert(rep() == ROW);
199  x = (*theCoPvec)[i] - SPxLP::lower(i);
200  if (x < 0)
201  return x;
202  // no break: next is else case
203  //lint -fallthrough
205  assert(rep() == ROW);
206  return SPxLP::upper(i) - (*theCoPvec)[i];
208  assert(rep() == ROW);
209  return (*theCoPvec)[i] - SPxLP::lower(i);
210 
212  assert(rep() == COLUMN);
213  return (*theCoPvec)[i] - maxRowObj(i); // slacks !
215  assert(rep() == COLUMN);
216  return maxRowObj(i) - (*theCoPvec)[i]; // slacks !
217 
218  default:
219  return 0;
220  }
221 }
222 
224 {
225  int i;
226  Real pricingTol = leavetol();
228  int ninfeasibilities = 0;
229  int sparsitythreshold = (int) (sparsePricingFactor * dim());
230  const SPxBasis::Desc& ds = desc();
231 
232  for (i = dim() - 1; i >= 0; --i)
233  {
234  SPxBasis::Desc::Status stat = ds.coStatus(i);
235  if (isBasic(stat))
236  {
237  theCoTest[i] = 0;
238  if( remainingRoundsEnter == 0 )
240  }
241  else
242  {
243  theCoTest[i] = coTest(i, stat);
244  if( remainingRoundsEnter == 0 )
245  {
246  if( theCoTest[i] < -pricingTol )
247  {
248  assert(infeasibilities.size() < infeasibilities.max());
251  ++ninfeasibilities;
252  }
253  else
255  if( ninfeasibilities > sparsitythreshold )
256  {
257  MSG_INFO2( (*spxout), (*spxout) << " --- using dense pricing"
258  << std::endl; )
260  sparsePricingEnter = false;
261  ninfeasibilities = 0;
262  }
263  }
264  }
265  }
266  if( ninfeasibilities == 0 && !sparsePricingEnter )
268  else if( ninfeasibilities <= sparsitythreshold && !sparsePricingEnter )
269  {
270  MSG_INFO2( (*spxout),
271  std::streamsize prec = spxout->precision();
272  if( hyperPricingEnter )
273  (*spxout) << " --- using hypersparse pricing, ";
274  else
275  (*spxout) << " --- using sparse pricing, ";
276  (*spxout) << "sparsity: "
277  << std::setw(6) << std::fixed << std::setprecision(4)
278  << (Real) ninfeasibilities/dim()
279  << std::scientific << std::setprecision(int(prec))
280  << std::endl;
281  )
282  sparsePricingEnter = true;
283  }
284 }
285 
286 
287 /*
288  The following methods require propersy initialized vectors |fVec| and
289  #coPvec#.
290  */
292 {
293  thePvec->delta().setup();
294 
295  const IdxSet& idx = thePvec->idx();
296  const SPxBasis::Desc& ds = desc();
297  Real pricingTol = leavetol();
298 
299  int i;
301  for (i = idx.size() - 1; i >= 0; --i)
302  {
303  int j = idx.index(i);
304  SPxBasis::Desc::Status stat = ds.status(j);
305  if (!isBasic(stat))
306  {
307  theTest[j] = test(j, stat);
308 
310  {
311  if( theTest[j] < -pricingTol )
312  {
313  assert(remainingRoundsEnterCo == 0);
315  {
318  }
319  if( hyperPricingEnter )
321  }
322  else
323  {
325  }
326  }
327  }
328  else
329  {
331  theTest[j] = 0;
332  }
333  }
334 }
335 
337 {
338  theCoPvec->delta().setup();
339 
340  const IdxSet& idx = theCoPvec->idx();
341  const SPxBasis::Desc& ds = desc();
342  Real pricingTol = leavetol();
343 
344  int i;
345  updateViols.clear();
346  for (i = idx.size() - 1; i >= 0; --i)
347  {
348  int j = idx.index(i);
349  SPxBasis::Desc::Status stat = ds.coStatus(j);
350  if (!isBasic(stat))
351  {
352  theCoTest[j] = coTest(j, stat);
353 
354  if( sparsePricingEnter )
355  {
356  if( theCoTest[j] < -pricingTol )
357  {
358  assert(remainingRoundsEnter == 0);
360  {
361  // if( !hyperPricingEnter )
364  }
365  if( hyperPricingEnter )
366  updateViols.addIdx(j);
367  }
368  else
369  {
370  // @todo do we need to remove index j from infeasibilitiesCo?
372  }
373  }
374  }
375  else
376  {
378  theCoTest[j] = 0;
379  }
380  }
381 }
382 
383 
384 
385 /* \Section{Compute statistics on entering variable}
386  Here is a list of variables relevant when including |Id| to the basis.
387  They are computed by |computeEnterStats()|.
388  */
390 (
391  SPxId enterId,
392  Real& enterTest,
393  Real& enterUB,
394  Real& enterLB,
395  Real& enterVal,
396  Real& enterMax,
397  Real& enterPric,
398  SPxBasis::Desc::Status& enterStat,
399  Real& enterRO,
400  Real& objChange
401 )
402 {
403  int enterIdx;
404  SPxBasis::Desc& ds = desc();
405 
406  if (enterId.isSPxColId())
407  {
408  enterIdx = number(SPxColId(enterId));
409  enterStat = ds.colStatus(enterIdx);
410  assert(!isBasic(enterStat));
411 
412  /* For an #Id# to enter the basis we better recompute the Test value.
413  */
414  if (rep() == COLUMN)
415  {
416  computePvec(enterIdx);
417  enterTest = computeTest(enterIdx);
418  theTest[enterIdx] = 0;
419  }
420  else
421  {
422  enterTest = coTest()[enterIdx];
423  theCoTest[enterIdx] = 0;
424  }
425 
426  switch (enterStat)
427  {
428  // primal/columnwise cases:
430  assert( rep() == COLUMN );
431  enterUB = theUCbound[enterIdx];
432  enterLB = theLCbound[enterIdx];
433  enterVal = enterUB;
434  enterMax = enterLB - enterUB;
435  enterPric = (*thePvec)[enterIdx];
436  enterRO = maxObj(enterIdx);
437  objChange -= enterVal * enterRO;
438  if( enterLB <= -infinity )
439  ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_LOWER;
440  else if( EQ( enterLB, enterUB ) )
441  ds.colStatus(enterIdx) = SPxBasis::Desc::D_FREE;
442  else
443  ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
444  break;
446  assert( rep() == COLUMN );
447  enterUB = theUCbound[enterIdx];
448  enterLB = theLCbound[enterIdx];
449  enterVal = enterLB;
450  enterMax = enterUB - enterLB;
451  enterPric = (*thePvec)[enterIdx];
452  enterRO = maxObj(enterIdx);
453  objChange -= enterVal * enterRO;
454  if( enterUB >= infinity )
455  ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_UPPER;
456  else if( EQ( enterLB, enterUB ) )
457  ds.colStatus(enterIdx) = SPxBasis::Desc::D_FREE;
458  else
459  ds.colStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
460  break;
462  assert( rep() == COLUMN );
463  enterUB = theUCbound[enterIdx];
464  enterLB = theLCbound[enterIdx];
465  enterVal = 0;
466  enterPric = (*thePvec)[enterIdx];
467  enterRO = maxObj(enterIdx);
468  ds.colStatus(enterIdx) = SPxBasis::Desc::D_UNDEFINED;
469  enterMax = (enterRO - enterPric > 0) ? infinity : -infinity;
470  break;
471 
472  // dual/rowwise cases:
474  assert( rep() == ROW );
475  assert(theUCbound[enterIdx] < infinity);
476  enterUB = theUCbound[enterIdx];
477  enterLB = -infinity;
478  enterMax = -infinity;
479  enterVal = enterUB;
480  enterPric = (*theCoPvec)[enterIdx];
481  enterRO = SPxLP::lower(enterIdx);
482  objChange -= enterRO * enterVal;
483  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
484  break;
486  assert( rep() == ROW );
487  assert(theLCbound[enterIdx] > -infinity);
488  enterLB = theLCbound[enterIdx];
489  enterUB = infinity;
490  enterMax = infinity;
491  enterVal = enterLB;
492  enterPric = (*theCoPvec)[enterIdx];
493  enterRO = SPxLP::upper(enterIdx);
494  objChange -= enterRO * enterVal;
495  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
496  break;
498  assert( rep() == ROW );
499  assert(SPxLP::lower(enterIdx) == SPxLP::upper(enterIdx));
500  enterUB = infinity;
501  enterLB = -infinity;
502  enterVal = 0;
503  enterRO = SPxLP::upper(enterIdx);
504  enterPric = (*theCoPvec)[enterIdx];
505  if (enterPric > enterRO)
506  enterMax = infinity;
507  else
508  enterMax = -infinity;
509  ds.colStatus(enterIdx) = SPxBasis::Desc::P_FIXED;
510  break;
512  assert( rep() == ROW );
513  enterPric = (*theCoPvec)[enterIdx];
514  if (enterPric > SPxLP::upper(enterIdx))
515  {
516  enterLB = theLCbound[enterIdx];
517  enterUB = infinity;
518  enterMax = infinity;
519  enterVal = enterLB;
520  enterRO = SPxLP::upper(enterIdx);
521  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
522  }
523  else
524  {
525  enterUB = theUCbound[enterIdx];
526  enterVal = enterUB;
527  enterRO = SPxLP::lower(enterIdx);
528  enterLB = -infinity;
529  enterMax = -infinity;
530  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
531  }
532  objChange -= theLCbound[enterIdx] * SPxLP::upper(enterIdx);
533  objChange -= theUCbound[enterIdx] * SPxLP::lower(enterIdx);
534  break;
535  default:
536  throw SPxInternalCodeException("XENTER01 This should never happen.");
537  }
538  MSG_DEBUG( std::cout << "DENTER03 SPxSolver::getEnterVals() : col " << enterIdx
539  << ": " << enterStat
540  << " -> " << ds.colStatus(enterIdx)
541  << " objChange: " << objChange
542  << std::endl; )
543  }
544 
545  else
546  {
547  assert(enterId.isSPxRowId());
548  enterIdx = number(SPxRowId(enterId));
549  enterStat = ds.rowStatus(enterIdx);
550  assert(!isBasic(enterStat));
551 
552  /* For an #Id# to enter the basis we better recompute the Test value.
553  */
554  if (rep() == ROW)
555  {
556  computePvec(enterIdx);
557  enterTest = computeTest(enterIdx);
558  theTest[enterIdx] = 0;
559  }
560  else
561  {
562  enterTest = coTest()[enterIdx];
563  theCoTest[enterIdx] = 0;
564  }
565 
566  switch (enterStat)
567  {
568  // primal/columnwise cases:
570  assert( rep() == COLUMN );
571  enterUB = theURbound[enterIdx];
572  enterLB = theLRbound[enterIdx];
573  enterVal = enterLB;
574  enterMax = enterUB - enterLB;
575  enterPric = (*theCoPvec)[enterIdx];
576  enterRO = maxRowObj(enterIdx);
577  objChange -= enterRO * enterVal;
578  if( enterUB >= infinity )
579  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_LOWER;
580  else if( EQ( enterLB, enterUB ) )
581  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_FREE;
582  else
583  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
584  break;
586  assert( rep() == COLUMN );
587  enterUB = theURbound[enterIdx];
588  enterLB = theLRbound[enterIdx];
589  enterVal = enterUB;
590  enterMax = enterLB - enterUB;
591  enterPric = (*theCoPvec)[enterIdx];
592  enterRO = maxRowObj(enterIdx);
593  objChange -= enterRO * enterVal;
594  if( enterLB <= -infinity )
595  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_UPPER;
596  else if( EQ( enterLB, enterUB ) )
597  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_FREE;
598  else
599  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_ON_BOTH;
600  break;
602  assert( rep() == COLUMN );
603 #if 1
604  throw SPxInternalCodeException("XENTER02 This should never happen.");
605 #else
606  MSG_ERROR( std::cerr << "EENTER99 ERROR: not yet debugged!" << std::endl; )
607  enterPric = (*theCoPvec)[enterIdx];
608  enterRO = maxRowObj(enterIdx);
609  ds.rowStatus(enterIdx) = SPxBasis::Desc::D_UNDEFINED;
610 #endif
611  break;
612 
613  // dual/rowwise cases:
615  assert( rep() == ROW );
616  assert(theURbound[enterIdx] < infinity);
617  enterUB = theURbound[enterIdx];
618  enterLB = -infinity;
619  enterVal = enterUB;
620  enterMax = -infinity;
621  enterPric = (*thePvec)[enterIdx];
622  enterRO = lhs(enterIdx);
623  objChange -= enterRO * enterVal;
624  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
625  break;
627  assert( rep() == ROW );
628  assert(theLRbound[enterIdx] > -infinity);
629  enterLB = theLRbound[enterIdx];
630  enterUB = infinity;
631  enterVal = enterLB;
632  enterMax = infinity;
633  enterPric = (*thePvec)[enterIdx];
634  enterRO = rhs(enterIdx);
635  objChange -= enterRO * enterVal;
636  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
637  break;
639  assert( rep() == ROW );
640  assert(rhs(enterIdx) == lhs(enterIdx));
641  enterUB = infinity;
642  enterLB = -infinity;
643  enterVal = 0;
644  enterPric = (*thePvec)[enterIdx];
645  enterRO = rhs(enterIdx);
646  enterMax = (enterPric > enterRO) ? infinity : -infinity;
647  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_FIXED;
648  break;
650  assert( rep() == ROW );
651  enterPric = (*thePvec)[enterIdx];
652  if (enterPric > rhs(enterIdx))
653  {
654  enterLB = theLRbound[enterIdx];
655  enterVal = enterLB;
656  enterUB = infinity;
657  enterMax = infinity;
658  enterRO = rhs(enterIdx);
659  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
660  }
661  else
662  {
663  enterUB = theURbound[enterIdx];
664  enterVal = enterUB;
665  enterLB = -infinity;
666  enterMax = -infinity;
667  enterRO = lhs(enterIdx);
668  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
669  }
670  objChange -= theLRbound[enterIdx] * rhs(enterIdx);
671  objChange -= theURbound[enterIdx] * lhs(enterIdx);
672  break;
673 
674  default:
675  throw SPxInternalCodeException("XENTER03 This should never happen.");
676  }
677  MSG_DEBUG( std::cout << "DENTER05 SPxSolver::getEnterVals() : row "
678  << enterIdx << ": " << enterStat
679  << " -> " << ds.rowStatus(enterIdx)
680  << " objChange: " << objChange
681  << std::endl; )
682  }
683 }
684 
685 /* process leaving variable
686  */
688 (
689  int leaveIdx,
690  Real enterMax,
691  Real& leavebound,
692  Real& objChange
693 )
694 {
695  int idx;
696  SPxBasis::Desc& ds = desc();
697  SPxId leftId = baseId(leaveIdx);
698 
699  if (leftId.isSPxRowId())
700  {
701  idx = number(SPxRowId(leftId));
702  SPxBasis::Desc::Status leaveStat = ds.rowStatus(idx);
703 
704  switch (leaveStat)
705  {
707  assert(rep() == ROW);
708  throw SPxInternalCodeException("XENTER04 This should never happen.");
709  break;
711  assert(rep() == ROW);
712  leavebound = theLBbound[leaveIdx];
713  theLRbound[idx] = leavebound;
714  ds.rowStatus(idx) = dualRowStatus(idx);
715  switch (ds.rowStatus(idx))
716  {
718  objChange += theURbound[idx] * lhs(idx);
719  break;
721  objChange += theLRbound[idx] * rhs(idx);
722  break;
724  objChange += theURbound[idx] * lhs(idx);
725  objChange += theLRbound[idx] * rhs(idx);
726  break;
727  default:
728  break;
729  }
730  break;
732  assert(rep() == ROW);
733  leavebound = theUBbound[leaveIdx];
734  theURbound[idx] = leavebound;
735  ds.rowStatus(idx) = dualRowStatus(idx);
736  switch (ds.rowStatus(idx))
737  {
739  objChange += theURbound[idx] * lhs(idx);
740  break;
742  objChange += theLRbound[idx] * rhs(idx);
743  break;
745  objChange += theURbound[idx] * lhs(idx);
746  objChange += theLRbound[idx] * rhs(idx);
747  break;
748  default:
749  break;
750  }
751  break;
753  assert(rep() == ROW);
754 #if 1
755  throw SPxInternalCodeException("XENTER05 This should never happen.");
756 #else
757  MSG_ERROR( std::cerr << "EENTER98 ERROR: not yet debugged!" << std::endl; )
758  if ((*theCoPvec)[leaveIdx] - theLBbound[leaveIdx] <
759  theUBbound[leaveIdx] - (*theCoPvec)[leaveIdx])
760  {
761  leavebound = theLBbound[leaveIdx];
762  theLRbound[idx] = leavebound;
763  }
764  else
765  {
766  leavebound = theUBbound[leaveIdx];
767  theURbound[idx] = leavebound;
768  }
770 #endif
771  break;
772  // primal/columnwise cases:
774  assert(rep() == COLUMN);
775  throw SPxInternalCodeException("XENTER06 This should never happen.");
776  break;
778  assert(rep() == COLUMN);
779  if (theFvec->delta()[leaveIdx] * enterMax < 0)
780  leavebound = theUBbound[leaveIdx];
781  else
782  leavebound = theLBbound[leaveIdx];
783  theLRbound[idx] = leavebound;
784  theURbound[idx] = leavebound;
785  objChange += leavebound * maxRowObj(leaveIdx);
787  break;
789  assert(rep() == COLUMN);
790  leavebound = theUBbound[leaveIdx];
791  theURbound[idx] = leavebound;
792  objChange += leavebound * maxRowObj(leaveIdx);
794  break;
796  assert(rep() == COLUMN);
797  leavebound = theLBbound[leaveIdx];
798  theLRbound[idx] = leavebound;
799  objChange += leavebound * maxRowObj(leaveIdx);
801  break;
803  assert(rep() == COLUMN);
804  if (enterMax * theFvec->delta()[leaveIdx] < 0)
805  {
806  leavebound = theUBbound[leaveIdx];
807  theURbound[idx] = leavebound;
808  objChange += leavebound * maxRowObj(leaveIdx);
810  }
811  else
812  {
813  leavebound = theLBbound[leaveIdx];
814  theLRbound[idx] = leavebound;
815  objChange += leavebound * maxRowObj(leaveIdx);
817  }
818  break;
819 
820  default:
821  throw SPxInternalCodeException("XENTER07 This should never happen.");
822  }
823  MSG_DEBUG( std::cout << "DENTER06 SPxSolver::getEnterVals2(): row "
824  << idx << ": " << leaveStat
825  << " -> " << ds.rowStatus(idx)
826  << " objChange: " << objChange
827  << std::endl; )
828  }
829 
830  else
831  {
832  assert(leftId.isSPxColId());
833  idx = number(SPxColId(leftId));
834  SPxBasis::Desc::Status leaveStat = ds.colStatus(idx);
835 
836  switch (leaveStat)
837  {
839  assert(rep() == ROW);
840  leavebound = theLBbound[leaveIdx];
841  theLCbound[idx] = leavebound;
842  ds.colStatus(idx) = dualColStatus(idx);
843  switch (ds.colStatus(idx))
844  {
846  objChange += theUCbound[idx] * lower(idx);
847  break;
849  objChange += theLCbound[idx] * upper(idx);
850  break;
852  objChange += theLCbound[idx] * upper(idx);
853  objChange += theUCbound[idx] * lower(idx);
854  break;
855  default:
856  break;
857  }
858  break;
860  assert(rep() == ROW);
861  leavebound = theUBbound[leaveIdx];
862  theUCbound[idx] = leavebound;
863  ds.colStatus(idx) = dualColStatus(idx);
864  switch (ds.colStatus(idx))
865  {
867  objChange += theUCbound[idx] * lower(idx);
868  break;
870  objChange += theLCbound[idx] * upper(idx);
871  break;
873  objChange += theLCbound[idx] * upper(idx);
874  objChange += theUCbound[idx] * lower(idx);
875  break;
876  default:
877  break;
878  }
879  break;
881  assert(rep() == ROW);
882  if (theFvec->delta()[leaveIdx] * enterMax > 0)
883  {
884  leavebound = theLBbound[leaveIdx];
885  theLCbound[idx] = leavebound;
886  }
887  else
888  {
889  leavebound = theUBbound[leaveIdx];
890  theUCbound[idx] = leavebound;
891  }
893  break;
895  assert(rep() == ROW);
896  throw SPxInternalCodeException("XENTER08 This should never happen.");
897  break;
898  // primal/columnwise cases:
900  assert(rep() == COLUMN);
901  if (theFvec->delta()[leaveIdx] * enterMax > 0)
902  leavebound = theLBbound[leaveIdx];
903  else
904  leavebound = theUBbound[leaveIdx];
905  theUCbound[idx] =
906  theLCbound[idx] = leavebound;
907  objChange += maxObj(idx) * leavebound;
909  break;
911  assert(rep() == COLUMN);
912  leavebound = theLBbound[leaveIdx];
913  theLCbound[idx] = leavebound;
914  objChange += maxObj(idx) * leavebound;
916  break;
918  assert(rep() == COLUMN);
919  leavebound = theUBbound[leaveIdx];
920  theUCbound[idx] = leavebound;
921  objChange += maxObj(idx) * leavebound;
923  break;
926  assert(rep() == COLUMN);
927  if (enterMax * theFvec->delta()[leaveIdx] < 0)
928  {
929  leavebound = theUBbound[leaveIdx];
930  theUCbound[idx] = leavebound;
931  objChange += maxObj(idx) * leavebound;
933  }
934  else
935  {
936  leavebound = theLBbound[leaveIdx];
937  theLCbound[idx] = leavebound;
938  objChange += maxObj(idx) * leavebound;
940  }
941  break;
942 
943  default:
944  throw SPxInternalCodeException("XENTER09 This should never happen.");
945  }
946  MSG_DEBUG( std::cout << "DENTER07 SPxSolver::getEnterVals2(): col "
947  << idx << ": " << leaveStat
948  << " -> " << ds.colStatus(idx)
949  << " objChange: " << objChange
950  << std::endl; )
951  }
952 }
953 
954 
955 void
957  SPxId enterId,
958  SPxBasis::Desc::Status enterStat,
959  Real leaveVal,
960  const SVector& vec,
961  Real& objChange
962 )
963 {
964  assert(rep() == COLUMN);
965  int enterIdx;
966  SPxBasis::Desc& ds = desc();
967 
968  if (enterId.isSPxColId())
969  {
970  enterIdx = number(SPxColId(enterId));
971  if (enterStat == SPxBasis::Desc::P_ON_UPPER)
972  {
973  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
974  objChange += theLCbound[enterIdx] * maxObj(enterIdx);
975  }
976  else
977  {
978  ds.colStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
979  objChange += theUCbound[enterIdx] * maxObj(enterIdx);
980  }
981  theFrhs->multAdd(leaveVal, vec);
982  }
983  else
984  {
985  enterIdx = number(SPxRowId(enterId));
986  assert(enterId.isSPxRowId());
987  if (enterStat == SPxBasis::Desc::P_ON_UPPER)
988  {
989  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_LOWER;
990  objChange += (theURbound[enterIdx]) * maxRowObj(enterIdx);
991  }
992  else
993  {
994  ds.rowStatus(enterIdx) = SPxBasis::Desc::P_ON_UPPER;
995  objChange += (theLRbound[enterIdx]) * maxRowObj(enterIdx);
996  }
997  (*theFrhs)[enterIdx] += leaveVal;
998  }
999  if (isId(enterId))
1000  {
1001  theTest[enterIdx] = 0;
1002  isInfeasibleCo[enterIdx] = SPxPricer::NOT_VIOLATED;
1003  }
1004  else
1005  {
1006  theCoTest[enterIdx] = 0;
1007  isInfeasible[enterIdx] = SPxPricer::NOT_VIOLATED;
1008  }
1009 }
1010 
1012  SPxId enterId,
1013  Real enterTest,
1014  SPxBasis::Desc::Status enterStat
1015 )
1016 {
1017  int enterIdx = number(enterId);
1018  if (isId(enterId))
1019  {
1020  theTest[enterIdx] = enterTest;
1021  desc().status(enterIdx) = enterStat;
1022  }
1023  else
1024  {
1025  theCoTest[enterIdx] = enterTest;
1026  desc().coStatus(enterIdx) = enterStat;
1027  }
1028 }
1029 
1030 bool SPxSolver::enter(SPxId& enterId)
1031 {
1032  assert(enterId.isValid());
1033  assert(type() == ENTER);
1034  assert(initialized);
1035 
1036  SPxId none; // invalid id used when enter fails
1037  Real enterTest; // correct test value of entering var
1038  Real enterUB; // upper bound of entering variable
1039  Real enterLB; // lower bound of entering variable
1040  Real enterVal; // current value of entering variable
1041  Real enterMax; // maximum value for entering shift
1042  Real enterPric; // priced value of entering variable
1043  SPxBasis::Desc::Status enterStat; // status of entering variable
1044  Real enterRO; // rhs/obj of entering variable
1045  Real objChange = 0.0;
1046  const SVector* enterVec = enterVector(enterId);
1047 
1048  bool instable = instableEnter;
1049  assert(!instable || instableEnterId.isValid());
1050 
1051  getEnterVals(enterId, enterTest, enterUB, enterLB,
1052  enterVal, enterMax, enterPric, enterStat, enterRO, objChange);
1053 
1054  if (enterTest > -epsilon())
1055  {
1056  rejectEnter(enterId, enterTest, enterStat);
1057  change(-1, none, 0);
1058  objChange = 0.0; // the nonbasicValue is not supposed to be updated in this case
1059 
1060  MSG_DEBUG( std::cout << "DENTER08 rejecting false enter pivot" << std::endl; )
1061 
1062  return true;
1063  }
1064 
1065  /* Before performing the actual basis update, we must determine, how this
1066  is to be accomplished.
1067  */
1068  // BH 2005-11-15: Obviously solve4update() is only called if theFvec.delta()
1069  // is setup (i.e. the indices of the NZEs are stored within it) and there are
1070  // 0 NZEs (???).
1071  // In that case theFvec->delta() is set such that
1072  // Base * theFvec->delta() = enterVec
1073  if (theFvec->delta().isSetup() && theFvec->delta().size() == 0)
1074  SPxBasis::solve4update(theFvec->delta(), *enterVec);
1075 #ifdef ENABLE_ADDITIONAL_CHECKS
1076  else
1077  {
1078  // BH 2005-11-29: This code block seems to check the assertion
1079  // || Base * theFvec->delta() - enterVec ||_2 <= entertol()
1080  DVector tmp(dim());
1081  // BH 2005-11-15: This cast is necessary since SSVector inherits protected from DVector.
1082  tmp = reinterpret_cast<DVector&>(theFvec->delta());
1083  multBaseWith(tmp);
1084  tmp -= *enterVec;
1085  if (tmp.length() > entertol()) {
1086  // This happens frequently and does usually not hurt, so print these
1087  // warnings only with verbose level INFO2 and higher.
1088  MSG_INFO2( (*spxout), (*spxout) << "WENTER09 fVec updated error = "
1089  << tmp.length() << std::endl; )
1090  }
1091  }
1092 #endif // ENABLE_ADDITIONAL_CHECKS
1093 
1094  if (m_numCycle > m_maxCycle)
1095  {
1096  if (-enterMax > 0)
1097  perturbMaxEnter();
1098  else
1099  perturbMinEnter();
1100  }
1101 
1102  Real leaveVal = -enterMax;
1103 
1104  boundflips = 0;
1105  int leaveIdx = theratiotester->selectLeave(leaveVal, enterTest);
1106 
1107  /* in row representation, fixed columns and rows should not leave the basis */
1108  assert(leaveIdx < 0 || !baseId(leaveIdx).isSPxColId() || desc().colStatus(number(SPxColId(baseId(leaveIdx)))) != SPxBasis::Desc::P_FIXED);
1109  assert(leaveIdx < 0 || !baseId(leaveIdx).isSPxRowId() || desc().rowStatus(number(SPxRowId(baseId(leaveIdx)))) != SPxBasis::Desc::P_FIXED);
1110 
1111  instableEnterVal = 0;
1112  instableEnterId = SPxId();
1113  instableEnter = false;
1114 
1115  /*
1116  We now tried to find a variable to leave the basis. If one has been
1117  found, a regular basis update is to be performed.
1118  */
1119  if (leaveIdx >= 0)
1120  {
1121  if (spxAbs(leaveVal) < entertol())
1122  {
1123  if (theUBbound[leaveIdx] != theLBbound[leaveIdx]
1124  && enterStat != Desc::P_FREE && enterStat != Desc::D_FREE)
1125  m_numCycle++;
1126  }
1127  else
1128  m_numCycle /= 2;
1129 
1130  // setup for updating the copricing vector
1132  {
1133  assert(coSolveVector2->isConsistent());
1134  assert(coSolveVector2rhs->isSetup());
1135  assert(coSolveVector3->isConsistent());
1136  assert(coSolveVector3rhs->isSetup());
1137  assert(boundflips > 0);
1139  , unitVecs[leaveIdx], *coSolveVector2rhs, *coSolveVector3rhs);
1140  (*theCoPvec) -= (*coSolveVector3);
1141  }
1142  else if (coSolveVector3)
1143  {
1144  assert(coSolveVector3->isConsistent());
1145  assert(coSolveVector3rhs->isSetup());
1146  assert(boundflips > 0);
1148  (*theCoPvec) -= (*coSolveVector3);
1149  }
1150  else if (coSolveVector2)
1152  else
1153  SPxBasis::coSolve(theCoPvec->delta(), unitVecs[leaveIdx]);
1154 
1155  if( boundflips > 0 )
1156  {
1157  for( int i = coSolveVector3->dim()-1; i >= 0; --i)
1158  {
1159  if( spxAbs((*coSolveVector3)[i]) > epsilon() )
1160  (*thePvec).multAdd(-(*coSolveVector3)[i],(*thecovectors)[i]);
1161  }
1162  // we need to update enterPric in case it was changed by bound flips
1163  if( enterId.isSPxColId() )
1164  enterPric = (*theCoPvec)[number(SPxColId(enterId))];
1165  else
1166  enterPric = (*thePvec)[number(SPxRowId(enterId))];
1167  MSG_INFO3( (*spxout), (*spxout) << "IEBFRT02 "
1168  << "breakpoints passed / bounds flipped = " << boundflips
1169  << std::endl; )
1171  }
1172 
1173  (*theCoPrhs)[leaveIdx] = enterRO;
1174  theCoPvec->value() = (enterRO - enterPric) / theFvec->delta()[leaveIdx];
1175 
1176  if (theCoPvec->value() > epsilon() || theCoPvec->value() < -epsilon())
1177  {
1178  if (pricing() == FULL)
1179  {
1180  thePvec->value() = theCoPvec->value();
1181  setupPupdate();
1182  }
1183  doPupdate();
1184  }
1185  assert(thePvec->isConsistent());
1186  assert(theCoPvec->isConsistent());
1187 
1188  assert(!baseId(leaveIdx).isSPxRowId() || desc().rowStatus(number(SPxRowId(baseId(leaveIdx)))) != SPxBasis::Desc::P_FIXED);
1189  assert(!baseId(leaveIdx).isSPxColId() || desc().colStatus(number(SPxColId(baseId(leaveIdx)))) != SPxBasis::Desc::P_FIXED);
1190 
1191  Real leavebound; // bound on which leaving variable moves
1192  try
1193  {
1194  getEnterVals2(leaveIdx, enterMax, leavebound, objChange);
1195  }
1196  catch( const SPxException& F )
1197  {
1198  rejectEnter(enterId, enterTest, enterStat);
1199  change(-1, none, 0);
1200  throw F;
1201  }
1202 
1203  // process entering variable
1204  theUBbound[leaveIdx] = enterUB;
1205  theLBbound[leaveIdx] = enterLB;
1206 
1207  // compute tests:
1208  updateCoTest();
1209  if (pricing() == FULL)
1210  updateTest();
1211 
1212  // update feasibility vectors
1213  theFvec->value() = leaveVal;
1214  theFvec->update();
1215  (*theFvec)[leaveIdx] = enterVal - leaveVal;
1216 
1217  if (leavebound > epsilon() || leavebound < -epsilon())
1218  theFrhs->multAdd(-leavebound, baseVec(leaveIdx));
1219 
1220  if (enterVal > epsilon() || enterVal < -epsilon())
1221  theFrhs->multAdd(enterVal, *enterVec);
1222 
1223  // update objective funtion value
1224  updateNonbasicValue(objChange);
1225 
1226  // change basis matrix
1227  change(leaveIdx, enterId, enterVec, &(theFvec->delta()));
1228  }
1229  /* No leaving vector could be found that would yield a stable pivot step.
1230  */
1231  else if (leaveVal != -enterMax)
1232  {
1233  /* In the ENTER algorithm, when for a selected entering variable we find only
1234  an instable leaving variable, then the basis change is not conducted.
1235  Instead, we save the entering variable's id in instableEnterId and set
1236  the test value to zero, hoping to find a different leaving
1237  variable with a stable leavingvariable.
1238  If this fails, however, and no more entering variable is found, we have to
1239  perform the instable basis change using instableEnterId. In this (and only
1240  in this) case, the flag instableEnter is set to true.
1241 
1242  leaveVal != enterMax is the case that selectLeave has found only an instable leaving
1243  variable. We store this leaving variable for later if we are not already in the
1244  instable case */
1245 
1246  objChange = 0.0; // the nonbasicValue is not supposed to be updated in this case
1247 
1248  if (!instable)
1249  {
1250  instableEnterId = enterId;
1251  instableEnterVal = enterTest;
1252 
1253  rejectEnter(enterId, 0.0, enterStat);
1254  change(-1, none, 0);
1255 
1256 
1257  return true;
1258  }
1259  else
1260  {
1261  rejectEnter(enterId, enterTest, enterStat);
1262  change(-1, none, 0);
1263  }
1264  }
1265  /* No leaving vector has been selected from the basis. However, if the
1266  shift amount for |fVec| is bounded, we are in the case, that the
1267  entering variable is moved from one bound to its other, before any of
1268  the basis feasibility variables reaches their bound. This may only
1269  happen in primal/columnwise case with upper and lower bounds on
1270  variables.
1271  */
1272  else if (leaveVal < infinity && leaveVal > -infinity)
1273  {
1274  assert(rep() == COLUMN);
1275  assert(leaveVal == -enterMax);
1276 
1277  change(-1, enterId, enterVec);
1278 
1279  theFvec->value() = leaveVal;
1280  theFvec->update();
1281 
1282  ungetEnterVal(enterId, enterStat, leaveVal, *enterVec, objChange);
1283 
1284  // update objective funtion value
1285  updateNonbasicValue(objChange);
1286  }
1287  /* No variable could be selected to leave the basis and even the entering
1288  variable is unbounded --- this is a failure.
1289  */
1290  else
1291  {
1292  /* The following line originally was in the "lastUpdate() > 1" case;
1293  we need it in the INFEASIBLE/UNBOUNDED case, too, to have the
1294  basis descriptor at the correct size.
1295  */
1296  rejectEnter(enterId, enterTest, enterStat);
1297  change(-1, none, 0);
1298 
1299  objChange = 0.0; // the nonbasicValue is not supposed to be updated in this case
1300 
1301  if (lastUpdate() > 1)
1302  {
1303  MSG_INFO3( (*spxout), (*spxout) << "IENTER01 factorization triggered in "
1304  << "enter() for feasibility test" << std::endl; )
1305  factorize();
1306 
1307  /* after a factorization, the entering column/row might not be infeasible or suboptimal anymore, hence we do
1308  * not try to call leave(leaveIdx), but rather return to the main solving loop and call the pricer again
1309  */
1310  return true;
1311  }
1312 
1313  /* do not exit with status infeasible or unbounded if there is only a very small violation
1314  * ROW: recompute the primal variables and activities for another, more precise, round of pricing
1315  */
1316  else if (spxAbs(enterTest) < entertol())
1317  {
1318  MSG_INFO3( (*spxout), (*spxout) << "IENTER11 clean up step to reduce numerical errors" << std::endl; )
1319 
1321  computePvec();
1322  computeCoTest();
1323  computeTest();
1324 
1325  return true;
1326  }
1327 
1328  MSG_INFO3( (*spxout), (*spxout) << "IENTER02 unboundedness/infeasibility found in "
1329  << "enter()" << std::endl; )
1330 
1331  if (rep() == ROW)
1332  {
1333  Real sign;
1334 
1335  dualFarkas.clear();
1336  dualFarkas.setMax(fVec().delta().size() + 1);
1337  sign = (leaveVal > 0 ? -1.0 : 1.0);
1338 
1339  for( int j = 0; j < fVec().delta().size(); ++j )
1340  {
1341  SPxId spxid = baseId(fVec().idx().index(j));
1342 
1343  if( spxid.isSPxRowId() )
1344  dualFarkas.add(number(SPxRowId(spxid)), sign * fVec().delta().value(j));
1345  }
1346 
1347  if( enterId.isSPxRowId() )
1348  dualFarkas.add(number(SPxRowId(enterId)), -sign);
1349 
1351  }
1352  /**@todo if shift() is not zero, we must not conclude primal unboundedness */
1353  else
1354  {
1355  Real sign;
1356 
1357  primalRay.clear();
1358  primalRay.setMax(fVec().delta().size() + 1);
1359  sign = leaveVal > 0 ? 1.0 : -1.0;
1360 
1361  for( int j = 0; j < fVec().delta().size(); ++j )
1362  {
1363  SPxId i = baseId(fVec().idx().index(j));
1364 
1365  if( i.isSPxColId() )
1366  primalRay.add(number(SPxColId(i)), sign*fVec().delta().value(j));
1367  }
1368 
1369  if( enterId.isSPxColId() )
1370  primalRay.add(number(SPxColId(enterId)), -sign);
1371 
1373  }
1374  return false;
1375  }
1376  return true;
1377 }
1378 } // namespace soplex