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