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 /* Updating the Basis for Leaving Variables
17 */
18 #include <assert.h>
19 #include <stdio.h>
20
21 #include "soplex/spxdefines.h"
22 #include "soplex/spxpricer.h"
23 #include "soplex/spxsolver.h"
24 #include "soplex/spxratiotester.h"
25 #include "soplex/spxout.h"
26 #include "soplex/exceptions.h"
27
28 namespace soplex
29 {
30 static const Real reject_leave_tol = 1e-10; // = LOWSTAB as defined in spxfastrt.hpp
31
32 /*
33 VectorBase<R> |fTest| gives the feasibility test of all basic variables. For its
34 computation |fVec|, |theUBbound| and |theLBbound| must be setup correctly.
35 Values of |fTest| $<0$ represent infeasible variables, which are eligible
36 for leaving the basis in the simplex loop.
37 */
38 template <class R>
39 void SPxSolverBase<R>::computeFtest()
40 {
41
42 assert(type() == LEAVE);
43
44 R theeps = entertol();
45 m_pricingViolUpToDate = true;
46 m_pricingViolCoUpToDate = true;
47 m_pricingViol = 0;
48 m_pricingViolCo = 0;
49 m_numViol = 0;
50 infeasibilities.clear();
51 int sparsitythreshold = (int)(sparsePricingFactor * dim());
52
53 for(int i = 0; i < dim(); ++i)
54 {
55 theCoTest[i] = ((*theFvec)[i] > theUBbound[i])
56 ? theUBbound[i] - (*theFvec)[i]
57 : (*theFvec)[i] - theLBbound[i];
58
59 if(remainingRoundsLeave == 0)
60 {
61 if(theCoTest[i] < -theeps)
62 {
63 m_pricingViol -= theCoTest[i];
64 infeasibilities.addIdx(i);
65 isInfeasible[i] = SPxPricer<R>::VIOLATED;
66 ++m_numViol;
67 }
68 else
69 isInfeasible[i] = SPxPricer<R>::NOT_VIOLATED;
70
71 if(infeasibilities.size() > sparsitythreshold)
72 {
73 MSG_INFO2((*this->spxout), (*this->spxout) << " --- using dense pricing"
74 << std::endl;)
75 remainingRoundsLeave = DENSEROUNDS;
76 sparsePricingLeave = false;
77 infeasibilities.clear();
78 }
79 }
80 else if(theCoTest[i] < -theeps)
81 {
82 m_pricingViol -= theCoTest[i];
83 m_numViol++;
84 }
85 }
86
87 if(infeasibilities.size() == 0 && !sparsePricingLeave)
88 {
89 --remainingRoundsLeave;
90 }
91 else if(infeasibilities.size() <= sparsitythreshold && !sparsePricingLeave)
92 {
93 MSG_INFO2((*this->spxout),
94 std::streamsize prec = spxout->precision();
95
96 if(hyperPricingLeave)
97 (*this->spxout) << " --- using hypersparse pricing, ";
98 else
99 (*this->spxout) << " --- using sparse pricing, ";
100 (*this->spxout) << "sparsity: "
101 << std::setw(6) << std::fixed << std::setprecision(4)
102 << (R) m_numViol / dim()
103 << std::scientific << std::setprecision(int(prec))
104 << std::endl;
105 )
106 sparsePricingLeave = true;
107 }
108 }
109
110 template <class R>
111 void SPxSolverBase<R>::updateFtest()
112 {
113 const IdxSet& idx = theFvec->idx();
114 VectorBase<R>& ftest = theCoTest; // |== fTest()|
115 assert(&ftest == &fTest());
116
117 assert(type() == LEAVE);
118
119 updateViols.clear();
120 R theeps = entertol();
121
122 for(int j = idx.size() - 1; j >= 0; --j)
123 {
124 int i = idx.index(j);
125
126 if(m_pricingViolUpToDate && ftest[i] < -theeps)
127 // violation was present before this iteration
128 m_pricingViol += ftest[i];
129
130 ftest[i] = ((*theFvec)[i] > theUBbound[i])
131 ? theUBbound[i] - (*theFvec)[i]
132 : (*theFvec)[i] - theLBbound[i];
133
134 if(sparsePricingLeave && ftest[i] < -theeps)
135 {
136 assert(remainingRoundsLeave == 0);
137
138 if(m_pricingViolUpToDate)
139 m_pricingViol -= ftest[i];
140
141 if(isInfeasible[i] == SPxPricer<R>::NOT_VIOLATED)
142 {
143 // this can cause problems - we cannot keep on adding indeces to infeasibilities,
144 // because they are not deleted in hyper mode...
145 // if( !hyperPricingLeave )
146 infeasibilities.addIdx(i);
147 isInfeasible[i] = SPxPricer<R>::VIOLATED;
148 }
149
150 if(hyperPricingLeave)
151 updateViols.addIdx(i);
152 }
153 else if(m_pricingViolUpToDate && ftest[i] < -theeps)
154 m_pricingViol -= ftest[i];
155 }
156
157 // if boundflips were performed, we need to update these indices as well
158 if(boundflips > 0)
159 {
160 R eps = epsilon();
161
162 for(int j = 0; j < solveVector3->size(); ++j)
163 {
164 if(spxAbs(solveVector3->value(j)) > eps)
165 {
166 int i = solveVector3->index(j);
167
168 if(m_pricingViolUpToDate && ftest[i] < -theeps)
169 m_pricingViol += ftest[i];
170
171 ftest[i] = ((*theFvec)[i] > theUBbound[i]) ? theUBbound[i] - (*theFvec)[i] :
172 (*theFvec)[i] - theLBbound[i];
173
174 if(sparsePricingLeave && ftest[i] < -theeps)
175 {
176 assert(remainingRoundsLeave == 0);
177
178 if(m_pricingViolUpToDate)
179 m_pricingViol -= ftest[i];
180
181 if(!isInfeasible[i])
182 {
183 infeasibilities.addIdx(i);
184 isInfeasible[i] = true;
185 }
186 }
187 else if(m_pricingViolUpToDate && ftest[i] < -theeps)
188 m_pricingViol -= ftest[i];
189 }
190 }
191 }
192 }
193
194
195 /* compute statistics on leaving variable
196 Compute a set of statistical values on the variable selected for leaving the
197 basis.
198 */
199 template <class R>
200 void SPxSolverBase<R>::getLeaveVals(
201 int leaveIdx,
202 typename SPxBasisBase<R>::Desc::Status& leaveStat,
203 SPxId& leaveId,
204 R& leaveMax,
205 R& leavebound,
206 int& leaveNum,
207 StableSum<R>& objChange)
208 {
209 typename SPxBasisBase<R>::Desc& ds = this->desc();
210 leaveId = this->baseId(leaveIdx);
211
212 if(leaveId.isSPxRowId())
213 {
214 leaveNum = this->number(SPxRowId(leaveId));
215 leaveStat = ds.rowStatus(leaveNum);
216
217 assert(isBasic(leaveStat));
218
219 switch(leaveStat)
220 {
221 case SPxBasisBase<R>::Desc::P_ON_UPPER :
222 assert(rep() == ROW);
223 ds.rowStatus(leaveNum) = this->dualRowStatus(leaveNum);
224 leavebound = 0;
225 leaveMax = R(-infinity);
226 break;
227
228 case SPxBasisBase<R>::Desc::P_ON_LOWER :
229 assert(rep() == ROW);
230 ds.rowStatus(leaveNum) = this->dualRowStatus(leaveNum);
231 leavebound = 0;
232 leaveMax = R(infinity);
233 break;
234
235 case SPxBasisBase<R>::Desc::P_FREE :
236 assert(rep() == ROW);
237 throw SPxInternalCodeException("XLEAVE01 This should never happen.");
238
239 case SPxBasisBase<R>::Desc::D_FREE :
240 assert(rep() == COLUMN);
241 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_FIXED;
242 assert(this->lhs(leaveNum) == this->rhs(leaveNum));
243 leavebound = -this->rhs(leaveNum);
244
245 if((*theFvec)[leaveIdx] < theLBbound[leaveIdx])
246 leaveMax = R(infinity);
247 else
248 leaveMax = R(-infinity);
249
250 break;
251
252 case SPxBasisBase<R>::Desc::D_ON_LOWER :
253 assert(rep() == COLUMN);
254 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER;
255 leavebound = -this->rhs(leaveNum); // slack !!
256 leaveMax = R(infinity);
257 objChange += theLRbound[leaveNum] * this->rhs(leaveNum);
258 break;
259
260 case SPxBasisBase<R>::Desc::D_ON_UPPER :
261 assert(rep() == COLUMN);
262 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER;
263 leavebound = -this->lhs(leaveNum); // slack !!
264 leaveMax = R(-infinity);
265 objChange += theURbound[leaveNum] * this->lhs(leaveNum);
266 break;
267
268 case SPxBasisBase<R>::Desc::D_ON_BOTH :
269 assert(rep() == COLUMN);
270
271 if((*theFvec)[leaveIdx] > theLBbound[leaveIdx])
272 {
273 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER;
274 theLRbound[leaveNum] = R(-infinity);
275 leavebound = -this->lhs(leaveNum); // slack !!
276 leaveMax = R(-infinity);
277 objChange += theURbound[leaveNum] * this->lhs(leaveNum);
278 }
279 else
280 {
281 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER;
282 theURbound[leaveNum] = R(infinity);
283 leavebound = -this->rhs(leaveNum); // slack !!
284 leaveMax = R(infinity);
285 objChange += theLRbound[leaveNum] * this->rhs(leaveNum);
286 }
287
288 break;
289
290 default:
291 throw SPxInternalCodeException("XLEAVE02 This should never happen.");
292 }
293
294 MSG_DEBUG(std::cout << "DLEAVE51 SPxSolverBase<R>::getLeaveVals() : row " << leaveNum
295 << ": " << leaveStat
296 << " -> " << ds.rowStatus(leaveNum)
297 << " objChange: " << objChange
298 << std::endl;)
299 }
300
301 else
302 {
303 assert(leaveId.isSPxColId());
304 leaveNum = this->number(SPxColId(leaveId));
305 leaveStat = ds.colStatus(leaveNum);
306
307 assert(isBasic(leaveStat));
308
309 switch(leaveStat)
310 {
311 case SPxBasisBase<R>::Desc::P_ON_UPPER :
312 assert(rep() == ROW);
313 ds.colStatus(leaveNum) = this->dualColStatus(leaveNum);
314 leavebound = 0;
315 leaveMax = R(-infinity);
316 break;
317
318 case SPxBasisBase<R>::Desc::P_ON_LOWER :
319 assert(rep() == ROW);
320 ds.colStatus(leaveNum) = this->dualColStatus(leaveNum);
321 leavebound = 0;
322 leaveMax = R(infinity);
323 break;
324
325 case SPxBasisBase<R>::Desc::P_FREE :
326 assert(rep() == ROW);
327 ds.colStatus(leaveNum) = this->dualColStatus(leaveNum);
328
329 if((*theFvec)[leaveIdx] < theLBbound[leaveIdx])
330 {
331 leavebound = theLBbound[leaveIdx];
332 leaveMax = R(-infinity);
333 }
334 else
335 {
336 leavebound = theUBbound[leaveIdx];
337 leaveMax = R(infinity);
338 }
339
340 break;
341
342 case SPxBasisBase<R>::Desc::D_FREE :
343 assert(rep() == COLUMN);
344 assert(SPxLPBase<R>::upper(leaveNum) == SPxLPBase<R>::lower(leaveNum));
345 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_FIXED;
346 leavebound = SPxLPBase<R>::upper(leaveNum);
347 objChange += this->maxObj(leaveNum) * leavebound;
348
349 if((*theFvec)[leaveIdx] < theLBbound[leaveIdx])
350 leaveMax = R(infinity);
351 else
352 leaveMax = R(-infinity);
353
354 break;
355
356 case SPxBasisBase<R>::Desc::D_ON_LOWER :
357 assert(rep() == COLUMN);
358 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER;
359 leavebound = SPxLPBase<R>::upper(leaveNum);
360 objChange += theUCbound[leaveNum] * leavebound;
361 leaveMax = R(-infinity);
362 break;
363
364 case SPxBasisBase<R>::Desc::D_ON_UPPER :
365 assert(rep() == COLUMN);
366 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER;
367 leavebound = SPxLPBase<R>::lower(leaveNum);
368 objChange += theLCbound[leaveNum] * leavebound;
369 leaveMax = R(infinity);
370 break;
371
372 case SPxBasisBase<R>::Desc::D_ON_BOTH :
373 assert(rep() == COLUMN);
374
375 if((*theFvec)[leaveIdx] > theUBbound[leaveIdx])
376 {
377 leaveMax = R(-infinity);
378 leavebound = SPxLPBase<R>::upper(leaveNum);
379 objChange += theUCbound[leaveNum] * leavebound;
380 theLCbound[leaveNum] = R(-infinity);
381 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER;
382 }
383 else
384 {
385 leaveMax = R(infinity);
386 leavebound = SPxLPBase<R>::lower(leaveNum);
387 objChange += theLCbound[leaveNum] * leavebound;
388 theUCbound[leaveNum] = R(infinity);
389 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER;
390 }
391
392 break;
393
394 default:
395 throw SPxInternalCodeException("XLEAVE03 This should never happen.");
396 }
397
398 MSG_DEBUG(std::cout << "DLEAVE52 SPxSolverBase<R>::getLeaveVals() : col " << leaveNum
399 << ": " << leaveStat
400 << " -> " << ds.colStatus(leaveNum)
401 << " objChange: " << objChange
402 << std::endl;)
403 }
404 }
405
406 template <class R>
407 void SPxSolverBase<R>::getLeaveVals2(
408 R leaveMax,
409 SPxId enterId,
410 R& enterBound,
411 R& newUBbound,
412 R& newLBbound,
413 R& newCoPrhs,
414 StableSum<R>& objChange
415 )
416 {
417 typename SPxBasisBase<R>::Desc& ds = this->desc();
418
419 enterBound = 0;
420
421 if(enterId.isSPxRowId())
422 {
423 int idx = this->number(SPxRowId(enterId));
424 typename SPxBasisBase<R>::Desc::Status enterStat = ds.rowStatus(idx);
425
426 switch(enterStat)
427 {
428 case SPxBasisBase<R>::Desc::D_FREE :
429 assert(rep() == ROW);
430
431 if(thePvec->delta()[idx] * leaveMax < 0)
432 newCoPrhs = theLRbound[idx];
433 else
434 newCoPrhs = theURbound[idx];
435
436 newUBbound = R(infinity);
437 newLBbound = R(-infinity);
438 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED;
439 break;
440
441 case SPxBasisBase<R>::Desc::D_ON_UPPER :
442 assert(rep() == ROW);
443 newUBbound = 0;
444 newLBbound = R(-infinity);
445 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
446 newCoPrhs = theLRbound[idx];
447 break;
448
449 case SPxBasisBase<R>::Desc::D_ON_LOWER :
450 assert(rep() == ROW);
451 newUBbound = R(infinity);
452 newLBbound = 0;
453 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
454 newCoPrhs = theURbound[idx];
455 break;
456
457 case SPxBasisBase<R>::Desc::D_ON_BOTH :
458 assert(rep() == ROW);
459
460 if(leaveMax * thePvec->delta()[idx] < 0)
461 {
462 newUBbound = 0;
463 newLBbound = R(-infinity);
464 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
465 newCoPrhs = theLRbound[idx];
466 }
467 else
468 {
469 newUBbound = R(infinity);
470 newLBbound = 0;
471 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
472 newCoPrhs = theURbound[idx];
473 }
474
475 break;
476
477 case SPxBasisBase<R>::Desc::P_ON_UPPER :
478 assert(rep() == COLUMN);
479 ds.rowStatus(idx) = this->dualRowStatus(idx);
480
481 if(this->lhs(idx) > R(-infinity))
482 theURbound[idx] = theLRbound[idx];
483
484 newCoPrhs = theLRbound[idx]; // slack !!
485 newUBbound = -this->lhs(idx);
486 newLBbound = -this->rhs(idx);
487 enterBound = -this->rhs(idx);
488 objChange -= newCoPrhs * this->rhs(idx);
489 break;
490
491 case SPxBasisBase<R>::Desc::P_ON_LOWER :
492 assert(rep() == COLUMN);
493 ds.rowStatus(idx) = this->dualRowStatus(idx);
494
495 if(this->rhs(idx) < R(infinity))
496 theLRbound[idx] = theURbound[idx];
497
498 newCoPrhs = theURbound[idx]; // slack !!
499 newLBbound = -this->rhs(idx);
500 newUBbound = -this->lhs(idx);
501 enterBound = -this->lhs(idx);
502 objChange -= newCoPrhs * this->lhs(idx);
503 break;
504
505 case SPxBasisBase<R>::Desc::P_FREE :
506 assert(rep() == COLUMN);
507 #if 1
508 throw SPxInternalCodeException("XLEAVE04 This should never happen.");
509 #else
510 MSG_ERROR(std::cerr << "ELEAVE53 ERROR: not yet debugged!" << std::endl;)
511 ds.rowStatus(idx) = this->dualRowStatus(idx);
512 newCoPrhs = theURbound[idx]; // slack !!
513 newUBbound = R(infinity);
514 newLBbound = R(-infinity);
515 enterBound = 0;
516 #endif
517 break;
518
519 case SPxBasisBase<R>::Desc::P_FIXED :
520 assert(rep() == COLUMN);
521 MSG_ERROR(std::cerr << "ELEAVE54 "
522 << "ERROR! Tried to put a fixed row variable into the basis: "
523 << "idx=" << idx
524 << ", lhs=" << this->lhs(idx)
525 << ", rhs=" << this->rhs(idx) << std::endl;)
526 throw SPxInternalCodeException("XLEAVE05 This should never happen.");
527
528 default:
529 throw SPxInternalCodeException("XLEAVE06 This should never happen.");
530 }
531
532 MSG_DEBUG(std::cout << "DLEAVE55 SPxSolverBase<R>::getLeaveVals2(): row " << idx
533 << ": " << enterStat
534 << " -> " << ds.rowStatus(idx)
535 << " objChange: " << objChange
536 << std::endl;)
537 }
538
539 else
540 {
541 assert(enterId.isSPxColId());
542 int idx = this->number(SPxColId(enterId));
543 typename SPxBasisBase<R>::Desc::Status enterStat = ds.colStatus(idx);
544
545 switch(enterStat)
546 {
547 case SPxBasisBase<R>::Desc::D_ON_UPPER :
548 assert(rep() == ROW);
549 newUBbound = 0;
550 newLBbound = R(-infinity);
551 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
552 newCoPrhs = theLCbound[idx];
553 break;
554
555 case SPxBasisBase<R>::Desc::D_ON_LOWER :
556 assert(rep() == ROW);
557 newUBbound = R(infinity);
558 newLBbound = 0;
559 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
560 newCoPrhs = theUCbound[idx];
561 break;
562
563 case SPxBasisBase<R>::Desc::D_FREE :
564 assert(rep() == ROW);
565 newUBbound = R(infinity);
566 newLBbound = R(-infinity);
567 newCoPrhs = theLCbound[idx];
568 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED;
569 break;
570
571 case SPxBasisBase<R>::Desc::D_ON_BOTH :
572 assert(rep() == ROW);
573
574 if(leaveMax * theCoPvec->delta()[idx] < 0)
575 {
576 newUBbound = 0;
577 newLBbound = R(-infinity);
578 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
579 newCoPrhs = theLCbound[idx];
580 }
581 else
582 {
583 newUBbound = R(infinity);
584 newLBbound = 0;
585 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
586 newCoPrhs = theUCbound[idx];
587 }
588
589 break;
590
591 case SPxBasisBase<R>::Desc::P_ON_UPPER :
592 assert(rep() == COLUMN);
593 ds.colStatus(idx) = this->dualColStatus(idx);
594
595 if(SPxLPBase<R>::lower(idx) > R(-infinity))
596 theLCbound[idx] = theUCbound[idx];
597
598 newCoPrhs = theUCbound[idx];
599 newUBbound = SPxLPBase<R>::upper(idx);
600 newLBbound = SPxLPBase<R>::lower(idx);
601 enterBound = SPxLPBase<R>::upper(idx);
602 objChange -= newCoPrhs * enterBound;
603 break;
604
605 case SPxBasisBase<R>::Desc::P_ON_LOWER :
606 assert(rep() == COLUMN);
607 ds.colStatus(idx) = this->dualColStatus(idx);
608
609 if(SPxLPBase<R>::upper(idx) < R(infinity))
610 theUCbound[idx] = theLCbound[idx];
611
612 newCoPrhs = theLCbound[idx];
613 newUBbound = SPxLPBase<R>::upper(idx);
614 newLBbound = SPxLPBase<R>::lower(idx);
615 enterBound = SPxLPBase<R>::lower(idx);
616 objChange -= newCoPrhs * enterBound;
617 break;
618
619 case SPxBasisBase<R>::Desc::P_FREE :
620 assert(rep() == COLUMN);
621 ds.colStatus(idx) = this->dualColStatus(idx);
622
623 if(thePvec->delta()[idx] * leaveMax > 0)
624 newCoPrhs = theUCbound[idx];
625 else
626 newCoPrhs = theLCbound[idx];
627
628 newUBbound = SPxLPBase<R>::upper(idx);
629 newLBbound = SPxLPBase<R>::lower(idx);
630 enterBound = 0;
631 break;
632
633 case SPxBasisBase<R>::Desc::P_FIXED :
634 assert(rep() == COLUMN);
635 MSG_ERROR(std::cerr << "ELEAVE56 "
636 << "ERROR! Tried to put a fixed column variable into the basis. "
637 << "idx=" << idx
638 << ", lower=" << this->lower(idx)
639 << ", upper=" << this->upper(idx) << std::endl;)
640 throw SPxInternalCodeException("XLEAVE07 This should never happen.");
641
642 default:
643 throw SPxInternalCodeException("XLEAVE08 This should never happen.");
644 }
645
646 MSG_DEBUG(std::cout << "DLEAVE57 SPxSolverBase<R>::getLeaveVals2(): col " << idx
647 << ": " << enterStat
648 << " -> " << ds.colStatus(idx)
649 << " objChange: " << objChange
650 << std::endl;)
651 }
652
653 }
654
655 template <class R>
656 void SPxSolverBase<R>::rejectLeave(
657 int leaveNum,
658 SPxId leaveId,
659 typename SPxBasisBase<R>::Desc::Status leaveStat,
660 const SVectorBase<R>* //newVec
661 )
662 {
663 typename SPxBasisBase<R>::Desc& ds = this->desc();
664
665 if(leaveId.isSPxRowId())
666 {
667 MSG_DEBUG(std::cout << "DLEAVE58 rejectLeave() : row " << leaveNum
668 << ": " << ds.rowStatus(leaveNum)
669 << " -> " << leaveStat << std::endl;)
670
671 if(leaveStat == SPxBasisBase<R>::Desc::D_ON_BOTH)
672 {
673 if(ds.rowStatus(leaveNum) == SPxBasisBase<R>::Desc::P_ON_LOWER)
674 theLRbound[leaveNum] = theURbound[leaveNum];
675 else
676 theURbound[leaveNum] = theLRbound[leaveNum];
677 }
678
679 ds.rowStatus(leaveNum) = leaveStat;
680 }
681 else
682 {
683 MSG_DEBUG(std::cout << "DLEAVE59 rejectLeave() : col " << leaveNum
684 << ": " << ds.colStatus(leaveNum)
685 << " -> " << leaveStat << std::endl;)
686
687 if(leaveStat == SPxBasisBase<R>::Desc::D_ON_BOTH)
688 {
689 if(ds.colStatus(leaveNum) == SPxBasisBase<R>::Desc::P_ON_UPPER)
690 theLCbound[leaveNum] = theUCbound[leaveNum];
691 else
692 theUCbound[leaveNum] = theLCbound[leaveNum];
693 }
694
695 ds.colStatus(leaveNum) = leaveStat;
696 }
697 }
698
699
700 template <class R>
701 void SPxSolverBase<R>::computePrimalray4Row(R direction)
702 {
703 R sign = (direction > 0 ? 1.0 : -1.0);
704
705 primalRay.clear();
706 primalRay.setMax(coPvec().delta().size());
707
708 for(int i = 0; i < coPvec().delta().size(); ++i)
709 primalRay.add(coPvec().delta().index(i), sign * coPvec().delta().value(i));
710 }
711
712 template <class R>
713 void SPxSolverBase<R>::computeDualfarkas4Col(R direction)
714 {
715 R sign = (direction > 0 ? -1.0 : 1.0);
716
717 dualFarkas.clear();
718 dualFarkas.setMax(coPvec().delta().size());
719
720 for(int i = 0; i < coPvec().delta().size(); ++i)
721 dualFarkas.add(coPvec().delta().index(i), sign * coPvec().delta().value(i));
722 }
723
724 template <class R>
725 bool SPxSolverBase<R>::leave(int leaveIdx, bool polish)
726 {
727 assert(leaveIdx < dim() && leaveIdx >= 0);
728 assert(type() == LEAVE);
729 assert(initialized);
730
731 bool instable = instableLeave;
732 assert(!instable || instableLeaveNum >= 0);
733
734 /*
735 Before performing the actual basis update, we must determine, how this
736 is to be accomplished.
737 When using steepest edge pricing this solve is already performed by the pricer
738 */
739 if(theCoPvec->delta().isSetup() && theCoPvec->delta().size() == 0)
740 {
741 this->coSolve(theCoPvec->delta(), unitVecs[leaveIdx]);
742 }
743
744 #ifdef ENABLE_ADDITIONAL_CHECKS
745 else
746 {
747 SSVectorBase<R> tmp(dim(), epsilon());
748 tmp.clear();
749 this->coSolve(tmp, unitVecs[leaveIdx]);
750 tmp -= theCoPvec->delta();
751
752 if(tmp.length() > leavetol())
753 {
754 // This happens very frequently and does usually not hurt, so print
755 // these warnings only with verbose level INFO2 and higher.
756 MSG_INFO2((*this->spxout), (*this->spxout) << "WLEAVE60 iteration=" << basis().iteration()
757 << ": coPvec.delta error = " << tmp.length()
758 << std::endl;)
759 }
760 }
761
762 #endif // ENABLE_ADDITIONAL_CHECKS
763
764 setupPupdate();
765
766 assert(thePvec->isConsistent());
767 assert(theCoPvec->isConsistent());
768
769 typename SPxBasisBase<R>::Desc::Status leaveStat; // status of leaving var
770 SPxId leaveId; // id of leaving var
771 SPxId none; // invalid id used if leave fails
772 R leaveMax; // maximium lambda of leaving var
773 R leavebound; // current fVec value of leaving var
774 int leaveNum; // number of leaveId in bounds
775 StableSum<R> objChange; // amount of change in the objective function
776
777 getLeaveVals(leaveIdx, leaveStat, leaveId, leaveMax, leavebound, leaveNum, objChange);
778
779 if(!polish && m_numCycle > m_maxCycle)
780 {
781 if(leaveMax > 0)
782 perturbMaxLeave();
783 else
784 perturbMinLeave();
785
786 //@ m_numCycle /= 2;
787 // perturbation invalidates the currently stored nonbasic value
788 forceRecompNonbasicValue();
789 }
790
791 //@ testBounds();
792
793 R enterVal = leaveMax;
794 boundflips = 0;
795 R oldShift = theShift;
796 SPxId enterId = theratiotester->selectEnter(enterVal, leaveIdx, polish);
797
798 if(NE(theShift, oldShift))
799 {
800 MSG_DEBUG(std::cout << "DLEAVE71 trigger recomputation of nonbasic value due to shifts in ratiotest"
801 << std::endl;)
802 forceRecompNonbasicValue();
803 }
804
805 assert(!enterId.isValid() || !isBasic(enterId));
806
807 instableLeaveNum = -1;
808 instableLeave = false;
809
810 /*
811 No variable could be selected to enter the basis and even the leaving
812 variable is unbounded.
813 */
814 if(!enterId.isValid())
815 {
816 /* the following line originally was below in "rejecting leave" case;
817 we need it in the unbounded/infeasible case, too, to have the
818 correct basis size */
819 rejectLeave(leaveNum, leaveId, leaveStat);
820 this->change(-1, none, 0);
821 objChange = R(0.0); // the nonbasicValue is not supposed to be updated in this case
822
823 if(polish)
824 return false;
825
826 if(NE(enterVal, leaveMax))
827 {
828 MSG_DEBUG(std::cout << "DLEAVE61 rejecting leave A (leaveIdx=" << leaveIdx
829 << ", theCoTest=" << theCoTest[leaveIdx] << ")"
830 << std::endl;)
831
832 /* In the LEAVE algorithm, when for a selected leaving variable we find only
833 an instable entering variable, then the basis change is not conducted.
834 Instead, we save the leaving variable's index in instableLeaveNum and scale
835 theCoTest[leaveIdx] down by some factor, hoping to find a different leaving
836 variable with a stable entering variable.
837 If this fails, however, and no more leaving variable is found, we have to
838 perform the instable basis change using instableLeaveNum. In this (and only
839 in this) case, the flag instableLeave is set to true.
840
841 enterVal != leaveMax is the case that selectEnter has found only an instable entering
842 variable. We store this leaving variable for later -- if we are not already in the
843 instable case: then we continue and conclude unboundedness/infeasibility */
844 if(!instable)
845 {
846 instableLeaveNum = leaveIdx;
847
848 // Note: These changes do not survive a refactorization
849 instableLeaveVal = theCoTest[leaveIdx];
850 theCoTest[leaveIdx] = instableLeaveVal / 10.0;
851
852 return true;
853 }
854 }
855
856 if(this->lastUpdate() > 1)
857 {
858 MSG_INFO3((*this->spxout), (*this->spxout) << "ILEAVE01 factorization triggered in "
859 << "leave() for feasibility test" << std::endl;)
860
861 try
862 {
863 factorize();
864 }
865 catch(const SPxStatusException& E)
866 {
867 // don't exit immediately but handle the singularity correctly
868 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::SINGULAR);
869 MSG_INFO3((*this->spxout), (*this->spxout) << "Caught exception in factorization: " << E.what() <<
870 std::endl;)
871 }
872
873 /* after a factorization, the leaving column/row might not be infeasible or suboptimal anymore, hence we do
874 * not try to call leave(leaveIdx), but rather return to the main solving loop and call the pricer again
875 */
876 return true;
877 }
878
879 /* do not exit with status infeasible or unbounded if there is only a very small violation */
880 if(!recomputedVectors && spxAbs(enterVal) < leavetol())
881 {
882 MSG_INFO3((*this->spxout), (*this->spxout) << "ILEAVE11 clean up step to reduce numerical errors" <<
883 std::endl;)
884
885 computeFrhs();
886 SPxBasisBase<R>::solve(*theFvec, *theFrhs);
887 computeFtest();
888
889 /* only do this once per solve */
890 recomputedVectors = true;
891
892 return true;
893 }
894
895 MSG_INFO3((*this->spxout), (*this->spxout) << "ILEAVE02 unboundedness/infeasibility found "
896 << "in leave()" << std::endl;)
897
898 if(rep() != COLUMN)
899 {
900 computePrimalray4Row(enterVal);
901 setBasisStatus(SPxBasisBase<R>::UNBOUNDED);
902 }
903 else
904 {
905 computeDualfarkas4Col(enterVal);
906 setBasisStatus(SPxBasisBase<R>::INFEASIBLE);
907 }
908
909 return false;
910 }
911 else
912 {
913 /*
914 If an entering variable has been found, a regular basis update is to
915 be performed.
916 */
917 if(enterId != this->baseId((leaveIdx)))
918 {
919 const SVectorBase<R>& newVector = *enterVector(enterId);
920
921 // update feasibility vectors
922 if(solveVector2 != NULL && solveVector3 != NULL)
923 {
924 assert(solveVector2->isConsistent());
925 assert(solveVector2rhs->isSetup());
926 assert(solveVector3->isConsistent());
927 assert(solveVector3rhs->isSetup());
928 assert(boundflips > 0);
929 SPxBasisBase<R>::solve4update(theFvec->delta(),
930 *solveVector2,
931 *solveVector3,
932 newVector,
933 *solveVector2rhs,
934 *solveVector3rhs);
935
936 // perform update of basic solution
937 primVec -= (*solveVector3);
938 MSG_DEBUG(std::cout << "ILBFRT02 breakpoints passed / bounds flipped = " << boundflips << std::endl;
939 )
940 totalboundflips += boundflips;
941 }
942 else if(solveVector2 != NULL)
943 {
944 assert(solveVector2->isConsistent());
945 assert(solveVector2rhs->isSetup());
946
947 SPxBasisBase<R>::solve4update(theFvec->delta(),
948 *solveVector2,
949 newVector,
950 *solveVector2rhs);
951 }
952 else if(solveVector3 != NULL)
953 {
954 assert(solveVector3->isConsistent());
955 assert(solveVector3rhs->isSetup());
956 assert(boundflips > 0);
957 SPxBasisBase<R>::solve4update(theFvec->delta(),
958 *solveVector3,
959 newVector,
960 *solveVector3rhs);
961
962 // perform update of basic solution
963 primVec -= (*solveVector3);
964 MSG_DEBUG(std::cout << "ILBFRT02 breakpoints passed / bounds flipped = " << boundflips << std::endl;
965 )
966 totalboundflips += boundflips;
967 }
968 else
969 SPxBasisBase<R>::solve4update(theFvec->delta(), newVector);
970
971 #ifdef ENABLE_ADDITIONAL_CHECKS
972 {
973 SSVectorBase<R> tmp(dim(), epsilon());
974 SPxBasisBase<R>::solve(tmp, newVector);
975 tmp -= fVec().delta();
976
977 if(tmp.length() > entertol())
978 {
979 // This happens very frequently and does usually not hurt, so print
980 // these warnings only with verbose level INFO2 and higher.
981 MSG_INFO2((*this->spxout), (*this->spxout) << "WLEAVE62\t(" << tmp.length() << ")\n";)
982 }
983 }
984 #endif // ENABLE_ADDITIONAL_CHECKS
985
986
987 if(spxAbs(theFvec->delta()[leaveIdx]) < reject_leave_tol)
988 {
989 if(instable)
990 {
991 /* We are in the case that for all leaving variables only instable entering
992 variables were found: Thus, above we already accepted such an instable
993 entering variable. Now even this seems to be impossible, thus we conclude
994 unboundedness/infeasibility. */
995 MSG_INFO3((*this->spxout), (*this->spxout) << "ILEAVE03 unboundedness/infeasibility found "
996 << "in leave()" << std::endl;)
997
998 rejectLeave(leaveNum, leaveId, leaveStat);
999 this->change(-1, none, 0);
1000 objChange = R(0.0); // the nonbasicValue is not supposed to be updated in this case
1001
1002 /**@todo if shift() is not zero we must not conclude unboundedness */
1003 if(rep() == ROW)
1004 {
1005 computePrimalray4Row(enterVal);
1006 setBasisStatus(SPxBasisBase<R>::UNBOUNDED);
1007 }
1008 else
1009 {
1010 computeDualfarkas4Col(enterVal);
1011 setBasisStatus(SPxBasisBase<R>::INFEASIBLE);
1012 }
1013
1014 return false;
1015 }
1016 else
1017 {
1018 theFvec->delta().clear();
1019 rejectLeave(leaveNum, leaveId, leaveStat, &newVector);
1020 this->change(-1, none, 0);
1021 objChange = R(0.0); // the nonbasicValue is not supposed to be updated in this case
1022
1023 MSG_DEBUG(std::cout << "DLEAVE63 rejecting leave B (leaveIdx=" << leaveIdx
1024 << ", theCoTest=" << theCoTest[leaveIdx]
1025 << ")" << std::endl;)
1026
1027 // Note: These changes do not survive a refactorization
1028 theCoTest[leaveIdx] *= 0.01;
1029
1030 return true;
1031 }
1032 }
1033
1034 // process leaving variable
1035 if(leavebound > epsilon() || leavebound < -epsilon())
1036 theFrhs->multAdd(-leavebound, this->baseVec(leaveIdx));
1037
1038 // process entering variable
1039 R enterBound;
1040 R newUBbound;
1041 R newLBbound;
1042 R newCoPrhs;
1043
1044 try
1045 {
1046 getLeaveVals2(leaveMax, enterId, enterBound, newUBbound, newLBbound, newCoPrhs, objChange);
1047 }
1048 catch(const SPxException& F)
1049 {
1050 rejectLeave(leaveNum, leaveId, leaveStat);
1051 this->change(-1, none, 0);
1052 objChange = R(0.0); // the nonbasicValue is not supposed to be updated in this case
1053 throw F;
1054 }
1055
1056 theUBbound[leaveIdx] = newUBbound;
1057 theLBbound[leaveIdx] = newLBbound;
1058 (*theCoPrhs)[leaveIdx] = newCoPrhs;
1059
1060 if(enterBound > epsilon() || enterBound < -epsilon())
1061 theFrhs->multAdd(enterBound, newVector);
1062
1063 // update pricing vectors
1064 theCoPvec->value() = enterVal;
1065 thePvec->value() = enterVal;
1066
1067 if(enterVal > epsilon() || enterVal < -epsilon())
1068 doPupdate();
1069
1070 // update feasibility vector
1071 theFvec->value() = -((*theFvec)[leaveIdx] - leavebound)
1072 / theFvec->delta()[leaveIdx];
1073 theFvec->update();
1074 (*theFvec)[leaveIdx] = enterBound - theFvec->value();
1075 updateFtest();
1076
1077 // update objective funtion value
1078 updateNonbasicValue(objChange);
1079
1080 // change basis matrix
1081 this->change(leaveIdx, enterId, &newVector, &(theFvec->delta()));
1082 }
1083
1084 /*
1085 No entering vector has been selected from the basis. However, if the
1086 shift amount for |coPvec| is bounded, we are in the case, that the
1087 entering variable is moved from one bound to its other, before any of
1088 the basis feasibility variables reaches their bound. This may only
1089 happen in primal/columnwise case with upper and lower bounds on
1090 variables.
1091 */
1092 else
1093 {
1094 // @todo update obj function value here!!!
1095 assert(rep() == ROW);
1096 typename SPxBasisBase<R>::Desc& ds = this->desc();
1097
1098 this->change(leaveIdx, none, 0);
1099
1100 if(leaveStat == SPxBasisBase<R>::Desc::P_ON_UPPER)
1101 {
1102 if(leaveId.isSPxRowId())
1103 {
1104 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1105 (*theCoPrhs)[leaveIdx] = theLRbound[leaveNum];
1106 }
1107 else
1108 {
1109 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1110 (*theCoPrhs)[leaveIdx] = theLCbound[leaveNum];
1111 }
1112
1113 theUBbound[leaveIdx] = 0;
1114 theLBbound[leaveIdx] = R(-infinity);
1115 }
1116 else
1117 {
1118 assert(leaveStat == SPxBasisBase<R>::Desc::P_ON_LOWER);
1119
1120 if(leaveId.isSPxRowId())
1121 {
1122 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1123 (*theCoPrhs)[leaveIdx] = theURbound[leaveNum];
1124 }
1125 else
1126 {
1127 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1128 (*theCoPrhs)[leaveIdx] = theUCbound[leaveNum];
1129 }
1130
1131 theUBbound[leaveIdx] = R(infinity);
1132 theLBbound[leaveIdx] = 0;
1133 }
1134
1135 // update copricing vector
1136 theCoPvec->value() = enterVal;
1137 thePvec->value() = enterVal;
1138
1139 if(enterVal > epsilon() || enterVal < -epsilon())
1140 doPupdate();
1141
1142 // update feasibility vectors
1143 theFvec->value() = 0;
1144 assert(theCoTest[leaveIdx] < 0.0);
1145 m_pricingViol += theCoTest[leaveIdx];
1146 theCoTest[leaveIdx] *= -1;
1147 }
1148
1149 if((leaveMax > entertol() && enterVal <= entertol()) || (leaveMax < -entertol()
1150 && enterVal >= -entertol()))
1151 {
1152 if((theUBbound[leaveIdx] < R(infinity) || theLBbound[leaveIdx] > R(-infinity))
1153 && leaveStat != SPxBasisBase<R>::Desc::P_FREE
1154 && leaveStat != SPxBasisBase<R>::Desc::D_FREE)
1155 {
1156 m_numCycle++;
1157 leaveCycles++;
1158 }
1159 }
1160 else
1161 m_numCycle /= 2;
1162
1163 #ifdef ENABLE_ADDITIONAL_CHECKS
1164 {
1165 VectorBase<R> tmp = fVec();
1166 this->multBaseWith(tmp);
1167 tmp -= fRhs();
1168
1169 if(tmp.length() > entertol())
1170 {
1171 // This happens very frequently and does usually not hurt, so print
1172 // these warnings only with verbose level INFO2 and higher.
1173 MSG_INFO2((*this->spxout), (*this->spxout) << "WLEAVE64\t" << basis().iteration()
1174 << ": fVec error = " << tmp.length() << std::endl;)
1175 SPxBasisBase<R>::solve(tmp, fRhs());
1176 tmp -= fVec();
1177 MSG_INFO2((*this->spxout), (*this->spxout) << "WLEAVE65\t(" << tmp.length() << ")\n";)
1178 }
1179 }
1180 #endif // ENABLE_ADDITIONAL_CHECKS
1181
1182 return true;
1183 }
1184 }
1185 } // namespace soplex
1186