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 /* \SubSection{Updating the Basis for Entering Variables}
17 */
18 #include <assert.h>
19
20 #include "soplex/spxdefines.h"
21 #include "soplex/spxratiotester.h"
22 #include "soplex/spxpricer.h"
23 #include "soplex/spxout.h"
24 #include "soplex/exceptions.h"
25
26 namespace soplex
27 {
28
29 /*
30 In the entering simplex algorithms (i.e. iteratively a vector is selected to
31 \em enter the simplex basis as in the dual rowwise and primal columnwise case)
32 let $A$ denote the current basis, $x$ and entering vector and $f$ the
33 feasibility vector. For a feasible basis $l \le f \le u$ holds.
34 For the rowwise case $f$ is obtained by solving $f^T = c^T A^{-1}$,
35 wherease in columnwisecase $f = A^{-1} b$.
36
37 Let us further consider the rowwise case. Exchanging $x$ with the $i$-th
38 vector of $A$ yields
39
40 \begin{equation}\label{update.eq}
41 A^{(i)} = E_i A \hbox{, with } E_i = I + e_i (x^T A^{-1} - e_i^T).
42 \end{equation}
43
44 With $E_i^{-1} = I + e_i \frac{e_i^T - \delta^T}{\delta}$,
45 $\delta^T = x^T A^{-1}$ one gets the new feasibility vector
46
47 \begin{eqnarray*}
48 (f^{(i)})^T
49 &=& c^T (A^{(i)})^{-1} \\
50 &=& c^T A^{-1} + c^T A^{-1} e_i \frac{e_i^T - \delta^T}{\delta_i} \\
51 &=& f^T + \frac{f_i}{\delta_i} e_i^T - \frac{f_i}{\delta_i} \delta^T. \\
52 \end{eqnarray*}
53
54 The selection of the leaving vector $i^*$ for the basis must ensure, that for
55 all $j \ne i^*$ $f^{(i^*)}_j$ remains within its bounds $l_j$ and $u_j$.
56 */
57
58
59 /*
60 Testing all values of |pVec| against its bounds. If $i$, say, is violated
61 the violation is saved as negative value in |theTest[i]|.
62 */
63
64 template <class R>
65 R SPxSolverBase<R>::test(int i, typename SPxBasisBase<R>::Desc::Status stat) const
66 {
67 assert(type() == ENTER);
68 assert(!isBasic(stat));
69
70 R x;
71
72 switch(stat)
73 {
74 case SPxBasisBase<R>::Desc::D_FREE:
75 case SPxBasisBase<R>::Desc::D_ON_BOTH:
76 assert(rep() == ROW);
77 x = (*thePvec)[i] - this->lhs(i);
78
79 if(x < 0)
80 return x;
81
82 // no break: next is else case
83 //lint -fallthrough
84 case SPxBasisBase<R>::Desc::D_ON_LOWER:
85 assert(rep() == ROW);
86 return this->rhs(i) - (*thePvec)[i];
87
88 case SPxBasisBase<R>::Desc::D_ON_UPPER:
89 assert(rep() == ROW);
90 return (*thePvec)[i] - this->lhs(i);
91
92 case SPxBasisBase<R>::Desc::P_ON_UPPER:
93 assert(rep() == COLUMN);
94 return this->maxObj(i) - (*thePvec)[i];
95
96 case SPxBasisBase<R>::Desc::P_ON_LOWER:
97 assert(rep() == COLUMN);
98 return (*thePvec)[i] - this->maxObj(i);
99
100 case SPxBasisBase<R>::Desc::P_FREE :
101 x = this->maxObj(i) - (*thePvec)[i];
102 return (x < 0) ? x : -x;
103
104 default:
105 return 0;
106 }
107 }
108
109 template <class R>
110 void SPxSolverBase<R>::computeTest()
111 {
112
113 const typename SPxBasisBase<R>::Desc& ds = this->desc();
114 R pricingTol = leavetol();
115 m_pricingViolCoUpToDate = true;
116 m_pricingViolCo = 0;
117
118 infeasibilitiesCo.clear();
119 int sparsitythreshold = (int)(sparsePricingFactor * coDim());
120
121 for(int i = 0; i < coDim(); ++i)
122 {
123 typename SPxBasisBase<R>::Desc::Status stat = ds.status(i);
124
125 if(isBasic(stat))
126 {
127 theTest[i] = 0.0;
128
129 if(remainingRoundsEnterCo == 0)
130 isInfeasibleCo[i] = SPxPricer<R>::NOT_VIOLATED;
131 }
132 else
133 {
134 assert(!isBasic(stat));
135 theTest[i] = test(i, stat);
136
137 if(remainingRoundsEnterCo == 0)
138 {
139 if(theTest[i] < -pricingTol)
140 {
141 assert(infeasibilitiesCo.size() < infeasibilitiesCo.max());
142 m_pricingViolCo -= theTest[i];
143 infeasibilitiesCo.addIdx(i);
144 isInfeasibleCo[i] = SPxPricer<R>::VIOLATED;
145 ++m_numViol;
146 }
147 else
148 isInfeasibleCo[i] = SPxPricer<R>::NOT_VIOLATED;
149
150 if(infeasibilitiesCo.size() > sparsitythreshold)
151 {
152 MSG_INFO2((*this->spxout), (*this->spxout) << " --- using dense pricing"
153 << std::endl;)
154 remainingRoundsEnterCo = DENSEROUNDS;
155 sparsePricingEnterCo = false;
156 infeasibilitiesCo.clear();
157 }
158 }
159 else if(theTest[i] < -pricingTol)
160 {
161 m_pricingViolCo -= theTest[i];
162 ++m_numViol;
163 }
164 }
165 }
166
167 if(infeasibilitiesCo.size() == 0 && !sparsePricingEnterCo)
168 --remainingRoundsEnterCo;
169 else if(infeasibilitiesCo.size() <= sparsitythreshold && !sparsePricingEnterCo)
170 {
171 MSG_INFO2((*this->spxout),
172 std::streamsize prec = spxout->precision();
173
174 if(hyperPricingEnter)
175 (*this->spxout) << " --- using hypersparse pricing, ";
176 else
177 (*this->spxout) << " --- using sparse pricing, ";
178 (*this->spxout) << "sparsity: "
179 << std::setw(6) << std::fixed << std::setprecision(4)
180 << (R) infeasibilitiesCo.size() / coDim()
181 << std::scientific << std::setprecision(int(prec))
182 << std::endl;
183 )
184 sparsePricingEnterCo = true;
185 }
186 }
187
188 template <class R>
189 R SPxSolverBase<R>::computePvec(int i)
190 {
191
192 return (*thePvec)[i] = vector(i) * (*theCoPvec);
193 }
194
195 template <class R>
196 R SPxSolverBase<R>::computeTest(int i)
197 {
198 typename SPxBasisBase<R>::Desc::Status stat = this->desc().status(i);
199
200 if(isBasic(stat))
201 return theTest[i] = 0;
202 else
203 return theTest[i] = test(i, stat);
204 }
205
206 /*
207 Testing all values of #coPvec# against its bounds. If $i$, say, is violated
208 the violation is saved as negative value in |theCoTest[i]|.
209 */
210 template <class R>
211 R SPxSolverBase<R>::coTest(int i, typename SPxBasisBase<R>::Desc::Status stat) const
212 {
213 assert(type() == ENTER);
214 assert(!isBasic(stat));
215
216 R x;
217
218 switch(stat)
219 {
220 case SPxBasisBase<R>::Desc::D_FREE:
221 case SPxBasisBase<R>::Desc::D_ON_BOTH :
222 assert(rep() == ROW);
223 x = (*theCoPvec)[i] - SPxLPBase<R>::lower(i);
224
225 if(x < 0)
226 return x;
227
228 // no break: next is else case
229 //lint -fallthrough
230 case SPxBasisBase<R>::Desc::D_ON_LOWER:
231 assert(rep() == ROW);
232 return SPxLPBase<R>::upper(i) - (*theCoPvec)[i];
233
234 case SPxBasisBase<R>::Desc::D_ON_UPPER:
235 assert(rep() == ROW);
236 return (*theCoPvec)[i] - SPxLPBase<R>::lower(i);
237
238 case SPxBasisBase<R>::Desc::P_ON_UPPER:
239 assert(rep() == COLUMN);
240 return (*theCoPvec)[i] - this->maxRowObj(i); // slacks !
241
242 case SPxBasisBase<R>::Desc::P_ON_LOWER:
243 assert(rep() == COLUMN);
244 return this->maxRowObj(i) - (*theCoPvec)[i]; // slacks !
245
246 default:
247 return 0;
248 }
249 }
250
251 template <class R>
252 void SPxSolverBase<R>::computeCoTest()
253 {
254 int i;
255 R pricingTol = leavetol();
256 m_pricingViolUpToDate = true;
257 m_pricingViol = 0;
258 m_numViol = 0;
259 infeasibilities.clear();
260 int sparsitythreshold = (int)(sparsePricingFactor * dim());
261 const typename SPxBasisBase<R>::Desc& ds = this->desc();
262
263 for(i = dim() - 1; i >= 0; --i)
264 {
265 typename SPxBasisBase<R>::Desc::Status stat = ds.coStatus(i);
266
267 if(isBasic(stat))
268 {
269 theCoTest[i] = 0;
270
271 if(remainingRoundsEnter == 0)
272 isInfeasible[i] = SPxPricer<R>::NOT_VIOLATED;
273 }
274 else
275 {
276 theCoTest[i] = coTest(i, stat);
277
278 if(remainingRoundsEnter == 0)
279 {
280 if(theCoTest[i] < -pricingTol)
281 {
282 assert(infeasibilities.size() < infeasibilities.max());
283 m_pricingViol -= theCoTest[i];
284 infeasibilities.addIdx(i);
285 isInfeasible[i] = SPxPricer<R>::VIOLATED;
286 ++m_numViol;
287 }
288 else
289 isInfeasible[i] = SPxPricer<R>::NOT_VIOLATED;
290
291 if(infeasibilities.size() > sparsitythreshold)
292 {
293 MSG_INFO2((*this->spxout), (*this->spxout) << " --- using dense pricing"
294 << std::endl;)
295 remainingRoundsEnter = DENSEROUNDS;
296 sparsePricingEnter = false;
297 infeasibilities.clear();
298 }
299 }
300 else if(theCoTest[i] < -pricingTol)
301 {
302 m_pricingViol -= theCoTest[i];
303 ++m_numViol;
304 }
305 }
306 }
307
308 if(infeasibilities.size() == 0 && !sparsePricingEnter)
309 --remainingRoundsEnter;
310 else if(infeasibilities.size() <= sparsitythreshold && !sparsePricingEnter)
311 {
312 MSG_INFO2((*this->spxout),
313 std::streamsize prec = spxout->precision();
314
315 if(hyperPricingEnter)
316 (*this->spxout) << " --- using hypersparse pricing, ";
317 else
318 (*this->spxout) << " --- using sparse pricing, ";
319 (*this->spxout) << "sparsity: "
320 << std::setw(6) << std::fixed << std::setprecision(4)
321 << (R) infeasibilities.size() / dim()
322 << std::scientific << std::setprecision(int(prec))
323 << std::endl;
324 )
325 sparsePricingEnter = true;
326 }
327 }
328
329 /*
330 The following methods require propersy initialized vectors |fVec| and
331 #coPvec#.
332 */
333 template <class R>
334 void SPxSolverBase<R>::updateTest()
335 {
336 thePvec->delta().setup();
337
338 const IdxSet& idx = thePvec->idx();
339 const typename SPxBasisBase<R>::Desc& ds = this->desc();
340 R pricingTol = leavetol();
341
342 int i;
343 updateViolsCo.clear();
344
345 for(i = idx.size() - 1; i >= 0; --i)
346 {
347 int j = idx.index(i);
348 typename SPxBasisBase<R>::Desc::Status stat = ds.status(j);
349
350 if(!isBasic(stat))
351 {
352 if(m_pricingViolCoUpToDate && theTest[j] < -pricingTol)
353 m_pricingViolCo += theTest[j];
354
355 theTest[j] = test(j, stat);
356
357 if(sparsePricingEnterCo)
358 {
359 if(theTest[j] < -pricingTol)
360 {
361 assert(remainingRoundsEnterCo == 0);
362 m_pricingViolCo -= theTest[j];
363
364 if(isInfeasibleCo[j] == SPxPricer<R>::NOT_VIOLATED)
365 {
366 infeasibilitiesCo.addIdx(j);
367 isInfeasibleCo[j] = SPxPricer<R>::VIOLATED;
368 }
369
370 if(hyperPricingEnter)
371 updateViolsCo.addIdx(j);
372 }
373 else
374 {
375 isInfeasibleCo[j] = SPxPricer<R>::NOT_VIOLATED;
376 }
377 }
378 else if(theTest[j] < -pricingTol)
379 m_pricingViolCo -= theTest[j];
380 }
381 else
382 {
383 isInfeasibleCo[j] = SPxPricer<R>::NOT_VIOLATED;
384 theTest[j] = 0;
385 }
386 }
387 }
388
389 template <class R>
390 void SPxSolverBase<R>::updateCoTest()
391 {
392 theCoPvec->delta().setup();
393
394 const IdxSet& idx = theCoPvec->idx();
395 const typename SPxBasisBase<R>::Desc& ds = this->desc();
396 R pricingTol = leavetol();
397
398 int i;
399 updateViols.clear();
400
401 for(i = idx.size() - 1; i >= 0; --i)
402 {
403 int j = idx.index(i);
404 typename SPxBasisBase<R>::Desc::Status stat = ds.coStatus(j);
405
406 if(!isBasic(stat))
407 {
408 if(m_pricingViolUpToDate && theCoTest[j] < -pricingTol)
409 m_pricingViol += theCoTest[j];
410
411 theCoTest[j] = coTest(j, stat);
412
413 if(sparsePricingEnter)
414 {
415 if(theCoTest[j] < -pricingTol)
416 {
417 assert(remainingRoundsEnter == 0);
418 m_pricingViol -= theCoTest[j];
419
420 if(isInfeasible[j] == SPxPricer<R>::NOT_VIOLATED)
421 {
422 // if( !hyperPricingEnter )
423 infeasibilities.addIdx(j);
424 isInfeasible[j] = SPxPricer<R>::VIOLATED;
425 }
426
427 if(hyperPricingEnter)
428 updateViols.addIdx(j);
429 }
430 else
431 {
432 // @todo do we need to remove index j from infeasibilitiesCo?
433 isInfeasible[j] = SPxPricer<R>::NOT_VIOLATED;
434 }
435 }
436 else if(theCoTest[j] < -pricingTol)
437 m_pricingViol -= theCoTest[j];
438 }
439 else
440 {
441 isInfeasible[j] = SPxPricer<R>::NOT_VIOLATED;
442 theCoTest[j] = 0;
443 }
444 }
445 }
446
447
448
449 /* \Section{Compute statistics on entering variable}
450 Here is a list of variables relevant when including |Id| to the basis.
451 They are computed by |computeEnterStats()|.
452 */
453 template <class R>
454 void SPxSolverBase<R>::getEnterVals
455 (
456 SPxId enterId,
457 R& enterTest,
458 R& enterUB,
459 R& enterLB,
460 R& enterVal,
461 R& enterMax,
462 R& enterPric,
463 typename SPxBasisBase<R>::Desc::Status& enterStat,
464 R& enterRO,
465 StableSum<R>& objChange
466 )
467 {
468 int enterIdx;
469 typename SPxBasisBase<R>::Desc& ds = this->desc();
470
471 if(enterId.isSPxColId())
472 {
473 enterIdx = this->number(SPxColId(enterId));
474 enterStat = ds.colStatus(enterIdx);
475 assert(!isBasic(enterStat));
476
477 /* For an #Id# to enter the basis we better recompute the Test value.
478 */
479 if(rep() == COLUMN)
480 {
481 computePvec(enterIdx);
482 enterTest = computeTest(enterIdx);
483 theTest[enterIdx] = 0;
484 }
485 else
486 {
487 enterTest = coTest()[enterIdx];
488 theCoTest[enterIdx] = 0;
489 }
490
491 switch(enterStat)
492 {
493 // primal/columnwise cases:
494 case SPxBasisBase<R>::Desc::P_ON_UPPER :
495 assert(rep() == COLUMN);
496 enterUB = theUCbound[enterIdx];
497 enterLB = theLCbound[enterIdx];
498 enterVal = enterUB;
499 enterMax = enterLB - enterUB;
500 enterPric = (*thePvec)[enterIdx];
501 enterRO = this->maxObj(enterIdx);
502 objChange -= enterVal * enterRO;
503
504 if(enterLB <= R(-infinity))
505 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_LOWER;
506 else if(EQ(enterLB, enterUB))
507 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE;
508 else
509 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH;
510
511 break;
512
513 case SPxBasisBase<R>::Desc::P_ON_LOWER :
514 assert(rep() == COLUMN);
515 enterUB = theUCbound[enterIdx];
516 enterLB = theLCbound[enterIdx];
517 enterVal = enterLB;
518 enterMax = enterUB - enterLB;
519 enterPric = (*thePvec)[enterIdx];
520 enterRO = this->maxObj(enterIdx);
521 objChange -= enterVal * enterRO;
522
523 if(enterUB >= R(infinity))
524 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_UPPER;
525 else if(EQ(enterLB, enterUB))
526 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE;
527 else
528 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH;
529
530 break;
531
532 case SPxBasisBase<R>::Desc::P_FREE :
533 assert(rep() == COLUMN);
534 enterUB = theUCbound[enterIdx];
535 enterLB = theLCbound[enterIdx];
536 enterVal = 0;
537 enterPric = (*thePvec)[enterIdx];
538 enterRO = this->maxObj(enterIdx);
539 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_UNDEFINED;
540 enterMax = (enterRO - enterPric > 0) ? R(infinity) : R(-infinity);
541 break;
542
543 // dual/rowwise cases:
544 case SPxBasisBase<R>::Desc::D_ON_UPPER :
545 assert(rep() == ROW);
546 assert(theUCbound[enterIdx] < R(infinity));
547 enterUB = theUCbound[enterIdx];
548 enterLB = R(-infinity);
549 enterMax = R(-infinity);
550 enterVal = enterUB;
551 enterPric = (*theCoPvec)[enterIdx];
552 enterRO = SPxLPBase<R>::lower(enterIdx);
553 objChange -= enterRO * enterVal;
554 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
555 break;
556
557 case SPxBasisBase<R>::Desc::D_ON_LOWER :
558 assert(rep() == ROW);
559 assert(theLCbound[enterIdx] > R(-infinity));
560 enterLB = theLCbound[enterIdx];
561 enterUB = R(infinity);
562 enterMax = R(infinity);
563 enterVal = enterLB;
564 enterPric = (*theCoPvec)[enterIdx];
565 enterRO = SPxLPBase<R>::upper(enterIdx);
566 objChange -= enterRO * enterVal;
567 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
568 break;
569
570 case SPxBasisBase<R>::Desc::D_FREE:
571 assert(rep() == ROW);
572 assert(SPxLPBase<R>::lower(enterIdx) == SPxLPBase<R>::upper(enterIdx));
573 enterUB = R(infinity);
574 enterLB = R(-infinity);
575 enterVal = 0;
576 enterRO = SPxLPBase<R>::upper(enterIdx);
577 enterPric = (*theCoPvec)[enterIdx];
578
579 if(enterPric > enterRO)
580 enterMax = R(infinity);
581 else
582 enterMax = R(-infinity);
583
584 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_FIXED;
585 break;
586
587 case SPxBasisBase<R>::Desc::D_ON_BOTH :
588 assert(rep() == ROW);
589 enterPric = (*theCoPvec)[enterIdx];
590
591 if(enterPric > SPxLPBase<R>::upper(enterIdx))
592 {
593 enterLB = theLCbound[enterIdx];
594 enterUB = R(infinity);
595 enterMax = R(infinity);
596 enterVal = enterLB;
597 enterRO = SPxLPBase<R>::upper(enterIdx);
598 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
599 }
600 else
601 {
602 enterUB = theUCbound[enterIdx];
603 enterVal = enterUB;
604 enterRO = SPxLPBase<R>::lower(enterIdx);
605 enterLB = R(-infinity);
606 enterMax = R(-infinity);
607 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
608 }
609
610 objChange -= theLCbound[enterIdx] * SPxLPBase<R>::upper(enterIdx);
611 objChange -= theUCbound[enterIdx] * SPxLPBase<R>::lower(enterIdx);
612 break;
613
614 default:
615 throw SPxInternalCodeException("XENTER01 This should never happen.");
616 }
617
618 MSG_DEBUG(std::cout << "DENTER03 SPxSolverBase::getEnterVals() : col " << enterIdx
619 << ": " << enterStat
620 << " -> " << ds.colStatus(enterIdx)
621 << " objChange: " << objChange
622 << std::endl;)
623 }
624
625 else
626 {
627 assert(enterId.isSPxRowId());
628 enterIdx = this->number(SPxRowId(enterId));
629 enterStat = ds.rowStatus(enterIdx);
630 assert(!isBasic(enterStat));
631
632 /* For an #Id# to enter the basis we better recompute the Test value.
633 */
634 if(rep() == ROW)
635 {
636 computePvec(enterIdx);
637 enterTest = computeTest(enterIdx);
638 theTest[enterIdx] = 0;
639 }
640 else
641 {
642 enterTest = coTest()[enterIdx];
643 theCoTest[enterIdx] = 0;
644 }
645
646 switch(enterStat)
647 {
648 // primal/columnwise cases:
649 case SPxBasisBase<R>::Desc::P_ON_UPPER :
650 assert(rep() == COLUMN);
651 enterUB = theURbound[enterIdx];
652 enterLB = theLRbound[enterIdx];
653 enterVal = enterLB;
654 enterMax = enterUB - enterLB;
655 enterPric = (*theCoPvec)[enterIdx];
656 enterRO = this->maxRowObj(enterIdx);
657 objChange -= enterRO * enterVal;
658
659 if(enterUB >= R(infinity))
660 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_LOWER;
661 else if(EQ(enterLB, enterUB))
662 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE;
663 else
664 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH;
665
666 break;
667
668 case SPxBasisBase<R>::Desc::P_ON_LOWER :
669 assert(rep() == COLUMN);
670 enterUB = theURbound[enterIdx];
671 enterLB = theLRbound[enterIdx];
672 enterVal = enterUB;
673 enterMax = enterLB - enterUB;
674 enterPric = (*theCoPvec)[enterIdx];
675 enterRO = this->maxRowObj(enterIdx);
676 objChange -= enterRO * enterVal;
677
678 if(enterLB <= R(-infinity))
679 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_UPPER;
680 else if(EQ(enterLB, enterUB))
681 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE;
682 else
683 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH;
684
685 break;
686
687 case SPxBasisBase<R>::Desc::P_FREE :
688 assert(rep() == COLUMN);
689 #if 1
690 throw SPxInternalCodeException("XENTER02 This should never happen.");
691 #else
692 MSG_ERROR(std::cerr << "EENTER99 ERROR: not yet debugged!" << std::endl;)
693 enterPric = (*theCoPvec)[enterIdx];
694 enterRO = this->maxRowObj(enterIdx);
695 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_UNDEFINED;
696 #endif
697 break;
698
699 // dual/rowwise cases:
700 case SPxBasisBase<R>::Desc::D_ON_UPPER :
701 assert(rep() == ROW);
702 assert(theURbound[enterIdx] < R(infinity));
703 enterUB = theURbound[enterIdx];
704 enterLB = R(-infinity);
705 enterVal = enterUB;
706 enterMax = R(-infinity);
707 enterPric = (*thePvec)[enterIdx];
708 enterRO = this->lhs(enterIdx);
709 objChange -= enterRO * enterVal;
710 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
711 break;
712
713 case SPxBasisBase<R>::Desc::D_ON_LOWER :
714 assert(rep() == ROW);
715 assert(theLRbound[enterIdx] > R(-infinity));
716 enterLB = theLRbound[enterIdx];
717 enterUB = R(infinity);
718 enterVal = enterLB;
719 enterMax = R(infinity);
720 enterPric = (*thePvec)[enterIdx];
721 enterRO = this->rhs(enterIdx);
722 objChange -= enterRO * enterVal;
723 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
724 break;
725
726 case SPxBasisBase<R>::Desc::D_FREE:
727 assert(rep() == ROW);
728 assert(this->rhs(enterIdx) == this->lhs(enterIdx));
729 enterUB = R(infinity);
730 enterLB = R(-infinity);
731 enterVal = 0;
732 enterPric = (*thePvec)[enterIdx];
733 enterRO = this->rhs(enterIdx);
734 enterMax = (enterPric > enterRO) ? R(infinity) : R(-infinity);
735 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_FIXED;
736 break;
737
738 case SPxBasisBase<R>::Desc::D_ON_BOTH :
739 assert(rep() == ROW);
740 enterPric = (*thePvec)[enterIdx];
741
742 if(enterPric > this->rhs(enterIdx))
743 {
744 enterLB = theLRbound[enterIdx];
745 enterVal = enterLB;
746 enterUB = R(infinity);
747 enterMax = R(infinity);
748 enterRO = this->rhs(enterIdx);
749 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
750 }
751 else
752 {
753 enterUB = theURbound[enterIdx];
754 enterVal = enterUB;
755 enterLB = R(-infinity);
756 enterMax = R(-infinity);
757 enterRO = this->lhs(enterIdx);
758 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
759 }
760
761 objChange -= theLRbound[enterIdx] * this->rhs(enterIdx);
762 objChange -= theURbound[enterIdx] * this->lhs(enterIdx);
763 break;
764
765 default:
766 throw SPxInternalCodeException("XENTER03 This should never happen.");
767 }
768
769 MSG_DEBUG(std::cout << "DENTER05 SPxSolverBase::getEnterVals() : row "
770 << enterIdx << ": " << enterStat
771 << " -> " << ds.rowStatus(enterIdx)
772 << " objChange: " << objChange
773 << std::endl;)
774 }
775 }
776
777 /* process leaving variable
778 */
779 template <class R>
780 void SPxSolverBase<R>::getEnterVals2
781 (
782 int leaveIdx,
783 R enterMax,
784 R& leavebound,
785 StableSum<R>& objChange
786 )
787 {
788 int idx;
789 typename SPxBasisBase<R>::Desc& ds = this->desc();
790 SPxId leftId = this->baseId(leaveIdx);
791
792 if(leftId.isSPxRowId())
793 {
794 idx = this->number(SPxRowId(leftId));
795 typename SPxBasisBase<R>::Desc::Status leaveStat = ds.rowStatus(idx);
796
797 switch(leaveStat)
798 {
799 case SPxBasisBase<R>::Desc::P_FIXED :
800 assert(rep() == ROW);
801 throw SPxInternalCodeException("XENTER04 This should never happen.");
802 break;
803
804 case SPxBasisBase<R>::Desc::P_ON_UPPER :
805 assert(rep() == ROW);
806 leavebound = theLBbound[leaveIdx];
807 theLRbound[idx] = leavebound;
808 ds.rowStatus(idx) = this->dualRowStatus(idx);
809
810 switch(ds.rowStatus(idx))
811 {
812 case SPxBasisBase<R>::Desc::D_ON_UPPER :
813 objChange += theURbound[idx] * this->lhs(idx);
814 break;
815
816 case SPxBasisBase<R>::Desc::D_ON_LOWER :
817 objChange += theLRbound[idx] * this->rhs(idx);
818 break;
819
820 case SPxBasisBase<R>::Desc::D_ON_BOTH :
821 objChange += theURbound[idx] * this->lhs(idx);
822 objChange += theLRbound[idx] * this->rhs(idx);
823 break;
824
825 default:
826 break;
827 }
828
829 break;
830
831 case SPxBasisBase<R>::Desc::P_ON_LOWER :
832 assert(rep() == ROW);
833 leavebound = theUBbound[leaveIdx];
834 theURbound[idx] = leavebound;
835 ds.rowStatus(idx) = this->dualRowStatus(idx);
836
837 switch(ds.rowStatus(idx))
838 {
839 case SPxBasisBase<R>::Desc::D_ON_UPPER :
840 objChange += theURbound[idx] * this->lhs(idx);
841 break;
842
843 case SPxBasisBase<R>::Desc::D_ON_LOWER :
844 objChange += theLRbound[idx] * this->rhs(idx);
845 break;
846
847 case SPxBasisBase<R>::Desc::D_ON_BOTH :
848 objChange += theURbound[idx] * this->lhs(idx);
849 objChange += theLRbound[idx] * this->rhs(idx);
850 break;
851
852 default:
853 break;
854 }
855
856 break;
857
858 case SPxBasisBase<R>::Desc::P_FREE :
859 assert(rep() == ROW);
860 #if 1
861 throw SPxInternalCodeException("XENTER05 This should never happen.");
862 #else
863 MSG_ERROR(std::cerr << "EENTER98 ERROR: not yet debugged!" << std::endl;)
864
865 if((*theCoPvec)[leaveIdx] - theLBbound[leaveIdx] <
866 theUBbound[leaveIdx] - (*theCoPvec)[leaveIdx])
867 {
868 leavebound = theLBbound[leaveIdx];
869 theLRbound[idx] = leavebound;
870 }
871 else
872 {
873 leavebound = theUBbound[leaveIdx];
874 theURbound[idx] = leavebound;
875 }
876
877 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::D_UNDEFINED;
878 #endif
879 break;
880
881 // primal/columnwise cases:
882 case SPxBasisBase<R>::Desc::D_UNDEFINED :
883 assert(rep() == COLUMN);
884 throw SPxInternalCodeException("XENTER06 This should never happen.");
885 break;
886
887 case SPxBasisBase<R>::Desc::D_FREE :
888 assert(rep() == COLUMN);
889
890 if(theFvec->delta()[leaveIdx] * enterMax < 0)
891 leavebound = theUBbound[leaveIdx];
892 else
893 leavebound = theLBbound[leaveIdx];
894
895 theLRbound[idx] = leavebound;
896 theURbound[idx] = leavebound;
897 objChange += leavebound * this->maxRowObj(leaveIdx);
898 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED;
899 break;
900
901 case SPxBasisBase<R>::Desc::D_ON_UPPER :
902 assert(rep() == COLUMN);
903 leavebound = theUBbound[leaveIdx];
904 theURbound[idx] = leavebound;
905 objChange += leavebound * this->maxRowObj(leaveIdx);
906 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
907 break;
908
909 case SPxBasisBase<R>::Desc::D_ON_LOWER :
910 assert(rep() == COLUMN);
911 leavebound = theLBbound[leaveIdx];
912 theLRbound[idx] = leavebound;
913 objChange += leavebound * this->maxRowObj(leaveIdx);
914 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
915 break;
916
917 case SPxBasisBase<R>::Desc::D_ON_BOTH :
918 assert(rep() == COLUMN);
919
920 if(enterMax * theFvec->delta()[leaveIdx] < 0)
921 {
922 leavebound = theUBbound[leaveIdx];
923 theURbound[idx] = leavebound;
924 objChange += leavebound * this->maxRowObj(leaveIdx);
925 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
926 }
927 else
928 {
929 leavebound = theLBbound[leaveIdx];
930 theLRbound[idx] = leavebound;
931 objChange += leavebound * this->maxRowObj(leaveIdx);
932 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
933 }
934
935 break;
936
937 default:
938 throw SPxInternalCodeException("XENTER07 This should never happen.");
939 }
940
941 MSG_DEBUG(std::cout << "DENTER06 SPxSolverBase::getEnterVals2(): row "
942 << idx << ": " << leaveStat
943 << " -> " << ds.rowStatus(idx)
944 << " objChange: " << objChange
945 << std::endl;)
946 }
947
948 else
949 {
950 assert(leftId.isSPxColId());
951 idx = this->number(SPxColId(leftId));
952 typename SPxBasisBase<R>::Desc::Status leaveStat = ds.colStatus(idx);
953
954 switch(leaveStat)
955 {
956 case SPxBasisBase<R>::Desc::P_ON_UPPER :
957 assert(rep() == ROW);
958 leavebound = theLBbound[leaveIdx];
959 theLCbound[idx] = leavebound;
960 ds.colStatus(idx) = this->dualColStatus(idx);
961
962 switch(ds.colStatus(idx))
963 {
964 case SPxBasisBase<R>::Desc::D_ON_UPPER :
965 objChange += theUCbound[idx] * this->lower(idx);
966 break;
967
968 case SPxBasisBase<R>::Desc::D_ON_LOWER :
969 objChange += theLCbound[idx] * this->upper(idx);
970 break;
971
972 case SPxBasisBase<R>::Desc::D_ON_BOTH :
973 objChange += theLCbound[idx] * this->upper(idx);
974 objChange += theUCbound[idx] * this->lower(idx);
975 break;
976
977 default:
978 break;
979 }
980
981 break;
982
983 case SPxBasisBase<R>::Desc::P_ON_LOWER :
984 assert(rep() == ROW);
985 leavebound = theUBbound[leaveIdx];
986 theUCbound[idx] = leavebound;
987 ds.colStatus(idx) = this->dualColStatus(idx);
988
989 switch(ds.colStatus(idx))
990 {
991 case SPxBasisBase<R>::Desc::D_ON_UPPER :
992 objChange += theUCbound[idx] * this->lower(idx);
993 break;
994
995 case SPxBasisBase<R>::Desc::D_ON_LOWER :
996 objChange += theLCbound[idx] * this->upper(idx);
997 break;
998
999 case SPxBasisBase<R>::Desc::D_ON_BOTH :
1000 objChange += theLCbound[idx] * this->upper(idx);
1001 objChange += theUCbound[idx] * this->lower(idx);
1002 break;
1003
1004 default:
1005 break;
1006 }
1007
1008 break;
1009
1010 case SPxBasisBase<R>::Desc::P_FREE :
1011 assert(rep() == ROW);
1012
1013 if(theFvec->delta()[leaveIdx] * enterMax > 0)
1014 {
1015 leavebound = theLBbound[leaveIdx];
1016 theLCbound[idx] = leavebound;
1017 }
1018 else
1019 {
1020 leavebound = theUBbound[leaveIdx];
1021 theUCbound[idx] = leavebound;
1022 }
1023
1024 ds.colStatus(idx) = SPxBasisBase<R>::Desc::D_UNDEFINED;
1025 break;
1026
1027 case SPxBasisBase<R>::Desc::P_FIXED:
1028 assert(rep() == ROW);
1029 throw SPxInternalCodeException("XENTER08 This should never happen.");
1030 break;
1031
1032 // primal/columnwise cases:
1033 case SPxBasisBase<R>::Desc::D_FREE :
1034 assert(rep() == COLUMN);
1035
1036 if(theFvec->delta()[leaveIdx] * enterMax > 0)
1037 leavebound = theLBbound[leaveIdx];
1038 else
1039 leavebound = theUBbound[leaveIdx];
1040
1041 theUCbound[idx] =
1042 theLCbound[idx] = leavebound;
1043 objChange += this->maxObj(idx) * leavebound;
1044 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED;
1045 break;
1046
1047 case SPxBasisBase<R>::Desc::D_ON_UPPER :
1048 assert(rep() == COLUMN);
1049 leavebound = theLBbound[leaveIdx];
1050 theLCbound[idx] = leavebound;
1051 objChange += this->maxObj(idx) * leavebound;
1052 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1053 break;
1054
1055 case SPxBasisBase<R>::Desc::D_ON_LOWER :
1056 assert(rep() == COLUMN);
1057 leavebound = theUBbound[leaveIdx];
1058 theUCbound[idx] = leavebound;
1059 objChange += this->maxObj(idx) * leavebound;
1060 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1061 break;
1062
1063 case SPxBasisBase<R>::Desc::D_ON_BOTH :
1064 case SPxBasisBase<R>::Desc::D_UNDEFINED :
1065 assert(rep() == COLUMN);
1066
1067 if(enterMax * theFvec->delta()[leaveIdx] < 0)
1068 {
1069 leavebound = theUBbound[leaveIdx];
1070 theUCbound[idx] = leavebound;
1071 objChange += this->maxObj(idx) * leavebound;
1072 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1073 }
1074 else
1075 {
1076 leavebound = theLBbound[leaveIdx];
1077 theLCbound[idx] = leavebound;
1078 objChange += this->maxObj(idx) * leavebound;
1079 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1080 }
1081
1082 break;
1083
1084 default:
1085 throw SPxInternalCodeException("XENTER09 This should never happen.");
1086 }
1087
1088 MSG_DEBUG(std::cout << "DENTER07 SPxSolverBase::getEnterVals2(): col "
1089 << idx << ": " << leaveStat
1090 << " -> " << ds.colStatus(idx)
1091 << " objChange: " << objChange
1092 << std::endl;)
1093 }
1094 }
1095
1096 template <class R>
1097 void
1098 SPxSolverBase<R>::ungetEnterVal(
1099 SPxId enterId,
1100 typename SPxBasisBase<R>::Desc::Status enterStat,
1101 R leaveVal,
1102 const SVectorBase<R>& vec,
1103 StableSum<R>& objChange
1104 )
1105 {
1106 assert(rep() == COLUMN);
1107 int enterIdx;
1108 typename SPxBasisBase<R>::Desc& ds = this->desc();
1109
1110 if(enterId.isSPxColId())
1111 {
1112 enterIdx = this->number(SPxColId(enterId));
1113
1114 if(enterStat == SPxBasisBase<R>::Desc::P_ON_UPPER)
1115 {
1116 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1117 objChange += theLCbound[enterIdx] * this->maxObj(enterIdx);
1118 }
1119 else
1120 {
1121 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1122 objChange += theUCbound[enterIdx] * this->maxObj(enterIdx);
1123 }
1124
1125 theFrhs->multAdd(leaveVal, vec);
1126 }
1127 else
1128 {
1129 enterIdx = this->number(SPxRowId(enterId));
1130 assert(enterId.isSPxRowId());
1131
1132 if(enterStat == SPxBasisBase<R>::Desc::P_ON_UPPER)
1133 {
1134 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER;
1135 objChange += (theURbound[enterIdx]) * this->maxRowObj(enterIdx);
1136 }
1137 else
1138 {
1139 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER;
1140 objChange += (theLRbound[enterIdx]) * this->maxRowObj(enterIdx);
1141 }
1142
1143 (*theFrhs)[enterIdx] += leaveVal;
1144 }
1145
1146 if(isId(enterId))
1147 {
1148 theTest[enterIdx] = 0;
1149 isInfeasibleCo[enterIdx] = SPxPricer<R>::NOT_VIOLATED;
1150 }
1151 else
1152 {
1153 theCoTest[enterIdx] = 0;
1154 isInfeasible[enterIdx] = SPxPricer<R>::NOT_VIOLATED;
1155 }
1156 }
1157
1158 template <class R>
1159 void SPxSolverBase<R>::rejectEnter(
1160 SPxId enterId,
1161 R enterTest,
1162 typename SPxBasisBase<R>::Desc::Status enterStat
1163 )
1164 {
1165 int enterIdx = this->number(enterId);
1166
1167 if(isId(enterId))
1168 {
1169 theTest[enterIdx] = enterTest;
1170 this->desc().status(enterIdx) = enterStat;
1171 }
1172 else
1173 {
1174 theCoTest[enterIdx] = enterTest;
1175 this->desc().coStatus(enterIdx) = enterStat;
1176 }
1177 }
1178
1179 template <class R>
1180 void SPxSolverBase<R>::computePrimalray4Col(R direction, SPxId enterId)
1181 {
1182 R sign = (direction > 0 ? 1.0 : -1.0);
1183
1184 primalRay.clear();
1185 primalRay.setMax(fVec().delta().size() + 1);
1186
1187 for(int j = 0; j < fVec().delta().size(); ++j)
1188 {
1189 SPxId i = this->baseId(fVec().idx().index(j));
1190
1191 if(i.isSPxColId())
1192 primalRay.add(this->number(SPxColId(i)), sign * fVec().delta().value(j));
1193 }
1194
1195 if(enterId.isSPxColId())
1196 primalRay.add(this->number(SPxColId(enterId)), -sign);
1197 }
1198
1199 template <class R>
1200 void SPxSolverBase<R>::computeDualfarkas4Row(R direction, SPxId enterId)
1201 {
1202 R sign = (direction > 0 ? -1.0 : 1.0);
1203
1204 dualFarkas.clear();
1205 dualFarkas.setMax(fVec().delta().size() + 1);
1206
1207 for(int j = 0; j < fVec().delta().size(); ++j)
1208 {
1209 SPxId spxid = this->baseId(fVec().idx().index(j));
1210
1211 if(spxid.isSPxRowId())
1212 dualFarkas.add(this->number(SPxRowId(spxid)), sign * fVec().delta().value(j));
1213 }
1214
1215 if(enterId.isSPxRowId())
1216 dualFarkas.add(this->number(SPxRowId(enterId)), -sign);
1217 }
1218
1219 template <class R>
1220 bool SPxSolverBase<R>::enter(SPxId& enterId, bool polish)
1221 {
1222 assert(enterId.isValid());
1223 assert(type() == ENTER);
1224 assert(initialized);
1225
1226 SPxId none; // invalid id used when enter fails
1227 R enterTest; // correct test value of entering var
1228 R enterUB; // upper bound of entering variable
1229 R enterLB; // lower bound of entering variable
1230 R enterVal; // current value of entering variable
1231 R enterMax; // maximum value for entering shift
1232 R enterPric; // priced value of entering variable
1233 typename SPxBasisBase<R>::Desc::Status enterStat; // status of entering variable
1234 R enterRO; // rhs/obj of entering variable
1235 StableSum<R> objChange;
1236 const SVectorBase<R>* enterVec = enterVector(enterId);
1237
1238 bool instable = instableEnter;
1239 assert(!instable || instableEnterId.isValid());
1240
1241 getEnterVals(enterId, enterTest, enterUB, enterLB,
1242 enterVal, enterMax, enterPric, enterStat, enterRO, objChange);
1243
1244 if(!polish && enterTest > -epsilon())
1245 {
1246 rejectEnter(enterId, enterTest, enterStat);
1247 this->change(-1, none, 0);
1248
1249 MSG_DEBUG(std::cout << "DENTER08 rejecting false enter pivot" << std::endl;)
1250
1251 return false;
1252 }
1253
1254 /* Before performing the actual basis update, we must determine, how this
1255 is to be accomplished.
1256 */
1257 // BH 2005-11-15: Obviously solve4update() is only called if theFvec.delta()
1258 // is setup (i.e. the indices of the NZEs are stored within it) and there are
1259 // 0 NZEs (???).
1260 // In that case theFvec->delta() is set such that
1261 // Base * theFvec->delta() = enterVec
1262 if(theFvec->delta().isSetup() && theFvec->delta().size() == 0)
1263 SPxBasisBase<R>::solve4update(theFvec->delta(), *enterVec);
1264
1265 #ifdef ENABLE_ADDITIONAL_CHECKS
1266 else
1267 {
1268 // BH 2005-11-29: This code block seems to check the assertion
1269 // || Base * theFvec->delta() - enterVec ||_2 <= entertol()
1270 VectorBase<R> tmp(dim());
1271 // BH 2005-11-15: This cast is necessary since SSVectorBase<R> inherits protected from VectorBase<R>.
1272 tmp = reinterpret_cast<VectorBase<R>&>(theFvec->delta());
1273 this->multBaseWith(tmp);
1274 tmp -= *enterVec;
1275
1276 if(tmp.length() > entertol())
1277 {
1278 // This happens frequently and does usually not hurt, so print these
1279 // warnings only with verbose level INFO2 and higher.
1280 MSG_INFO2((*this->spxout), (*this->spxout) << "WENTER09 fVec updated error = "
1281 << tmp.length() << std::endl;)
1282 }
1283 }
1284
1285 #endif // ENABLE_ADDITIONAL_CHECKS
1286
1287 if(!polish && m_numCycle > m_maxCycle)
1288 {
1289 if(-enterMax > 0)
1290 perturbMaxEnter();
1291 else
1292 perturbMinEnter();
1293 }
1294
1295 R leaveVal = -enterMax;
1296
1297 boundflips = 0;
1298 int leaveIdx = theratiotester->selectLeave(leaveVal, enterTest, polish);
1299
1300 /* in row representation, fixed columns and rows should not leave the basis */
1301 assert(leaveIdx < 0 || !this->baseId(leaveIdx).isSPxColId()
1302 || this->desc().colStatus(this->number(SPxColId(this->baseId(leaveIdx)))) !=
1303 SPxBasisBase<R>::Desc::P_FIXED);
1304 assert(leaveIdx < 0 || !this->baseId(leaveIdx).isSPxRowId()
1305 || this->desc().rowStatus(this->number(SPxRowId(this->baseId(leaveIdx)))) !=
1306 SPxBasisBase<R>::Desc::P_FIXED);
1307
1308 instableEnterVal = 0;
1309 instableEnterId = SPxId();
1310 instableEnter = false;
1311
1312 /*
1313 We now tried to find a variable to leave the basis. If one has been
1314 found, a regular basis update is to be performed.
1315 */
1316 if(leaveIdx >= 0)
1317 {
1318 if(spxAbs(leaveVal) < entertol())
1319 {
1320 if(NE(theUBbound[leaveIdx], theLBbound[leaveIdx])
1321 && enterStat != SPxBasisBase<R>::Desc::P_FREE && enterStat != SPxBasisBase<R>::Desc::D_FREE)
1322 {
1323 m_numCycle++;
1324 enterCycles++;
1325 }
1326 }
1327 else
1328 m_numCycle /= 2;
1329
1330 // setup for updating the copricing vector
1331 if(coSolveVector3 && coSolveVector2)
1332 {
1333 assert(coSolveVector2->isConsistent());
1334 assert(coSolveVector2rhs->isSetup());
1335 assert(coSolveVector3->isConsistent());
1336 assert(coSolveVector3rhs->isSetup());
1337 assert(boundflips > 0);
1338 SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector2, *coSolveVector3
1339 , unitVecs[leaveIdx], *coSolveVector2rhs, *coSolveVector3rhs);
1340 (*theCoPvec) -= (*coSolveVector3);
1341 }
1342 else if(coSolveVector3)
1343 {
1344 assert(coSolveVector3->isConsistent());
1345 assert(coSolveVector3rhs->isSetup());
1346 assert(boundflips > 0);
1347 SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector3, unitVecs[leaveIdx],
1348 *coSolveVector3rhs);
1349 (*theCoPvec) -= (*coSolveVector3);
1350 }
1351 else if(coSolveVector2)
1352 SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector2, unitVecs[leaveIdx],
1353 *coSolveVector2rhs);
1354 else
1355 SPxBasisBase<R>::coSolve(theCoPvec->delta(), unitVecs[leaveIdx]);
1356
1357 if(boundflips > 0)
1358 {
1359 for(int i = coSolveVector3->dim() - 1; i >= 0; --i)
1360 {
1361 if(spxAbs((*coSolveVector3)[i]) > epsilon())
1362 (*thePvec).multAdd(-(*coSolveVector3)[i], (*thecovectors)[i]);
1363 }
1364
1365 // we need to update enterPric in case it was changed by bound flips
1366 if(enterId.isSPxColId())
1367 enterPric = (*theCoPvec)[this->number(SPxColId(enterId))];
1368 else
1369 enterPric = (*thePvec)[this->number(SPxRowId(enterId))];
1370
1371 MSG_DEBUG(std::cout << "IEBFRT02 breakpoints passed / bounds flipped = " << boundflips << std::endl;
1372 )
1373 totalboundflips += boundflips;
1374 }
1375
1376 (*theCoPrhs)[leaveIdx] = enterRO;
1377 theCoPvec->value() = (enterRO - enterPric) / theFvec->delta()[leaveIdx];
1378
1379 if(theCoPvec->value() > epsilon() || theCoPvec->value() < -epsilon())
1380 {
1381 if(pricing() == FULL)
1382 {
1383 thePvec->value() = theCoPvec->value();
1384 setupPupdate();
1385 }
1386
1387 doPupdate();
1388 }
1389
1390 assert(thePvec->isConsistent());
1391 assert(theCoPvec->isConsistent());
1392
1393 assert(!this->baseId(leaveIdx).isSPxRowId()
1394 || this->desc().rowStatus(this->number(SPxRowId(this->baseId(leaveIdx)))) !=
1395 SPxBasisBase<R>::Desc::P_FIXED);
1396 assert(!this->baseId(leaveIdx).isSPxColId()
1397 || this->desc().colStatus(this->number(SPxColId(this->baseId(leaveIdx)))) !=
1398 SPxBasisBase<R>::Desc::P_FIXED);
1399
1400 R leavebound; // bound on which leaving variable moves
1401
1402 try
1403 {
1404 getEnterVals2(leaveIdx, enterMax, leavebound, objChange);
1405 }
1406 catch(const SPxException& F)
1407 {
1408 rejectEnter(enterId, enterTest, enterStat);
1409 this->change(-1, none, 0);
1410 throw F;
1411 }
1412
1413 // process entering variable
1414 theUBbound[leaveIdx] = enterUB;
1415 theLBbound[leaveIdx] = enterLB;
1416
1417 // compute tests:
1418 updateCoTest();
1419
1420 if(pricing() == FULL)
1421 updateTest();
1422
1423 // update feasibility vectors
1424 theFvec->value() = leaveVal;
1425 theFvec->update();
1426 (*theFvec)[leaveIdx] = enterVal - leaveVal;
1427
1428 if(leavebound > epsilon() || leavebound < -epsilon())
1429 theFrhs->multAdd(-leavebound, this->baseVec(leaveIdx));
1430
1431 if(enterVal > epsilon() || enterVal < -epsilon())
1432 theFrhs->multAdd(enterVal, *enterVec);
1433
1434 // update objective funtion value
1435 updateNonbasicValue(objChange);
1436
1437 // change basis matrix
1438 this->change(leaveIdx, enterId, enterVec, &(theFvec->delta()));
1439
1440 return true;
1441 }
1442 /* No leaving vector could be found that would yield a stable pivot step.
1443 */
1444 else if(NE(leaveVal, -enterMax))
1445 {
1446 /* In the ENTER algorithm, when for a selected entering variable we find only
1447 an instable leaving variable, then the basis change is not conducted.
1448 Instead, we save the entering variable's id in instableEnterId and set
1449 the test value to zero, hoping to find a different leaving
1450 variable with a stable leavingvariable.
1451 If this fails, however, and no more entering variable is found, we have to
1452 perform the instable basis change using instableEnterId. In this (and only
1453 in this) case, the flag instableEnter is set to true.
1454
1455 leaveVal != enterMax is the case that selectLeave has found only an instable leaving
1456 variable. We store this leaving variable for later if we are not already in the
1457 instable case */
1458
1459 if(!instable)
1460 {
1461 instableEnterId = enterId;
1462 instableEnterVal = enterTest;
1463
1464 MSG_DEBUG(std::cout << "DENTER09 rejecting enter pivot and looking for others" << std::endl;)
1465
1466 rejectEnter(enterId, enterTest / 10.0, enterStat);
1467 this->change(-1, none, 0);
1468 }
1469 else
1470 {
1471 MSG_DEBUG(std::cout << "DENTER10 rejecting enter pivot in instable state, resetting values" <<
1472 std::endl;)
1473 rejectEnter(enterId, enterTest, enterStat);
1474 this->change(-1, none, 0);
1475 }
1476
1477 return false;
1478 }
1479 /* No leaving vector has been selected from the basis. However, if the
1480 shift amount for |fVec| is bounded, we are in the case, that the
1481 entering variable is moved from one bound to its other, before any of
1482 the basis feasibility variables reaches their bound. This may only
1483 happen in primal/columnwise case with upper and lower bounds on
1484 variables.
1485 */
1486 else if(!polish && leaveVal < R(infinity) && leaveVal > R(-infinity))
1487 {
1488 assert(rep() == COLUMN);
1489 assert(EQ(leaveVal, -enterMax));
1490
1491 this->change(-1, enterId, enterVec);
1492
1493 theFvec->value() = leaveVal;
1494 theFvec->update();
1495
1496 ungetEnterVal(enterId, enterStat, leaveVal, *enterVec, objChange);
1497
1498 // update objective funtion value
1499 updateNonbasicValue(objChange);
1500
1501 MSG_DEBUG(std::cout << "DENTER11 moving entering variable from one bound to the other" << std::endl;
1502 )
1503
1504 return false;
1505 }
1506 /* No variable could be selected to leave the basis and even the entering
1507 variable is unbounded --- this is a failure.
1508 */
1509 else
1510 {
1511 /* The following line originally was in the "lastUpdate() > 1" case;
1512 we need it in the INFEASIBLE/UNBOUNDED case, too, to have the
1513 basis descriptor at the correct size.
1514 */
1515 rejectEnter(enterId, enterTest, enterStat);
1516 this->change(-1, none, 0);
1517
1518 if(polish)
1519 return false;
1520
1521 else if(this->lastUpdate() > 1)
1522 {
1523 MSG_INFO3((*this->spxout), (*this->spxout) << "IENTER01 factorization triggered in "
1524 << "enter() for feasibility test" << std::endl;)
1525
1526 try
1527 {
1528 factorize();
1529 }
1530 catch(const SPxStatusException& E)
1531 {
1532 // don't exit immediately but handle the singularity correctly
1533 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::SINGULAR);
1534 MSG_INFO3((*this->spxout), (*this->spxout) << "Caught exception in factorization: " << E.what() <<
1535 std::endl;)
1536 }
1537
1538 /* after a factorization, the entering column/row might not be infeasible or suboptimal anymore, hence we do
1539 * not try to call leave(leaveIdx), but rather return to the main solving loop and call the pricer again
1540 */
1541 return false;
1542 }
1543
1544 /* do not exit with status infeasible or unbounded if there is only a very small violation
1545 * ROW: recompute the primal variables and activities for another, more precise, round of pricing
1546 */
1547 else if(spxAbs(enterTest) < entertol())
1548 {
1549 MSG_INFO3((*this->spxout), (*this->spxout) << "IENTER11 clean up step to reduce numerical errors" <<
1550 std::endl;)
1551
1552 SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs);
1553 computePvec();
1554 computeCoTest();
1555 computeTest();
1556
1557 return false;
1558 }
1559
1560 MSG_INFO3((*this->spxout), (*this->spxout) << "IENTER02 unboundedness/infeasibility found in "
1561 << "enter()" << std::endl;)
1562
1563 if(rep() == ROW)
1564 {
1565 computeDualfarkas4Row(leaveVal, enterId);
1566 setBasisStatus(SPxBasisBase<R>::INFEASIBLE);
1567 }
1568 /**@todo if shift() is not zero, we must not conclude primal unboundedness */
1569 else
1570 {
1571 computePrimalray4Col(leaveVal, enterId);
1572 setBasisStatus(SPxBasisBase<R>::UNBOUNDED);
1573 }
1574
1575 return false;
1576 }
1577 }
1578 } // namespace soplex
1579