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   	
30   	namespace soplex
31   	{
32   	/**
33   	 * Here comes the ratio test for selecting a variable to leave the basis.
34   	 * It is assumed that Vec.delta() and fVec.idx() have been setup
35   	 * correctly!
36   	 *
37   	 * The leaving variable is selected such that the update of fVec() (using
38   	 * fVec.value() * fVec.delta()) keeps the basis feasible within
39   	 * solver()->entertol(). Hence, fVec.value() must be chosen such that one
40   	 * updated value of theFvec just reaches its bound and no other one exceeds
41   	 * them by more than solver()->entertol(). Further, fVec.value() must have the
42   	 * same sign as argument \p val.
43   	 *
44   	 * The return value of selectLeave() is the number of a variable in the
45   	 * basis selected to leave the basis. -1 indicates that no variable could be
46   	 * selected. Otherwise, parameter \p val contains the chosen fVec.value().
47   	 */
48   	template <class R>
49   	int SPxDefaultRT<R>::selectLeave(R& val, R, bool)
50   	{
51   	   this->solver()->fVec().delta().setup();
52   	
53   	   const R*   vec = this->solver()->fVec().get_const_ptr();
54   	   const R*   upd = this->solver()->fVec().delta().values();
55   	   const IdxSet& idx = this->solver()->fVec().idx();
56   	   const R*   ub  = this->solver()->ubBound().get_const_ptr();
57   	   const R*   lb  = this->solver()->lbBound().get_const_ptr();
58   	
59   	   R epsilon = this->solver()->epsilon();
60   	   int  leave   = -1;
61   	
62   	   R x;
63   	   int  i;
64   	   int  j;
65   	
66   	   // PARALLEL the j loop could be parallelized
67   	   if(val > 0)
68   	   {
69   	      // Loop over NZEs of delta vector.
70   	      for(j = 0; j < idx.size(); ++j)
71   	      {
72   	         i = idx.index(j);
73   	         x = upd[i];
74   	
75   	         if(x > epsilon)
76   	         {
77   	            if(ub[i] < R(infinity))
78   	            {
79   	               R y = (ub[i] - vec[i] + this->delta) / x;
80   	
81   	               if(y < val)
82   	               {
83   	                  leave = i;
84   	                  val   = y;
85   	               }
86   	            }
87   	         }
88   	         else if(x < -epsilon)
89   	         {
90   	            if(lb[i] > R(-infinity))
91   	            {
92   	               R y = (lb[i] - vec[i] - this->delta) / x;
93   	
94   	               if(y < val)
95   	               {
96   	                  leave = i;
97   	                  val   = y;
98   	               }
99   	            }
100  	         }
101  	      }
102  	
103  	      if(leave >= 0)
104  	      {
105  	         x   = upd[leave];
106  	
107  	         // BH 2005-11-30: It may well happen that the basis is degenerate and the
108  	         // selected leaving variable is (at most this->delta) beyond its bound. (This
109  	         // happens for instance on LP/netlib/adlittle.mps with setting -r -t0.)
110  	         // In that case we do a pivot step with length zero to avoid difficulties.
111  	         if((x > epsilon  && vec[leave] >= ub[leave]) ||
112  	               (x < -epsilon && vec[leave] <= lb[leave]))
113  	         {
114  	            val = 0.0;
115  	         }
116  	         else
117  	         {
118  	            val = (x > epsilon) ? ub[leave] : lb[leave];
119  	            val = (val - vec[leave]) / x;
120  	         }
121  	      }
122  	
123  	      SOPLEX_ASSERT_WARN("WDEFRT01", val > -epsilon);
124  	   }
125  	   else
126  	   {
127  	      for(j = 0; j < idx.size(); ++j)
128  	      {
129  	         i = idx.index(j);
130  	         x = upd[i];
131  	
132  	         if(x < -epsilon)
133  	         {
134  	            if(ub[i] < R(infinity))
135  	            {
136  	               R y = (ub[i] - vec[i] + this->delta) / x;
137  	
138  	               if(y > val)
139  	               {
140  	                  leave = i;
141  	                  val   = y;
142  	               }
143  	            }
144  	         }
145  	         else if(x > epsilon)
146  	         {
147  	            if(lb[i] > R(-infinity))
148  	            {
149  	               R y = (lb[i] - vec[i] - this->delta) / x;
150  	
151  	               if(y > val)
152  	               {
153  	                  leave = i;
154  	                  val   = y;
155  	               }
156  	            }
157  	         }
158  	      }
159  	
160  	      if(leave >= 0)
161  	      {
162  	         x   = upd[leave];
163  	
164  	         // See comment above.
165  	         if((x < -epsilon && vec[leave] >= ub[leave]) ||
166  	               (x > epsilon  && vec[leave] <= lb[leave]))
167  	         {
168  	            val = 0.0;
169  	         }
170  	         else
171  	         {
172  	            val = (x < epsilon) ? ub[leave] : lb[leave];
173  	            val = (val - vec[leave]) / x;
174  	         }
175  	      }
176  	
177  	      SOPLEX_ASSERT_WARN("WDEFRT02", val < epsilon);
178  	   }
179  	
180  	   return leave;
181  	}
182  	
183  	/**
184  	   Here comes the ratio test. It is assumed that theCoPvec.this->delta() and
185  	   theCoPvec.idx() have been setup correctly!
186  	*/
187  	template <class R>
188  	SPxId SPxDefaultRT<R>::selectEnter(R& max, int, bool)
189  	{
190  	   this->solver()->coPvec().delta().setup();
191  	   this->solver()->pVec().delta().setup();
192  	
193  	   const R*   pvec = this->solver()->pVec().get_const_ptr();
194  	   const R*   pupd = this->solver()->pVec().delta().values();
195  	   const IdxSet& pidx = this->solver()->pVec().idx();
196  	   const R*   lpb  = this->solver()->lpBound().get_const_ptr();
197  	   const R*   upb  = this->solver()->upBound().get_const_ptr();
198  	
199  	   const R*   cvec = this->solver()->coPvec().get_const_ptr();
200  	   const R*   cupd = this->solver()->coPvec().delta().values();
201  	   const IdxSet& cidx = this->solver()->coPvec().idx();
202  	   const R*   lcb  = this->solver()->lcBound().get_const_ptr();
203  	   const R*   ucb  = this->solver()->ucBound().get_const_ptr();
204  	
205  	   R epsilon = this->solver()->epsilon();
206  	   R val     = max;
207  	   int  pnum    = -1;
208  	   int  cnum    = -1;
209  	
210  	   SPxId enterId;
211  	   int   i;
212  	   int   j;
213  	   R  x;
214  	
215  	   // PARALLEL the j loops could be parallelized
216  	   if(val > 0)
217  	   {
218  	      for(j = 0; j < pidx.size(); ++j)
219  	      {
220  	         i = pidx.index(j);
221  	         x = pupd[i];
222  	
223  	         if(x > epsilon)
224  	         {
225  	            if(upb[i] < R(infinity))
226  	            {
227  	               R y = (upb[i] - pvec[i] + this->delta) / x;
228  	
229  	               if(y < val)
230  	               {
231  	                  enterId = this->solver()->id(i);
232  	                  val     = y;
233  	                  pnum    = j;
234  	               }
235  	            }
236  	         }
237  	         else if(x < -epsilon)
238  	         {
239  	            if(lpb[i] > R(-infinity))
240  	            {
241  	               R y = (lpb[i] - pvec[i] - this->delta) / x;
242  	
243  	               if(y < val)
244  	               {
245  	                  enterId = this->solver()->id(i);
246  	                  val     = y;
247  	                  pnum    = j;
248  	               }
249  	            }
250  	         }
251  	      }
252  	
253  	      for(j = 0; j < cidx.size(); ++j)
254  	      {
255  	         i = cidx.index(j);
256  	         x = cupd[i];
257  	
258  	         if(x > epsilon)
259  	         {
260  	            if(ucb[i] < R(infinity))
261  	            {
262  	               R y = (ucb[i] - cvec[i] + this->delta) / x;
263  	
264  	               if(y < val)
265  	               {
266  	                  enterId = this->solver()->coId(i);
267  	                  val     = y;
268  	                  cnum    = j;
269  	               }
270  	            }
271  	         }
272  	         else if(x < -epsilon)
273  	         {
274  	            if(lcb[i] > R(-infinity))
275  	            {
276  	               R y = (lcb[i] - cvec[i] - this->delta) / x;
277  	
278  	               if(y < val)
279  	               {
280  	                  enterId = this->solver()->coId(i);
281  	                  val     = y;
282  	                  cnum    = j;
283  	               }
284  	            }
285  	         }
286  	      }
287  	
288  	      if(cnum >= 0)
289  	      {
290  	         i   = cidx.index(cnum);
291  	         x   = cupd[i];
292  	         val = (x > epsilon) ? ucb[i] : lcb[i];
293  	         val = (val - cvec[i]) / x;
294  	      }
295  	      else if(pnum >= 0)
296  	      {
297  	         i   = pidx.index(pnum);
298  	         x   = pupd[i];
299  	         val = (x > epsilon) ? upb[i] : lpb[i];
300  	         val = (val - pvec[i]) / x;
301  	      }
302  	   }
303  	   else
304  	   {
305  	      for(j = 0; j < pidx.size(); ++j)
306  	      {
307  	         i = pidx.index(j);
308  	         x = pupd[i];
309  	
310  	         if(x > epsilon)
311  	         {
312  	            if(lpb[i] > R(-infinity))
313  	            {
314  	               R y = (lpb[i] - pvec[i] - this->delta) / x;
315  	
316  	               if(y > val)
317  	               {
318  	                  enterId = this->solver()->id(i);
319  	                  val     = y;
320  	                  pnum    = j;
321  	               }
322  	            }
323  	         }
324  	         else if(x < -epsilon)
325  	         {
326  	            if(upb[i] < R(infinity))
327  	            {
328  	               R y = (upb[i] - pvec[i] + this->delta) / x;
329  	
330  	               if(y > val)
331  	               {
332  	                  enterId = this->solver()->id(i);
333  	                  val     = y;
334  	                  pnum    = j;
335  	               }
336  	            }
337  	         }
338  	      }
339  	
340  	      for(j = 0; j < cidx.size(); ++j)
341  	      {
342  	         i = cidx.index(j);
343  	         x = cupd[i];
344  	
345  	         if(x > epsilon)
346  	         {
347  	            if(lcb[i] > R(-infinity))
348  	            {
349  	               R y = (lcb[i] - cvec[i] - this->delta) / x;
350  	
351  	               if(y > val)
352  	               {
353  	                  enterId = this->solver()->coId(i);
354  	                  val     = y;
355  	                  cnum    = j;
356  	               }
357  	            }
358  	         }
359  	         else if(x < -epsilon)
360  	         {
361  	            if(ucb[i] < R(infinity))
362  	            {
363  	               R y = (ucb[i] - cvec[i] + this->delta) / x;
364  	
365  	               if(y > val)
366  	               {
367  	                  enterId = this->solver()->coId(i);
368  	                  val     = y;
369  	                  cnum    = j;
370  	               }
371  	            }
372  	         }
373  	      }
374  	
375  	      if(cnum >= 0)
376  	      {
377  	         i   = cidx.index(cnum);
378  	         x   = cupd[i];
379  	         val = (x < epsilon) ? ucb[i] : lcb[i];
380  	         val = (val - cvec[i]) / x;
381  	      }
382  	      else if(pnum >= 0)
383  	      {
384  	         i   = pidx.index(pnum);
385  	         x   = pupd[i];
386  	         val = (x < epsilon) ? upb[i] : lpb[i];
387  	         val = (val - pvec[i]) / x;
388  	      }
389  	   }
390  	
391  	   if(enterId.isValid() && this->solver()->isBasic(enterId))
392  	   {
393  	      SPxOut::debug(this, "DDEFRT01 isValid() && isBasic(): max={}\n", max);
394  	
395  	      if(cnum >= 0)
396  	         this->solver()->coPvec().delta().clearNum(cnum);
397  	      else if(pnum >= 0)
398  	         this->solver()->pVec().delta().clearNum(pnum);
399  	
400  	      return SPxDefaultRT<R>::selectEnter(max, 0, false);
401  	   }
402  	
403  	   SPX_DEBUG(
404  	
405  	      if(!enterId.isValid())
406  	      std::cout << "DDEFRT02 !isValid(): max=" << max << ", x=" << x << std::endl;
407  	   )
408  	      max = val;
409  	
410  	   return enterId;
411  	}
412  	
413  	} // namespace soplex
414