1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the class library                   */
4    	/*       SoPlex --- the Sequential object-oriented simPlex.                  */
5    	/*                                                                           */
6    	/*  Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SoPlex; see the file LICENSE. If not email to soplex@zib.de.  */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/*      \SubSection{Updating the Basis for Entering Variables}
26   	 */
27   	#include <assert.h>
28   	
29   	#include "soplex/spxdefines.h"
30   	#include "soplex/spxratiotester.h"
31   	#include "soplex/spxpricer.h"
32   	#include "soplex/spxout.h"
33   	#include "soplex/exceptions.h"
34   	
35   	namespace soplex
36   	{
37   	
38   	/*
39   	  In the entering simplex algorithms (i.e. iteratively a vector is selected to
40   	  \em enter the simplex basis as in the dual rowwise and primal columnwise case)
41   	  let $A$ denote the current basis, $x$ and entering vector and $f$ the
42   	  feasibility vector. For a feasible basis $l \le f \le u$ holds.
43   	  For the rowwise case $f$ is obtained by solving $f^T = c^T A^{-1}$,
44   	  wherease in columnwisecase $f = A^{-1} b$.
45   	
46   	  Let us further consider the rowwise case. Exchanging $x$ with the $i$-th
47   	  vector of $A$ yields
48   	
49   	  \begin{equation}\label{update.eq}
50   	  A^{(i)} = E_i A \hbox{, with } E_i = I + e_i (x^T A^{-1} - e_i^T).
51   	  \end{equation}
52   	
53   	  With $E_i^{-1} = I + e_i \frac{e_i^T - \delta^T}{\delta}$,
54   	  $\delta^T = x^T A^{-1}$ one gets the new feasibility vector
55   	
56   	  \begin{eqnarray*}
57   	  (f^{(i)})^T
58   	  &=& c^T (A^{(i)})^{-1}      \\
59   	  &=& c^T A^{-1} + c^T A^{-1} e_i \frac{e_i^T - \delta^T}{\delta_i} \\
60   	  &=& f^T + \frac{f_i}{\delta_i} e_i^T - \frac{f_i}{\delta_i} \delta^T. \\
61   	  \end{eqnarray*}
62   	
63   	  The selection of the leaving vector $i^*$ for the basis must ensure, that for
64   	  all $j \ne i^*$ $f^{(i^*)}_j$ remains within its bounds $l_j$ and $u_j$.
65   	*/
66   	
67   	
68   	/*
69   	  Testing all values of |pVec| against its bounds. If $i$, say, is violated
70   	  the violation is saved as negative value in |theTest[i]|.
71   	*/
72   	
73   	template <class R>
74   	R SPxSolverBase<R>::test(int i, typename SPxBasisBase<R>::Desc::Status stat) const
75   	{
76   	   assert(type() == ENTER);
77   	   assert(!isBasic(stat));
78   	
79   	   R x;
80   	
81   	   switch(stat)
82   	   {
83   	   case SPxBasisBase<R>::Desc::D_FREE:
84   	   case SPxBasisBase<R>::Desc::D_ON_BOTH:
85   	      assert(rep() == ROW);
86   	      x = (*thePvec)[i] - this->lhs(i);
87   	
88   	      if(x < 0)
89   	         return x;
90   	
91   	   // no break: next is else case
92   	   //lint -fallthrough
93   	   case SPxBasisBase<R>::Desc::D_ON_LOWER:
94   	      assert(rep() == ROW);
95   	      return this->rhs(i) - (*thePvec)[i];
96   	
97   	   case SPxBasisBase<R>::Desc::D_ON_UPPER:
98   	      assert(rep() == ROW);
99   	      return (*thePvec)[i] - this->lhs(i);
100  	
101  	   case SPxBasisBase<R>::Desc::P_ON_UPPER:
102  	      assert(rep() == COLUMN);
103  	      return this->maxObj(i) - (*thePvec)[i];
104  	
105  	   case SPxBasisBase<R>::Desc::P_ON_LOWER:
106  	      assert(rep() == COLUMN);
107  	      return (*thePvec)[i] - this->maxObj(i);
108  	
109  	   case SPxBasisBase<R>::Desc::P_FREE :
110  	      x = this->maxObj(i) - (*thePvec)[i];
111  	      return (x < 0) ? x : -x;
112  	
113  	   default:
114  	      return 0;
115  	   }
116  	}
117  	
118  	template <class R>
119  	void SPxSolverBase<R>::computeTest()
120  	{
121  	
122  	   const typename SPxBasisBase<R>::Desc& ds = this->desc();
123  	   R pricingTol = leavetol();
124  	   m_pricingViolCoUpToDate = true;
125  	   m_pricingViolCo = 0;
126  	
127  	   infeasibilitiesCo.clear();
128  	   int sparsitythreshold = (int)(sparsePricingFactor * coDim());
129  	
130  	   for(int i = 0; i < coDim(); ++i)
131  	   {
132  	      typename SPxBasisBase<R>::Desc::Status stat = ds.status(i);
133  	
134  	      if(isBasic(stat))
135  	      {
136  	         theTest[i] = 0.0;
137  	
138  	         if(remainingRoundsEnterCo == 0)
139  	            isInfeasibleCo[i] = SPxPricer<R>::NOT_VIOLATED;
140  	      }
141  	      else
142  	      {
143  	         assert(!isBasic(stat));
144  	         theTest[i] = test(i, stat);
145  	
146  	         if(remainingRoundsEnterCo == 0)
147  	         {
148  	            if(theTest[i] < -pricingTol)
149  	            {
150  	               assert(infeasibilitiesCo.size() < infeasibilitiesCo.max());
151  	               m_pricingViolCo -= theTest[i];
152  	               infeasibilitiesCo.addIdx(i);
153  	               isInfeasibleCo[i] = SPxPricer<R>::VIOLATED;
154  	               ++m_numViol;
155  	            }
156  	            else
157  	               isInfeasibleCo[i] = SPxPricer<R>::NOT_VIOLATED;
158  	
159  	            if(infeasibilitiesCo.size() > sparsitythreshold)
160  	            {
161  	               SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- using dense pricing"
162  	                             << std::endl;)
163  	               remainingRoundsEnterCo = SOPLEX_DENSEROUNDS;
164  	               sparsePricingEnterCo = false;
165  	               infeasibilitiesCo.clear();
166  	            }
167  	         }
168  	         else if(theTest[i] < -pricingTol)
169  	         {
170  	            m_pricingViolCo -= theTest[i];
171  	            ++m_numViol;
172  	         }
173  	      }
174  	   }
175  	
176  	   if(infeasibilitiesCo.size() == 0 && !sparsePricingEnterCo)
177  	      --remainingRoundsEnterCo;
178  	   else if(infeasibilitiesCo.size() <= sparsitythreshold && !sparsePricingEnterCo)
179  	   {
180  	      SPX_MSG_INFO2((*this->spxout),
181  	                    std::streamsize prec = spxout->precision();
182  	
183  	                    if(hyperPricingEnter)
184  	                    (*this->spxout) << " --- using hypersparse pricing, ";
185  	                    else
186  	                       (*this->spxout) << " --- using sparse pricing, ";
187  	                       (*this->spxout) << "sparsity: "
188  	                       << std::setw(6) << std::fixed << std::setprecision(4)
189  	                       << (R) infeasibilitiesCo.size() / coDim()
190  	                       << std::scientific << std::setprecision(int(prec))
191  	                       << std::endl;)
192  	            sparsePricingEnterCo = true;
193  	   }
194  	}
195  	
196  	template <class R>
197  	R SPxSolverBase<R>::computePvec(int i)
198  	{
199  	
200  	   return (*thePvec)[i] = vector(i) * (*theCoPvec);
201  	}
202  	
203  	template <class R>
204  	R SPxSolverBase<R>::computeTest(int i)
205  	{
206  	   typename SPxBasisBase<R>::Desc::Status stat = this->desc().status(i);
207  	
208  	   if(isBasic(stat))
209  	      return theTest[i] = 0;
210  	   else
211  	      return theTest[i] = test(i, stat);
212  	}
213  	
214  	/*
215  	  Testing all values of #coPvec# against its bounds. If $i$, say, is violated
216  	  the violation is saved as negative value in |theCoTest[i]|.
217  	*/
218  	template <class R>
219  	R SPxSolverBase<R>::coTest(int i, typename SPxBasisBase<R>::Desc::Status stat) const
220  	{
221  	   assert(type() == ENTER);
222  	   assert(!isBasic(stat));
223  	
224  	   R x;
225  	
226  	   switch(stat)
227  	   {
228  	   case SPxBasisBase<R>::Desc::D_FREE:
229  	   case SPxBasisBase<R>::Desc::D_ON_BOTH :
230  	      assert(rep() == ROW);
231  	      x = (*theCoPvec)[i] - SPxLPBase<R>::lower(i);
232  	
233  	      if(x < 0)
234  	         return x;
235  	
236  	   // no break: next is else case
237  	   //lint -fallthrough
238  	   case SPxBasisBase<R>::Desc::D_ON_LOWER:
239  	      assert(rep() == ROW);
240  	      return SPxLPBase<R>::upper(i) - (*theCoPvec)[i];
241  	
242  	   case SPxBasisBase<R>::Desc::D_ON_UPPER:
243  	      assert(rep() == ROW);
244  	      return (*theCoPvec)[i] - SPxLPBase<R>::lower(i);
245  	
246  	   case SPxBasisBase<R>::Desc::P_ON_UPPER:
247  	      assert(rep() == COLUMN);
248  	      return (*theCoPvec)[i] - this->maxRowObj(i);             // slacks !
249  	
250  	   case SPxBasisBase<R>::Desc::P_ON_LOWER:
251  	      assert(rep() == COLUMN);
252  	      return this->maxRowObj(i) - (*theCoPvec)[i];             // slacks !
253  	
254  	   default:
255  	      return 0;
256  	   }
257  	}
258  	
259  	template <class R>
260  	void SPxSolverBase<R>::computeCoTest()
261  	{
262  	   int i;
263  	   R pricingTol = leavetol();
264  	   m_pricingViolUpToDate = true;
265  	   m_pricingViol = 0;
266  	   m_numViol = 0;
267  	   infeasibilities.clear();
268  	   int sparsitythreshold = (int)(sparsePricingFactor * dim());
269  	   const typename SPxBasisBase<R>::Desc& ds = this->desc();
270  	
271  	   for(i = dim() - 1; i >= 0; --i)
272  	   {
273  	      typename SPxBasisBase<R>::Desc::Status stat = ds.coStatus(i);
274  	
275  	      if(isBasic(stat))
276  	      {
277  	         theCoTest[i] = 0;
278  	
279  	         if(remainingRoundsEnter == 0)
280  	            isInfeasible[i] = SPxPricer<R>::NOT_VIOLATED;
281  	      }
282  	      else
283  	      {
284  	         theCoTest[i] = coTest(i, stat);
285  	
286  	         if(remainingRoundsEnter == 0)
287  	         {
288  	            if(theCoTest[i] < -pricingTol)
289  	            {
290  	               assert(infeasibilities.size() < infeasibilities.max());
291  	               m_pricingViol -= theCoTest[i];
292  	               infeasibilities.addIdx(i);
293  	               isInfeasible[i] = SPxPricer<R>::VIOLATED;
294  	               ++m_numViol;
295  	            }
296  	            else
297  	               isInfeasible[i] = SPxPricer<R>::NOT_VIOLATED;
298  	
299  	            if(infeasibilities.size() > sparsitythreshold)
300  	            {
301  	               SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- using dense pricing"
302  	                             << std::endl;)
303  	               remainingRoundsEnter = SOPLEX_DENSEROUNDS;
304  	               sparsePricingEnter = false;
305  	               infeasibilities.clear();
306  	            }
307  	         }
308  	         else if(theCoTest[i] < -pricingTol)
309  	         {
310  	            m_pricingViol -= theCoTest[i];
311  	            ++m_numViol;
312  	         }
313  	      }
314  	   }
315  	
316  	   if(infeasibilities.size() == 0 && !sparsePricingEnter)
317  	      --remainingRoundsEnter;
318  	   else if(infeasibilities.size() <= sparsitythreshold && !sparsePricingEnter)
319  	   {
320  	      SPX_MSG_INFO2((*this->spxout),
321  	                    std::streamsize prec = spxout->precision();
322  	
323  	                    if(hyperPricingEnter)
324  	                    (*this->spxout) << " --- using hypersparse pricing, ";
325  	                    else
326  	                       (*this->spxout) << " --- using sparse pricing, ";
327  	                       (*this->spxout) << "sparsity: "
328  	                       << std::setw(6) << std::fixed << std::setprecision(4)
329  	                       << (R) infeasibilities.size() / dim()
330  	                       << std::scientific << std::setprecision(int(prec))
331  	                       << std::endl;)
332  	            sparsePricingEnter = true;
333  	   }
334  	}
335  	
336  	/*
337  	  The following methods require propersy initialized vectors |fVec| and
338  	  #coPvec#.
339  	*/
340  	template <class R>
341  	void SPxSolverBase<R>::updateTest()
342  	{
343  	   thePvec->delta().setup();
344  	
345  	   const IdxSet& idx = thePvec->idx();
346  	   const typename SPxBasisBase<R>::Desc& ds = this->desc();
347  	   R pricingTol = leavetol();
348  	
349  	   int i;
350  	   updateViolsCo.clear();
351  	
352  	   for(i = idx.size() - 1; i >= 0; --i)
353  	   {
354  	      int j = idx.index(i);
355  	      typename SPxBasisBase<R>::Desc::Status stat = ds.status(j);
356  	
357  	      if(!isBasic(stat))
358  	      {
359  	         if(m_pricingViolCoUpToDate && theTest[j] < -pricingTol)
360  	            m_pricingViolCo += theTest[j];
361  	
362  	         theTest[j] = test(j, stat);
363  	
364  	         if(sparsePricingEnterCo)
365  	         {
366  	            if(theTest[j] < -pricingTol)
367  	            {
368  	               assert(remainingRoundsEnterCo == 0);
369  	               m_pricingViolCo -= theTest[j];
370  	
371  	               if(isInfeasibleCo[j] == SPxPricer<R>::NOT_VIOLATED)
372  	               {
373  	                  infeasibilitiesCo.addIdx(j);
374  	                  isInfeasibleCo[j] = SPxPricer<R>::VIOLATED;
375  	               }
376  	
377  	               if(hyperPricingEnter)
378  	                  updateViolsCo.addIdx(j);
379  	            }
380  	            else
381  	            {
382  	               isInfeasibleCo[j] = SPxPricer<R>::NOT_VIOLATED;
383  	            }
384  	         }
385  	         else if(theTest[j] < -pricingTol)
386  	            m_pricingViolCo -= theTest[j];
387  	      }
388  	      else
389  	      {
390  	         isInfeasibleCo[j] = SPxPricer<R>::NOT_VIOLATED;
391  	         theTest[j] = 0;
392  	      }
393  	   }
394  	}
395  	
396  	template <class R>
397  	void SPxSolverBase<R>::updateCoTest()
398  	{
399  	   theCoPvec->delta().setup();
400  	
401  	   const IdxSet& idx = theCoPvec->idx();
402  	   const typename SPxBasisBase<R>::Desc& ds = this->desc();
403  	   R pricingTol = leavetol();
404  	
405  	   int i;
406  	   updateViols.clear();
407  	
408  	   for(i = idx.size() - 1; i >= 0; --i)
409  	   {
410  	      int j = idx.index(i);
411  	      typename SPxBasisBase<R>::Desc::Status stat = ds.coStatus(j);
412  	
413  	      if(!isBasic(stat))
414  	      {
415  	         if(m_pricingViolUpToDate && theCoTest[j] < -pricingTol)
416  	            m_pricingViol += theCoTest[j];
417  	
418  	         theCoTest[j] = coTest(j, stat);
419  	
420  	         if(sparsePricingEnter)
421  	         {
422  	            if(theCoTest[j] < -pricingTol)
423  	            {
424  	               assert(remainingRoundsEnter == 0);
425  	               m_pricingViol -= theCoTest[j];
426  	
427  	               if(isInfeasible[j] == SPxPricer<R>::NOT_VIOLATED)
428  	               {
429  	                  //                if( !hyperPricingEnter )
430  	                  infeasibilities.addIdx(j);
431  	                  isInfeasible[j] = SPxPricer<R>::VIOLATED;
432  	               }
433  	
434  	               if(hyperPricingEnter)
435  	                  updateViols.addIdx(j);
436  	            }
437  	            else
438  	            {
439  	               // @todo do we need to remove index j from infeasibilitiesCo?
440  	               isInfeasible[j] = SPxPricer<R>::NOT_VIOLATED;
441  	            }
442  	         }
443  	         else if(theCoTest[j] < -pricingTol)
444  	            m_pricingViol -= theCoTest[j];
445  	      }
446  	      else
447  	      {
448  	         isInfeasible[j] = SPxPricer<R>::NOT_VIOLATED;
449  	         theCoTest[j] = 0;
450  	      }
451  	   }
452  	}
453  	
454  	
455  	
456  	/*  \Section{Compute statistics on entering variable}
457  	    Here is a list of variables relevant when including |Id| to the basis.
458  	    They are computed by |computeEnterStats()|.
459  	*/
460  	template <class R>
461  	void SPxSolverBase<R>::getEnterVals
462  	(
463  	   SPxId enterId,
464  	   R& enterTest,
465  	   R& enterUB,
466  	   R& enterLB,
467  	   R& enterVal,
468  	   R& enterMax,
469  	   R& enterPric,
470  	   typename SPxBasisBase<R>::Desc::Status& enterStat,
471  	   R& enterRO,
472  	   StableSum<R>& objChange
473  	)
474  	{
475  	   int enterIdx;
476  	   typename SPxBasisBase<R>::Desc& ds = this->desc();
477  	
478  	   if(enterId.isSPxColId())
479  	   {
480  	      enterIdx = this->number(SPxColId(enterId));
481  	      enterStat = ds.colStatus(enterIdx);
482  	      assert(!isBasic(enterStat));
483  	
484  	      /*      For an #Id# to enter the basis we better recompute the Test value.
485  	       */
486  	      if(rep() == COLUMN)
487  	      {
488  	         computePvec(enterIdx);
489  	         enterTest = computeTest(enterIdx);
490  	         theTest[enterIdx] = 0;
491  	      }
492  	      else
493  	      {
494  	         enterTest = coTest()[enterIdx];
495  	         theCoTest[enterIdx] = 0;
496  	      }
497  	
498  	      switch(enterStat)
499  	      {
500  	      // primal/columnwise cases:
501  	      case SPxBasisBase<R>::Desc::P_ON_UPPER :
502  	         assert(rep() == COLUMN);
503  	         enterUB = theUCbound[enterIdx];
504  	         enterLB = theLCbound[enterIdx];
505  	         enterVal = enterUB;
506  	         enterMax = enterLB - enterUB;
507  	         enterPric = (*thePvec)[enterIdx];
508  	         enterRO = this->maxObj(enterIdx);
509  	         objChange -= enterVal * enterRO;
510  	
511  	         if(enterLB <= R(-infinity))
512  	            ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_LOWER;
513  	         else if(EQ(enterLB, enterUB, this->epsilon()))
514  	            ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE;
515  	         else
516  	            ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH;
517  	
518  	         break;
519  	
520  	      case SPxBasisBase<R>::Desc::P_ON_LOWER :
521  	         assert(rep() == COLUMN);
522  	         enterUB = theUCbound[enterIdx];
523  	         enterLB = theLCbound[enterIdx];
524  	         enterVal = enterLB;
525  	         enterMax = enterUB - enterLB;
526  	         enterPric = (*thePvec)[enterIdx];
527  	         enterRO = this->maxObj(enterIdx);
528  	         objChange -= enterVal * enterRO;
529  	
530  	         if(enterUB >= R(infinity))
531  	            ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_UPPER;
532  	         else if(EQ(enterLB, enterUB, this->epsilon()))
533  	            ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE;
534  	         else
535  	            ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH;
536  	
537  	         break;
538  	
539  	      case SPxBasisBase<R>::Desc::P_FREE :
540  	         assert(rep() == COLUMN);
541  	         enterUB = theUCbound[enterIdx];
542  	         enterLB = theLCbound[enterIdx];
543  	         enterVal = 0;
544  	         enterPric = (*thePvec)[enterIdx];
545  	         enterRO = this->maxObj(enterIdx);
546  	         ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_UNDEFINED;
547  	         enterMax = (enterRO - enterPric > 0) ? R(infinity) : R(-infinity);
548  	         break;
549  	
550  	      // dual/rowwise cases:
551  	      case SPxBasisBase<R>::Desc::D_ON_UPPER :
552  	         assert(rep() == ROW);
553  	         assert(theUCbound[enterIdx] < R(infinity));
554  	         enterUB = theUCbound[enterIdx];
555  	         enterLB = R(-infinity);
556  	         enterMax = R(-infinity);
557  	         enterVal = enterUB;
558  	         enterPric = (*theCoPvec)[enterIdx];
559  	         enterRO = SPxLPBase<R>::lower(enterIdx);
560  	         objChange -= enterRO * enterVal;
561  	         ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
562  	         break;
563  	
564  	      case SPxBasisBase<R>::Desc::D_ON_LOWER :
565  	         assert(rep() == ROW);
566  	         assert(theLCbound[enterIdx] > R(-infinity));
567  	         enterLB = theLCbound[enterIdx];
568  	         enterUB = R(infinity);
569  	         enterMax = R(infinity);
570  	         enterVal = enterLB;
571  	         enterPric = (*theCoPvec)[enterIdx];
572  	         enterRO = SPxLPBase<R>::upper(enterIdx);
573  	         objChange -= enterRO * enterVal;
574  	         ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
575  	         break;
576  	
577  	      case SPxBasisBase<R>::Desc::D_FREE:
578  	         assert(rep() == ROW);
579  	         assert(SPxLPBase<R>::lower(enterIdx) == SPxLPBase<R>::upper(enterIdx));
580  	         enterUB = R(infinity);
581  	         enterLB = R(-infinity);
582  	         enterVal = 0;
583  	         enterRO = SPxLPBase<R>::upper(enterIdx);
584  	         enterPric = (*theCoPvec)[enterIdx];
585  	
586  	         if(enterPric > enterRO)
587  	            enterMax = R(infinity);
588  	         else
589  	            enterMax = R(-infinity);
590  	
591  	         ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_FIXED;
592  	         break;
593  	
594  	      case SPxBasisBase<R>::Desc::D_ON_BOTH :
595  	         assert(rep() == ROW);
596  	         enterPric = (*theCoPvec)[enterIdx];
597  	
598  	         if(enterPric > SPxLPBase<R>::upper(enterIdx))
599  	         {
600  	            enterLB = theLCbound[enterIdx];
601  	            enterUB = R(infinity);
602  	            enterMax = R(infinity);
603  	            enterVal = enterLB;
604  	            enterRO = SPxLPBase<R>::upper(enterIdx);
605  	            ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
606  	         }
607  	         else
608  	         {
609  	            enterUB = theUCbound[enterIdx];
610  	            enterVal = enterUB;
611  	            enterRO = SPxLPBase<R>::lower(enterIdx);
612  	            enterLB = R(-infinity);
613  	            enterMax = R(-infinity);
614  	            ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
615  	         }
616  	
617  	         objChange -= theLCbound[enterIdx] * SPxLPBase<R>::upper(enterIdx);
618  	         objChange -= theUCbound[enterIdx] * SPxLPBase<R>::lower(enterIdx);
619  	         break;
620  	
621  	      default:
622  	         throw SPxInternalCodeException("XENTER01 This should never happen.");
623  	      }
624  	
625  	      SPxOut::debug(this, "DENTER03 SPxSolverBase::getEnterVals() : col {}: {} -> {} objChange: {}\n",
626  	                    enterIdx, enterStat, ds.colStatus(enterIdx), objChange);
627  	   }
628  	
629  	   else
630  	   {
631  	      assert(enterId.isSPxRowId());
632  	      enterIdx = this->number(SPxRowId(enterId));
633  	      enterStat = ds.rowStatus(enterIdx);
634  	      assert(!isBasic(enterStat));
635  	
636  	      /*      For an #Id# to enter the basis we better recompute the Test value.
637  	       */
638  	      if(rep() == ROW)
639  	      {
640  	         computePvec(enterIdx);
641  	         enterTest = computeTest(enterIdx);
642  	         theTest[enterIdx] = 0;
643  	      }
644  	      else
645  	      {
646  	         enterTest = coTest()[enterIdx];
647  	         theCoTest[enterIdx] = 0;
648  	      }
649  	
650  	      switch(enterStat)
651  	      {
652  	      // primal/columnwise cases:
653  	      case SPxBasisBase<R>::Desc::P_ON_UPPER :
654  	         assert(rep() == COLUMN);
655  	         enterUB = theURbound[enterIdx];
656  	         enterLB = theLRbound[enterIdx];
657  	         enterVal = enterLB;
658  	         enterMax = enterUB - enterLB;
659  	         enterPric = (*theCoPvec)[enterIdx];
660  	         enterRO = this->maxRowObj(enterIdx);
661  	         objChange -= enterRO * enterVal;
662  	
663  	         if(enterUB >= R(infinity))
664  	            ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_LOWER;
665  	         else if(EQ(enterLB, enterUB, this->epsilon()))
666  	            ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE;
667  	         else
668  	            ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH;
669  	
670  	         break;
671  	
672  	      case SPxBasisBase<R>::Desc::P_ON_LOWER :
673  	         assert(rep() == COLUMN);
674  	         enterUB = theURbound[enterIdx];
675  	         enterLB = theLRbound[enterIdx];
676  	         enterVal = enterUB;
677  	         enterMax = enterLB - enterUB;
678  	         enterPric = (*theCoPvec)[enterIdx];
679  	         enterRO = this->maxRowObj(enterIdx);
680  	         objChange -= enterRO * enterVal;
681  	
682  	         if(enterLB <= R(-infinity))
683  	            ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_UPPER;
684  	         else if(EQ(enterLB, enterUB, this->epsilon()))
685  	            ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE;
686  	         else
687  	            ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH;
688  	
689  	         break;
690  	
691  	      case SPxBasisBase<R>::Desc::P_FREE :
692  	         assert(rep() == COLUMN);
693  	#if 1
694  	         throw SPxInternalCodeException("XENTER02 This should never happen.");
695  	#else
696  	         SPX_MSG_ERROR(std::cerr << "EENTER99 ERROR: not yet debugged!" << std::endl;)
697  	         enterPric = (*theCoPvec)[enterIdx];
698  	         enterRO = this->maxRowObj(enterIdx);
699  	         ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_UNDEFINED;
700  	#endif
701  	         break;
702  	
703  	      // dual/rowwise cases:
704  	      case SPxBasisBase<R>::Desc::D_ON_UPPER :
705  	         assert(rep() == ROW);
706  	         assert(theURbound[enterIdx] < R(infinity));
707  	         enterUB = theURbound[enterIdx];
708  	         enterLB = R(-infinity);
709  	         enterVal = enterUB;
710  	         enterMax = R(-infinity);
711  	         enterPric = (*thePvec)[enterIdx];
712  	         enterRO = this->lhs(enterIdx);
713  	         objChange -= enterRO * enterVal;
714  	         ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
715  	         break;
716  	
717  	      case SPxBasisBase<R>::Desc::D_ON_LOWER :
718  	         assert(rep() == ROW);
719  	         assert(theLRbound[enterIdx] > R(-infinity));
720  	         enterLB = theLRbound[enterIdx];
721  	         enterUB = R(infinity);
722  	         enterVal = enterLB;
723  	         enterMax = R(infinity);
724  	         enterPric = (*thePvec)[enterIdx];
725  	         enterRO = this->rhs(enterIdx);
726  	         objChange -= enterRO * enterVal;
727  	         ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
728  	         break;
729  	
730  	      case SPxBasisBase<R>::Desc::D_FREE:
731  	         assert(rep() == ROW);
732  	         assert(this->rhs(enterIdx) == this->lhs(enterIdx));
733  	         enterUB = R(infinity);
734  	         enterLB = R(-infinity);
735  	         enterVal = 0;
736  	         enterPric = (*thePvec)[enterIdx];
737  	         enterRO = this->rhs(enterIdx);
738  	         enterMax = (enterPric > enterRO) ? R(infinity) : R(-infinity);
739  	         ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_FIXED;
740  	         break;
741  	
742  	      case SPxBasisBase<R>::Desc::D_ON_BOTH :
743  	         assert(rep() == ROW);
744  	         enterPric = (*thePvec)[enterIdx];
745  	
746  	         if(enterPric > this->rhs(enterIdx))
747  	         {
748  	            enterLB = theLRbound[enterIdx];
749  	            enterVal = enterLB;
750  	            enterUB = R(infinity);
751  	            enterMax = R(infinity);
752  	            enterRO = this->rhs(enterIdx);
753  	            ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
754  	         }
755  	         else
756  	         {
757  	            enterUB = theURbound[enterIdx];
758  	            enterVal = enterUB;
759  	            enterLB = R(-infinity);
760  	            enterMax = R(-infinity);
761  	            enterRO = this->lhs(enterIdx);
762  	            ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
763  	         }
764  	
765  	         objChange -= theLRbound[enterIdx] * this->rhs(enterIdx);
766  	         objChange -= theURbound[enterIdx] * this->lhs(enterIdx);
767  	         break;
768  	
769  	      default:
770  	         throw SPxInternalCodeException("XENTER03 This should never happen.");
771  	      }
772  	
773  	      SPxOut::debug(this, "DENTER05 SPxSolverBase::getEnterVals() : row {}: {} -> {} objChange: {}\n",
774  	                    enterIdx, enterStat, ds.rowStatus(enterIdx), objChange);
775  	   }
776  	}
777  	
778  	/*      process leaving variable
779  	 */
780  	template <class R>
781  	void SPxSolverBase<R>::getEnterVals2
782  	(
783  	   int leaveIdx,
784  	   R enterMax,
785  	   R& leavebound,
786  	   StableSum<R>& objChange
787  	)
788  	{
789  	   int idx;
790  	   typename SPxBasisBase<R>::Desc& ds = this->desc();
791  	   SPxId leftId = this->baseId(leaveIdx);
792  	
793  	   if(leftId.isSPxRowId())
794  	   {
795  	      idx = this->number(SPxRowId(leftId));
796  	      typename SPxBasisBase<R>::Desc::Status leaveStat = ds.rowStatus(idx);
797  	
798  	      // coverity[switch_selector_expr_is_constant]
799  	      switch(leaveStat)
800  	      {
801  	      case SPxBasisBase<R>::Desc::P_FIXED :
802  	         assert(rep() == ROW);
803  	         throw SPxInternalCodeException("XENTER04 This should never happen.");
804  	         break;
805  	
806  	      case SPxBasisBase<R>::Desc::P_ON_UPPER :
807  	         assert(rep() == ROW);
808  	         leavebound = theLBbound[leaveIdx];
809  	         theLRbound[idx] = leavebound;
810  	         ds.rowStatus(idx) = this->dualRowStatus(idx);
811  	
812  	         switch(ds.rowStatus(idx))
813  	         {
814  	         case SPxBasisBase<R>::Desc::D_ON_UPPER :
815  	            objChange += theURbound[idx] * this->lhs(idx);
816  	            break;
817  	
818  	         case SPxBasisBase<R>::Desc::D_ON_LOWER :
819  	            objChange += theLRbound[idx] * this->rhs(idx);
820  	            break;
821  	
822  	         case SPxBasisBase<R>::Desc::D_ON_BOTH :
823  	            objChange += theURbound[idx] * this->lhs(idx);
824  	            objChange += theLRbound[idx] * this->rhs(idx);
825  	            break;
826  	
827  	         default:
828  	            break;
829  	         }
830  	
831  	         break;
832  	
833  	      case SPxBasisBase<R>::Desc::P_ON_LOWER :
834  	         assert(rep() == ROW);
835  	         leavebound = theUBbound[leaveIdx];
836  	         theURbound[idx] = leavebound;
837  	         ds.rowStatus(idx) = this->dualRowStatus(idx);
838  	
839  	         switch(ds.rowStatus(idx))
840  	         {
841  	         case SPxBasisBase<R>::Desc::D_ON_UPPER :
842  	            objChange += theURbound[idx] * this->lhs(idx);
843  	            break;
844  	
845  	         case SPxBasisBase<R>::Desc::D_ON_LOWER :
846  	            objChange += theLRbound[idx] * this->rhs(idx);
847  	            break;
848  	
849  	         case SPxBasisBase<R>::Desc::D_ON_BOTH :
850  	            objChange += theURbound[idx] * this->lhs(idx);
851  	            objChange += theLRbound[idx] * this->rhs(idx);
852  	            break;
853  	
854  	         default:
855  	            break;
856  	         }
857  	
858  	         break;
859  	
860  	      case SPxBasisBase<R>::Desc::P_FREE :
861  	         assert(rep() == ROW);
862  	#if 1
863  	         throw SPxInternalCodeException("XENTER05 This should never happen.");
864  	#else
865  	         SPX_MSG_ERROR(std::cerr << "EENTER98 ERROR: not yet debugged!" << std::endl;)
866  	
867  	         if((*theCoPvec)[leaveIdx] - theLBbound[leaveIdx] <
868  	               theUBbound[leaveIdx] - (*theCoPvec)[leaveIdx])
869  	         {
870  	            leavebound = theLBbound[leaveIdx];
871  	            theLRbound[idx] = leavebound;
872  	         }
873  	         else
874  	         {
875  	            leavebound = theUBbound[leaveIdx];
876  	            theURbound[idx] = leavebound;
877  	         }
878  	
879  	         ds.rowStatus(idx) = SPxBasisBase<R>::Desc::D_UNDEFINED;
880  	#endif
881  	         break;
882  	
883  	      // primal/columnwise cases:
884  	      case SPxBasisBase<R>::Desc::D_UNDEFINED :
885  	         assert(rep() == COLUMN);
886  	         throw SPxInternalCodeException("XENTER06 This should never happen.");
887  	         break;
888  	
889  	      case SPxBasisBase<R>::Desc::D_FREE :
890  	         assert(rep() == COLUMN);
891  	
892  	         if(theFvec->delta()[leaveIdx] * enterMax < 0)
893  	            leavebound = theUBbound[leaveIdx];
894  	         else
895  	            leavebound = theLBbound[leaveIdx];
896  	
897  	         theLRbound[idx] = leavebound;
898  	         theURbound[idx] = leavebound;
899  	         objChange += leavebound * this->maxRowObj(leaveIdx);
900  	         ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED;
901  	         break;
902  	
903  	      case SPxBasisBase<R>::Desc::D_ON_UPPER :
904  	         assert(rep() == COLUMN);
905  	         leavebound = theUBbound[leaveIdx];
906  	         theURbound[idx] = leavebound;
907  	         objChange += leavebound * this->maxRowObj(leaveIdx);
908  	         ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
909  	         break;
910  	
911  	      case SPxBasisBase<R>::Desc::D_ON_LOWER :
912  	         assert(rep() == COLUMN);
913  	         leavebound = theLBbound[leaveIdx];
914  	         theLRbound[idx] = leavebound;
915  	         objChange += leavebound * this->maxRowObj(leaveIdx);
916  	         ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
917  	         break;
918  	
919  	      case SPxBasisBase<R>::Desc::D_ON_BOTH :
920  	         assert(rep() == COLUMN);
921  	
922  	         if(enterMax * theFvec->delta()[leaveIdx] < 0)
923  	         {
924  	            leavebound = theUBbound[leaveIdx];
925  	            theURbound[idx] = leavebound;
926  	            objChange += leavebound * this->maxRowObj(leaveIdx);
927  	            ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
928  	         }
929  	         else
930  	         {
931  	            leavebound = theLBbound[leaveIdx];
932  	            theLRbound[idx] = leavebound;
933  	            objChange += leavebound * this->maxRowObj(leaveIdx);
934  	            ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
935  	         }
936  	
937  	         break;
938  	
939  	      default:
940  	         throw SPxInternalCodeException("XENTER07 This should never happen.");
941  	      }
942  	
943  	      SPxOut::debug(this, "DENTER06 SPxSolverBase::getEnterVals2(): row {}: {} -> {} objChange: {}\n",
944  	                    idx, leaveStat, ds.rowStatus(idx), objChange);
945  	   }
946  	
947  	   else
948  	   {
949  	      assert(leftId.isSPxColId());
950  	      idx = this->number(SPxColId(leftId));
951  	      typename SPxBasisBase<R>::Desc::Status leaveStat = ds.colStatus(idx);
952  	
953  	      // coverity[switch_selector_expr_is_constant]
954  	      switch(leaveStat)
955  	      {
956  	      case SPxBasisBase<R>::Desc::P_ON_UPPER :
957  	         assert(rep() == ROW);
958  	         leavebound = theLBbound[leaveIdx];
959  	         theLCbound[idx] = leavebound;
960  	         ds.colStatus(idx) = this->dualColStatus(idx);
961  	
962  	         switch(ds.colStatus(idx))
963  	         {
964  	         case SPxBasisBase<R>::Desc::D_ON_UPPER :
965  	            objChange += theUCbound[idx] * this->lower(idx);
966  	            break;
967  	
968  	         case SPxBasisBase<R>::Desc::D_ON_LOWER :
969  	            objChange += theLCbound[idx] * this->upper(idx);
970  	            break;
971  	
972  	         case SPxBasisBase<R>::Desc::D_ON_BOTH :
973  	            objChange += theLCbound[idx] * this->upper(idx);
974  	            objChange += theUCbound[idx] * this->lower(idx);
975  	            break;
976  	
977  	         default:
978  	            break;
979  	         }
980  	
981  	         break;
982  	
983  	      case SPxBasisBase<R>::Desc::P_ON_LOWER :
984  	         assert(rep() == ROW);
985  	         leavebound = theUBbound[leaveIdx];
986  	         theUCbound[idx] = leavebound;
987  	         ds.colStatus(idx) = this->dualColStatus(idx);
988  	
989  	         switch(ds.colStatus(idx))
990  	         {
991  	         case SPxBasisBase<R>::Desc::D_ON_UPPER :
992  	            objChange += theUCbound[idx] * this->lower(idx);
993  	            break;
994  	
995  	         case SPxBasisBase<R>::Desc::D_ON_LOWER :
996  	            objChange += theLCbound[idx] * this->upper(idx);
997  	            break;
998  	
999  	         case SPxBasisBase<R>::Desc::D_ON_BOTH :
1000 	            objChange += theLCbound[idx] * this->upper(idx);
1001 	            objChange += theUCbound[idx] * this->lower(idx);
1002 	            break;
1003 	
1004 	         default:
1005 	            break;
1006 	         }
1007 	
1008 	         break;
1009 	
1010 	      case SPxBasisBase<R>::Desc::P_FREE :
1011 	         assert(rep() == ROW);
1012 	
1013 	         if(theFvec->delta()[leaveIdx] * enterMax > 0)
1014 	         {
1015 	            leavebound = theLBbound[leaveIdx];
1016 	            theLCbound[idx] = leavebound;
1017 	         }
1018 	         else
1019 	         {
1020 	            leavebound = theUBbound[leaveIdx];
1021 	            theUCbound[idx] = leavebound;
1022 	         }
1023 	
1024 	         ds.colStatus(idx) = SPxBasisBase<R>::Desc::D_UNDEFINED;
1025 	         break;
1026 	
1027 	      case SPxBasisBase<R>::Desc::P_FIXED:
1028 	         assert(rep() == ROW);
1029 	         throw SPxInternalCodeException("XENTER08 This should never happen.");
1030 	         break;
1031 	
1032 	      // primal/columnwise cases:
1033 	      case SPxBasisBase<R>::Desc::D_FREE :
1034 	         assert(rep() == COLUMN);
1035 	
1036 	         if(theFvec->delta()[leaveIdx] * enterMax > 0)
1037 	            leavebound = theLBbound[leaveIdx];
1038 	         else
1039 	            leavebound = theUBbound[leaveIdx];
1040 	
1041 	         theUCbound[idx] =
1042 	            theLCbound[idx] = leavebound;
1043 	         objChange += this->maxObj(idx) * leavebound;
1044 	         ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED;
1045 	         break;
1046 	
1047 	      case SPxBasisBase<R>::Desc::D_ON_UPPER :
1048 	         assert(rep() == COLUMN);
1049 	         leavebound = theLBbound[leaveIdx];
1050 	         theLCbound[idx] = leavebound;
1051 	         objChange += this->maxObj(idx) * leavebound;
1052 	         ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1053 	         break;
1054 	
1055 	      case SPxBasisBase<R>::Desc::D_ON_LOWER :
1056 	         assert(rep() == COLUMN);
1057 	         leavebound = theUBbound[leaveIdx];
1058 	         theUCbound[idx] = leavebound;
1059 	         objChange += this->maxObj(idx) * leavebound;
1060 	         ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1061 	         break;
1062 	
1063 	      case SPxBasisBase<R>::Desc::D_ON_BOTH :
1064 	      case SPxBasisBase<R>::Desc::D_UNDEFINED :
1065 	         assert(rep() == COLUMN);
1066 	
1067 	         if(enterMax * theFvec->delta()[leaveIdx] < 0)
1068 	         {
1069 	            leavebound = theUBbound[leaveIdx];
1070 	            theUCbound[idx] = leavebound;
1071 	            objChange += this->maxObj(idx) * leavebound;
1072 	            ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1073 	         }
1074 	         else
1075 	         {
1076 	            leavebound = theLBbound[leaveIdx];
1077 	            theLCbound[idx] = leavebound;
1078 	            objChange += this->maxObj(idx) * leavebound;
1079 	            ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1080 	         }
1081 	
1082 	         break;
1083 	
1084 	      default:
1085 	         throw SPxInternalCodeException("XENTER09 This should never happen.");
1086 	      }
1087 	
1088 	      SPxOut::debug(this, "DENTER07 SPxSolverBase::getEnterVals2(): col {}: {} -> {} objChange: {}\n",
1089 	                    idx, leaveStat, ds.colStatus(idx), objChange);
1090 	   }
1091 	}
1092 	
1093 	template <class R>
1094 	void
1095 	SPxSolverBase<R>::ungetEnterVal(
1096 	   SPxId enterId,
1097 	   typename SPxBasisBase<R>::Desc::Status enterStat,
1098 	   R leaveVal,
1099 	   const SVectorBase<R>& vec,
1100 	   StableSum<R>& objChange
1101 	)
1102 	{
1103 	   assert(rep() == COLUMN);
1104 	   int enterIdx;
1105 	   typename SPxBasisBase<R>::Desc& ds = this->desc();
1106 	
1107 	   if(enterId.isSPxColId())
1108 	   {
1109 	      enterIdx = this->number(SPxColId(enterId));
1110 	
1111 	      if(enterStat == SPxBasisBase<R>::Desc::P_ON_UPPER)
1112 	      {
1113 	         ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1114 	         objChange += theLCbound[enterIdx] * this->maxObj(enterIdx);
1115 	      }
1116 	      else
1117 	      {
1118 	         ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1119 	         objChange += theUCbound[enterIdx] * this->maxObj(enterIdx);
1120 	      }
1121 	
1122 	      theFrhs->multAdd(leaveVal, vec);
1123 	   }
1124 	   else
1125 	   {
1126 	      enterIdx = this->number(SPxRowId(enterId));
1127 	      assert(enterId.isSPxRowId());
1128 	
1129 	      if(enterStat == SPxBasisBase<R>::Desc::P_ON_UPPER)
1130 	      {
1131 	         ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1132 	         objChange += (theURbound[enterIdx]) * this->maxRowObj(enterIdx);
1133 	      }
1134 	      else
1135 	      {
1136 	         ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1137 	         objChange += (theLRbound[enterIdx]) * this->maxRowObj(enterIdx);
1138 	      }
1139 	
1140 	      (*theFrhs)[enterIdx] += leaveVal;
1141 	   }
1142 	
1143 	   if(isId(enterId))
1144 	   {
1145 	      theTest[enterIdx] = 0;
1146 	      isInfeasibleCo[enterIdx] = SPxPricer<R>::NOT_VIOLATED;
1147 	   }
1148 	   else
1149 	   {
1150 	      theCoTest[enterIdx] = 0;
1151 	      isInfeasible[enterIdx] = SPxPricer<R>::NOT_VIOLATED;
1152 	   }
1153 	}
1154 	
1155 	template <class R>
1156 	void SPxSolverBase<R>::rejectEnter(
1157 	   SPxId enterId,
1158 	   R enterTest,
1159 	   typename SPxBasisBase<R>::Desc::Status enterStat
1160 	)
1161 	{
1162 	   int enterIdx = this->number(enterId);
1163 	
1164 	   if(isId(enterId))
1165 	   {
1166 	      theTest[enterIdx] = enterTest;
1167 	      this->desc().status(enterIdx) = enterStat;
1168 	   }
1169 	   else
1170 	   {
1171 	      theCoTest[enterIdx] = enterTest;
1172 	      this->desc().coStatus(enterIdx) = enterStat;
1173 	   }
1174 	}
1175 	
1176 	template <class R>
1177 	void SPxSolverBase<R>::computePrimalray4Col(R direction, SPxId enterId)
1178 	{
1179 	   R sign = (direction > 0 ? 1.0 : -1.0);
1180 	
1181 	   primalRay.clear();
1182 	   primalRay.setMax(fVec().delta().size() + 1);
1183 	
1184 	   for(int j = 0; j < fVec().delta().size(); ++j)
1185 	   {
1186 	      SPxId i = this->baseId(fVec().idx().index(j));
1187 	
1188 	      if(i.isSPxColId())
1189 	         primalRay.add(this->number(SPxColId(i)), sign * fVec().delta().value(j));
1190 	   }
1191 	
1192 	   if(enterId.isSPxColId())
1193 	      primalRay.add(this->number(SPxColId(enterId)), -sign);
1194 	}
1195 	
1196 	template <class R>
1197 	void SPxSolverBase<R>::computeDualfarkas4Row(R direction, SPxId enterId)
1198 	{
1199 	   R sign = (direction > 0 ? -1.0 : 1.0);
1200 	
1201 	   dualFarkas.clear();
1202 	   dualFarkas.setMax(fVec().delta().size() + 1);
1203 	
1204 	   for(int j = 0; j < fVec().delta().size(); ++j)
1205 	   {
1206 	      SPxId spxid = this->baseId(fVec().idx().index(j));
1207 	
1208 	      if(spxid.isSPxRowId())
1209 	         dualFarkas.add(this->number(SPxRowId(spxid)), sign * fVec().delta().value(j));
1210 	   }
1211 	
1212 	   if(enterId.isSPxRowId())
1213 	      dualFarkas.add(this->number(SPxRowId(enterId)), -sign);
1214 	}
1215 	
1216 	template <class R>
1217 	bool SPxSolverBase<R>::enter(SPxId& enterId, bool polish)
1218 	{
1219 	   assert(enterId.isValid());
1220 	   assert(type() == ENTER);
1221 	   assert(initialized);
1222 	
1223 	   SPxId none;          // invalid id used when enter fails
1224 	   R enterTest;      // correct test value of entering var
1225 	   R enterUB;        // upper bound of entering variable
1226 	   R enterLB;        // lower bound of entering variable
1227 	   R enterVal;       // current value of entering variable
1228 	   R enterMax;       // maximum value for entering shift
1229 	   R enterPric;      // priced value of entering variable
1230 	   typename SPxBasisBase<R>::Desc::Status enterStat;      // status of entering variable
1231 	   R enterRO;        // rhs/obj of entering variable
1232 	   StableSum<R> objChange;
1233 	   const SVectorBase<R>* enterVec = enterVector(enterId);
1234 	
1235 	   bool instable = instableEnter;
1236 	   assert(!instable || instableEnterId.isValid());
1237 	
1238 	   getEnterVals(enterId, enterTest, enterUB, enterLB,
1239 	                enterVal, enterMax, enterPric, enterStat, enterRO, objChange);
1240 	
1241 	   if(!polish && enterTest > -epsilon())
1242 	   {
1243 	      rejectEnter(enterId, enterTest, enterStat);
1244 	      this->change(-1, none, 0);
1245 	
1246 	      SPxOut::debug(this, "DENTER08 rejecting false enter pivot\n");
1247 	
1248 	      return false;
1249 	   }
1250 	
1251 	   /*  Before performing the actual basis update, we must determine, how this
1252 	       is to be accomplished.
1253 	   */
1254 	   // BH 2005-11-15: Obviously solve4update() is only called if theFvec.delta()
1255 	   // is setup (i.e. the indices of the NZEs are stored within it) and there are
1256 	   // 0 NZEs (???).
1257 	   // In that case theFvec->delta() is set such that
1258 	   //   Base * theFvec->delta() = enterVec
1259 	   if(theFvec->delta().isSetup() && theFvec->delta().size() == 0)
1260 	      SPxBasisBase<R>::solve4update(theFvec->delta(), *enterVec);
1261 	
1262 	#ifdef ENABLE_ADDITIONAL_CHECKS
1263 	   else
1264 	   {
1265 	      // BH 2005-11-29: This code block seems to check the assertion
1266 	      //   || Base * theFvec->delta() - enterVec ||_2 <= entertol()
1267 	      VectorBase<R> tmp(dim());
1268 	      // BH 2005-11-15: This cast is necessary since SSVectorBase<R>  inherits protected from VectorBase<R>.
1269 	      tmp = reinterpret_cast<VectorBase<R>&>(theFvec->delta());
1270 	      this->multBaseWith(tmp);
1271 	      tmp -= *enterVec;
1272 	
1273 	      if(tmp.length() > entertol())
1274 	      {
1275 	         // This happens frequently and does usually not hurt, so print these
1276 	         // warnings only with verbose level INFO2 and higher.
1277 	         SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "WENTER09 fVec updated error = "
1278 	                       << tmp.length() << std::endl;)
1279 	      }
1280 	   }
1281 	
1282 	#endif  // ENABLE_ADDITIONAL_CHECKS
1283 	
1284 	   if(!polish && m_numCycle > m_maxCycle)
1285 	   {
1286 	      if(-enterMax > 0)
1287 	         perturbMaxEnter();
1288 	      else
1289 	         perturbMinEnter();
1290 	   }
1291 	
1292 	   R leaveVal = -enterMax;
1293 	
1294 	   boundflips = 0;
1295 	   int leaveIdx = theratiotester->selectLeave(leaveVal, enterTest, polish);
1296 	
1297 	   /* in row representation, fixed columns and rows should not leave the basis */
1298 	   assert(leaveIdx < 0 || !this->baseId(leaveIdx).isSPxColId()
1299 	          || this->desc().colStatus(this->number(SPxColId(this->baseId(leaveIdx)))) !=
1300 	          SPxBasisBase<R>::Desc::P_FIXED);
1301 	   assert(leaveIdx < 0 || !this->baseId(leaveIdx).isSPxRowId()
1302 	          || this->desc().rowStatus(this->number(SPxRowId(this->baseId(leaveIdx)))) !=
1303 	          SPxBasisBase<R>::Desc::P_FIXED);
1304 	
1305 	   instableEnterVal = 0;
1306 	   instableEnterId = SPxId();
1307 	   instableEnter = false;
1308 	
1309 	   /*
1310 	     We now tried to find a variable to leave the basis. If one has been
1311 	     found, a regular basis update is to be performed.
1312 	   */
1313 	   if(leaveIdx >= 0)
1314 	   {
1315 	      if(spxAbs(leaveVal) < entertol())
1316 	      {
1317 	         if(NE(theUBbound[leaveIdx], theLBbound[leaveIdx], this->tolerances()->epsilon())
1318 	               && enterStat != SPxBasisBase<R>::Desc::P_FREE && enterStat != SPxBasisBase<R>::Desc::D_FREE)
1319 	         {
1320 	            m_numCycle++;
1321 	            enterCycles++;
1322 	         }
1323 	      }
1324 	      else
1325 	         m_numCycle /= 2;
1326 	
1327 	      // setup for updating the copricing vector
1328 	      if(coSolveVector3 && coSolveVector2)
1329 	      {
1330 	         assert(coSolveVector2->isConsistent());
1331 	         assert(coSolveVector2rhs->isSetup());
1332 	         assert(coSolveVector3->isConsistent());
1333 	         assert(coSolveVector3rhs->isSetup());
1334 	         assert(boundflips > 0);
1335 	         SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector2, *coSolveVector3
1336 	                                  , unitVecs[leaveIdx], *coSolveVector2rhs, *coSolveVector3rhs);
1337 	         (*theCoPvec) -= (*coSolveVector3);
1338 	      }
1339 	      else if(coSolveVector3)
1340 	      {
1341 	         assert(coSolveVector3->isConsistent());
1342 	         assert(coSolveVector3rhs->isSetup());
1343 	         assert(boundflips > 0);
1344 	         SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector3, unitVecs[leaveIdx],
1345 	                                  *coSolveVector3rhs);
1346 	         (*theCoPvec) -= (*coSolveVector3);
1347 	      }
1348 	      else if(coSolveVector2)
1349 	         SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector2, unitVecs[leaveIdx],
1350 	                                  *coSolveVector2rhs);
1351 	      else
1352 	         SPxBasisBase<R>::coSolve(theCoPvec->delta(), unitVecs[leaveIdx]);
1353 	
1354 	      if(boundflips > 0)
1355 	      {
1356 	         for(int i = coSolveVector3->dim() - 1; i >= 0; --i)
1357 	         {
1358 	            if(spxAbs((*coSolveVector3)[i]) > epsilon())
1359 	               (*thePvec).multAdd(-(*coSolveVector3)[i], (*thecovectors)[i]);
1360 	         }
1361 	
1362 	         // we need to update enterPric in case it was changed by bound flips
1363 	         if(enterId.isSPxColId())
1364 	            enterPric = (*theCoPvec)[this->number(SPxColId(enterId))];
1365 	         else
1366 	            enterPric = (*thePvec)[this->number(SPxRowId(enterId))];
1367 	
1368 	         SPxOut::debug(this, "IEBFRT02 breakpoints passed / bounds flipped = {}\n", boundflips);
1369 	         totalboundflips += boundflips;
1370 	      }
1371 	
1372 	      (*theCoPrhs)[leaveIdx] = enterRO;
1373 	      theCoPvec->value() = (enterRO - enterPric) / theFvec->delta()[leaveIdx];
1374 	
1375 	      if(theCoPvec->value() > epsilon() || theCoPvec->value() < -epsilon())
1376 	      {
1377 	         if(pricing() == FULL)
1378 	         {
1379 	            thePvec->value() = theCoPvec->value();
1380 	            setupPupdate();
1381 	         }
1382 	
1383 	         doPupdate();
1384 	      }
1385 	
1386 	      assert(thePvec->isConsistent());
1387 	      assert(theCoPvec->isConsistent());
1388 	
1389 	      assert(!this->baseId(leaveIdx).isSPxRowId()
1390 	             || this->desc().rowStatus(this->number(SPxRowId(this->baseId(leaveIdx)))) !=
1391 	             SPxBasisBase<R>::Desc::P_FIXED);
1392 	      assert(!this->baseId(leaveIdx).isSPxColId()
1393 	             || this->desc().colStatus(this->number(SPxColId(this->baseId(leaveIdx)))) !=
1394 	             SPxBasisBase<R>::Desc::P_FIXED);
1395 	
1396 	      R leavebound;             // bound on which leaving variable moves
1397 	
1398 	      try
1399 	      {
1400 	         getEnterVals2(leaveIdx, enterMax, leavebound, objChange);
1401 	      }
1402 	      catch(const SPxException& F)
1403 	      {
1404 	         rejectEnter(enterId, enterTest, enterStat);
1405 	         this->change(-1, none, 0);
1406 	         throw F;
1407 	      }
1408 	
1409 	      //  process entering variable
1410 	      theUBbound[leaveIdx] = enterUB;
1411 	      theLBbound[leaveIdx] = enterLB;
1412 	
1413 	      //  compute tests:
1414 	      updateCoTest();
1415 	
1416 	      if(pricing() == FULL)
1417 	         updateTest();
1418 	
1419 	      // update feasibility vectors
1420 	      theFvec->value() = leaveVal;
1421 	      theFvec->update();
1422 	      (*theFvec)[leaveIdx] = enterVal - leaveVal;
1423 	
1424 	      if(leavebound > epsilon() || leavebound < -epsilon())
1425 	         theFrhs->multAdd(-leavebound, this->baseVec(leaveIdx));
1426 	
1427 	      if(enterVal > epsilon() || enterVal < -epsilon())
1428 	         theFrhs->multAdd(enterVal, *enterVec);
1429 	
1430 	      // update objective funtion value
1431 	      updateNonbasicValue(objChange);
1432 	
1433 	      //  change basis matrix
1434 	      this->change(leaveIdx, enterId, enterVec, &(theFvec->delta()));
1435 	
1436 	      return true;
1437 	   }
1438 	   /*  No leaving vector could be found that would yield a stable pivot step.
1439 	    */
1440 	   else if(NE(leaveVal, -enterMax, this->tolerances()->epsilon()))
1441 	   {
1442 	      /* In the ENTER algorithm, when for a selected entering variable we find only
1443 	         an instable leaving variable, then the basis change is not conducted.
1444 	         Instead, we save the entering variable's id in instableEnterId and set
1445 	         the test value to zero, hoping to find a different leaving
1446 	         variable with a stable leavingvariable.
1447 	         If this fails, however, and no more entering variable is found, we have to
1448 	         perform the instable basis change using instableEnterId. In this (and only
1449 	         in this) case, the flag instableEnter is set to true.
1450 	
1451 	         leaveVal != enterMax is the case that selectLeave has found only an instable leaving
1452 	         variable. We store this leaving variable for later if we are not already in the
1453 	         instable case */
1454 	
1455 	      if(!instable)
1456 	      {
1457 	         instableEnterId = enterId;
1458 	         instableEnterVal = enterTest;
1459 	
1460 	         SPxOut::debug(this, "DENTER09 rejecting enter pivot and looking for others\n");
1461 	
1462 	         rejectEnter(enterId, enterTest / 10.0, enterStat);
1463 	         this->change(-1, none, 0);
1464 	      }
1465 	      else
1466 	      {
1467 	         SPxOut::debug(this, "DENTER10 rejecting enter pivot in instable state, resetting values\n");
1468 	         rejectEnter(enterId, enterTest, enterStat);
1469 	         this->change(-1, none, 0);
1470 	      }
1471 	
1472 	      return false;
1473 	   }
1474 	   /*  No leaving vector has been selected from the basis. However, if the
1475 	       shift amount for |fVec| is bounded, we are in the case, that the
1476 	       entering variable is moved from one bound to its other, before any of
1477 	       the basis feasibility variables reaches their bound. This may only
1478 	       happen in primal/columnwise case with upper and lower bounds on
1479 	       variables.
1480 	   */
1481 	   else if(!polish && leaveVal < R(infinity) && leaveVal > R(-infinity))
1482 	   {
1483 	      assert(rep() == COLUMN);
1484 	      assert(EQ(leaveVal, -enterMax, this->epsilon()));
1485 	
1486 	      this->change(-1, enterId, enterVec);
1487 	
1488 	      theFvec->value() = leaveVal;
1489 	      theFvec->update();
1490 	
1491 	      ungetEnterVal(enterId, enterStat, leaveVal, *enterVec, objChange);
1492 	
1493 	      // update objective funtion value
1494 	      updateNonbasicValue(objChange);
1495 	
1496 	      SPxOut::debug(this, "DENTER11 moving entering variable from one bound to the other\n");
1497 	
1498 	      return false;
1499 	   }
1500 	   /*  No variable could be selected to leave the basis and even the entering
1501 	       variable is unbounded --- this is a failure.
1502 	   */
1503 	   else
1504 	   {
1505 	      /* The following line originally was in the "lastUpdate() > 1" case;
1506 	         we need it in the INFEASIBLE/UNBOUNDED case, too, to have the
1507 	         basis descriptor at the correct size.
1508 	      */
1509 	      rejectEnter(enterId, enterTest, enterStat);
1510 	      this->change(-1, none, 0);
1511 	
1512 	      if(polish)
1513 	         return false;
1514 	
1515 	      else if(this->lastUpdate() > 1)
1516 	      {
1517 	         SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IENTER01 factorization triggered in "
1518 	                       << "enter() for feasibility test" << std::endl;)
1519 	
1520 	         try
1521 	         {
1522 	            factorize();
1523 	         }
1524 	         catch(const SPxStatusException& E)
1525 	         {
1526 	            // don't exit immediately but handle the singularity correctly
1527 	            assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::SINGULAR);
1528 	            SPX_MSG_INFO3((*this->spxout),
1529 	                          (*this->spxout) << "Caught exception in factorization: " << E.what() <<
1530 	                          std::endl;)
1531 	         }
1532 	
1533 	         /* after a factorization, the entering column/row might not be infeasible or suboptimal anymore, hence we do
1534 	          * not try to call leave(leaveIdx), but rather return to the main solving loop and call the pricer again
1535 	          */
1536 	         return false;
1537 	      }
1538 	
1539 	      /* do not exit with status infeasible or unbounded if there is only a very small violation
1540 	       * ROW: recompute the primal variables and activities for another, more precise, round of pricing
1541 	       */
1542 	      else if(spxAbs(enterTest) < entertol())
1543 	      {
1544 	         SPX_MSG_INFO3((*this->spxout),
1545 	                       (*this->spxout) << "IENTER11 clean up step to reduce numerical errors" <<
1546 	                       std::endl;)
1547 	
1548 	         SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs);
1549 	         computePvec();
1550 	         computeCoTest();
1551 	         computeTest();
1552 	
1553 	         return false;
1554 	      }
1555 	
1556 	      SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IENTER02 unboundedness/infeasibility found in "
1557 	                    << "enter()" << std::endl;)
1558 	
1559 	      if(rep() == ROW)
1560 	      {
1561 	         computeDualfarkas4Row(leaveVal, enterId);
1562 	         setBasisStatus(SPxBasisBase<R>::INFEASIBLE);
1563 	      }
1564 	      /**@todo if shift() is not zero, we must not conclude primal unboundedness */
1565 	      else
1566 	      {
1567 	         computePrimalray4Col(leaveVal, enterId);
1568 	         setBasisStatus(SPxBasisBase<R>::UNBOUNDED);
1569 	      }
1570 	
1571 	      return false;
1572 	   }
1573 	}
1574 	} // namespace soplex
1575