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   	namespace soplex
26   	{
27   	/**@todo suspicious: *max is not set, but it is used
28   	 * (with the default setting *max=1) in selectLeave and selectEnter
29   	 * The question might be if max shouldn't be updated with themax?
30   	 *
31   	 * numCycle and maxCycle are integers. So degeneps will be
32   	 * exactly delta until numCycle >= maxCycle. Then it will be
33   	 * 0 until numCycle >= 2 * maxCycle, after wich it becomes
34   	 * negative. This does not look ok.
35   	 */
36   	template <class R>
37   	R SPxHarrisRT<R>::degenerateEps() const
38   	{
39   	   return this->solver()->delta()
40   	          * (1.0 - this->solver()->numCycle() / this->solver()->maxCycle());
41   	}
42   	
43   	template <class R>
44   	int SPxHarrisRT<R>::maxDelta(
45   	   R* /*max*/,             /* max abs value in upd */
46   	   R* val,             /* initial and chosen value */
47   	   int num,             /* # of indices in idx */
48   	   const int* idx,             /* nonzero indices in upd */
49   	   const R* upd,             /* update VectorBase<R> for vec */
50   	   const R* vec,             /* current VectorBase<R> */
51   	   const R* low,             /* lower bounds for vec */
52   	   const R* up               /* upper bounds for vec */
53   	)  const
54   	{
55   	   R x;
56   	   R theval;
57   	   /**@todo patch suggests using *max instead of themax */
58   	   R themax;
59   	   int sel;
60   	   int i;
61   	   R epsilon = this->tolerances()->epsilon();
62   	
63   	   assert(*val >= 0);
64   	
65   	   theval = *val;
66   	   themax = 0;
67   	   sel = -1;
68   	
69   	   while(num--)
70   	   {
71   	      i = idx[num];
72   	      x = upd[i];
73   	
74   	      if(x > epsilon)
75   	      {
76   	         themax = (x > themax) ? x : themax;
77   	         x = (up[i] - vec[i] + this->delta) / x;
78   	
79   	         if(x < theval && up[i] < R(infinity))
80   	            theval = x;
81   	      }
82   	      else if(x < -epsilon)
83   	      {
84   	         themax = (-x > themax) ? -x : themax;
85   	         x = (low[i] - vec[i] - this->delta) / x;
86   	
87   	         if(x < theval && low[i] > R(-infinity))
88   	            theval = x;
89   	      }
90   	   }
91   	
92   	   *val = theval;
93   	   return sel;
94   	}
95   	
96   	/**@todo suspicious: *max is not set, but it is used
97   	   (with the default setting *max=1)
98   	   in selectLeave and selectEnter
99   	*/
100  	template <class R>
101  	int SPxHarrisRT<R>::minDelta(
102  	   R* /*max*/,             /* max abs value in upd */
103  	   R* val,             /* initial and chosen value */
104  	   int num,             /* # of indices in idx */
105  	   const int* idx,             /* nonzero indices in upd */
106  	   const R* upd,             /* update VectorBase<R> for vec */
107  	   const R* vec,             /* current VectorBase<R> */
108  	   const R* low,             /* lower bounds for vec */
109  	   const R* up               /* upper bounds for vec */
110  	) const
111  	{
112  	   R x;
113  	   R theval;
114  	   /**@todo patch suggests using *max instead of themax */
115  	   R themax;
116  	   int sel;
117  	   int i;
118  	   R epsilon = this->tolerances()->epsilon();
119  	
120  	   assert(*val < 0);
121  	
122  	   theval = *val;
123  	   themax = 0;
124  	   sel = -1;
125  	
126  	   while(num--)
127  	   {
128  	      i = idx[num];
129  	      x = upd[i];
130  	
131  	      if(x > epsilon)
132  	      {
133  	         themax = (x > themax) ? x : themax;
134  	         x = (low[i] - vec[i] - this->delta) / x;
135  	
136  	         if(x > theval && low[i] > R(-infinity))
137  	            theval = x;
138  	      }
139  	      else if(x < -epsilon)
140  	      {
141  	         themax = (-x > themax) ? -x : themax;
142  	         x = (up[i] - vec[i] + this->delta) / x;
143  	
144  	         if(x > theval && up[i] < R(infinity))
145  	            theval = x;
146  	      }
147  	   }
148  	
149  	   *val = theval;
150  	   return sel;
151  	}
152  	
153  	/**
154  	   Here comes our implementation of the Harris procedure improved by shifting
155  	   bounds. The basic idea is to used the tolerated infeasibility within
156  	   solver()->entertol() for searching numerically stable pivots.
157  	
158  	   The algorithms operates in two phases. In a first phase, the maximum \p val
159  	   is determined, when infeasibility within solver()->entertol() is allowed. In the second
160  	   phase, between all variables with values < \p val the one is selected which
161  	   gives the best step forward in the simplex iteration. However, this may not
162  	   allways yield an improvement. In that case, we shift the variable toward
163  	   infeasibility and retry. This avoids cycling in the shifted LP.
164  	*/
165  	template <class R>
166  	int SPxHarrisRT<R>::selectLeave(R& val, R, bool)
167  	{
168  	   int i, j;
169  	   R stab, x, y;
170  	   R max;
171  	   R sel;
172  	   R lastshift;
173  	   R useeps;
174  	   int leave = -1;
175  	   R maxabs = 1;
176  	
177  	   R epsilon  = this->solver()->epsilon();
178  	   R degeneps = degenerateEps();
179  	
180  	   SSVectorBase<R>& upd = this->solver()->fVec().delta();
181  	   VectorBase<R>& vec = this->solver()->fVec();
182  	
183  	   const VectorBase<R>& up = this->solver()->ubBound();
184  	   const VectorBase<R>& low = this->solver()->lbBound();
185  	
186  	   assert(this->delta > epsilon);
187  	   assert(epsilon > 0);
188  	   assert(this->solver()->maxCycle() > 0);
189  	
190  	   max = val;
191  	   lastshift = this->solver()->shift();
192  	
193  	   this->solver()->fVec().delta().setup();
194  	
195  	   if(max > epsilon)
196  	   {
197  	      // phase 1:
198  	      maxDelta(
199  	         &maxabs,             /* max abs value in upd */
200  	         &max,                /* initial and chosen value */
201  	         upd.size(),          /* # of indices in upd */
202  	         upd.indexMem(),      /* nonzero indices in upd */
203  	         upd.values(),        /* update VectorBase<R> for vec */
204  	         vec.get_const_ptr(),         /* current VectorBase<R> */
205  	         low.get_const_ptr(),                 /* lower bounds for vec */
206  	         up.get_const_ptr()                   /* upper bounds for vec */
207  	      );
208  	
209  	      if(max == val)
210  	         return -1;
211  	
212  	
213  	      // phase 2:
214  	      stab = 0;
215  	      sel = R(-infinity);
216  	      useeps = maxabs * epsilon * 0.001;
217  	
218  	      if(useeps < epsilon)
219  	         useeps = epsilon;
220  	
221  	      for(j = upd.size() - 1; j >= 0; --j)
222  	      {
223  	         i = upd.index(j);
224  	         x = upd[i];
225  	
226  	         if(x > useeps)
227  	         {
228  	            y = up[i] - vec[i];
229  	
230  	            if(y < -degeneps)
231  	               this->solver()->shiftUBbound(i, vec[i]); // ensure simplex improvement
232  	            else
233  	            {
234  	               y /= x;
235  	
236  	               if(y <= max && y > sel - epsilon && x > stab)
237  	               {
238  	                  sel = y;
239  	                  leave = i;
240  	                  stab = x;
241  	               }
242  	            }
243  	         }
244  	         else if(x < -useeps)
245  	         {
246  	            y = low[i] - vec[i];
247  	
248  	            if(y > degeneps)
249  	               this->solver()->shiftLBbound(i, vec[i]); // ensure simplex improvement
250  	            else
251  	            {
252  	               y /= x;
253  	
254  	               if(y <= max && y > sel - epsilon && -x > stab)
255  	               {
256  	                  sel = y;
257  	                  leave = i;
258  	                  stab = -x;
259  	               }
260  	            }
261  	         }
262  	         else
263  	            upd.clearNum(j);
264  	      }
265  	   }
266  	
267  	
268  	   else if(max < -epsilon)
269  	   {
270  	      // phase 1:
271  	      minDelta(
272  	         &maxabs,             /* max abs value in upd */
273  	         &max,                /* initial and chosen value */
274  	         upd.size(),          /* # of indices in upd */
275  	         upd.indexMem(),      /* nonzero indices in upd */
276  	         upd.values(),        /* update VectorBase<R> for vec */
277  	         vec.get_const_ptr(),                 /* current VectorBase<R> */
278  	         low.get_const_ptr(),                 /* lower bounds for vec */
279  	         up.get_const_ptr()                   /* upper bounds for vec */
280  	      );
281  	
282  	      if(max == val)
283  	         return -1;
284  	
285  	      // phase 2:
286  	      stab = 0;
287  	      sel = R(infinity);
288  	      useeps = maxabs * epsilon * 0.001;
289  	
290  	      if(useeps < epsilon)
291  	         useeps = epsilon;
292  	
293  	      for(j = upd.size() - 1; j >= 0; --j)
294  	      {
295  	         i = upd.index(j);
296  	         x = upd[i];
297  	
298  	         if(x < -useeps)
299  	         {
300  	            y = up[i] - vec[i];
301  	
302  	            if(y < -degeneps)
303  	               this->solver()->shiftUBbound(i, vec[i]);   // ensure simplex improvement
304  	            else
305  	            {
306  	               y /= x;
307  	
308  	               if(y >= max && y < sel + epsilon && -x > stab)
309  	               {
310  	                  sel = y;
311  	                  leave = i;
312  	                  stab = -x;
313  	               }
314  	            }
315  	         }
316  	         else if(x > useeps)
317  	         {
318  	            y = low[i] - vec[i];
319  	
320  	            if(y > degeneps)
321  	               this->solver()->shiftLBbound(i, vec[i]); // ensure simplex improvement
322  	            else
323  	            {
324  	               y /= x;
325  	
326  	               if(y >= max && y < sel + epsilon && x > stab)
327  	               {
328  	                  sel = y;
329  	                  leave = i;
330  	                  stab = x;
331  	               }
332  	            }
333  	         }
334  	         else
335  	            upd.clearNum(j);
336  	      }
337  	   }
338  	
339  	   else
340  	      return -1;
341  	
342  	
343  	   if(lastshift != this->solver()->shift())
344  	      return selectLeave(val, 0, false);
345  	
346  	   assert(leave >= 0);
347  	
348  	   val = sel;
349  	   return leave;
350  	}
351  	
352  	template <class R>
353  	SPxId SPxHarrisRT<R>::selectEnter(R& val, int, bool)
354  	{
355  	   int i, j;
356  	   SPxId enterId;
357  	   R stab, x, y;
358  	   R max = 0.0;
359  	   R sel = 0.0;
360  	   R lastshift;
361  	   R cuseeps;
362  	   R ruseeps;
363  	   R cmaxabs = 1;
364  	   R rmaxabs = 1;
365  	   int pnr, cnr;
366  	
367  	   R minStability = this->tolerances()->scaleAccordingToEpsilon(1e-4);
368  	   R epsilon      = this->solver()->epsilon();
369  	   R degeneps     = degenerateEps();
370  	
371  	   VectorBase<R>& pvec = this->solver()->pVec();
372  	   SSVectorBase<R>& pupd = this->solver()->pVec().delta();
373  	
374  	   VectorBase<R>& cvec = this->solver()->coPvec();
375  	   SSVectorBase<R>& cupd = this->solver()->coPvec().delta();
376  	
377  	   const VectorBase<R>& upb = this->solver()->upBound();
378  	   const VectorBase<R>& lpb = this->solver()->lpBound();
379  	   const VectorBase<R>& ucb = this->solver()->ucBound();
380  	   const VectorBase<R>& lcb = this->solver()->lcBound();
381  	
382  	   assert(this->delta > epsilon);
383  	   assert(epsilon > 0);
384  	   assert(this->solver()->maxCycle() > 0);
385  	
386  	   this->solver()->coPvec().delta().setup();
387  	   this->solver()->pVec().delta().setup();
388  	
389  	   if(val > epsilon)
390  	   {
391  	      for(;;)
392  	      {
393  	         pnr = -1;
394  	         cnr = -1;
395  	         max = val;
396  	         lastshift = this->solver()->shift();
397  	         assert(this->delta > epsilon);
398  	
399  	         // phase 1:
400  	         maxDelta(
401  	            &rmaxabs,            /* max abs value in upd */
402  	            &max,                /* initial and chosen value */
403  	            pupd.size(),         /* # of indices in pupd */
404  	            pupd.indexMem(),     /* nonzero indices in pupd */
405  	            pupd.values(),       /* update VectorBase<R> for vec */
406  	            pvec.get_const_ptr(),                /* current VectorBase<R> */
407  	            lpb.get_const_ptr(),                 /* lower bounds for vec */
408  	            upb.get_const_ptr()                  /* upper bounds for vec */
409  	         );
410  	
411  	         maxDelta(
412  	            &cmaxabs,            /* max abs value in upd */
413  	            &max,                /* initial and chosen value */
414  	            cupd.size(),         /* # of indices in cupd */
415  	            cupd.indexMem(),     /* nonzero indices in cupd */
416  	            cupd.values(),       /* update VectorBase<R> for vec */
417  	            cvec.get_const_ptr(),                /* current VectorBase<R> */
418  	            lcb.get_const_ptr(),                 /* lower bounds for vec */
419  	            ucb.get_const_ptr()                  /* upper bounds for vec */
420  	         );
421  	
422  	         if(max == val)
423  	            return enterId;
424  	
425  	
426  	         // phase 2:
427  	         stab = 0;
428  	         sel = R(-infinity);
429  	         ruseeps = rmaxabs * 0.001 * epsilon;
430  	
431  	         if(ruseeps < epsilon)
432  	            ruseeps = epsilon;
433  	
434  	         cuseeps = cmaxabs * 0.001 * epsilon;
435  	
436  	         if(cuseeps < epsilon)
437  	            cuseeps = epsilon;
438  	
439  	         for(j = pupd.size() - 1; j >= 0; --j)
440  	         {
441  	            i = pupd.index(j);
442  	            x = pupd[i];
443  	
444  	            if(x > ruseeps)
445  	            {
446  	               y = upb[i] - pvec[i];
447  	
448  	               if(y < -degeneps)
449  	                  this->solver()->shiftUPbound(i, pvec[i] - degeneps);
450  	               else
451  	               {
452  	                  y /= x;
453  	
454  	                  if(y <= max && x >= stab)        // &&  y > sel-epsilon
455  	                  {
456  	                     enterId = this->solver()->id(i);
457  	                     sel = y;
458  	                     pnr = i;
459  	                     stab = x;
460  	                  }
461  	               }
462  	            }
463  	            else if(x < -ruseeps)
464  	            {
465  	               y = lpb[i] - pvec[i];
466  	
467  	               if(y > degeneps)
468  	                  this->solver()->shiftLPbound(i, pvec[i] + degeneps);
469  	               else
470  	               {
471  	                  y /= x;
472  	
473  	                  if(y <= max && -x >= stab)       // &&  y > sel-epsilon
474  	                  {
475  	                     enterId = this->solver()->id(i);
476  	                     sel = y;
477  	                     pnr = i;
478  	                     stab = -x;
479  	                  }
480  	               }
481  	            }
482  	            else
483  	            {
484  	               SPxOut::debug(this, "DHARRI01 removing value {}\n", pupd[i]);
485  	               pupd.clearNum(j);
486  	            }
487  	         }
488  	
489  	         for(j = cupd.size() - 1; j >= 0; --j)
490  	         {
491  	            i = cupd.index(j);
492  	            x = cupd[i];
493  	
494  	            if(x > cuseeps)
495  	            {
496  	               y = ucb[i] - cvec[i];
497  	
498  	               if(y < -degeneps)
499  	                  this->solver()->shiftUCbound(i, cvec[i] - degeneps);
500  	               else
501  	               {
502  	                  y /= x;
503  	
504  	                  if(y <= max && x >= stab)        // &&  y > sel-epsilon
505  	                  {
506  	                     enterId = this->solver()->coId(i);
507  	                     sel = y;
508  	                     cnr = j;
509  	                     stab = x;
510  	                  }
511  	               }
512  	            }
513  	            else if(x < -cuseeps)
514  	            {
515  	               y = lcb[i] - cvec[i];
516  	
517  	               if(y > degeneps)
518  	                  this->solver()->shiftLCbound(i, cvec[i] + degeneps);
519  	               else
520  	               {
521  	                  y /= x;
522  	
523  	                  if(y <= max && -x >= stab)       // &&  y > sel-epsilon
524  	                  {
525  	                     enterId = this->solver()->coId(i);
526  	                     sel = y;
527  	                     cnr = j;
528  	                     stab = -x;
529  	                  }
530  	               }
531  	            }
532  	            else
533  	            {
534  	               SPxOut::debug(this, "DHARRI02 removing value {}\n", cupd[i]);
535  	               cupd.clearNum(j);
536  	            }
537  	         }
538  	
539  	         if(lastshift == this->solver()->shift())
540  	         {
541  	            if(cnr >= 0)
542  	            {
543  	               if(this->solver()->isBasic(enterId))
544  	               {
545  	                  cupd.clearNum(cnr);
546  	                  continue;
547  	               }
548  	               else
549  	                  break;
550  	            }
551  	            else if(pnr >= 0)
552  	            {
553  	               pvec[pnr] = this->solver()->vector(pnr) * cvec;
554  	
555  	               if(this->solver()->isBasic(enterId))
556  	               {
557  	                  pupd.setValue(pnr, 0.0);
558  	                  continue;
559  	               }
560  	               else
561  	               {
562  	                  x = pupd[pnr];
563  	
564  	                  if(x > 0)
565  	                  {
566  	                     sel = upb[pnr] - pvec[pnr];
567  	
568  	                     if(x < minStability && sel < this->delta)
569  	                     {
570  	                        minStability /= 2.0;
571  	                        this->solver()->shiftUPbound(pnr, pvec[pnr]);
572  	                        continue;
573  	                     }
574  	                  }
575  	                  else
576  	                  {
577  	                     sel = lpb[pnr] - pvec[pnr];
578  	
579  	                     if(-x < minStability && -sel < this->delta)
580  	                     {
581  	                        minStability /= 2.0;
582  	                        this->solver()->shiftLPbound(pnr, pvec[pnr]);
583  	                        continue;
584  	                     }
585  	                  }
586  	
587  	                  sel /= x;
588  	               }
589  	            }
590  	            else
591  	            {
592  	               val = 0;
593  	               enterId.inValidate();
594  	               return enterId;
595  	            }
596  	
597  	            if(sel > max)              // instability detected => recompute
598  	               continue;               // ratio test with corrected value
599  	
600  	            break;
601  	         }
602  	      }
603  	   }
604  	   else if(val < -epsilon)
605  	   {
606  	      for(;;)
607  	      {
608  	         pnr = -1;
609  	         cnr = -1;
610  	         max = val;
611  	         lastshift = this->solver()->shift();
612  	         assert(this->delta > epsilon);
613  	
614  	
615  	         // phase 1:
616  	         minDelta
617  	         (
618  	            &rmaxabs,            /* max abs value in upd */
619  	            &max,                /* initial and chosen value */
620  	            pupd.size(),         /* # of indices in pupd */
621  	            pupd.indexMem(),     /* nonzero indices in pupd */
622  	            pupd.values(),       /* update VectorBase<R> for vec */
623  	            pvec.get_const_ptr(),                /* current VectorBase<R> */
624  	            lpb.get_const_ptr(),                 /* lower bounds for vec */
625  	            upb.get_const_ptr()                  /* upper bounds for vec */
626  	         );
627  	
628  	         minDelta
629  	         (
630  	            &cmaxabs,            /* max abs value in upd */
631  	            &max,                /* initial and chosen value */
632  	            cupd.size(),         /* # of indices in cupd */
633  	            cupd.indexMem(),     /* nonzero indices in cupd */
634  	            cupd.values(),       /* update VectorBase<R> for vec */
635  	            cvec.get_const_ptr(),                /* current VectorBase<R> */
636  	            lcb.get_const_ptr(),                 /* lower bounds for vec */
637  	            ucb.get_const_ptr()                  /* upper bounds for vec */
638  	         );
639  	
640  	         if(max == val)
641  	            return enterId;
642  	
643  	
644  	         // phase 2:
645  	         stab = 0;
646  	         sel = R(infinity);
647  	         ruseeps = rmaxabs * epsilon * 0.001;
648  	         cuseeps = cmaxabs * epsilon * 0.001;
649  	
650  	         for(j = pupd.size() - 1; j >= 0; --j)
651  	         {
652  	            i = pupd.index(j);
653  	            x = pupd[i];
654  	
655  	            if(x > ruseeps)
656  	            {
657  	               y = lpb[i] - pvec[i];
658  	
659  	               if(y > degeneps)
660  	                  this->solver()->shiftLPbound(i, pvec[i]);  // ensure simplex improvement
661  	               else
662  	               {
663  	                  y /= x;
664  	
665  	                  if(y >= max && x > stab)         // &&  y < sel+epsilon
666  	                  {
667  	                     enterId = this->solver()->id(i);
668  	                     sel = y;
669  	                     pnr = i;
670  	                     stab = x;
671  	                  }
672  	               }
673  	            }
674  	            else if(x < -ruseeps)
675  	            {
676  	               y = upb[i] - pvec[i];
677  	
678  	               if(y < -degeneps)
679  	                  this->solver()->shiftUPbound(i, pvec[i]);  // ensure simplex improvement
680  	               else
681  	               {
682  	                  y /= x;
683  	
684  	                  if(y >= max && -x > stab)        // &&  y < sel+epsilon
685  	                  {
686  	                     enterId = this->solver()->id(i);
687  	                     sel = y;
688  	                     pnr = i;
689  	                     stab = -x;
690  	                  }
691  	               }
692  	            }
693  	            else
694  	            {
695  	               SPxOut::debug(this, "DHARRI03 removing value {}\n", pupd[i]);
696  	               pupd.clearNum(j);
697  	            }
698  	         }
699  	
700  	         for(j = cupd.size() - 1; j >= 0; --j)
701  	         {
702  	            i = cupd.index(j);
703  	            x = cupd[i];
704  	
705  	            if(x > cuseeps)
706  	            {
707  	               y = lcb[i] - cvec[i];
708  	
709  	               if(y > degeneps)
710  	                  this->solver()->shiftLCbound(i, cvec[i]);  // ensure simplex improvement
711  	               else
712  	               {
713  	                  y /= x;
714  	
715  	                  if(y >= max && x > stab)         // &&  y < sel+epsilon
716  	                  {
717  	                     enterId = this->solver()->coId(i);
718  	                     sel = y;
719  	                     cnr = j;
720  	                     stab = x;
721  	                  }
722  	               }
723  	            }
724  	            else if(x < -cuseeps)
725  	            {
726  	               y = ucb[i] - cvec[i];
727  	
728  	               if(y < -degeneps)
729  	                  this->solver()->shiftUCbound(i, cvec[i]);  // ensure simplex improvement
730  	               else
731  	               {
732  	                  y /= x;
733  	
734  	                  if(y >= max && -x > stab)        // &&  y < sel+epsilon
735  	                  {
736  	                     enterId = this->solver()->coId(i);
737  	                     sel = y;
738  	                     cnr = j;
739  	                     stab = -x;
740  	                  }
741  	               }
742  	            }
743  	            else
744  	            {
745  	               SPxOut::debug(this, "DHARRI04 removing value {}\n", x);
746  	               cupd.clearNum(j);
747  	            }
748  	         }
749  	
750  	         if(lastshift == this->solver()->shift())
751  	         {
752  	            if(cnr >= 0)
753  	            {
754  	               if(this->solver()->isBasic(enterId))
755  	               {
756  	                  cupd.clearNum(cnr);
757  	                  continue;
758  	               }
759  	               else
760  	                  break;
761  	            }
762  	            else if(pnr >= 0)
763  	            {
764  	               pvec[pnr] = this->solver()->vector(pnr) * cvec;
765  	
766  	               if(this->solver()->isBasic(enterId))
767  	               {
768  	                  pupd.setValue(pnr, 0.0);
769  	                  continue;
770  	               }
771  	               else
772  	               {
773  	                  x = pupd[pnr];
774  	
775  	                  if(x > 0)
776  	                  {
777  	                     sel = lpb[pnr] - pvec[pnr];
778  	
779  	                     if(x < minStability && -sel < this->delta)
780  	                     {
781  	                        minStability /= 2;
782  	                        this->solver()->shiftLPbound(pnr, pvec[pnr]);
783  	                        continue;
784  	                     }
785  	                  }
786  	                  else
787  	                  {
788  	                     sel = upb[pnr] - pvec[pnr];
789  	
790  	                     if(-x < minStability && sel < this->delta)
791  	                     {
792  	                        minStability /= 2;
793  	                        this->solver()->shiftUPbound(pnr, pvec[pnr]);
794  	                        continue;
795  	                     }
796  	                  }
797  	
798  	                  sel /= x;
799  	               }
800  	            }
801  	            else
802  	            {
803  	               val = 0;
804  	               enterId.inValidate();
805  	               return enterId;
806  	            }
807  	
808  	            if(sel < max)              // instability detected => recompute
809  	               continue;               // ratio test with corrected value
810  	
811  	            break;
812  	         }
813  	      }
814  	   }
815  	
816  	   assert(max * val >= 0);
817  	   assert(enterId.type() != SPxId::INVALID);
818  	
819  	   val = sel;
820  	
821  	   return enterId;
822  	}
823  	} // namespace soplex
824