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 <iostream>
27   	
28   	#include "soplex/spxdefines.h"
29   	#include "soplex/spxsolver.h"
30   	#include "soplex/spxout.h"
31   	
32   	namespace soplex
33   	{
34   	template <class R>
35   	void SPxSolverBase<R>::shiftFvec()
36   	{
37   	
38   	   /* the allowed tolerance is (rep() == COLUMN) ? tolerances()->floatingPointFeastol() : tolerances()->floatingPointOpttol() because theFvec is the primal VectorBase<R> in COLUMN
39   	    * and the dual VectorBase<R> in ROW representation; this is equivalent to entertol()
40   	    */
41   	   R minrandom = 10.0 * entertol();
42   	   R maxrandom = 100.0 * entertol();
43   	   R allow = entertol() - epsilon();
44   	
45   	   assert(type() == ENTER);
46   	   assert(allow > 0);
47   	
48   	   for(int i = dim() - 1; i >= 0; --i)
49   	   {
50   	      if(theUBbound[i] + allow < (*theFvec)[i])
51   	      {
52   	         SPxOut::debug(this, "DSHIFT08 theUBbound[{}] violated by {}", i,
53   	                       (*theFvec)[i] - theUBbound[i] - allow);
54   	
55   	         if(theUBbound[i] != theLBbound[i])
56   	         {
57   	            // since minrandom and maxrandom are of the order 10 different,
58   	            // we currently doesn't care about higher precision random
59   	            // numbers. Hence the cast to double.
60   	            shiftUBbound(i, (*theFvec)[i] + random.next((double)minrandom, (double)maxrandom));
61   	         }
62   	         else
63   	         {
64   	            shiftUBbound(i, (*theFvec)[i]);
65   	            theLBbound[i] = theUBbound[i];
66   	         }
67   	      }
68   	      else if((*theFvec)[i] < theLBbound[i] - allow)
69   	      {
70   	         SPxOut::debug(this, "DSHIFT08 theLBbound[{}] violated by {}", i,
71   	                       theLBbound[i] - (*theFvec)[i] - allow);
72   	
73   	         if(theUBbound[i] != theLBbound[i])
74   	            shiftLBbound(i, (*theFvec)[i] - random.next((double)minrandom, (double)maxrandom));
75   	         else
76   	         {
77   	            shiftLBbound(i, (*theFvec)[i]);
78   	            theUBbound[i] = theLBbound[i];
79   	         }
80   	      }
81   	   }
82   	
83   	#ifndef NDEBUG
84   	   testBounds();
85   	   SPxOut::debug(this, "DSHIFT01 shiftFvec: OK\n");
86   	#endif
87   	}
88   	
89   	// -----------------------------------------------------------------
90   	
91   	/*
92   	  This methods assumes correctly setup vectors |pVec| and |coPvec| and bound
93   	  vectors for leaving simplex. Then it checks all values of |pVec| and
94   	  |coPvec| to obey these bounds and enlarges them if neccessary.
95   	*/
96   	template <class R>
97   	void SPxSolverBase<R>::shiftPvec()
98   	{
99   	
100  	   /* the allowed tolerance is (rep() == ROW) ? tolerances()->floatingPointFeastol() : tolerances()->floatingPointOpttol() because thePvec is the primal VectorBase<R> in ROW and the
101  	    * dual VectorBase<R> in COLUMN representation; this is equivalent to leavetol()
102  	    */
103  	   R minrandom = 10.0 * leavetol();
104  	   R maxrandom = 100.0 * leavetol();
105  	   R allow = leavetol() - epsilon();
106  	   bool tmp;
107  	   int i;
108  	
109  	   assert(type() == LEAVE);
110  	   assert(allow > 0.0);
111  	
112  	   for(i = dim() - 1; i >= 0; --i)
113  	   {
114  	      tmp = !isBasic(coId(i));
115  	
116  	      if((*theCoUbound)[i] + allow <= (*theCoPvec)[i] && tmp)
117  	      {
118  	         if((*theCoUbound)[i] != (*theCoLbound)[i])
119  	            shiftUCbound(i, (*theCoPvec)[i] + random.next((double)minrandom, (double)maxrandom));
120  	         else
121  	         {
122  	            shiftUCbound(i, (*theCoPvec)[i]);
123  	            (*theCoLbound)[i] = (*theCoUbound)[i];
124  	         }
125  	      }
126  	      else if((*theCoLbound)[i] - allow >= (*theCoPvec)[i] && tmp)
127  	      {
128  	         if((*theCoUbound)[i] != (*theCoLbound)[i])
129  	            shiftLCbound(i, (*theCoPvec)[i] - random.next((double)minrandom, (double)maxrandom));
130  	         else
131  	         {
132  	            shiftLCbound(i, (*theCoPvec)[i]);
133  	            (*theCoUbound)[i] = (*theCoLbound)[i];
134  	         }
135  	      }
136  	   }
137  	
138  	   for(i = coDim() - 1; i >= 0; --i)
139  	   {
140  	      tmp = !isBasic(id(i));
141  	
142  	      if((*theUbound)[i] + allow <= (*thePvec)[i] && tmp)
143  	      {
144  	         if((*theUbound)[i] != (*theLbound)[i])
145  	            shiftUPbound(i, (*thePvec)[i] + random.next((double)minrandom, (double)maxrandom));
146  	         else
147  	         {
148  	            shiftUPbound(i, (*thePvec)[i]);
149  	            (*theLbound)[i] = (*theUbound)[i];
150  	         }
151  	      }
152  	      else if((*theLbound)[i] - allow >= (*thePvec)[i] && tmp)
153  	      {
154  	         if((*theUbound)[i] != (*theLbound)[i])
155  	            shiftLPbound(i, (*thePvec)[i] - random.next((double)minrandom, (double)maxrandom));
156  	         else
157  	         {
158  	            shiftLPbound(i, (*thePvec)[i]);
159  	            (*theUbound)[i] = (*theLbound)[i];
160  	         }
161  	      }
162  	   }
163  	
164  	#ifndef NDEBUG
165  	   testBounds();
166  	   SPxOut::debug(this, "DSHIFT02 shiftPvec: OK\n");
167  	#endif
168  	}
169  	// -----------------------------------------------------------------
170  	template <class R>
171  	void SPxSolverBase<R>::perturbMin(
172  	   const UpdateVector<R>& uvec,
173  	   VectorBase<R>& p_low,
174  	   VectorBase<R>& p_up,
175  	   R eps,
176  	   R p_delta,
177  	   int start,
178  	   int incr)
179  	{
180  	   assert(uvec.dim() == p_low.dim());
181  	   assert(uvec.dim() == p_up.dim());
182  	
183  	   const R* vec = uvec.get_const_ptr();
184  	   R minrandom = 10.0 * p_delta;
185  	   R maxrandom = 100.0 * p_delta;
186  	   R x, l, u;
187  	   int i;
188  	
189  	   if(fullPerturbation)
190  	   {
191  	      eps = p_delta;
192  	
193  	      for(i = uvec.dim() - start - 1; i >= 0; i -= incr)
194  	      {
195  	         u = p_up[i];
196  	         l = p_low[i];
197  	         x = vec[i];
198  	
199  	         if(LT(u, R(infinity), eps) && NE(l, u, eps) && u <= x + eps)
200  	         {
201  	            p_up[i] = x + random.next((double) minrandom, (double)maxrandom);
202  	            theShift += p_up[i] - u;
203  	         }
204  	
205  	         if(GT(l, R(-infinity), eps) && NE(l, u, eps) && l >= x - eps)
206  	         {
207  	            p_low[i] = x - random.next((double)minrandom, (double)maxrandom);
208  	            theShift -= p_low[i] - l;
209  	         }
210  	      }
211  	   }
212  	   else
213  	   {
214  	      const R* upd = uvec.delta().values();
215  	      const IdxSet& idx = uvec.delta().indices();
216  	
217  	      for(int j = uvec.delta().size() - start - 1; j >= 0; j -= incr)
218  	      {
219  	         i = idx.index(j);
220  	         x = upd[i];
221  	         u = p_up[i];
222  	         l = p_low[i];
223  	
224  	         // do not permute these bounds! c.f. with computeFrhs2() in spxvecs.cpp
225  	         if(this->dualStatus(this->baseId(i)) == SPxBasisBase<R>::Desc::D_ON_BOTH)
226  	         {
227  	            continue;
228  	         }
229  	
230  	         if(x < -eps)
231  	         {
232  	            if(LT(u, R(infinity), eps) && NE(l, u, eps) && vec[i] >= u - eps)
233  	            {
234  	               p_up[i] = vec[i] + random.next((double)minrandom, (double)maxrandom);
235  	               theShift += p_up[i] - u;
236  	            }
237  	         }
238  	         else if(x > eps)
239  	         {
240  	            if(GT(l, R(-infinity), eps) && NE(l, u, eps) && vec[i] <= l + eps)
241  	            {
242  	               p_low[i] = vec[i] - random.next((double)minrandom, (double)maxrandom);
243  	               theShift -= p_low[i] - l;
244  	            }
245  	         }
246  	      }
247  	   }
248  	}
249  	// -----------------------------------------------------------------
250  	template <class R>
251  	void SPxSolverBase<R>::perturbMax(
252  	   const UpdateVector<R>& uvec,
253  	   VectorBase<R>& p_low,
254  	   VectorBase<R>& p_up,
255  	   R eps,
256  	   R p_delta,
257  	   int start,
258  	   int incr)
259  	{
260  	   assert(uvec.dim() == p_low.dim());
261  	   assert(uvec.dim() == p_up.dim());
262  	
263  	   const R* vec = uvec.get_const_ptr();
264  	   R minrandom = 10.0 * p_delta;
265  	   R maxrandom = 100.0 * p_delta;
266  	   R x, l, u;
267  	   int i;
268  	
269  	   if(fullPerturbation)
270  	   {
271  	      eps = p_delta;
272  	
273  	      for(i = uvec.dim() - start - 1; i >= 0; i -= incr)
274  	      {
275  	         u = p_up[i];
276  	         l = p_low[i];
277  	         x = vec[i];
278  	
279  	         if(LT(u, R(infinity), eps) && NE(l, u, eps) && u <= x + eps)
280  	         {
281  	            p_up[i] = x + random.next((double)minrandom, (double)maxrandom);
282  	            theShift += p_up[i] - u;
283  	         }
284  	
285  	         if(GT(l, R(-infinity), eps) && NE(l, u, eps) && l >= x - eps)
286  	         {
287  	            p_low[i] = x - random.next((double)minrandom, (double)maxrandom);
288  	            theShift -= p_low[i] - l;
289  	         }
290  	      }
291  	   }
292  	   else
293  	   {
294  	      const R* upd = uvec.delta().values();
295  	      const IdxSet& idx = uvec.delta().indices();
296  	
297  	      for(int j = uvec.delta().size() - start - 1; j >= 0; j -= incr)
298  	      {
299  	         i = idx.index(j);
300  	         x = upd[i];
301  	         u = p_up[i];
302  	         l = p_low[i];
303  	
304  	         // do not permute these bounds! c.f. computeFrhs2() in spxvecs.cpp
305  	         if(this->dualStatus(this->baseId(i)) == SPxBasisBase<R>::Desc::D_ON_BOTH)
306  	         {
307  	            continue;
308  	         }
309  	
310  	         if(x > eps)
311  	         {
312  	            if(LT(u, R(infinity), eps) && NE(l, u, eps) && vec[i] >= u - eps)
313  	            {
314  	               p_up[i] = vec[i] + random.next((double)minrandom, (double)maxrandom);
315  	               theShift += p_up[i] - u;
316  	            }
317  	         }
318  	         else if(x < -eps)
319  	         {
320  	            if(GT(l, R(-infinity), eps) && NE(l, u, eps) && vec[i] <= l + eps)
321  	            {
322  	               p_low[i] = vec[i] - random.next((double)minrandom, (double)maxrandom);
323  	               theShift -= p_low[i] - l;
324  	            }
325  	         }
326  	      }
327  	   }
328  	}
329  	
330  	template <class R>
331  	void SPxSolverBase<R>::perturbMinEnter(void)
332  	{
333  	   SPxOut::debug(this, "DSHIFT03 iteration= {} perturbing {}", this->iteration(), shift());
334  	   fVec().delta().setup();
335  	   perturbMin(fVec(), lbBound(), ubBound(), epsilon(), entertol());
336  	   SPxOut::debug(this, "\t->{}\n", shift());
337  	}
338  	
339  	
340  	template <class R>
341  	void SPxSolverBase<R>::perturbMaxEnter(void)
342  	{
343  	   SPxOut::debug(this, "DSHIFT04 iteration= {} perturbing {}", this->iteration(), shift());
344  	   fVec().delta().setup();
345  	   perturbMax(fVec(), lbBound(), ubBound(), epsilon(), entertol());
346  	   SPxOut::debug(this, "\t->{}\n", shift());
347  	}
348  	
349  	
350  	template <class R>
351  	R SPxSolverBase<R>::perturbMin(
352  	   const UpdateVector<R>& uvec,
353  	   VectorBase<R>& p_low,
354  	   VectorBase<R>& p_up,
355  	   R eps,
356  	   R p_delta,
357  	   const typename SPxBasisBase<R>::Desc::Status* stat,
358  	   int start,
359  	   int incr)
360  	{
361  	   assert(uvec.dim() == p_low.dim());
362  	   assert(uvec.dim() == p_up.dim());
363  	
364  	   const R* vec = uvec.get_const_ptr();
365  	   R minrandom = 10.0 * p_delta;
366  	   R maxrandom = 100.0 * p_delta;
367  	   R x, l, u;
368  	   int i;
369  	   R l_theShift = 0;
370  	
371  	   if(fullPerturbation)
372  	   {
373  	      eps = p_delta;
374  	
375  	      for(i = uvec.dim() - start - 1; i >= 0; i -= incr)
376  	      {
377  	         u = p_up[i];
378  	         l = p_low[i];
379  	         x = vec[i];
380  	
381  	         if(LT(u, R(infinity), eps) && NE(l, u, eps) && u <= x + eps && rep() * stat[i] < 0)
382  	         {
383  	            p_up[i] = vec[i] + random.next((double)minrandom, (double)maxrandom);
384  	            l_theShift += p_up[i] - u;
385  	         }
386  	
387  	         if(GT(l, R(-infinity), eps) && NE(l, u, eps) && l >= x - eps && rep() * stat[i] < 0)
388  	         {
389  	            p_low[i] = vec[i] - random.next((double)minrandom, (double)maxrandom);
390  	            l_theShift -= p_low[i] - l;
391  	         }
392  	      }
393  	   }
394  	   else
395  	   {
396  	      const R* upd = uvec.delta().values();
397  	      const IdxSet& idx = uvec.delta().indices();
398  	
399  	      for(int j = uvec.delta().size() - start - 1; j >= 0; j -= incr)
400  	      {
401  	         i = idx.index(j);
402  	         x = upd[i];
403  	         u = p_up[i];
404  	         l = p_low[i];
405  	
406  	         if(x < -eps)
407  	         {
408  	            if(LT(u, R(infinity), eps) && NE(l, u, eps) && vec[i] >= u - eps && rep() * stat[i] < 0)
409  	            {
410  	               p_up[i] = vec[i] + random.next((double)minrandom, (double)maxrandom);
411  	               l_theShift += p_up[i] - u;
412  	            }
413  	         }
414  	         else if(x > eps)
415  	         {
416  	            if(GT(l, R(-infinity), eps) && NE(l, u, eps) && vec[i] <= l + eps && rep() * stat[i] < 0)
417  	            {
418  	               p_low[i] = vec[i] - random.next((double)minrandom, (double)maxrandom);
419  	               l_theShift -= p_low[i] - l;
420  	            }
421  	         }
422  	      }
423  	   }
424  	
425  	   return l_theShift;
426  	}
427  	
428  	template <class R>
429  	R SPxSolverBase<R>::perturbMax(
430  	   const UpdateVector<R>& uvec,
431  	   VectorBase<R>& p_low,
432  	   VectorBase<R>& p_up,
433  	   R eps,
434  	   R p_delta,
435  	   const typename SPxBasisBase<R>::Desc::Status* stat,
436  	   int start,
437  	   int incr)
438  	{
439  	   assert(uvec.dim() == p_low.dim());
440  	   assert(uvec.dim() == p_up.dim());
441  	
442  	   const R* vec = uvec.get_const_ptr();
443  	   R minrandom = 10.0 * p_delta;
444  	   R maxrandom = 100.0 * p_delta;
445  	   R x, l, u;
446  	   int i;
447  	   R l_theShift = 0;
448  	
449  	   if(fullPerturbation)
450  	   {
451  	      eps = p_delta;
452  	
453  	      for(i = uvec.dim() - start - 1; i >= 0; i -= incr)
454  	      {
455  	         u = p_up[i];
456  	         l = p_low[i];
457  	         x = vec[i];
458  	
459  	         if(LT(u, R(infinity), eps) && NE(l, u, eps) && u <= x + eps && rep() * stat[i] < 0)
460  	         {
461  	            p_up[i] = vec[i] + random.next((double)minrandom, (double)maxrandom);
462  	            l_theShift += p_up[i] - u;
463  	         }
464  	
465  	         if(GT(l, R(-infinity), eps) && NE(l, u, eps) && l >= x - eps && rep() * stat[i] < 0)
466  	         {
467  	            p_low[i] = vec[i] - random.next((double)minrandom, (double)maxrandom);
468  	            l_theShift -= p_low[i] - l;
469  	         }
470  	      }
471  	   }
472  	   else
473  	   {
474  	      const R* upd = uvec.delta().values();
475  	      const IdxSet& idx = uvec.delta().indices();
476  	
477  	      for(int j = uvec.delta().size() - start - 1; j >= 0; j -= incr)
478  	      {
479  	         i = idx.index(j);
480  	         x = upd[i];
481  	         u = p_up[i];
482  	         l = p_low[i];
483  	
484  	         if(x > eps)
485  	         {
486  	            if(LT(u, R(infinity), eps) && NE(l, u, eps) && vec[i] >= u - eps && rep() * stat[i] < 0)
487  	            {
488  	               p_up[i] = vec[i] + random.next((double)minrandom, (double)maxrandom);
489  	               l_theShift += p_up[i] - u;
490  	            }
491  	         }
492  	         else if(x < -eps)
493  	         {
494  	            if(GT(l, R(-infinity), eps) && NE(l, u, eps) && vec[i] <= l + eps && rep() * stat[i] < 0)
495  	            {
496  	               p_low[i] = vec[i] - random.next((double)minrandom, (double)maxrandom);
497  	               l_theShift -= p_low[i] - l;
498  	            }
499  	         }
500  	      }
501  	   }
502  	
503  	   return l_theShift;
504  	}
505  	
506  	
507  	template <class R>
508  	void SPxSolverBase<R>::perturbMinLeave(void)
509  	{
510  	   SPxOut::debug(this, "DSHIFT05 iteration= {} perturbing {}", this->iteration(), shift());
511  	   pVec().delta().setup();
512  	   coPvec().delta().setup();
513  	   theShift += perturbMin(pVec(), lpBound(), upBound(), epsilon(), leavetol(),
514  	                          this->desc().status(), 0, 1);
515  	   theShift += perturbMin(coPvec(), lcBound(), ucBound(), epsilon(), leavetol(),
516  	                          this->desc().coStatus(), 0, 1);
517  	   SPxOut::debug(this, "\t->{}\n", shift());
518  	}
519  	
520  	
521  	template <class R>
522  	void SPxSolverBase<R>::perturbMaxLeave(void)
523  	{
524  	   SPxOut::debug(this, "DSHIFT06 iteration= {} perturbing {}", this->iteration(), shift());
525  	   pVec().delta().setup();
526  	   coPvec().delta().setup();
527  	   theShift += perturbMax(pVec(), lpBound(), upBound(), epsilon(), leavetol(),
528  	                          this->desc().status(), 0, 1);
529  	   theShift += perturbMax(coPvec(), lcBound(), ucBound(), epsilon(), leavetol(),
530  	                          this->desc().coStatus(), 0, 1);
531  	   SPxOut::debug(this, "\t->{}\n", shift());
532  	}
533  	
534  	
535  	template <class R>
536  	void SPxSolverBase<R>::unShift(void)
537  	{
538  	   SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "DSHIFT07 = " << "unshifting ..." << std::endl;);
539  	
540  	   if(isInitialized())
541  	   {
542  	      int i;
543  	      R t_up, t_low;
544  	      const typename SPxBasisBase<R>::Desc& ds = this->desc();
545  	
546  	      theShift = 0;
547  	
548  	      if(type() == ENTER)
549  	      {
550  	         R eps = entertol();
551  	
552  	         if(rep() == COLUMN)
553  	         {
554  	            for(i = dim(); i-- > 0;)
555  	            {
556  	               SPxId l_id = this->baseId(i);
557  	               int l_num = this->number(l_id);
558  	
559  	               if(l_id.type() == SPxId::ROW_ID)
560  	               {
561  	                  t_up = -this->lhs(l_num);
562  	                  t_low = -this->rhs(l_num);
563  	               }
564  	               else
565  	               {
566  	                  assert(l_id.type() == SPxId::COL_ID);
567  	                  t_up = this->upper(l_num);
568  	                  t_low = this->lower(l_num);
569  	               }
570  	
571  	               if(t_up != t_low)
572  	               {
573  	                  if((*theFvec)[i] < t_up + eps)  // check allowed violation
574  	                     theUBbound[i] = t_up; // reset shifted bound to original
575  	                  else if((*theFvec)[i] > t_up)  // shifted bound is required for feasibility
576  	                     theShift += theUBbound[i] - t_up;
577  	
578  	                  if((*theFvec)[i] > t_low - eps)  // check allowed violation
579  	                     theLBbound[i] = t_low; // reset shifted bound to original
580  	                  else if((*theFvec)[i] < t_low)  // shifted bound is required for feasibility
581  	                     theShift -= theLBbound[i] - t_low;
582  	               }
583  	               else
584  	               {
585  	                  if(theUBbound[i] > t_up)
586  	                     theShift += theUBbound[i] - t_up;
587  	                  else if(theLBbound[i] < t_low)
588  	                     theShift += t_low - theLBbound[i];
589  	               }
590  	            }
591  	
592  	            for(i = this->nRows(); i-- > 0;)
593  	            {
594  	               if(!isBasic(ds.rowStatus(i)))
595  	               {
596  	                  t_up = -this->lhs(i);
597  	                  t_low = -this->rhs(i);
598  	
599  	                  if(theURbound[i] > t_up)  // what about t_up == t_low ?
600  	                     theShift += theURbound[i] - t_up;
601  	
602  	                  if(t_low > theLRbound[i])  // what about t_up == t_low ?
603  	                     theShift += t_low - theLRbound[i];
604  	               }
605  	            }
606  	
607  	            for(i = this->nCols(); i-- > 0;)
608  	            {
609  	               if(!isBasic(ds.colStatus(i)))
610  	               {
611  	                  t_up = this->upper(i);
612  	                  t_low = this->lower(i);
613  	
614  	                  if(theUCbound[i] > t_up)  // what about t_up == t_low ?
615  	                     theShift += theUCbound[i] - t_up;
616  	
617  	                  if(t_low > theLCbound[i])  // what about t_up == t_low ?
618  	                     theShift += t_low - theLCbound[i];
619  	               }
620  	            }
621  	         }
622  	         else
623  	         {
624  	            assert(rep() == ROW);
625  	
626  	            for(i = dim(); i-- > 0;)
627  	            {
628  	               SPxId l_id = this->baseId(i);
629  	               int l_num = this->number(l_id);
630  	               t_up = t_low = 0;
631  	
632  	               if(l_id.type() == SPxId::ROW_ID)
633  	                  clearDualBounds(ds.rowStatus(l_num), t_up, t_low);
634  	               else
635  	                  clearDualBounds(ds.colStatus(l_num), t_up, t_low);
636  	
637  	               if(theUBbound[i] != theLBbound[i])
638  	               {
639  	                  if(theUBbound[i] > t_up)
640  	                     theShift += theUBbound[i] - t_up;
641  	                  else
642  	                     theShift -= theUBbound[i] - t_up;
643  	               }
644  	               else
645  	               {
646  	                  /* if the basic (primal or dual) variable is fixed (e.g., basis status P_FREE in row representation)
647  	                   * then shiftFvec() and shiftPvec() do not relax the bounds, but shift both, hence they may be outside
648  	                   * of [t_low,t_up] */
649  	                  assert(theLBbound[i] == theUBbound[i] || theUBbound[i] >= t_up);
650  	                  assert(theLBbound[i] == theUBbound[i] || theLBbound[i] <= t_low);
651  	
652  	                  if((*theFvec)[i] < t_up - eps)
653  	                     theUBbound[i] = t_up;
654  	                  else if((*theFvec)[i] > t_up)
655  	                     theShift += theUBbound[i] - t_up;
656  	
657  	                  if((*theFvec)[i] > t_low + eps)
658  	                     theLBbound[i] = t_low;
659  	                  else if((*theFvec)[i] < t_low)
660  	                     theShift -= theLBbound[i] - t_low;
661  	               }
662  	            }
663  	
664  	            for(i = this->nRows(); i-- > 0;)
665  	            {
666  	               if(!isBasic(ds.rowStatus(i)))
667  	               {
668  	                  t_up = t_low = 0;
669  	                  clearDualBounds(ds.rowStatus(i), t_up, t_low);
670  	
671  	                  if(theURbound[i] > t_up)  // what about t_up == t_low ?
672  	                     theShift += theURbound[i] - t_up;
673  	
674  	                  if(t_low > theLRbound[i])  // what about t_up == t_low ?
675  	                     theShift += t_low - theLRbound[i];
676  	               }
677  	            }
678  	
679  	            for(i = this->nCols(); i-- > 0;)
680  	            {
681  	               if(!isBasic(ds.colStatus(i)))
682  	               {
683  	                  t_up = t_low = 0;
684  	                  clearDualBounds(ds.colStatus(i), t_up, t_low);
685  	
686  	                  if(theUCbound[i] > t_up)  // what about t_up == t_low ?
687  	                     theShift += theUCbound[i] - t_up;
688  	
689  	                  if(t_low > theLCbound[i])  // what about t_up == t_low ?
690  	                     theShift += t_low - theLCbound[i];
691  	               }
692  	            }
693  	         }
694  	      }
695  	      else
696  	      {
697  	         assert(type() == LEAVE);
698  	
699  	         R eps = leavetol();
700  	
701  	         if(rep() == COLUMN)
702  	         {
703  	            for(i = this->nRows(); i-- > 0;)
704  	            {
705  	               t_up = t_low = this->maxRowObj(i);
706  	               clearDualBounds(ds.rowStatus(i), t_up, t_low);
707  	
708  	               if(!isBasic(ds.rowStatus(i)))
709  	               {
710  	                  if((*theCoPvec)[i] < t_up + eps)
711  	                  {
712  	                     theURbound[i] = t_up; // reset bound to original value
713  	
714  	                     if(t_up == t_low)
715  	                        theLRbound[i] = t_low; // for fixed rows we change both bounds
716  	                  }
717  	                  else
718  	                     theShift += theURbound[i] - t_up;
719  	
720  	                  if((*theCoPvec)[i] > t_low - eps)
721  	                  {
722  	                     theLRbound[i] = t_low; // reset bound to original value
723  	
724  	                     if(t_up == t_low)
725  	                        theURbound[i] = t_up; // for fixed rows we change both bounds
726  	                  }
727  	                  else
728  	                     theShift += t_low - theLRbound[i];
729  	               }
730  	               else if(theURbound[i] > t_up)
731  	                  theShift += theURbound[i] - t_up;
732  	               else if(theLRbound[i] < t_low)
733  	                  theShift += t_low - theLRbound[i];
734  	            }
735  	
736  	            for(i = this->nCols(); i-- > 0;)
737  	            {
738  	               t_up = t_low = -this->maxObj(i);
739  	               clearDualBounds(ds.colStatus(i), t_low, t_up);
740  	
741  	               if(!isBasic(ds.colStatus(i)))
742  	               {
743  	                  if((*thePvec)[i] < -t_up + eps)
744  	                  {
745  	                     theUCbound[i] = -t_up; // reset bound to original value
746  	
747  	                     if(t_up == t_low)
748  	                        theLCbound[i] = -t_low; // for fixed variables we change both bounds
749  	                  }
750  	                  else
751  	                     theShift += theUCbound[i] - (-t_up);
752  	
753  	                  if((*thePvec)[i] > -t_low - eps)
754  	                  {
755  	                     theLCbound[i] = -t_low; // reset bound to original value
756  	
757  	                     if(t_up == t_low)
758  	                        theUCbound[i] = -t_up; // for fixed variables we change both bounds
759  	                  }
760  	                  else
761  	                     theShift += (-t_low) - theLCbound[i];
762  	               }
763  	               else if(theUCbound[i] > -t_up)
764  	                  theShift += theUCbound[i] - (-t_up);
765  	               else if(theLCbound[i] < -t_low)
766  	                  theShift += (-t_low) - theLCbound[i];
767  	            }
768  	         }
769  	         else
770  	         {
771  	            assert(rep() == ROW);
772  	
773  	            for(i = this->nRows(); i-- > 0;)
774  	            {
775  	               t_up = this->rhs(i);
776  	               t_low = this->lhs(i);
777  	
778  	               if(t_up == t_low)
779  	               {
780  	                  if(theURbound[i] > t_up)
781  	                     theShift += theURbound[i] - t_up;
782  	                  else
783  	                     theShift += t_low - theLRbound[i];
784  	               }
785  	               else if(!isBasic(ds.rowStatus(i)))
786  	               {
787  	                  if((*thePvec)[i] < t_up + eps)
788  	                     theURbound[i] = t_up; // reset bound to original value
789  	                  else
790  	                     theShift += theURbound[i] - t_up;
791  	
792  	                  if((*thePvec)[i] > t_low - eps)
793  	                     theLRbound[i] = t_low; // reset bound to original value
794  	                  else
795  	                     theShift += t_low - theLRbound[i];
796  	               }
797  	               else if(theURbound[i] > t_up)
798  	                  theShift += theURbound[i] - t_up;
799  	               else if(theLRbound[i] < t_low)
800  	                  theShift += t_low - theLRbound[i];
801  	            }
802  	
803  	            for(i = this->nCols(); i-- > 0;)
804  	            {
805  	               t_up = this->upper(i);
806  	               t_low = this->lower(i);
807  	
808  	               if(t_up == t_low)
809  	               {
810  	                  if(theUCbound[i] > t_up)
811  	                     theShift += theUCbound[i] - t_up;
812  	                  else
813  	                     theShift += t_low - theLCbound[i];
814  	               }
815  	               else if(!isBasic(ds.colStatus(i)))
816  	               {
817  	                  if((*theCoPvec)[i] < t_up + eps)
818  	                     theUCbound[i] = t_up; // reset bound to original value
819  	                  else
820  	                     theShift += theUCbound[i] - t_up;
821  	
822  	                  if((*theCoPvec)[i] > t_low - eps)
823  	                     theLCbound[i] = t_low; // reset bound to original value
824  	                  else
825  	                     theShift += t_low - theLCbound[i];
826  	               }
827  	               else if(theUCbound[i] > t_up)
828  	                  theShift += theUCbound[i] - t_up;
829  	               else if(theLCbound[i] < t_low)
830  	                  theShift += t_low - theLCbound[i];
831  	            }
832  	         }
833  	      }
834  	   }
835  	}
836  	} // namespace soplex
837