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 "soplex/spxdefines.h"
18 #include "soplex/spxboundflippingrt.h"
19 #include "soplex/sorter.h"
20 #include "soplex/spxsolver.h"
21 #include "soplex/spxout.h"
22 #include "soplex/spxid.h"
23
24 namespace soplex
25 {
26
27 #define LOWSTAB 1e-10
28 #define MAX_RELAX_COUNT 2
29 #define LONGSTEP_FREQ 100
30
31
32 /** perform necessary bound flips to restore dual feasibility */
33 template <class R>
34 void SPxBoundFlippingRT<R>::flipAndUpdate(
35 int& nflips /**< number of bounds that should be flipped */
36 )
37 {
38 assert(nflips > 0);
39
40 // number of bound flips that are not performed
41 int skipped;
42
43 updPrimRhs.setup();
44 updPrimRhs.reDim(this->thesolver->dim());
45 updPrimVec.reDim(this->thesolver->dim());
46 updPrimRhs.clear();
47 updPrimVec.clear();
48
49 skipped = 0;
50
51 for(int i = 0; i < nflips; ++i)
52 {
53 int idx;
54 idx = breakpoints[i].idx;
55
56 if(idx < 0)
57 {
58 ++skipped;
59 continue;
60 }
61
62 R range;
63 R upper;
64 R lower;
65 R objChange = 0.0;
66 typename SPxBasisBase<R>::Desc::Status stat;
67 typename SPxBasisBase<R>::Desc& ds = this->thesolver->basis().desc();
68
69 range = 0;
70
71 if(breakpoints[i].src == PVEC)
72 {
73 assert(this->thesolver->rep() == SPxSolverBase<R>::COLUMN);
74 stat = ds.status(idx);
75 upper = this->thesolver->upper(idx);
76 lower = this->thesolver->lower(idx);
77
78 switch(stat)
79 {
80 case SPxBasisBase<R>::Desc::P_ON_UPPER :
81 ds.status(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
82 range = lower - upper;
83 assert((*this->thesolver->theLbound)[idx] == R(-infinity));
84 (*this->thesolver->theLbound)[idx] = (*this->thesolver->theUbound)[idx];
85 (*this->thesolver->theUbound)[idx] = R(infinity);
86 objChange = range * (*this->thesolver->theLbound)[idx];
87 break;
88
89 case SPxBasisBase<R>::Desc::P_ON_LOWER :
90 ds.status(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
91 range = upper - lower;
92 assert((*this->thesolver->theUbound)[idx] == R(infinity));
93 (*this->thesolver->theUbound)[idx] = (*this->thesolver->theLbound)[idx];
94 (*this->thesolver->theLbound)[idx] = R(-infinity);
95 objChange = range * (*this->thesolver->theUbound)[idx];
96 break;
97
98 default :
99 ++skipped;
100 MSG_WARNING((*this->thesolver->spxout),
101 (*this->thesolver->spxout) << "PVEC unexpected status: " << static_cast<int>(stat)
102 << " index: " << idx
103 << " val: " << this->thesolver->pVec()[idx]
104 << " upd: " << this->thesolver->pVec().delta()[idx]
105 << " lower: " << lower
106 << " upper: " << upper
107 << " bp.val: " << breakpoints[i].val
108 << std::endl;)
109 }
110
111 MSG_DEBUG(std::cout << "PVEC flipped from: " << stat
112 << " index: " << idx
113 << " val: " << this->thesolver->pVec()[idx]
114 << " upd: " << this->thesolver->pVec().delta()[idx]
115 << " lower: " << lower
116 << " upper: " << upper
117 << " bp.val: " << breakpoints[i].val
118 << " UCbound: " << this->thesolver->theUCbound[idx]
119 << " LCbound: " << this->thesolver->theLCbound[idx]
120 << std::endl;)
121 assert(spxAbs(range) < 1e20);
122 updPrimRhs.multAdd(range, this->thesolver->vector(idx));
123
124 if(objChange != 0.0)
125 this->thesolver->updateNonbasicValue(objChange);
126 }
127 else if(breakpoints[i].src == COPVEC)
128 {
129 assert(this->thesolver->rep() == SPxSolverBase<R>::COLUMN);
130 stat = ds.coStatus(idx);
131 upper = this->thesolver->rhs(idx);
132 lower = this->thesolver->lhs(idx);
133
134 switch(stat)
135 {
136 case SPxBasisBase<R>::Desc::P_ON_UPPER :
137 ds.coStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
138 range = lower - upper;
139 assert((*this->thesolver->theCoUbound)[idx] == R(infinity));
140 (*this->thesolver->theCoUbound)[idx] = -(*this->thesolver->theCoLbound)[idx];
141 (*this->thesolver->theCoLbound)[idx] = R(-infinity);
142 objChange = range * (*this->thesolver->theCoUbound)[idx];
143 break;
144
145 case SPxBasisBase<R>::Desc::P_ON_LOWER :
146 ds.coStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
147 range = upper - lower;
148 assert((*this->thesolver->theCoLbound)[idx] == R(-infinity));
149 (*this->thesolver->theCoLbound)[idx] = -(*this->thesolver->theCoUbound)[idx];
150 (*this->thesolver->theCoUbound)[idx] = R(infinity);
151 objChange = range * (*this->thesolver->theCoLbound)[idx];
152 break;
153
154 default :
155 ++skipped;
156 MSG_WARNING((*this->thesolver->spxout),
157 (*this->thesolver->spxout) << "COPVEC unexpected status: " << static_cast<int>(stat)
158 << " index: " << idx
159 << " val: " << this->thesolver->coPvec()[idx]
160 << " upd: " << this->thesolver->coPvec().delta()[idx]
161 << " lower: " << lower
162 << " upper: " << upper
163 << " bp.val: " << breakpoints[i].val
164 << std::endl;)
165 }
166
167 MSG_DEBUG(std::cout << "COPVEC flipped from: " << stat
168 << " index: " << idx
169 << " val: " << this->thesolver->coPvec()[idx]
170 << " upd: " << this->thesolver->coPvec().delta()[idx]
171 << " lower: " << lower
172 << " upper: " << upper
173 << " bp.val: " << breakpoints[i].val
174 << " URbound: " << this->thesolver->theURbound[idx]
175 << " LRbound: " << this->thesolver->theLRbound[idx]
176 << std::endl;)
177 assert(spxAbs(range) < 1e20);
178 updPrimRhs.setValue(idx, updPrimRhs[idx] - range);
179
180 if(objChange != 0.0)
181 this->thesolver->updateNonbasicValue(objChange);
182 }
183 else if(breakpoints[i].src == FVEC)
184 {
185 assert(this->thesolver->rep() == SPxSolverBase<R>::ROW);
186 SPxId baseId = this->thesolver->basis().baseId(idx);
187 int IdNumber;
188
189 if(baseId.isSPxRowId())
190 {
191 IdNumber = this->thesolver->number(SPxRowId(baseId));
192 stat = ds.rowStatus(IdNumber);
193 upper = this->thesolver->rhs(IdNumber);
194 lower = this->thesolver->lhs(IdNumber);
195
196 switch(stat)
197 {
198 case SPxBasisBase<R>::Desc::P_ON_UPPER :
199 ds.rowStatus(IdNumber) = SPxBasisBase<R>::Desc::P_ON_LOWER;
200 range = upper - lower;
201 assert(this->thesolver->theUBbound[idx] == R(infinity));
202 this->thesolver->theUBbound[idx] = -this->thesolver->theLBbound[idx];
203 this->thesolver->theLBbound[idx] = R(-infinity);
204 break;
205
206 case SPxBasisBase<R>::Desc::P_ON_LOWER :
207 ds.rowStatus(IdNumber) = SPxBasisBase<R>::Desc::P_ON_UPPER;
208 range = lower - upper;
209 assert(this->thesolver->theLBbound[idx] == R(-infinity));
210 this->thesolver->theLBbound[idx] = -this->thesolver->theUBbound[idx];
211 this->thesolver->theUBbound[idx] = R(infinity);
212 break;
213
214 default :
215 ++skipped;
216 MSG_WARNING((*this->thesolver->spxout),
217 (*this->thesolver->spxout) << "unexpected basis status: " << static_cast<int>(stat)
218 << " index: " << idx
219 << " val: " << this->thesolver->fVec()[idx]
220 << " upd: " << this->thesolver->fVec().delta()[idx]
221 << " lower: " << lower
222 << " upper: " << upper
223 << " bp.val: " << breakpoints[i].val
224 << std::endl;)
225 }
226 }
227 else
228 {
229 assert(baseId.isSPxColId());
230 IdNumber = this->thesolver->number(SPxColId(baseId));
231 stat = ds.colStatus(IdNumber);
232 upper = this->thesolver->upper(IdNumber);
233 lower = this->thesolver->lower(IdNumber);
234
235 switch(stat)
236 {
237 case SPxBasisBase<R>::Desc::P_ON_UPPER :
238 ds.colStatus(IdNumber) = SPxBasisBase<R>::Desc::P_ON_LOWER;
239 range = upper - lower;
240 assert(this->thesolver->theUBbound[idx] == R(infinity));
241 this->thesolver->theUBbound[idx] = -this->thesolver->theLBbound[idx];
242 this->thesolver->theLBbound[idx] = R(-infinity);
243 break;
244
245 case SPxBasisBase<R>::Desc::P_ON_LOWER :
246 ds.colStatus(IdNumber) = SPxBasisBase<R>::Desc::P_ON_UPPER;
247 range = lower - upper;
248 assert(this->thesolver->theLBbound[idx] == R(-infinity));
249 this->thesolver->theLBbound[idx] = -this->thesolver->theUBbound[idx];
250 this->thesolver->theUBbound[idx] = R(infinity);
251 break;
252
253 default :
254 ++skipped;
255 MSG_WARNING((*this->thesolver->spxout),
256 (*this->thesolver->spxout) << "FVEC unexpected status: " << static_cast<int>(stat)
257 << " index: " << idx
258 << " val: " << this->thesolver->fVec()[idx]
259 << " upd: " << this->thesolver->fVec().delta()[idx]
260 << " lower: " << lower
261 << " upper: " << upper
262 << " bp.val: " << breakpoints[i].val
263 << std::endl;)
264 }
265 }
266
267 MSG_DEBUG(std::cout << "basic row/col flipped from: " << stat
268 << " index: " << idx
269 << " val: " << this->thesolver->fVec()[idx]
270 << " upd: " << this->thesolver->fVec().delta()[idx]
271 << " lower: " << lower
272 << " upper: " << upper
273 << " bp.val: " << breakpoints[i].val
274 << std::endl;)
275 assert(spxAbs(range) < 1e20);
276 assert(updPrimRhs[idx] == 0);
277 updPrimRhs.add(idx, range);
278 }
279 }
280
281 nflips -= skipped;
282
283 if(nflips > 0)
284 {
285 if(this->thesolver->rep() == SPxSolverBase<R>::ROW)
286 {
287 assert(this->m_type == SPxSolverBase<R>::ENTER);
288 (*this->thesolver->theCoPrhs) -= updPrimRhs;
289 this->thesolver->setup4coSolve2(&updPrimVec, &updPrimRhs);
290 }
291 else
292 {
293 assert(this->thesolver->rep() == SPxSolverBase<R>::COLUMN);
294 assert(this->m_type == SPxSolverBase<R>::LEAVE);
295 (*this->thesolver->theFrhs) -= updPrimRhs;
296 this->thesolver->setup4solve2(&updPrimVec, &updPrimRhs);
297 }
298 }
299
300 return;
301 }
302
303 /** store all available pivots/breakpoints in an array (positive pivot search direction) */
304 template <class R>
305 void SPxBoundFlippingRT<R>::collectBreakpointsMax(
306 int& nBp, /**< number of found breakpoints so far */
307 int& minIdx, /**< index to current minimal breakpoint */
308 const int* idx, /**< pointer to indices of current VectorBase<R> */
309 int nnz, /**< number of nonzeros in current VectorBase<R> */
310 const R* upd, /**< pointer to update values of current VectorBase<R> */
311 const R* vec, /**< pointer to values of current VectorBase<R> */
312 const R* upp, /**< pointer to upper bound/rhs of current VectorBase<R> */
313 const R* low, /**< pointer to lower bound/lhs of current VectorBase<R> */
314 BreakpointSource src /**< type of VectorBase<R> (pVec, coPvec or fVec)*/
315 )
316 {
317 R minVal;
318 R curVal;
319 const int* last;
320
321 minVal = (nBp == 0) ? R(infinity) : breakpoints[minIdx].val;
322
323 last = idx + nnz;
324
325 for(; idx < last; ++idx)
326 {
327 int i = *idx;
328 R x = upd[i];
329
330 if(x > this->epsilon)
331 {
332 if(upp[i] < R(infinity))
333 {
334 R y = upp[i] - vec[i];
335 curVal = (y <= 0) ? this->fastDelta / x : (y + this->fastDelta) / x;
336 assert(curVal > 0);
337
338 breakpoints[nBp].idx = i;
339 breakpoints[nBp].src = src;
340 breakpoints[nBp].val = curVal;
341
342 if(curVal < minVal)
343 {
344 minVal = curVal;
345 minIdx = nBp;
346 }
347
348 nBp++;
349 }
350 }
351 else if(x < -this->epsilon)
352 {
353 if(low[i] > R(-infinity))
354 {
355 R y = low[i] - vec[i];
356 curVal = (y >= 0) ? -this->fastDelta / x : (y - this->fastDelta) / x;
357 assert(curVal > 0);
358
359 breakpoints[nBp].idx = i;
360 breakpoints[nBp].src = src;
361 breakpoints[nBp].val = curVal;
362
363 if(curVal < minVal)
364 {
365 minVal = curVal;
366 minIdx = nBp;
367 }
368
369 nBp++;
370 }
371 }
372
373 if(nBp >= breakpoints.size())
374 breakpoints.reSize(nBp * 2);
375 }
376
377 return;
378 }
379
380 /** store all available pivots/breakpoints in an array (negative pivot search direction) */
381 template <class R>
382 void SPxBoundFlippingRT<R>::collectBreakpointsMin(
383 int& nBp, /**< number of found breakpoints so far */
384 int& minIdx, /**< index to current minimal breakpoint */
385 const int* idx, /**< pointer to indices of current VectorBase<R> */
386 int nnz, /**< number of nonzeros in current VectorBase<R> */
387 const R* upd, /**< pointer to update values of current VectorBase<R> */
388 const R* vec, /**< pointer to values of current VectorBase<R> */
389 const R* upp, /**< pointer to upper bound/rhs of current VectorBase<R> */
390 const R* low, /**< pointer to lower bound/lhs of current VectorBase<R> */
391 BreakpointSource src /**< type of VectorBase<R> (pVec, coPvec or fVec)*/
392 )
393 {
394 R minVal;
395 R curVal;
396 const int* last;
397
|
(1) Event cond_false: |
Condition "nBp == 0", taking false branch. |
|
(2) Event neg_sink_parm_call: |
Passing "minIdx" to "operator []", which cannot accept a negative number. [details] |
398 minVal = (nBp == 0) ? R(infinity) : breakpoints[minIdx].val;
399
400 last = idx + nnz;
401
402 for(; idx < last; ++idx)
403 {
404 int i = *idx;
405 R x = upd[i];
406
407 if(x > this->epsilon)
408 {
409 if(low[i] > R(-infinity))
410 {
411 R y = low[i] - vec[i];
412
413 curVal = (y >= 0) ? this->fastDelta / x : (this->fastDelta - y) / x;
414 assert(curVal > 0);
415
416 breakpoints[nBp].idx = i;
417 breakpoints[nBp].src = src;
418 breakpoints[nBp].val = curVal;
419
420 if(curVal < minVal)
421 {
422 minVal = curVal;
423 minIdx = nBp;
424 }
425
426 nBp++;
427 }
428 }
429 else if(x < -this->epsilon)
430 {
431 if(upp[i] < R(infinity))
432 {
433 R y = upp[i] - vec[i];
434 curVal = (y <= 0) ? -this->fastDelta / x : -(y + this->fastDelta) / x;
435 assert(curVal > 0);
436
437 breakpoints[nBp].idx = i;
438 breakpoints[nBp].src = src;
439 breakpoints[nBp].val = curVal;
440
441 if(curVal < minVal)
442 {
443 minVal = curVal;
444 minIdx = nBp;
445 }
446
447 nBp++;
448 }
449 }
450
451 if(nBp >= breakpoints.size())
452 breakpoints.reSize(nBp * 2);
453 }
454
455 return;
456 }
457
458 /** get values for entering index and perform shifts if necessary */
459 template <class R>
460 bool SPxBoundFlippingRT<R>::getData(
461 R& val,
462 SPxId& enterId,
463 int idx,
464 R stab,
465 R degeneps,
466 const R* upd,
467 const R* vec,
468 const R* low,
469 const R* upp,
470 BreakpointSource src,
471 R max
472 )
473 {
474 if(src == PVEC)
475 {
476 this->thesolver->pVec()[idx] = this->thesolver->vector(idx) * this->thesolver->coPvec();
477 R x = upd[idx];
478
479 // skip breakpoint if it is too small
480 if(spxAbs(x) < stab)
481 {
482 return false;
483 }
484
485 enterId = this->thesolver->id(idx);
486 val = (max * x > 0) ? upp[idx] : low[idx];
487 val = (val - vec[idx]) / x;
488
489 if(upp[idx] == low[idx])
490 {
491 val = 0.0;
492
493 if(vec[idx] > upp[idx])
494 this->thesolver->theShift += vec[idx] - upp[idx];
495 else
496 this->thesolver->theShift += low[idx] - vec[idx];
497
498 this->thesolver->upBound()[idx] = this->thesolver->lpBound()[idx] = vec[idx];
499 }
500 else if((max > 0 && val < -degeneps) || (max < 0 && val > degeneps))
501 {
502 val = 0.0;
503
504 if(max * x > 0)
505 this->thesolver->shiftUPbound(idx, vec[idx]);
506 else
507 this->thesolver->shiftLPbound(idx, vec[idx]);
508 }
509 }
510 else // src == COPVEC
511 {
512 R x = upd[idx];
513
514 if(spxAbs(x) < stab)
515 {
516 return false;
517 }
518
519 enterId = this->thesolver->coId(idx);
520 val = (max * x > 0.0) ? upp[idx] : low[idx];
521 val = (val - vec[idx]) / x;
522
523 if(upp[idx] == low[idx])
524 {
525 val = 0.0;
526
527 if(vec[idx] > upp[idx])
528 this->thesolver->theShift += vec[idx] - upp[idx];
529 else
530 this->thesolver->theShift += low[idx] - vec[idx];
531
532 this->thesolver->ucBound()[idx] = this->thesolver->lcBound()[idx] = vec[idx];
533 }
534 else if((max > 0 && val < -degeneps) || (max < 0 && val > degeneps))
535 {
536 val = 0.0;
537
538 if(max * x > 0)
539 this->thesolver->shiftUCbound(idx, vec[idx]);
540 else
541 this->thesolver->shiftLCbound(idx, vec[idx]);
542 }
543 }
544
545 return true;
546 }
547
548 /** get values for leaving index and perform shifts if necessary */
549 template <class R>
550 bool SPxBoundFlippingRT<R>::getData(
551 R& val,
552 int& leaveIdx,
553 int idx,
554 R stab,
555 R degeneps,
556 const R* upd,
557 const R* vec,
558 const R* low,
559 const R* upp,
560 BreakpointSource src,
561 R max
562 )
563 {
564 assert(src == FVEC);
565
566 R x = upd[idx];
567
568 // skip breakpoint if it is too small
569 if(spxAbs(x) < stab)
570 {
571 return false;
572 }
573
574 leaveIdx = idx;
575 val = (max * x > 0) ? upp[idx] : low[idx];
576 val = (val - vec[idx]) / x;
577
578 if(upp[idx] == low[idx])
579 {
580 val = 0.0;
581 this->thesolver->shiftLBbound(idx, vec[idx]);
582 this->thesolver->shiftUBbound(idx, vec[idx]);
583 }
584 else if((max > 0 && val < -degeneps) || (max < 0 && val > degeneps))
585 {
586 val = 0.0;
587
588 if(this->thesolver->dualStatus(this->thesolver->baseId(idx)) != SPxBasisBase<R>::Desc::D_ON_BOTH)
589 {
590 if(max * x > 0)
591 this->thesolver->shiftUBbound(idx, vec[idx]);
592 else
593 this->thesolver->shiftLBbound(idx, vec[idx]);
594 }
595 }
596
597 return true;
598 }
599
600 /** determine entering row/column */
601 template <class R>
602 SPxId SPxBoundFlippingRT<R>::selectEnter(
603 R& val,
604 int leaveIdx,
605 bool polish
606 )
607 {
608 assert(this->m_type == SPxSolverBase<R>::LEAVE);
609 assert(this->thesolver->boundflips == 0);
610
611 // reset the history and try again to do some long steps
|
(1) Event cond_true: |
Condition "this->thesolver->leaveCount % 100 == 0", taking true branch. |
612 if(this->thesolver->leaveCount % LONGSTEP_FREQ == 0)
613 {
614 MSG_DEBUG(std::cout << "DLBFRT06 resetting long step history" << std::endl;)
615 flipPotential = 1;
616 }
617
|
(2) Event cond_false: |
Condition "!this->enableBoundFlips", taking false branch. |
|
(3) Event cond_false: |
Condition "polish", taking false branch. |
|
(4) Event cond_false: |
Condition "this->thesolver->rep() == soplex::SPxSolverBase<double>::ROW", taking false branch. |
|
(5) Event cond_false: |
Condition "this->flipPotential <= 0", taking false branch. |
618 if(!enableBoundFlips || polish || this->thesolver->rep() == SPxSolverBase<R>::ROW
619 || flipPotential <= 0)
620 {
621 MSG_DEBUG(std::cout << "DLBFRT07 switching to fast ratio test" << std::endl;)
622 return SPxFastRT<R>::selectEnter(val, leaveIdx, polish);
|
(6) Event if_end: |
End of if statement. |
623 }
624
625 const R* pvec = this->thesolver->pVec().get_const_ptr();
626 const R* pupd = this->thesolver->pVec().delta().values();
627 const int* pidx = this->thesolver->pVec().delta().indexMem();
628 int pupdnnz = this->thesolver->pVec().delta().size();
629 const R* lpb = this->thesolver->lpBound().get_const_ptr();
630 const R* upb = this->thesolver->upBound().get_const_ptr();
631
632 const R* cvec = this->thesolver->coPvec().get_const_ptr();
633 const R* cupd = this->thesolver->coPvec().delta().values();
634 const int* cidx = this->thesolver->coPvec().delta().indexMem();
635 int cupdnnz = this->thesolver->coPvec().delta().size();
636 const R* lcb = this->thesolver->lcBound().get_const_ptr();
637 const R* ucb = this->thesolver->ucBound().get_const_ptr();
638
639 this->resetTols();
640
641 R max;
642
643 // index in breakpoint array of minimal value (i.e. choice of normal RT)
644 int minIdx;
645
646 // temporary breakpoint data structure to make swaps possible
647 Breakpoint tmp;
648
649 // most stable pivot value in candidate set
650 R moststable;
651
652 // initialize invalid enterId
653 SPxId enterId;
654
655 // slope of objective function improvement
656 R slope;
657
658 // number of found breakpoints
659 int nBp;
660
661 // number of passed breakpoints
662 int npassedBp;
663
664 R degeneps;
665 R stab;
666 bool instable;
667
668 max = val;
669 val = 0.0;
670 moststable = 0.0;
671 nBp = 0;
|
(7) Event var_tested_neg: |
Assigning: "minIdx" = a negative value. |
| Also see events: |
[negative_returns] |
672 minIdx = -1;
673
674 // get breakpoints and and determine the index of the minimal value
|
(8) Event cond_false: |
Condition "max > 0", taking false branch. |
675 if(max > 0)
676 {
677 collectBreakpointsMax(nBp, minIdx, pidx, pupdnnz, pupd, pvec, upb, lpb, PVEC);
678 collectBreakpointsMax(nBp, minIdx, cidx, cupdnnz, cupd, cvec, ucb, lcb, COPVEC);
679 }
680 else
|
(9) Event else_branch: |
Reached else branch. |
681 {
682 collectBreakpointsMin(nBp, minIdx, pidx, pupdnnz, pupd, pvec, upb, lpb, PVEC);
|
(10) Event negative_returns: |
"minIdx" is passed to a parameter that cannot be negative. [details] |
| Also see events: |
[var_tested_neg] |
683 collectBreakpointsMin(nBp, minIdx, cidx, cupdnnz, cupd, cvec, ucb, lcb, COPVEC);
684 }
685
686 if(nBp == 0)
687 {
688 val = max;
689 return enterId;
690 }
691
692 assert(minIdx >= 0);
693
694 // swap smallest breakpoint to the front to skip the sorting phase if no bound flip is possible
695 tmp = breakpoints[minIdx];
696 breakpoints[minIdx] = breakpoints[0];
697 breakpoints[0] = tmp;
698
699 // get initial slope
700 slope = spxAbs(this->thesolver->fTest()[leaveIdx]);
701
702 if(slope == 0)
703 {
704 // this may only happen if SoPlex decides to make an instable pivot
705 assert(this->thesolver->instableLeaveNum >= 0);
706 // restore original slope
707 slope = spxAbs(this->thesolver->instableLeaveVal);
708 }
709
710 // set up structures for the quicksort implementation
711 BreakpointCompare compare;
712 compare.entry = breakpoints.get_const_ptr();
713
714 // pointer to end of sorted part of breakpoints
715 int sorted = 0;
716 // minimum number of entries that are supposed to be sorted by partial sort
717 int sortsize = 4;
718
719 // get all skipable breakpoints
720 for(npassedBp = 0; npassedBp < nBp && slope > 0; ++npassedBp)
721 {
722 // sort breakpoints only partially to save time
723 if(npassedBp > sorted)
724 {
725 sorted = SPxQuicksortPart(breakpoints.get_ptr(), compare, sorted + 1, nBp, sortsize);
726 }
727
728 int i = breakpoints[npassedBp].idx;
729
730 // compute new slope
731 if(breakpoints[npassedBp].src == PVEC)
732 {
733 if(this->thesolver->isBasic(i))
734 {
735 // mark basic indices
736 breakpoints[npassedBp].idx = -1;
737 this->thesolver->pVec().delta().clearIdx(i);
738 }
739 else
740 {
741 R absupd = spxAbs(pupd[i]);
742 slope -= (this->thesolver->upper(i) * absupd) - (this->thesolver->lower(i) * absupd);
743
744 // get most stable pivot
745 if(absupd > moststable)
746 moststable = absupd;
747 }
748 }
749 else
750 {
751 assert(breakpoints[npassedBp].src == COPVEC);
752
753 if(this->thesolver->isCoBasic(i))
754 {
755 // mark basic indices
756 breakpoints[npassedBp].idx = -1;
757 this->thesolver->coPvec().delta().clearIdx(i);
758 }
759 else
760 {
761 R absupd = spxAbs(cupd[i]);
762 slope -= (this->thesolver->rhs(i) * absupd) - (this->thesolver->lhs(i) * absupd);
763
764 if(absupd > moststable)
765 moststable = absupd;
766 }
767 }
768 }
769
770 --npassedBp;
771 assert(npassedBp >= 0);
772
773 // check for unboundedness/infeasibility
774 if(slope > this->delta && npassedBp >= nBp - 1)
775 {
776 MSG_DEBUG(std::cout << "DLBFRT02 " << this->thesolver->basis().iteration()
777 << ": unboundedness in ratio test" << std::endl;)
778 flipPotential -= 0.5;
779 val = max;
780 return SPxFastRT<R>::selectEnter(val, leaveIdx);
781 }
782
783 MSG_DEBUG(std::cout << "DLBFRT01 "
784 << this->thesolver->basis().iteration()
785 << ": number of flip candidates: "
786 << npassedBp
787 << std::endl;)
788
789 // try to get a more stable pivot by looking at those with similar step length
790 int stableBp; // index to walk over additional breakpoints (after slope change)
791 int bestBp = -1; // breakpoints index with best possible stability
792 R bestDelta = breakpoints[npassedBp].val; // best step length (after bound flips)
793
794 for(stableBp = npassedBp + 1; stableBp < nBp; ++stableBp)
795 {
796 R stableDelta = 0;
797
798 // get next breakpoints in increasing order
799 if(stableBp > sorted)
800 {
801 sorted = SPxQuicksortPart(breakpoints.get_ptr(), compare, sorted + 1, nBp, sortsize);
802 }
803
804 int idx = breakpoints[stableBp].idx;
805
806 if(breakpoints[stableBp].src == PVEC)
807 {
808 if(this->thesolver->isBasic(idx))
809 {
810 // mark basic indices
811 breakpoints[stableBp].idx = -1;
812 this->thesolver->pVec().delta().clearIdx(idx);
813 continue;
814 }
815
816 R x = pupd[idx];
817
818 if(spxAbs(x) > moststable)
819 {
820 this->thesolver->pVec()[idx] = this->thesolver->vector(idx) * this->thesolver->coPvec();
821 stableDelta = (x > 0.0) ? upb[idx] : lpb[idx];
822 stableDelta = (stableDelta - pvec[idx]) / x;
823
824 if(stableDelta <= bestDelta)
825 {
826 moststable = spxAbs(x);
827 bestBp = stableBp;
828 }
829 }
830 }
831 else
832 {
833 if(this->thesolver->isCoBasic(idx))
834 {
835 // mark basic indices
836 breakpoints[stableBp].idx = -1;
837 this->thesolver->coPvec().delta().clearIdx(idx);
838 continue;
839 }
840
841 R x = cupd[idx];
842
843 if(spxAbs(x) > moststable)
844 {
845 stableDelta = (x > 0.0) ? ucb[idx] : lcb[idx];
846 stableDelta = (stableDelta - cvec[idx]) / x;
847
848 if(stableDelta <= bestDelta)
849 {
850 moststable = spxAbs(x);
851 bestBp = stableBp;
852 }
853 }
854 }
855
856 // stop searching if the step length is too big
857 if(stableDelta > this->delta + bestDelta)
858 break;
859 }
860
861 degeneps = this->fastDelta / moststable; /* as in SPxFastRT */
862 // get stability requirements
863 instable = this->thesolver->instableLeave;
864 assert(!instable || this->thesolver->instableLeaveNum >= 0);
865 stab = instable ? LOWSTAB : SPxFastRT<R>::minStability(moststable);
866
867 bool foundStable = false;
868
869 if(bestBp >= 0)
870 {
871 // found a more stable pivot
872 if(moststable > stab)
873 {
874 // stability requirements are satisfied
875 int idx = breakpoints[bestBp].idx;
876 assert(idx >= 0);
877
878 if(breakpoints[bestBp].src == PVEC)
879 foundStable = getData(val, enterId, idx, stab, degeneps, pupd, pvec, lpb, upb, PVEC, max);
880 else
881 foundStable = getData(val, enterId, idx, stab, degeneps, cupd, cvec, lcb, ucb, COPVEC, max);
882 }
883 }
884
885 else
886 {
887 // scan passed breakpoints from back to front and stop as soon as a good one is found
888 while(!foundStable && npassedBp >= 0)
889 {
890 int idx = breakpoints[npassedBp].idx;
891
892 // only look for non-basic variables
893 if(idx >= 0)
894 {
895 if(breakpoints[npassedBp].src == PVEC)
896 foundStable = getData(val, enterId, idx, stab, degeneps, pupd, pvec, lpb, upb, PVEC, max);
897 else
898 foundStable = getData(val, enterId, idx, stab, degeneps, cupd, cvec, lcb, ucb, COPVEC, max);
899 }
900
901 --npassedBp;
902 }
903
904 ++npassedBp;
905 }
906
907 if(!foundStable)
908 {
909 assert(!enterId.isValid());
910
911 if(relax_count < MAX_RELAX_COUNT)
912 {
913 MSG_DEBUG(std::cout << "DLBFRT04 "
914 << this->thesolver->basis().iteration()
915 << ": no valid enterId found - relaxing..."
916 << std::endl;)
917 this->relax();
918 ++relax_count;
919 // restore original value
920 val = max;
921 // try again with relaxed delta
922 return SPxBoundFlippingRT<R>::selectEnter(val, leaveIdx);
923 }
924 else
925 {
926 MSG_DEBUG(std::cout << "DLBFRT05 "
927 << this->thesolver->basis().iteration()
928 << " no valid enterId found - breaking..."
929 << std::endl;)
930 return enterId;
931 }
932 }
933 else
934 {
935 relax_count = 0;
936 this->tighten();
937 }
938
939 // flip bounds of skipped breakpoints only if a nondegenerate step is to be performed
940 if(npassedBp > 0 && spxAbs(breakpoints[npassedBp].val) > this->fastDelta)
941 {
942 flipAndUpdate(npassedBp);
943 this->thesolver->boundflips = npassedBp;
944
945 if(npassedBp >= 10)
946 flipPotential = 1;
947 else
948 flipPotential -= 0.05;
949 }
950 else
951 {
952 this->thesolver->boundflips = 0;
953 flipPotential -= 0.1;
954 }
955
956 MSG_DEBUG(std::cout << "DLBFRT06 "
957 << this->thesolver->basis().iteration()
958 << ": selected Id: "
959 << enterId
960 << " number of candidates: "
961 << nBp
962 << std::endl;)
963 return enterId;
964 }
965
966 /** determine leaving row/column */
967 template <class R>
968 int SPxBoundFlippingRT<R>::selectLeave(
969 R& val,
970 R enterTest,
971 bool polish
972 )
973 {
974 assert(this->m_type == SPxSolverBase<R>::ENTER);
975 assert(this->thesolver->boundflips == 0);
976
977 // reset the history and try again to do some long steps
978 if(this->thesolver->enterCount % LONGSTEP_FREQ == 0)
979 {
980 MSG_DEBUG(std::cout << "DEBFRT06 resetting long step history" << std::endl;)
981 flipPotential = 1;
982 }
983
984 if(polish || !enableBoundFlips || !enableRowBoundFlips
985 || this->thesolver->rep() == SPxSolverBase<R>::COLUMN || flipPotential <= 0)
986 {
987 MSG_DEBUG(std::cout << "DEBFRT07 switching to fast ratio test" << std::endl;)
988 return SPxFastRT<R>::selectLeave(val, enterTest, polish);
989 }
990
991 const R* vec =
992 this->thesolver->fVec().get_const_ptr(); /**< pointer to values of current VectorBase<R> */
993 const R* upd =
994 this->thesolver->fVec().delta().values(); /**< pointer to update values of current VectorBase<R> */
995 const int* idx =
996 this->thesolver->fVec().delta().indexMem(); /**< pointer to indices of current VectorBase<R> */
997 int updnnz =
998 this->thesolver->fVec().delta().size(); /**< number of nonzeros in update VectorBase<R> */
999 const R* lb =
1000 this->thesolver->lbBound().get_const_ptr(); /**< pointer to lower bound/lhs of current VectorBase<R> */
1001 const R* ub =
1002 this->thesolver->ubBound().get_const_ptr(); /**< pointer to upper bound/rhs of current VectorBase<R> */
1003
1004 this->resetTols();
1005
1006 R max;
1007
1008 // index in breakpoint array of minimal value (i.e. choice of normal RT)
1009 int minIdx;
1010
1011 // temporary breakpoint data structure to make swaps possible
1012 Breakpoint tmp;
1013
1014 // most stable pivot value in candidate set
1015 R moststable;
1016
1017 // initialize invalid leaving index
1018 int leaveIdx = -1;
1019
1020 // slope of objective function improvement
1021 R slope;
1022
1023 // number of found breakpoints
1024 int nBp;
1025
1026 // number of passed breakpoints
1027 int npassedBp;
1028
1029 R degeneps;
1030 R stab;
1031 bool instable;
1032
1033 max = val;
1034 val = 0.0;
1035 moststable = 0.0;
1036 nBp = 0;
1037 minIdx = -1;
1038
1039 assert(this->thesolver->fVec().delta().isSetup());
1040
1041 // get breakpoints and and determine the index of the minimal value
1042 if(max > 0)
1043 {
1044 collectBreakpointsMax(nBp, minIdx, idx, updnnz, upd, vec, ub, lb, FVEC);
1045 }
1046 else
1047 {
1048 collectBreakpointsMin(nBp, minIdx, idx, updnnz, upd, vec, ub, lb, FVEC);
1049 }
1050
1051 // return -1 if no BP was found
1052 if(nBp == 0)
1053 {
1054 val = max;
1055 return leaveIdx;
1056 }
1057
1058 assert(minIdx >= 0);
1059
1060 // swap smallest breakpoint to the front to skip the sorting phase if no bound flip is possible
1061 tmp = breakpoints[minIdx];
1062 breakpoints[minIdx] = breakpoints[0];
1063 breakpoints[0] = tmp;
1064
1065 // get initial slope
1066 slope = spxAbs(enterTest);
1067
1068 if(slope == 0)
1069 {
1070 // this may only happen if SoPlex decides to make an instable pivot
1071 assert(this->thesolver->instableEnterId.isValid());
1072 // restore original slope
1073 slope = this->thesolver->instableEnterVal;
1074 }
1075
1076 // set up structures for the quicksort implementation
1077 BreakpointCompare compare;
1078 compare.entry = breakpoints.get_const_ptr();
1079
1080 // pointer to end of sorted part of breakpoints
1081 int sorted = 0;
1082 // minimum number of entries that are supposed to be sorted by partial sort
1083 int sortsize = 4;
1084
1085 // get all skipable breakpoints
1086 for(npassedBp = 0; npassedBp < nBp && slope > 0; ++npassedBp)
1087 {
1088 // sort breakpoints only partially to save time
1089 if(npassedBp > sorted)
1090 {
1091 sorted = SPxQuicksortPart(breakpoints.get_ptr(), compare, sorted + 1, nBp, sortsize);
1092 }
1093
1094 assert(breakpoints[npassedBp].src == FVEC);
1095 int breakpointidx = breakpoints[npassedBp].idx;
1096 // compute new slope
1097 R upper;
1098 R lower;
1099 R absupd = spxAbs(upd[breakpointidx]);
1100 SPxId baseId = this->thesolver->baseId(breakpointidx);
1101 int i = this->thesolver->number(baseId);
1102
1103 if(baseId.isSPxColId())
1104 {
1105 upper = this->thesolver->upper(i);
1106 lower = this->thesolver->lower(i);
1107 }
1108 else
1109 {
1110 assert(baseId.isSPxRowId());
1111 upper = this->thesolver->rhs(i);
1112 lower = this->thesolver->lhs(i);
1113 }
1114
1115 slope -= (upper * absupd) - (lower * absupd);
1116
1117 // get most stable pivot
1118 if(absupd > moststable)
1119 moststable = absupd;
1120 }
1121
1122 --npassedBp;
1123 assert(npassedBp >= 0);
1124
1125 // check for unboundedness/infeasibility
1126 if(slope > this->delta && npassedBp >= nBp - 1)
1127 {
1128 MSG_DEBUG(std::cout << "DEBFRT02 " << this->thesolver->basis().iteration()
1129 << ": unboundedness in ratio test" << std::endl;)
1130 flipPotential -= 0.5;
1131 val = max;
1132 return SPxFastRT<R>::selectLeave(val, enterTest);
1133 }
1134
1135 MSG_DEBUG(std::cout << "DEBFRT01 "
1136 << this->thesolver->basis().iteration()
1137 << ": number of flip candidates: "
1138 << npassedBp
1139 << std::endl;)
1140
1141 // try to get a more stable pivot by looking at those with similar step length
1142 int stableBp; // index to walk over additional breakpoints (after slope change)
1143 int bestBp = -1; // breakpoints index with best possible stability
1144 R bestDelta = breakpoints[npassedBp].val; // best step length (after bound flips)
1145
1146 for(stableBp = npassedBp + 1; stableBp < nBp; ++stableBp)
1147 {
1148 R stableDelta = 0;
1149
1150 // get next breakpoints in increasing order
1151 if(stableBp > sorted)
1152 {
1153 sorted = SPxQuicksortPart(breakpoints.get_ptr(), compare, sorted + 1, nBp, sortsize);
1154 }
1155
1156 int breakpointidx = breakpoints[stableBp].idx;
1157 assert(breakpoints[stableBp].src == FVEC);
1158 R x = upd[breakpointidx];
1159
1160 if(spxAbs(x) > moststable)
1161 {
1162 stableDelta = (x > 0.0) ? ub[breakpointidx] : lb[breakpointidx];
1163 stableDelta = (stableDelta - vec[breakpointidx]) / x;
1164
1165 if(stableDelta <= bestDelta)
1166 {
1167 moststable = spxAbs(x);
1168 bestBp = stableBp;
1169 }
1170 }
1171 // stop searching if the step length is too big
1172 else if(stableDelta > this->delta + bestDelta)
1173 break;
1174 }
1175
1176 degeneps = this->fastDelta / moststable; /* as in SPxFastRT */
1177 // get stability requirements
1178 instable = this->thesolver->instableEnter;
1179 assert(!instable || this->thesolver->instableEnterId.isValid());
1180 stab = instable ? LOWSTAB : SPxFastRT<R>::minStability(moststable);
1181
1182 bool foundStable = false;
1183
1184 if(bestBp >= 0)
1185 {
1186 // found a more stable pivot
1187 if(moststable > stab)
1188 {
1189 // stability requirements are satisfied
1190 int breakpointidx = breakpoints[bestBp].idx;
1191 assert(breakpointidx >= 0);
1192 foundStable = getData(val, leaveIdx, breakpointidx, moststable, degeneps, upd, vec, lb, ub, FVEC,
1193 max);
1194 }
1195 }
1196
1197 else
1198 {
1199 // scan passed breakpoints from back to front and stop as soon as a good one is found
1200 while(!foundStable && npassedBp >= 0)
1201 {
1202 int breakpointidx = breakpoints[npassedBp].idx;
1203
1204 // only look for non-basic variables
1205 if(breakpointidx >= 0)
1206 {
1207 foundStable = getData(val, leaveIdx, breakpointidx, moststable, degeneps, upd, vec, lb, ub, FVEC,
1208 max);
1209 }
1210
1211 --npassedBp;
1212 }
1213
1214 ++npassedBp;
1215 }
1216
1217 if(!foundStable)
1218 {
1219 assert(leaveIdx < 0);
1220
1221 if(relax_count < MAX_RELAX_COUNT)
1222 {
1223 MSG_DEBUG(std::cout << "DEBFRT04 "
1224 << this->thesolver->basis().iteration()
1225 << ": no valid leaveIdx found - relaxing..."
1226 << std::endl;)
1227 this->relax();
1228 ++relax_count;
1229 // restore original value
1230 val = max;
1231 // try again with relaxed delta
1232 return SPxBoundFlippingRT<R>::selectLeave(val, enterTest);
1233 }
1234 else
1235 {
1236 MSG_DEBUG(std::cout << "DEBFRT05 "
1237 << this->thesolver->basis().iteration()
1238 << " no valid leaveIdx found - breaking..."
1239 << std::endl;)
1240 return leaveIdx;
1241 }
1242 }
1243 else
1244 {
1245 relax_count = 0;
1246 this->tighten();
1247 }
1248
1249 // flip bounds of skipped breakpoints only if a nondegenerate step is to be performed
1250 if(npassedBp > 0 && spxAbs(breakpoints[npassedBp].val) > this->fastDelta)
1251 {
1252 flipAndUpdate(npassedBp);
1253 this->thesolver->boundflips = npassedBp;
1254
1255 if(npassedBp >= 10)
1256 flipPotential = 1;
1257 else
1258 flipPotential -= 0.05;
1259 }
1260 else
1261 {
1262 this->thesolver->boundflips = 0;
1263 flipPotential -= 0.1;
1264 }
1265
1266 MSG_DEBUG(std::cout << "DEBFRT06 "
1267 << this->thesolver->basis().iteration()
1268 << ": selected Index: "
1269 << leaveIdx
1270 << " number of candidates: "
1271 << nBp
1272 << std::endl;)
1273
1274 return leaveIdx;
1275 }
1276
1277
1278 } // namespace soplex
1279