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   	#include <assert.h>
26   	#include <stdio.h>
27   	
28   	#include "soplex/spxdefines.h"
29   	
30   	/*
31   	  Here comes our implementation of the Harris procedure improved by shifting
32   	  bounds. The basic idea is to allow a slight infeasibility within |delta| to
33   	  allow for more freedom when selecting the leaving variable. This freedom
34   	  may then be used for selecting numerically stable variables yielding great
35   	  improvements.
36   	
37   	  The algorithms operates in two phases. In a first phase, the maximum value
38   	  |val| is determined, when infeasibility within |delta| is allowed. In the second
39   	  phase, between all variables with value < |val| the one is selected which
40   	  has the largest update Vector component. However, this may not
41   	  always yield an improvement. In that case, we shift the variable towards
42   	  infeasibility and retry. This avoids cycling in the shifted LP.
43   	*/
44   	
45   	namespace soplex
46   	{
47   	
48   	
49   	#define SOPLEX_MINSTAB          1e-5
50   	#define SOPLEX_LOWSTAB          1e-10
51   	#define SOPLEX_TRIES                   2
52   	#define SOPLEX_SHORTVAL         1e-5
53   	#define SOPLEX_DELTA_SHIFT      1e-5
54   	#define SOPLEX_EPSILON          1e-10
55   	
56   	template <class R>
57   	void SPxFastRT<R>::resetTols()
58   	{
59   	   epsilon = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_EPSILON);
60   	}
61   	
62   	template <class R>
63   	void SPxFastRT<R>::tighten()
64   	{
65   	   R delta_shift = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_DELTA_SHIFT);
66   	
67   	   if(fastDelta >= this->delta + delta_shift)
68   	   {
69   	      fastDelta -= delta_shift;
70   	
71   	      if(fastDelta > this->tolerances()->scaleAccordingToEpsilon(1e-4))
72   	         fastDelta -= 2 * delta_shift;
73   	   }
74   	
75   	   if(minStab < this->tolerances()->scaleAccordingToEpsilon(SOPLEX_MINSTAB))
76   	   {
77   	      minStab /= 0.90;
78   	
79   	      if(minStab < this->tolerances()->floatingPointFeastol())
80   	         minStab /= 0.90;
81   	   }
82   	}
83   	
84   	template <class R>
85   	void SPxFastRT<R>::relax()
86   	{
87   	   R delta_shift = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_DELTA_SHIFT);
88   	   minStab *= 0.95;
89   	   fastDelta += 3 * delta_shift;
90   	}
91   	
92   	template <class R>
93   	R SPxFastRT<R>::minStability(R maxabs)
94   	{
95   	   if(maxabs < 1000.0)
96   	      return minStab;
97   	
98   	   return maxabs * minStab / 1000.0;
99   	}
100  	
101  	/* The code below implements phase 1 of the ratio test. It differs from the description in the
102  	 * Ph.D. thesis page 57 as follows: It uses \f$\delta_i = d_i - s_i - \delta\f$ if \f$d_i > s_i\f$.
103  	 *
104  	 * This change leads to the following behavior of the source code. Consider the first case (x >
105  	 * epsilon, u < infinity): If u - vec[i] <= 0, vec[i] violates the upper bound. In the Harris ratio
106  	 * test, we would compute (u - vec[i] + delta)/upd[i]. The code computes instead delta/upd[i].
107  	 */
108  	template <class R>
109  	int SPxFastRT<R>::maxDelta(
110  	   R& val,                                /* on return: maximum step length */
111  	   R& maxabs,                             /* on return: maximum absolute value in upd VectorBase<R> */
112  	   UpdateVector<R>& update,
113  	   const VectorBase<R>& lowBound,
114  	   const VectorBase<R>& upBound,
115  	   int start,
116  	   int incr) const
117  	{
118  	   int i, sel;
119  	   R x, y, max;
120  	   R u, l;
121  	   bool leaving = this->m_type == SPxSolverBase<R>::LEAVE;
122  	   bool enterrowrep = !leaving && this->thesolver->theRep == SPxSolverBase<R>::ROW;
123  	
124  	   R mabs = maxabs;
125  	
126  	   const R* up = upBound.get_const_ptr();
127  	   const R* low = lowBound.get_const_ptr();
128  	   const R* vec = update.get_const_ptr();
129  	   const R* upd = update.delta().values();
130  	   const int* idx = update.delta().indexMem();
131  	
132  	   sel = -1;
133  	   max = val;
134  	
135  	   if(update.delta().isSetup())
136  	   {
137  	      const int* last = idx + update.delta().size();
138  	
139  	      for(idx += start; idx < last; idx += incr)
140  	      {
141  	         i = *idx;
142  	
143  	         /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
144  	         if(leaving && ((iscoid && this->thesolver->isCoBasic(i)) || (!iscoid
145  	                        && this->thesolver->isBasic(i))))
146  	            continue;
147  	
148  	         if(enterrowrep && this->thesolver->baseId(i).isSPxColId()
149  	               && this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(i))))
150  	               == SPxBasisBase<R>::Desc::P_FIXED)
151  	            continue;
152  	
153  	         x = upd[i];
154  	
155  	         if(x > epsilonZero())
156  	         {
157  	            // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
158  	            mabs = (x > mabs) ? x : mabs;
159  	            u = up[i];
160  	
161  	            if(u < R(infinity))
162  	            {
163  	               y = u - vec[i];
164  	
165  	               if(y <= 0)
166  	                  x = fastDelta / x;
167  	               else
168  	                  x = (y + fastDelta) / x;
169  	
170  	               if(x < max)
171  	               {
172  	                  max = x;
173  	                  sel = i;
174  	               }
175  	            }
176  	         }
177  	         else if(x < -epsilonZero())
178  	         {
179  	            // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
180  	            mabs = (-x > mabs) ? -x : mabs;
181  	            l = low[i];
182  	
183  	            if(l > R(-infinity))
184  	            {
185  	               y = l - vec[i];
186  	
187  	               if(y >= 0)
188  	                  x = - fastDelta / x;
189  	               else
190  	                  x = (y - fastDelta) / x;
191  	
192  	               if(x < max)
193  	               {
194  	                  max = x;
195  	                  sel = i;
196  	               }
197  	            }
198  	         }
199  	      }
200  	   }
201  	   else
202  	   {
203  	      /* In this case, the indices of the semi-sparse Vector update.delta() are not set up and are filled below. */
204  	      int* l_idx = update.delta().altIndexMem();
205  	      R* uval = update.delta().altValues();
206  	      const R* uend = uval + update.dim();
207  	      assert(uval == upd);
208  	
209  	      for(i = 0; uval < uend; ++uval, ++i)
210  	      {
211  	         x = *uval;
212  	
213  	         if(x != 0.0)
214  	         {
215  	            if(x >= -epsilonZero() && x <= epsilonZero())
216  	            {
217  	               *uval = 0.0;
218  	               continue;
219  	            }
220  	            else
221  	               *l_idx++ = i;
222  	
223  	            /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
224  	            if(leaving && ((iscoid && this->thesolver->isCoBasic(i)) || (!iscoid
225  	                           && this->thesolver->isBasic(i))))
226  	               continue;
227  	
228  	            if(enterrowrep && this->thesolver->baseId(i).isSPxColId()
229  	                  && this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(i))))
230  	                  == SPxBasisBase<R>::Desc::P_FIXED)
231  	               continue;
232  	
233  	            if(x > epsilonZero())
234  	            {
235  	               mabs = (x > mabs) ? x : mabs;
236  	               u = up[i];
237  	
238  	               if(u < R(infinity))
239  	               {
240  	                  y = u - vec[i];
241  	
242  	                  if(y <= 0)
243  	                     x = fastDelta / x;
244  	                  else
245  	                     x = (y + fastDelta) / x;
246  	
247  	                  if(x < max)
248  	                  {
249  	                     max = x;
250  	                     sel = i;
251  	                  }
252  	               }
253  	            }
254  	            else if(x < -epsilonZero())
255  	            {
256  	               mabs = (-x > mabs) ? -x : mabs;
257  	               l = low[i];
258  	
259  	               if(l > R(-infinity))
260  	               {
261  	                  y = l - vec[i];
262  	
263  	                  if(y >= 0)
264  	                     x = - fastDelta / x;
265  	                  else
266  	                     x = (y - fastDelta) / x;
267  	
268  	                  if(x < max)
269  	                  {
270  	                     max = x;
271  	                     sel = i;
272  	                  }
273  	               }
274  	            }
275  	         }
276  	      }
277  	
278  	      update.delta().setSize(int(l_idx - update.delta().indexMem()));
279  	      update.delta().forceSetup();
280  	   }
281  	
282  	   val = max;
283  	   maxabs = mabs;
284  	   return sel;
285  	}
286  	
287  	/* See maxDelta() */
288  	template <class R>
289  	int SPxFastRT<R>::minDelta(
290  	   R& val,
291  	   R& maxabs,
292  	   UpdateVector<R>& update,
293  	   const VectorBase<R>& lowBound,
294  	   const VectorBase<R>& upBound,
295  	   int start,
296  	   int incr) const
297  	{
298  	   int i, sel;
299  	   R x, y, max;
300  	   R u, l;
301  	   bool leaving = (this->m_type == SPxSolverBase<R>::LEAVE);
302  	   bool enterrowrep = !leaving && this->thesolver->theRep == SPxSolverBase<R>::ROW;
303  	
304  	   R mabs = maxabs;
305  	
306  	   const R* up = upBound.get_const_ptr();
307  	   const R* low = lowBound.get_const_ptr();
308  	   const R* vec = update.get_const_ptr();
309  	   const R* upd = update.delta().values();
310  	   const int* idx = update.delta().indexMem();
311  	
312  	   sel = -1;
313  	   max = val;
314  	
315  	   if(update.delta().isSetup())
316  	   {
317  	      const int* last = idx + update.delta().size();
318  	
319  	      for(idx += start; idx < last; idx += incr)
320  	      {
321  	         i = *idx;
322  	         x = upd[i];
323  	
324  	         /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
325  	         if(leaving && ((iscoid && this->thesolver->isCoBasic(i)) || (!iscoid
326  	                        && this->thesolver->isBasic(i))))
327  	            continue;
328  	
329  	         if(enterrowrep && this->thesolver->baseId(i).isSPxColId()
330  	               && this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(i))))
331  	               == SPxBasisBase<R>::Desc::P_FIXED)
332  	            continue;
333  	
334  	         if(x > epsilonZero())
335  	         {
336  	            // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
337  	            mabs = (x > mabs) ? x : mabs;
338  	            l = low[i];
339  	
340  	            if(l > R(-infinity))
341  	            {
342  	               y = l - vec[i];
343  	
344  	               if(y >= 0)
345  	                  x = - fastDelta / x;
346  	               else
347  	                  x = (y - fastDelta) / x;
348  	
349  	               if(x > max)
350  	               {
351  	                  max = x;
352  	                  sel = i;
353  	               }
354  	            }
355  	         }
356  	         else if(x < -epsilonZero())
357  	         {
358  	            // @todo check wether mabs should be computed only over bounded vars, i.e., in the if block below
359  	            mabs = (-x > mabs) ? -x : mabs;
360  	            u = up[i];
361  	
362  	            if(u < R(infinity))
363  	            {
364  	               y = u - vec[i];
365  	
366  	               if(y <= 0)
367  	                  x = fastDelta / x;
368  	               else
369  	                  x = (y + fastDelta) / x;
370  	
371  	               if(x > max)
372  	               {
373  	                  max = x;
374  	                  sel = i;
375  	               }
376  	            }
377  	         }
378  	      }
379  	   }
380  	   else
381  	   {
382  	      /* In this case, the indices of the semi-sparse VectorBase<R> update.delta() are not set up and are filled below. */
383  	      int* l_idx = update.delta().altIndexMem();
384  	      R* uval = update.delta().altValues();
385  	      const R* uend = uval + update.dim();
386  	      assert(uval == upd);
387  	
388  	      for(i = 0; uval < uend; ++uval, ++i)
389  	      {
390  	         x = *uval;
391  	
392  	         if(x != 0.0)
393  	         {
394  	            if(x >= -epsilonZero() && x <= epsilonZero())
395  	            {
396  	               *uval = 0.0;
397  	               continue;
398  	            }
399  	            else
400  	               *l_idx++ = i;
401  	
402  	            /* in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables */
403  	            if(leaving && ((iscoid && this->thesolver->isCoBasic(i)) || (!iscoid
404  	                           && this->thesolver->isBasic(i))))
405  	               continue;
406  	
407  	            if(enterrowrep && this->thesolver->baseId(i).isSPxColId()
408  	                  && this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(i))))
409  	                  == SPxBasisBase<R>::Desc::P_FIXED)
410  	               continue;
411  	
412  	
413  	            if(x > epsilonZero())
414  	            {
415  	               mabs = (x > mabs) ? x : mabs;
416  	               l = low[i];
417  	
418  	               if(l > R(-infinity))
419  	               {
420  	                  y = l - vec[i];
421  	
422  	                  if(y >= 0)
423  	                     x = - fastDelta / x;
424  	                  else
425  	                     x = (y - fastDelta) / x;
426  	
427  	                  if(x > max)
428  	                  {
429  	                     max = x;
430  	                     sel = i;
431  	                  }
432  	               }
433  	            }
434  	            else if(x < -epsilonZero())
435  	            {
436  	               mabs = (-x > mabs) ? -x : mabs;
437  	               u = up[i];
438  	
439  	               if(u < R(infinity))
440  	               {
441  	                  y = u - vec[i];
442  	
443  	                  if(y <= 0)
444  	                     x = fastDelta / x;
445  	                  else
446  	                     x = (y + fastDelta) / x;
447  	
448  	                  if(x > max)
449  	                  {
450  	                     max = x;
451  	                     sel = i;
452  	                  }
453  	               }
454  	            }
455  	         }
456  	      }
457  	
458  	      update.delta().setSize(int(l_idx - update.delta().indexMem()));
459  	      update.delta().forceSetup();
460  	   }
461  	
462  	   val = max;
463  	   maxabs = mabs;
464  	
465  	   return sel;
466  	}
467  	
468  	template <class R>
469  	int SPxFastRT<R>::maxDelta(
470  	   R& val,
471  	   R& maxabs)
472  	{
473  	   assert(this->m_type == SPxSolverBase<R>::ENTER);
474  	   return maxDelta(val, maxabs,
475  	                   this->thesolver->fVec(), this->thesolver->lbBound(), this->thesolver->ubBound(), 0, 1);
476  	}
477  	
478  	template <class R>
479  	int SPxFastRT<R>::minDelta(
480  	   R& val,
481  	   R& maxabs)
482  	{
483  	   assert(this->m_type == SPxSolverBase<R>::ENTER);
484  	   return minDelta(val, maxabs,
485  	                   this->thesolver->fVec(), this->thesolver->lbBound(), this->thesolver->ubBound(), 0, 1);
486  	}
487  	
488  	template <class R>
489  	SPxId SPxFastRT<R>::maxDelta(
490  	   int& nr,
491  	   R& max,                                /* on return: maximum step length */
492  	   R& maxabs)                             /* on return: maximum absolute value in delta VectorBase<R> */
493  	{
494  	   /* The following cause side effects on coPvec and pVec - both changes may be needed later in
495  	      maxSelect(). We can therefore not move the first function after the (indp >= 0) check. */
496  	   iscoid = true;
497  	   int indc = maxDelta(max, maxabs,
498  	                       this->thesolver->coPvec(), this->thesolver->lcBound(), this->thesolver->ucBound(), 0, 1);
499  	   iscoid = false;
500  	   int indp = maxDelta(max, maxabs,
501  	                       this->thesolver->pVec(), this->thesolver->lpBound(), this->thesolver->upBound(), 0, 1);
502  	
503  	   if(indp >= 0)
504  	   {
505  	      nr = indp;
506  	      return this->thesolver->id(indp);
507  	   }
508  	
509  	   if(indc >= 0)
510  	   {
511  	      nr = indc;
512  	      return this->thesolver->coId(indc);
513  	   }
514  	
515  	   nr = -1;
516  	   return SPxId();
517  	}
518  	
519  	template <class R>
520  	SPxId SPxFastRT<R>::minDelta(
521  	   int& nr,
522  	   R& max,
523  	   R& maxabs)
524  	{
525  	   /* The following cause side effects on coPvec and pVec - both changes may be needed later in
526  	      minSelect(). We can therefore not move the first function after the (indp >= 0) check. */
527  	   iscoid = true;
528  	   const int indc = minDelta(max, maxabs,
529  	                             this->thesolver->coPvec(), this->thesolver->lcBound(), this->thesolver->ucBound(), 0, 1);
530  	   iscoid = false;
531  	   const int indp = minDelta(max, maxabs,
532  	                             this->thesolver->pVec(), this->thesolver->lpBound(), this->thesolver->upBound(), 0, 1);
533  	
534  	   if(indp >= 0)
535  	   {
536  	      nr = indp;
537  	      return this->thesolver->id(indp);
538  	   }
539  	
540  	   if(indc >= 0)
541  	   {
542  	      nr = indc;
543  	      return this->thesolver->coId(indc);
544  	   }
545  	
546  	   nr = -1;
547  	   return SPxId();
548  	}
549  	
550  	/* \p best returns the minimum update value such that the corresponding value of \p upd.delta() is
551  	 * at least \p stab and the update value is smaller than \p max. If no valid update value has been
552  	 * found \p bestDelta returns the slack to the bound corresponding to the index used for \p best. */
553  	template <class R>
554  	int SPxFastRT<R>::minSelect(
555  	   R& val,
556  	   R& stab,
557  	   R& best,
558  	   R& bestDelta,
559  	   R max,
560  	   const UpdateVector<R>& update,
561  	   const VectorBase<R>& lowBound,
562  	   const VectorBase<R>& upBound,
563  	   int start,
564  	   int incr) const
565  	{
566  	   int i;
567  	   R x, y;
568  	   bool leaving = this->m_type == SPxSolverBase<R>::LEAVE;
569  	   bool enterrowrep = !leaving && this->thesolver->theRep == SPxSolverBase<R>::ROW;
570  	
571  	   const R* up = upBound.get_const_ptr();
572  	   const R* low = lowBound.get_const_ptr();
573  	   const R* vec = update.get_const_ptr();
574  	   const R* upd = update.delta().values();
575  	   const int* idx = update.delta().indexMem();
576  	   const int* last = idx + update.delta().size();
577  	
578  	   int nr = -1;
579  	   int bestNr = -1;
580  	
581  	   for(idx += start; idx < last; idx += incr)
582  	   {
583  	      i = *idx;
584  	      x = upd[i];
585  	
586  	      // in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables
587  	      if(leaving && ((iscoid && this->thesolver->isCoBasic(i)) || (!iscoid
588  	                     && this->thesolver->isBasic(i))))
589  	         continue;
590  	
591  	      if(enterrowrep && this->thesolver->baseId(i).isSPxColId()
592  	            && this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(i))))
593  	            == SPxBasisBase<R>::Desc::P_FIXED)
594  	         continue;
595  	
596  	      if(x > stab)
597  	      {
598  	         y = (low[i] - vec[i]) / x;
599  	
600  	         if(y >= max)
601  	         {
602  	            val = y;
603  	            nr = i;
604  	            stab = x;
605  	         }
606  	         else if(y < best)
607  	         {
608  	            best = y;
609  	            bestNr = i;
610  	         }
611  	      }
612  	      else if(x < -stab)
613  	      {
614  	         y = (up[i] - vec[i]) / x;
615  	
616  	         if(y >= max)
617  	         {
618  	            val = y;
619  	            nr = i;
620  	            stab = -x;
621  	         }
622  	         else if(y < best)
623  	         {
624  	            best = y;
625  	            bestNr = i;
626  	         }
627  	      }
628  	   }
629  	
630  	   if(nr < 0 && bestNr > 0)
631  	   {
632  	      if(upd[bestNr] < 0)
633  	         bestDelta = up[bestNr] - vec[bestNr];
634  	      else
635  	         bestDelta = vec[bestNr] - low[bestNr];
636  	   }
637  	
638  	   return nr;
639  	}
640  	
641  	/* See minSelect() */
642  	template <class R>
643  	int SPxFastRT<R>::maxSelect(
644  	   R& val,
645  	   R& stab,
646  	   R& best,
647  	   R& bestDelta,
648  	   R max,
649  	   const UpdateVector<R>& update,
650  	   const VectorBase<R>& lowBound,
651  	   const VectorBase<R>& upBound,
652  	   int start,
653  	   int incr) const
654  	{
655  	   int i;
656  	   R x, y;
657  	   bool leaving = this->m_type == SPxSolverBase<R>::LEAVE;
658  	   bool enterrowrep = !leaving && this->thesolver->theRep == SPxSolverBase<R>::ROW;
659  	
660  	   const R* up = upBound.get_const_ptr();
661  	   const R* low = lowBound.get_const_ptr();
662  	   const R* vec = update.get_const_ptr();
663  	   const R* upd = update.delta().values();
664  	   const int* idx = update.delta().indexMem();
665  	   const int* last = idx + update.delta().size();
666  	
667  	   int nr = -1;
668  	   int bestNr = -1;
669  	
670  	   for(idx += start; idx < last; idx += incr)
671  	   {
672  	      i = *idx;
673  	      x = upd[i];
674  	
675  	      // in the dual algorithm, bound flips cannot happen, hence we only consider nonbasic variables
676  	      if(leaving && ((iscoid && this->thesolver->isCoBasic(i)) || (!iscoid
677  	                     && this->thesolver->isBasic(i))))
678  	         continue;
679  	
680  	      if(enterrowrep && this->thesolver->baseId(i).isSPxColId()
681  	            && this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(i))))
682  	            == SPxBasisBase<R>::Desc::P_FIXED)
683  	         continue;
684  	
685  	      if(x > stab)
686  	      {
687  	         y = (up[i] - vec[i]) / x;
688  	
689  	         if(y <= max)
690  	         {
691  	            val = y;
692  	            nr = i;
693  	            stab = x;
694  	         }
695  	         else if(y > best)
696  	         {
697  	            best = y;
698  	            bestNr = i;
699  	         }
700  	      }
701  	      else if(x < -stab)
702  	      {
703  	         y = (low[i] - vec[i]) / x;
704  	
705  	         if(y <= max)
706  	         {
707  	            val = y;
708  	            nr = i;
709  	            stab = -x;
710  	         }
711  	         else if(y > best)
712  	         {
713  	            best = y;
714  	            bestNr = i;
715  	         }
716  	      }
717  	   }
718  	
719  	   if(nr < 0 && bestNr > 0)
720  	   {
721  	      if(upd[bestNr] > 0)
722  	         bestDelta = up[bestNr] - vec[bestNr];
723  	      else
724  	         bestDelta = vec[bestNr] - low[bestNr];
725  	   }
726  	
727  	   return nr;
728  	}
729  	
730  	template <class R>
731  	int SPxFastRT<R>::maxSelect(
732  	   R& val,
733  	   R& stab,
734  	   R& bestDelta,
735  	   R max)
736  	{
737  	   R best = R(-infinity);
738  	   bestDelta = 0.0;
739  	   assert(this->m_type == SPxSolverBase<R>::ENTER);
740  	   return maxSelect(val, stab, best, bestDelta, max,
741  	                    this->thesolver->fVec(), this->thesolver->lbBound(), this->thesolver->ubBound(),  0, 1);
742  	}
743  	
744  	template <class R>
745  	SPxId SPxFastRT<R>::maxSelect(
746  	   int& nr,
747  	   R& val,
748  	   R& stab,
749  	   R& bestDelta,
750  	   R max
751  	)
752  	{
753  	   int indp, indc;
754  	   R best = R(-infinity);
755  	   bestDelta = 0.0;
756  	   iscoid = true;
757  	   indc = maxSelect(val, stab, best, bestDelta, max,
758  	                    this->thesolver->coPvec(), this->thesolver->lcBound(), this->thesolver->ucBound(), 0, 1);
759  	   iscoid = false;
760  	   indp = maxSelect(val, stab, best, bestDelta, max,
761  	                    this->thesolver->pVec(), this->thesolver->lpBound(), this->thesolver->upBound(), 0, 1);
762  	
763  	   if(indp >= 0)
764  	   {
765  	      nr = indp;
766  	      return this->thesolver->id(indp);
767  	   }
768  	
769  	   if(indc >= 0)
770  	   {
771  	      nr = indc;
772  	      return this->thesolver->coId(indc);
773  	   }
774  	
775  	   nr = -1;
776  	   return SPxId();
777  	}
778  	
779  	template <class R>
780  	int SPxFastRT<R>::minSelect(
781  	   R& val,
782  	   R& stab,
783  	   R& bestDelta,
784  	   R max)
785  	{
786  	   R best = R(infinity);
787  	   bestDelta = 0.0;
788  	   assert(this->m_type == SPxSolverBase<R>::ENTER);
789  	   return minSelect(val, stab, best, bestDelta, max,
790  	                    this->thesolver->fVec(), this->thesolver->lbBound(), this->thesolver->ubBound(), 0, 1);
791  	}
792  	
793  	template <class R>
794  	SPxId SPxFastRT<R>::minSelect(
795  	   int& nr,
796  	   R& val,
797  	   R& stab,
798  	   R& bestDelta,
799  	   R max)
800  	{
801  	   R best = R(infinity);
802  	   bestDelta = 0.0;
803  	   iscoid = true;
804  	   int indc = minSelect(val, stab, best, bestDelta, max,
805  	                        this->thesolver->coPvec(), this->thesolver->lcBound(), this->thesolver->ucBound(), 0, 1);
806  	   iscoid = false;
807  	   int indp = minSelect(val, stab, best, bestDelta, max,
808  	                        this->thesolver->pVec(), this->thesolver->lpBound(), this->thesolver->upBound(), 0, 1);
809  	
810  	   if(indp >= 0)
811  	   {
812  	      nr = indp;
813  	      return this->thesolver->id(indp);
814  	   }
815  	
816  	   if(indc >= 0)
817  	   {
818  	      nr = indc;
819  	      return this->thesolver->coId(indc);
820  	   }
821  	
822  	   nr = -1;
823  	   return SPxId();
824  	}
825  	
826  	template <class R>
827  	bool SPxFastRT<R>::maxShortLeave(R& sel, int leave, R maxabs)
828  	{
829  	   assert(leave >= 0);
830  	   assert(maxabs >= 0);
831  	
832  	   R shortval = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_SHORTVAL);
833  	
834  	   sel = this->thesolver->fVec().delta()[leave];
835  	
836  	   if(sel > maxabs * shortval)
837  	   {
838  	      sel = (this->thesolver->ubBound()[leave] - this->thesolver->fVec()[leave]) / sel;
839  	      return true;
840  	   }
841  	
842  	   if(sel < -maxabs * shortval)
843  	   {
844  	      sel = (this->thesolver->lbBound()[leave] - this->thesolver->fVec()[leave]) / sel;
845  	      return true;
846  	   }
847  	
848  	   return false;
849  	}
850  	
851  	template <class R>
852  	bool SPxFastRT<R>::minShortLeave(R& sel, int leave, R maxabs)
853  	{
854  	   assert(leave >= 0);
855  	   assert(maxabs >= 0);
856  	
857  	   R shortval = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_SHORTVAL);
858  	
859  	   sel = this->thesolver->fVec().delta()[leave];
860  	
861  	   if(sel > maxabs * shortval)
862  	   {
863  	      sel = (this->thesolver->lbBound()[leave] - this->thesolver->fVec()[leave]) / sel;
864  	      return true;
865  	   }
866  	
867  	   if(sel < -maxabs * shortval)
868  	   {
869  	      sel = (this->thesolver->ubBound()[leave] - this->thesolver->fVec()[leave]) / sel;
870  	      return true;
871  	   }
872  	
873  	   return false;
874  	}
875  	
876  	template <class R>
877  	bool SPxFastRT<R>::maxReLeave(R& sel, int leave, R maxabs, bool polish)
878  	{
879  	   UpdateVector<R>& vec = this->thesolver->fVec();
880  	   VectorBase<R>& low = this->thesolver->lbBound();
881  	   VectorBase<R>& up = this->thesolver->ubBound();
882  	
883  	   if(leave < 0)
884  	      return true;
885  	
886  	   if(up[leave] > low[leave])
887  	   {
888  	      R x = vec.delta()[leave];
889  	
890  	      if(sel < -fastDelta / maxabs)
891  	      {
892  	         sel = 0.0;
893  	
894  	         // prevent shifts in polishing mode to avoid a final cleanup step (i.e. simplex type switch)
895  	         if(!polish
896  	               && this->thesolver->dualStatus(this->thesolver->baseId(leave)) != SPxBasisBase<R>::Desc::D_ON_BOTH)
897  	         {
898  	            if(x < 0.0)
899  	               this->thesolver->shiftLBbound(leave, vec[leave]);
900  	            else
901  	               this->thesolver->shiftUBbound(leave, vec[leave]);
902  	         }
903  	      }
904  	   }
905  	   else
906  	   {
907  	      sel = 0.0;
908  	
909  	      // prevent shifts in polishing mode to avoid a final cleanup step (i.e. simplex type switch)
910  	      if(!polish)
911  	      {
912  	         this->thesolver->shiftLBbound(leave, vec[leave]);
913  	         this->thesolver->shiftUBbound(leave, vec[leave]);
914  	      }
915  	   }
916  	
917  	   return false;
918  	}
919  	
920  	template <class R>
921  	bool SPxFastRT<R>::minReLeave(R& sel, int leave, R maxabs, bool polish)
922  	{
923  	   UpdateVector<R>& vec = this->thesolver->fVec();
924  	   VectorBase<R>& low = this->thesolver->lbBound();
925  	   VectorBase<R>& up = this->thesolver->ubBound();
926  	
927  	   if(leave < 0)
928  	      return true;
929  	
930  	   if(up[leave] > low[leave])
931  	   {
932  	      R x = vec.delta()[leave];
933  	
934  	      if(sel > fastDelta / maxabs)
935  	      {
936  	         sel = 0.0;
937  	
938  	         if(!polish
939  	               && this->thesolver->dualStatus(this->thesolver->baseId(leave)) != SPxBasisBase<R>::Desc::D_ON_BOTH)
940  	            // prevent shifts in polishing mode to avoid a final cleanup step (i.e. simplex type switch)
941  	         {
942  	            if(x > 0.0)
943  	               this->thesolver->shiftLBbound(leave, vec[leave]);
944  	            else
945  	               this->thesolver->shiftUBbound(leave, vec[leave]);
946  	         }
947  	      }
948  	   }
949  	   else
950  	   {
951  	      sel = 0.0;
952  	
953  	      // prevent shifts in polishing mode to avoid a final cleanup step (i.e. simplex type switch)
954  	      if(!polish)
955  	      {
956  	         this->thesolver->shiftLBbound(leave, vec[leave]);
957  	         this->thesolver->shiftUBbound(leave, vec[leave]);
958  	      }
959  	   }
960  	
961  	   return false;
962  	}
963  	
964  	template <class R>
965  	int SPxFastRT<R>::selectLeave(R& val, R, bool polish)
966  	{
967  	   R maxabs, max, sel;
968  	   R delta_shift = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_DELTA_SHIFT);
969  	   int leave = -1;
970  	   int cnt = 0;
971  	
972  	   assert(this->m_type == SPxSolverBase<R>::ENTER);
973  	
974  	   // force instable pivot iff true (see explanation in enter.cpp and spxsolve.hpp)
975  	   bool instable = this->solver()->instableEnter;
976  	   R lowstab = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_LOWSTAB);
977  	   assert(!instable || this->solver()->instableEnterId.isValid());
978  	
979  	   resetTols();
980  	
981  	   if(val > epsilonZero())
982  	   {
983  	      do
984  	      {
985  	         // phase 1:
986  	         max = val;
987  	         maxabs = 0.0;
988  	         leave = maxDelta(max, maxabs);
989  	
990  	         assert(leave < 0 || !(this->thesolver->baseId(leave).isSPxColId()) ||
991  	                this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(
992  	                         leave)))) != SPxBasisBase<R>::Desc::P_FIXED);
993  	
994  	         if(max == val || leave == -1)
995  	         {
996  	            assert(max == val && leave == -1);
997  	            return -1;
998  	         }
999  	
1000 	         if(!maxShortLeave(sel, leave, maxabs))
1001 	         {
1002 	            // phase 2:
1003 	            R stab, bestDelta;
1004 	
1005 	            stab = 100.0 * minStability(maxabs);
1006 	
1007 	            // force instable pivot iff instable is true (see explanation in enter.hpp and spxsolve.hpp)
1008 	            if(instable)
1009 	               leave = maxSelect(sel, lowstab, bestDelta, max);
1010 	            else
1011 	               leave = maxSelect(sel, stab, bestDelta, max);
1012 	
1013 	            if(bestDelta < delta_shift * SOPLEX_TRIES)
1014 	               cnt++;
1015 	            else
1016 	               cnt += SOPLEX_TRIES;
1017 	         }
1018 	
1019 	         if(!maxReLeave(sel, leave, maxabs, polish))
1020 	            break;
1021 	
1022 	         relax();
1023 	      }
1024 	      while(cnt < SOPLEX_TRIES);
1025 	   }
1026 	   else if(val < -epsilonZero())
1027 	   {
1028 	      do
1029 	      {
1030 	         max = val;
1031 	         maxabs = 0;
1032 	         leave = minDelta(max, maxabs);
1033 	
1034 	         assert(leave < 0 || !(this->thesolver->baseId(leave).isSPxColId()) ||
1035 	                this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(
1036 	                         leave)))) != SPxBasisBase<R>::Desc::P_FIXED);
1037 	
1038 	         if(max == val || leave == -1)
1039 	         {
1040 	            assert(max == val && leave == -1);
1041 	            return -1;
1042 	         }
1043 	
1044 	         if(!minShortLeave(sel, leave, maxabs))
1045 	         {
1046 	            // phase 2:
1047 	            R stab, bestDelta;
1048 	
1049 	            stab = 100.0 * minStability(maxabs);
1050 	
1051 	            // force instable pivot iff instable is true (see explanation in enter.hpp and spxsolve.hpp)
1052 	            if(instable)
1053 	               leave = minSelect(sel, lowstab, bestDelta, max);
1054 	            else
1055 	               leave = minSelect(sel, stab, bestDelta, max);
1056 	
1057 	            assert(leave < 0 || !(this->thesolver->baseId(leave).isSPxColId())
1058 	                   || this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(
1059 	                            leave)))) != SPxBasisBase<R>::Desc::P_FIXED);
1060 	
1061 	            if(bestDelta < delta_shift * SOPLEX_TRIES)
1062 	               cnt++;
1063 	            else
1064 	               cnt += SOPLEX_TRIES;
1065 	         }
1066 	
1067 	         if(!minReLeave(sel, leave, maxabs, polish))
1068 	            break;
1069 	
1070 	         relax();
1071 	      }
1072 	      while(cnt < SOPLEX_TRIES);
1073 	   }
1074 	   else
1075 	      return -1;
1076 	
1077 	   SPX_DEBUG(
1078 	
1079 	      if(leave >= 0)
1080 	      std::cout
1081 	      << "DFSTRT01 "
1082 	      << this->thesolver->basis().iteration() << "("
1083 	      << std::setprecision(6) << this->thesolver->value() << ","
1084 	      << std::setprecision(2) << this->thesolver->basis().stability() << "):"
1085 	      << leave << "\t"
1086 	      << std::setprecision(4) << sel << " "
1087 	      << std::setprecision(4) << this->thesolver->fVec().delta()[leave] << " "
1088 	      << std::setprecision(6) << maxabs
1089 	      << std::endl;
1090 	      else
1091 	         std::cout << "DFSTRT02 " << this->thesolver->basis().iteration()
1092 	         << ": skipping instable pivot" << std::endl;
1093 	      )
1094 	
1095 	         if(polish && leave >= 0)
1096 	         {
1097 	            assert(this->thesolver->rep() == SPxSolverBase<R>::COLUMN);
1098 	            SPxId leaveId = this->thesolver->baseId(leave);
1099 	
1100 	            // decide whether the chosen leave index contributes to the polishing objective
1101 	            if(this->thesolver->polishObj == SPxSolverBase<R>::POLISH_INTEGRALITY)
1102 	            {
1103 	               // only allow (integer) variables to leave the basis
1104 	               if(leaveId.isSPxRowId())
1105 	                  return -1;
1106 	               else if(this->thesolver->integerVariables.size() == this->thesolver->nCols())
1107 	               {
1108 	                  if(leaveId.isSPxColId() && this->thesolver->integerVariables[this->thesolver->number(leaveId)] == 0)
1109 	                     return -1;
1110 	               }
1111 	            }
1112 	            else if(this->thesolver->polishObj == SPxSolverBase<R>::POLISH_FRACTIONALITY)
1113 	            {
1114 	               // only allow slacks and continuous variables to leave the basis
1115 	               if(this->thesolver->integerVariables.size() == this->thesolver->nCols())
1116 	               {
1117 	                  if(this->thesolver->baseId(leave).isSPxColId()
1118 	                        && this->thesolver->integerVariables[this->thesolver->number(leaveId)] == 1)
1119 	                     return -1;
1120 	               }
1121 	               else if(this->thesolver->baseId(leave).isSPxColId())
1122 	                  return -1;
1123 	            }
1124 	         }
1125 	
1126 	   if(leave >= 0 || minStab > 2 * this->solver()->epsilon())
1127 	   {
1128 	      val = sel;
1129 	
1130 	      if(leave >= 0)
1131 	         tighten();
1132 	   }
1133 	
1134 	   assert(leave < 0 || !(this->thesolver->baseId(leave).isSPxColId())
1135 	          || this->thesolver->desc().colStatus(this->thesolver->number(SPxColId(this->thesolver->baseId(
1136 	                   leave)))) != SPxBasisBase<R>::Desc::P_FIXED);
1137 	
1138 	   return leave;
1139 	}
1140 	
1141 	
1142 	template <class R>
1143 	bool SPxFastRT<R>::maxReEnter(R& sel,
1144 	                              R maxabs,
1145 	                              const SPxId& id,
1146 	                              int nr,
1147 	                              bool polish)
1148 	{
1149 	   R x, d;
1150 	   VectorBase<R>* up;
1151 	   VectorBase<R>* low;
1152 	
1153 	   UpdateVector<R>& pvec = this->thesolver->pVec();
1154 	   SSVectorBase<R>& pupd = this->thesolver->pVec().delta();
1155 	   VectorBase<R>& upb = this->thesolver->upBound();
1156 	   VectorBase<R>& lpb = this->thesolver->lpBound();
1157 	   UpdateVector<R>& cvec = this->thesolver->coPvec();
1158 	   SSVectorBase<R>& cupd = this->thesolver->coPvec().delta();
1159 	   VectorBase<R>& ucb = this->thesolver->ucBound();
1160 	   VectorBase<R>& lcb = this->thesolver->lcBound();
1161 	
1162 	   if(this->thesolver->isCoId(id))
1163 	   {
1164 	      if(this->thesolver->isCoBasic(nr))
1165 	      {
1166 	         cupd.clearIdx(nr);
1167 	         return true;
1168 	      }
1169 	
1170 	      x = cvec[nr];
1171 	      d = cupd[nr];
1172 	      up = &ucb;
1173 	      low = &lcb;
1174 	
1175 	      if(d < 0.0)
1176 	         sel = (lcb[nr] - cvec[nr]) / d;
1177 	      else
1178 	         sel = (ucb[nr] - cvec[nr]) / d;
1179 	   }
1180 	   else if(this->thesolver->isId(id))
1181 	   {
1182 	      pvec[nr] = this->thesolver->vector(nr) * cvec;
1183 	
1184 	      if(this->thesolver->isBasic(nr))
1185 	      {
1186 	         pupd.clearIdx(nr);
1187 	         return true;
1188 	      }
1189 	
1190 	      x = pvec[nr];
1191 	      d = pupd[nr];
1192 	      up = &upb;
1193 	      low = &lpb;
1194 	
1195 	      if(d < 0.0)
1196 	         sel = (lpb[nr] - pvec[nr]) / d;
1197 	      else
1198 	         sel = (upb[nr] - pvec[nr]) / d;
1199 	   }
1200 	   else
1201 	      return true;
1202 	
1203 	   if((*up)[nr] != (*low)[nr])
1204 	   {
1205 	      if(sel < -fastDelta / maxabs)
1206 	      {
1207 	         sel = 0.0;
1208 	
1209 	         // prevent shifts in polishing mode to avoid a final cleanup step (i.e. simplex type switch)
1210 	         if(!polish)
1211 	         {
1212 	            if(d > 0.0)
1213 	            {
1214 	               this->thesolver->theShift -= (*up)[nr];
1215 	               (*up)[nr] = x;
1216 	               this->thesolver->theShift += (*up)[nr];
1217 	            }
1218 	            else
1219 	            {
1220 	               this->thesolver->theShift += (*low)[nr];
1221 	               (*low)[nr] = x;
1222 	               this->thesolver->theShift -= (*low)[nr];
1223 	            }
1224 	         }
1225 	      }
1226 	   }
1227 	   else
1228 	   {
1229 	      sel = 0.0;
1230 	
1231 	      // prevent shifts in polishing mode to avoid a final cleanup step (i.e. simplex type switch)
1232 	      if(!polish)
1233 	      {
1234 	         if(x > (*up)[nr])
1235 	            this->thesolver->theShift += x - (*up)[nr];
1236 	         else
1237 	            this->thesolver->theShift += (*low)[nr] - x;
1238 	
1239 	         (*up)[nr] = (*low)[nr] = x;
1240 	      }
1241 	   }
1242 	
1243 	   return false;
1244 	}
1245 	
1246 	template <class R>
1247 	bool SPxFastRT<R>::minReEnter(
1248 	   R& sel,
1249 	   R maxabs,
1250 	   const SPxId& id,
1251 	   int nr,
1252 	   bool polish)
1253 	{
1254 	   R x, d;
1255 	   VectorBase<R>* up;
1256 	   VectorBase<R>* low;
1257 	
1258 	   UpdateVector<R>& pvec = this->thesolver->pVec();
1259 	   SSVectorBase<R>& pupd = this->thesolver->pVec().delta();
1260 	   VectorBase<R>& upb = this->thesolver->upBound();
1261 	   VectorBase<R>& lpb = this->thesolver->lpBound();
1262 	   UpdateVector<R>& cvec = this->thesolver->coPvec();
1263 	   SSVectorBase<R>& cupd = this->thesolver->coPvec().delta();
1264 	   VectorBase<R>& ucb = this->thesolver->ucBound();
1265 	   VectorBase<R>& lcb = this->thesolver->lcBound();
1266 	
1267 	   if(this->thesolver->isCoId(id))
1268 	   {
1269 	      if(this->thesolver->isCoBasic(nr))
1270 	      {
1271 	         cupd.clearIdx(nr);
1272 	         return true;
1273 	      }
1274 	
1275 	      x = cvec[nr];
1276 	      d = cupd[nr];
1277 	      up = &ucb;
1278 	      low = &lcb;
1279 	
1280 	      if(d > 0.0)
1281 	         sel = (this->thesolver->lcBound()[nr] - cvec[nr]) / d;
1282 	      else
1283 	         sel = (this->thesolver->ucBound()[nr] - cvec[nr]) / d;
1284 	   }
1285 	
1286 	   else if(this->thesolver->isId(id))
1287 	   {
1288 	      pvec[nr] = this->thesolver->vector(nr) * cvec;
1289 	
1290 	      if(this->thesolver->isBasic(nr))
1291 	      {
1292 	         pupd.clearIdx(nr);
1293 	         return true;
1294 	      }
1295 	
1296 	      x = pvec[nr];
1297 	      d = pupd[nr];
1298 	      up = &upb;
1299 	      low = &lpb;
1300 	
1301 	      if(d > 0.0)
1302 	         sel = (this->thesolver->lpBound()[nr] - pvec[nr]) / d;
1303 	      else
1304 	         sel = (this->thesolver->upBound()[nr] - pvec[nr]) / d;
1305 	   }
1306 	
1307 	   else
1308 	      return true;
1309 	
1310 	   if((*up)[nr] != (*low)[nr])
1311 	   {
1312 	      if(sel > fastDelta / maxabs)
1313 	      {
1314 	         sel = 0.0;
1315 	
1316 	         // prevent shifts in polishing mode to avoid a final cleanup step (i.e. simplex type switch)
1317 	         if(!polish)
1318 	         {
1319 	            if(d < 0.0)
1320 	            {
1321 	               this->thesolver->theShift -= (*up)[nr];
1322 	               (*up)[nr] = x;
1323 	               this->thesolver->theShift += (*up)[nr];
1324 	            }
1325 	            else
1326 	            {
1327 	               this->thesolver->theShift += (*low)[nr];
1328 	               (*low)[nr] = x;
1329 	               this->thesolver->theShift -= (*low)[nr];
1330 	            }
1331 	         }
1332 	      }
1333 	   }
1334 	   else
1335 	   {
1336 	      sel = 0.0;
1337 	
1338 	      // prevent shifts in polishing mode to avoid a final cleanup step (i.e. simplex type switch)
1339 	      if(!polish)
1340 	      {
1341 	         if(x > (*up)[nr])
1342 	            this->thesolver->theShift += x - (*up)[nr];
1343 	         else
1344 	            this->thesolver->theShift += (*low)[nr] - x;
1345 	
1346 	         (*up)[nr] = (*low)[nr] = x;
1347 	      }
1348 	   }
1349 	
1350 	   return false;
1351 	}
1352 	
1353 	template <class R>
1354 	bool SPxFastRT<R>::shortEnter(
1355 	   const SPxId& enterId,
1356 	   int nr,
1357 	   R max,
1358 	   R maxabs) const
1359 	{
1360 	   R shortval = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_SHORTVAL);
1361 	
1362 	   if(this->thesolver->isCoId(enterId))
1363 	   {
1364 	      if(max != 0.0)
1365 	      {
1366 	         R x = this->thesolver->coPvec().delta()[nr];
1367 	
1368 	         if(x < maxabs * shortval && -x < maxabs * shortval)
1369 	            return false;
1370 	      }
1371 	
1372 	      return true;
1373 	   }
1374 	   else if(this->thesolver->isId(enterId))
1375 	   {
1376 	      if(max != 0.0)
1377 	      {
1378 	         R x = this->thesolver->pVec().delta()[nr];
1379 	
1380 	         if(x < maxabs * shortval && -x < maxabs * shortval)
1381 	            return false;
1382 	      }
1383 	
1384 	      return true;
1385 	   }
1386 	
1387 	   return false;
1388 	}
1389 	
1390 	template <class R>
1391 	SPxId SPxFastRT<R>::selectEnter(R& val, int, bool polish)
1392 	{
1393 	   SPxId enterId;
1394 	   R max, sel;
1395 	   R maxabs = 0.0;
1396 	   R delta_shift = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_DELTA_SHIFT);
1397 	   int nr;
1398 	   int cnt = 0;
1399 	
1400 	   assert(this->m_type == SPxSolverBase<R>::LEAVE);
1401 	
1402 	   resetTols();
1403 	
1404 	   // force instable pivot iff true (see explanation in leave.hpp and spxsolve.hpp)
1405 	   bool instable = this->solver()->instableLeave;
1406 	   R lowstab = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_LOWSTAB);
1407 	   assert(!instable || this->solver()->instableLeaveNum >= 0);
1408 	
1409 	   sel = 0.0;
1410 	
1411 	   if(val > epsilonZero())
1412 	   {
1413 	      do
1414 	      {
1415 	         maxabs = 0.0;
1416 	         max = val;
1417 	
1418 	         enterId = maxDelta(nr, max, maxabs);
1419 	
1420 	         if(!enterId.isValid())
1421 	            return enterId;
1422 	
1423 	         assert(max >= 0.0);
1424 	         assert(!enterId.isValid() || !this->solver()->isBasic(enterId));
1425 	
1426 	         if(!shortEnter(enterId, nr, max, maxabs))
1427 	         {
1428 	            R bestDelta, stab;
1429 	
1430 	            stab = minStability(maxabs);
1431 	
1432 	            // force instable pivot iff instable is true (see explanation in leave.hpp and spxsolve.hpp)
1433 	            if(instable)
1434 	            {
1435 	               enterId = maxSelect(nr, sel, lowstab, bestDelta, max);
1436 	            }
1437 	            else
1438 	            {
1439 	               enterId = maxSelect(nr, sel, stab, bestDelta, max);
1440 	            }
1441 	
1442 	            if(bestDelta < delta_shift * SOPLEX_TRIES)
1443 	               cnt++;
1444 	            else
1445 	               cnt += SOPLEX_TRIES;
1446 	         }
1447 	
1448 	         if(!maxReEnter(sel, maxabs, enterId, nr, polish))
1449 	            break;
1450 	
1451 	         relax();
1452 	      }
1453 	      while(cnt < SOPLEX_TRIES);
1454 	   }
1455 	   else if(val < -epsilonZero())
1456 	   {
1457 	      do
1458 	      {
1459 	         maxabs = 0.0;
1460 	         max = val;
1461 	         enterId = minDelta(nr, max, maxabs);
1462 	
1463 	         if(!enterId.isValid())
1464 	            return enterId;
1465 	
1466 	         assert(max <= 0.0);
1467 	         assert(!enterId.isValid() || !this->solver()->isBasic(enterId));
1468 	
1469 	         if(!shortEnter(enterId, nr, max, maxabs))
1470 	         {
1471 	            R bestDelta, stab;
1472 	
1473 	            stab = minStability(maxabs);
1474 	
1475 	            // force instable pivot iff instable is true (see explanation in leave.hpp and spxsolve.hpp)
1476 	            if(instable)
1477 	            {
1478 	               enterId = minSelect(nr, sel, lowstab, bestDelta, max);
1479 	            }
1480 	            else
1481 	            {
1482 	               enterId = minSelect(nr, sel, stab, bestDelta, max);
1483 	            }
1484 	
1485 	            if(bestDelta < delta_shift * SOPLEX_TRIES)
1486 	               cnt++;
1487 	            else
1488 	               cnt += SOPLEX_TRIES;
1489 	         }
1490 	
1491 	         if(!minReEnter(sel, maxabs, enterId, nr, polish))
1492 	            break;
1493 	
1494 	         relax();
1495 	      }
1496 	      while(cnt < SOPLEX_TRIES);
1497 	   }
1498 	
1499 	   SPX_DEBUG(
1500 	
1501 	      if(enterId.isValid())
1502 	{
1503 	   assert(!enterId.isValid() || !this->solver()->isBasic(enterId));
1504 	
1505 	      R x;
1506 	
1507 	      if(this->thesolver->isCoId(enterId))
1508 	         x = this->thesolver->coPvec().delta()[ this->thesolver->number(enterId) ];
1509 	      else
1510 	         x = this->thesolver->pVec().delta()[ this->thesolver->number(enterId) ];
1511 	
1512 	      std::cout << "DFSTRT03 " << this->thesolver->basis().iteration() << ": "
1513 	                << sel << '\t' << x << " (" << maxabs << ")" << std::endl;
1514 	   }
1515 	   else
1516 	      std::cout << "DFSTRT04 " << this->thesolver->basis().iteration()
1517 	                << ": skipping instable pivot" << std::endl;
1518 	      )
1519 	
1520 	      if(polish && enterId.isValid())
1521 	      {
1522 	         assert(this->thesolver->rep() == SPxSolverBase<R>::ROW);
1523 	
1524 	         // decide whether the chosen entering index contributes to the polishing objective
1525 	         if(this->thesolver->polishObj == SPxSolverBase<R>::POLISH_INTEGRALITY)
1526 	         {
1527 	            // only allow (integer) variables to enter the basis
1528 	            if(enterId.isSPxRowId())
1529 	               return SPxId();
1530 	            else if(this->thesolver->integerVariables.size() == this->thesolver->nCols()
1531 	                    && this->thesolver->integerVariables[this->thesolver->number(enterId)] == 0)
1532 	               return SPxId();
1533 	         }
1534 	         else if(this->thesolver->polishObj == SPxSolverBase<R>::POLISH_FRACTIONALITY)
1535 	         {
1536 	            // only allow slacks and continuous variables to enter the basis
1537 	            if(this->thesolver->integerVariables.size() == this->thesolver->nCols())
1538 	            {
1539 	               if(enterId.isSPxColId() && this->thesolver->integerVariables[this->thesolver->number(enterId)] == 1)
1540 	                  return SPxId();
1541 	            }
1542 	            else if(enterId.isSPxColId())
1543 	               return SPxId();
1544 	         }
1545 	      }
1546 	
1547 	
1548 	   if(enterId.isValid() || minStab > 2 * epsilonZero())
1549 	   {
1550 	      val = sel;
1551 	
1552 	      if(enterId.isValid())
1553 	         tighten();
1554 	   }
1555 	
1556 	   assert(!enterId.isValid() || !this->solver()->isBasic(enterId));
1557 	
1558 	   return enterId;
1559 	}
1560 	
1561 	template <class R>
1562 	void SPxFastRT<R>::load(SPxSolverBase<R>* spx)
1563 	{
1564 	   this->thesolver = spx;
1565 	   setType(spx->type());
1566 	}
1567 	
1568 	template <class R>
1569 	void SPxFastRT<R>::setType(typename SPxSolverBase<R>::Type type)
1570 	{
1571 	   this->m_type = type;
1572 	
1573 	   minStab = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_MINSTAB);
1574 	   fastDelta = this->delta;
1575 	}
1576 	} // namespace soplex
1577