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 "soplex/spxdefines.h"
|
(1) Event include_recursion: |
#include file "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scipoptsuite-8.0.1/soplex/src/soplex/spxdevexpr.h" includes itself: spxdevexpr.h -> spxdevexpr.hpp -> spxdevexpr.h |
|
(2) Event caretline: |
^ |
17 #include "soplex/spxdevexpr.h"
18
19 #define DEVEX_REFINETOL 2.0
20
21 namespace soplex
22 {
23 // Definition of signature to avoid the specialization after instantiation error
24
25 template <class R>
26 void SPxDevexPR<R>::load(SPxSolverBase<R>* base)
27 {
28 this->thesolver = base;
29 setRep(base->rep());
30 assert(isConsistent());
31 }
32
33 template <class R>
34 bool SPxDevexPR<R>::isConsistent() const
35 {
36 #ifdef ENABLE_CONSISTENCY_CHECKS
37
38 if(this->thesolver != 0)
39 if(this->thesolver->weights.dim() != this->thesolver->coDim()
40 || this->thesolver->coWeights.dim() != this->thesolver->dim())
41 return MSGinconsistent("SPxDevexPR");
42
43 #endif
44
45 return true;
46 }
47
48 template <class R>
49 void SPxDevexPR<R>::setupWeights(typename SPxSolverBase<R>::Type tp)
50 {
51 int i;
52 int coWeightSize = 0;
53 int weightSize = 0;
54
55 VectorBase<R>& weights = this->thesolver->weights;
56 VectorBase<R>& coWeights = this->thesolver->coWeights;
57
58 if(tp == SPxSolverBase<R>::ENTER)
59 {
60 coWeights.reDim(this->thesolver->dim(), false);
61
62 for(i = this->thesolver->dim() - 1; i >= coWeightSize; --i)
63 coWeights[i] = 2.0;
64
65 weights.reDim(this->thesolver->coDim(), false);
66
67 for(i = this->thesolver->coDim() - 1; i >= weightSize; --i)
68 weights[i] = 2.0;
69 }
70 else
71 {
72 coWeights.reDim(this->thesolver->dim(), false);
73
74 for(i = this->thesolver->dim() - 1; i >= coWeightSize; --i)
75 coWeights[i] = 1.0;
76 }
77
78 this->thesolver->weightsAreSetup = true;
79 }
80
81 template <class R>
82 void SPxDevexPR<R>::setType(typename SPxSolverBase<R>::Type tp)
83 {
84 setupWeights(tp);
85 refined = false;
86
87 bestPrices.clear();
88 bestPrices.setMax(this->thesolver->dim());
89 prices.reSize(this->thesolver->dim());
90
91 if(tp == SPxSolverBase<R>::ENTER)
92 {
93 bestPricesCo.clear();
94 bestPricesCo.setMax(this->thesolver->coDim());
95 pricesCo.reSize(this->thesolver->coDim());
96 }
97
98 assert(isConsistent());
99 }
100
101 /**@todo suspicious: Shouldn't the relation between dim, coDim, Vecs,
102 * and CoVecs be influenced by the representation ?
103 */
104 template <class R>
105 void SPxDevexPR<R>::setRep(typename SPxSolverBase<R>::Representation)
106 {
107 if(this->thesolver != 0)
108 {
109 // resize weights and initialize new entries
110 addedVecs(this->thesolver->coDim());
111 addedCoVecs(this->thesolver->dim());
112 assert(isConsistent());
113 }
114 }
115
116 // A namespace to avoid collision with other pricers
117 namespace devexpr
118 {
119 template <class R>
120 R computePrice(R viol, R weight, R tol)
121 {
122 if(weight < tol)
123 return viol * viol / tol;
124 else
125 return viol * viol / weight;
126 }
127 }
128
129 template <class R>
130 int SPxDevexPR<R>::buildBestPriceVectorLeave(R feastol)
131 {
132 int idx;
133 int nsorted;
134 R fTesti;
135 const R* fTest = this->thesolver->fTest().get_const_ptr();
136 const R* cpen = this->thesolver->coWeights.get_const_ptr();
137 typename SPxPricer<R>::IdxElement price;
138 prices.clear();
139 bestPrices.clear();
140
141 // TODO we should check infeasiblities for duplicates or loop over dimension
142 // bestPrices may then also contain duplicates!
143 // construct vector of all prices
144 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i)
145 {
146 idx = this->thesolver->infeasibilities.index(i);
147 fTesti = fTest[idx];
148
149 if(fTesti < -feastol)
150 {
151 this->thesolver->isInfeasible[idx] = this->VIOLATED;
152 price.idx = idx;
153 price.val = devexpr::computePrice(fTesti, cpen[idx], feastol);
154 prices.append(price);
155 }
156 }
157
158 // set up structures for the quicksort implementation
159 this->compare.elements = prices.get_const_ptr();
160 // do a partial sort to move the best ones to the front
161 // TODO this can be done more efficiently, since we only need the indices
162 nsorted = SPxQuicksortPart(prices.get_ptr(), this->compare, 0, prices.size(), HYPERPRICINGSIZE);
163
164 // copy indices of best values to bestPrices
165 for(int i = 0; i < nsorted; ++i)
166 {
167 bestPrices.addIdx(prices[i].idx);
168 this->thesolver->isInfeasible[prices[i].idx] = this->VIOLATED_AND_CHECKED;
169 }
170
171 if(nsorted > 0)
172 return prices[0].idx;
173 else
174 return -1;
175 }
176
177 template <class R>
178 int SPxDevexPR<R>::selectLeave()
179 {
180 int retid;
181
182 if(this->thesolver->hyperPricingLeave && this->thesolver->sparsePricingLeave)
183 {
184 if(bestPrices.size() < 2 || this->thesolver->basis().lastUpdate() == 0)
185 {
186 // call init method to build up price-vector and return index of largest price
187 retid = buildBestPriceVectorLeave(this->theeps);
188 }
189 else
190 retid = selectLeaveHyper(this->theeps);
191 }
192 else if(this->thesolver->sparsePricingLeave)
193 retid = selectLeaveSparse(this->theeps);
194 else
195 retid = selectLeaveX(this->theeps);
196
197 if(retid < 0 && !refined)
198 {
199 refined = true;
200 MSG_INFO3((*this->thesolver->spxout),
201 (*this->thesolver->spxout) << "WDEVEX02 trying refinement step..\n";)
202 retid = selectLeaveX(this->theeps / DEVEX_REFINETOL);
203 }
204
205 assert(retid < this->thesolver->dim());
206
207 return retid;
208 }
209
210 template <class R>
211 int SPxDevexPR<R>::selectLeaveX(R feastol, int start, int incr)
212 {
213 R x;
214
215 const R* fTest = this->thesolver->fTest().get_const_ptr();
216 const R* cpen = this->thesolver->coWeights.get_const_ptr();
217 R best = 0;
218 int bstI = -1;
219 int end = this->thesolver->coWeights.dim();
220
221 for(; start < end; start += incr)
222 {
223 if(fTest[start] < -feastol)
224 {
225 x = devexpr::computePrice(fTest[start], cpen[start], feastol);
226
227 if(x > best)
228 {
229 best = x;
230 bstI = start;
231 last = cpen[start];
232 }
233 }
234 }
235
236 return bstI;
237 }
238
239 template <class R>
240 int SPxDevexPR<R>::selectLeaveSparse(R feastol)
241 {
242 R x;
243
244 const R* fTest = this->thesolver->fTest().get_const_ptr();
245 const R* cpen = this->thesolver->coWeights.get_const_ptr();
246 R best = 0;
247 int bstI = -1;
248 int idx = -1;
249
250 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i)
251 {
252 idx = this->thesolver->infeasibilities.index(i);
253 x = fTest[idx];
254
255 if(x < -feastol)
256 {
257 x = devexpr::computePrice(x, cpen[idx], feastol);
258
259 if(x > best)
260 {
261 best = x;
262 bstI = idx;
263 last = cpen[idx];
264 }
265 }
266 else
267 {
268 this->thesolver->infeasibilities.remove(i);
269 assert(this->thesolver->isInfeasible[idx] == this->VIOLATED
270 || this->thesolver->isInfeasible[idx] == this->VIOLATED_AND_CHECKED);
271 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
272 }
273 }
274
275 return bstI;
276 }
277
278 template <class R>
279 int SPxDevexPR<R>::selectLeaveHyper(R feastol)
280 {
281 R x;
282
283 const R* fTest = this->thesolver->fTest().get_const_ptr();
284 const R* cpen = this->thesolver->coWeights.get_const_ptr();
285 R best = 0;
286 R leastBest = -1;
287 int bstI = -1;
288 int idx = -1;
289
290 // find the best price from the short candidate list
291 for(int i = bestPrices.size() - 1; i >= 0; --i)
292 {
293 idx = bestPrices.index(i);
294 x = fTest[idx];
295
296 if(x < -feastol)
297 {
298 x = devexpr::computePrice(x, cpen[idx], feastol);
299
300 assert(x >= 0);
301
302 // update the best price of candidate list
303 if(x > best)
304 {
305 best = x;
306 bstI = idx;
307 last = cpen[idx];
308 }
309
310 // update the smallest price of candidate list
311 if(x < leastBest || leastBest < 0)
312 leastBest = x;
313 }
314 else
315 {
316 bestPrices.remove(i);
317 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
318 }
319 }
320
321 // scan the updated indices for a better price
322 for(int i = this->thesolver->updateViols.size() - 1; i >= 0; --i)
323 {
324 idx = this->thesolver->updateViols.index(i);
325
326 // only look at indeces that were not checked already
327 if(this->thesolver->isInfeasible[idx] == this->VIOLATED)
328 {
329 x = fTest[idx];
330 assert(x < -feastol);
331 x = devexpr::computePrice(x, cpen[idx], feastol);
332
333 if(x > leastBest)
334 {
335 if(x > best)
336 {
337 best = x;
338 bstI = idx;
339 last = cpen[idx];
340 }
341
342 // put index into candidate list
343 this->thesolver->isInfeasible[idx] = this->VIOLATED_AND_CHECKED;
344 bestPrices.addIdx(idx);
345 }
346 }
347 }
348
349 return bstI;
350 }
351
352 template <class R>
353 void SPxDevexPR<R>::left4(int n, SPxId id)
354 {
355 VectorBase<R>& coWeights = this->thesolver->coWeights;
356
357 if(id.isValid())
358 {
359 int i, j;
360 R x;
361 const R* rhoVec = this->thesolver->fVec().delta().values();
362 R rhov_1 = 1 / rhoVec[n];
363 R beta_q = this->thesolver->coPvec().delta().length2() * rhov_1 * rhov_1;
364
365 #ifndef NDEBUG
366
367 if(spxAbs(rhoVec[n]) < this->theeps)
368 {
369 MSG_INFO3((*this->thesolver->spxout), (*this->thesolver->spxout) << "WDEVEX01: rhoVec = "
370 << rhoVec[n] << " with smaller absolute value than this->theeps = " << this->theeps << std::endl;)
371 }
372
373 #endif // NDEBUG
374
375 // Update #coPenalty# vector
376 const IdxSet& rhoIdx = this->thesolver->fVec().idx();
377 int len = this->thesolver->fVec().idx().size();
378
379 for(i = len - 1; i >= 0; --i)
380 {
381 j = rhoIdx.index(i);
382 x = rhoVec[j] * rhoVec[j] * beta_q;
383 // if(x > coPenalty[j])
384 coWeights[j] += x;
385 }
386
387 coWeights[n] = beta_q;
388 }
389 }
390
391 template <class R>
392 SPxId SPxDevexPR<R>::buildBestPriceVectorEnterDim(R& best, R feastol)
393 {
394 int idx;
395 int nsorted;
396 R x;
397 const R* coTest = this->thesolver->coTest().get_const_ptr();
398 const R* cpen = this->thesolver->coWeights.get_const_ptr();
399 typename SPxPricer<R>::IdxElement price;
400 prices.clear();
401 bestPrices.clear();
402
403 // construct vector of all prices
404 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i)
405 {
406 idx = this->thesolver->infeasibilities.index(i);
407 x = coTest[idx];
408
409 if(x < -feastol)
410 {
411 this->thesolver->isInfeasible[idx] = this->VIOLATED;
412 price.idx = idx;
413 price.val = devexpr::computePrice(x, cpen[idx], feastol);
414 prices.append(price);
415 }
416 else
417 {
418 this->thesolver->infeasibilities.remove(i);
419 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
420 }
421 }
422
423 // set up structures for the quicksort implementation
424 this->compare.elements = prices.get_const_ptr();
425 // do a partial sort to move the best ones to the front
426 // TODO this can be done more efficiently, since we only need the indices
427 nsorted = SPxQuicksortPart(prices.get_ptr(), this->compare, 0, prices.size(), HYPERPRICINGSIZE);
428
429 // copy indices of best values to bestPrices
430 for(int i = 0; i < nsorted; ++i)
431 {
432 bestPrices.addIdx(prices[i].idx);
433 this->thesolver->isInfeasible[prices[i].idx] = this->VIOLATED_AND_CHECKED;
434 }
435
436 if(nsorted > 0)
437 {
438 best = prices[0].val;
439 return this->thesolver->coId(prices[0].idx);
440 }
441 else
442 return SPxId();
443 }
444
445 template <class R>
446 SPxId SPxDevexPR<R>::buildBestPriceVectorEnterCoDim(R& best, R feastol)
447 {
448 int idx;
449 int nsorted;
450 R x;
451 const R* test = this->thesolver->test().get_const_ptr();
452 const R* pen = this->thesolver->weights.get_const_ptr();
453 typename SPxPricer<R>::IdxElement price;
454 pricesCo.clear();
455 bestPricesCo.clear();
456
457 // construct vector of all prices
458 for(int i = this->thesolver->infeasibilitiesCo.size() - 1; i >= 0; --i)
459 {
460 idx = this->thesolver->infeasibilitiesCo.index(i);
461 x = test[idx];
462
463 if(x < -feastol)
464 {
465 this->thesolver->isInfeasibleCo[idx] = this->VIOLATED;
466 price.idx = idx;
467 price.val = devexpr::computePrice(x, pen[idx], feastol);
468 pricesCo.append(price);
469 }
470 else
471 {
472 this->thesolver->infeasibilitiesCo.remove(i);
473 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED;
474 }
475 }
476
477 // set up structures for the quicksort implementation
478 this->compare.elements = pricesCo.get_const_ptr();
479 // do a partial sort to move the best ones to the front
480 // TODO this can be done more efficiently, since we only need the indices
481 nsorted = SPxQuicksortPart(pricesCo.get_ptr(), this->compare, 0, pricesCo.size(), HYPERPRICINGSIZE);
482
483 // copy indices of best values to bestPrices
484 for(int i = 0; i < nsorted; ++i)
485 {
486 bestPricesCo.addIdx(pricesCo[i].idx);
487 this->thesolver->isInfeasibleCo[pricesCo[i].idx] = this->VIOLATED_AND_CHECKED;
488 }
489
490 if(nsorted > 0)
491 {
492 best = pricesCo[0].val;
493 return this->thesolver->id(pricesCo[0].idx);
494 }
495 else
496 return SPxId();
497 }
498
499 template <class R>
500 SPxId SPxDevexPR<R>::selectEnter()
501 {
502 assert(this->thesolver != 0);
503
504 SPxId enterId;
505
506 enterId = selectEnterX(this->theeps);
507
508 if(enterId.isSPxColId() && this->thesolver->isBasic(SPxColId(enterId)))
509 enterId.info = 0;
510
511 if(enterId.isSPxRowId() && this->thesolver->isBasic(SPxRowId(enterId)))
512 enterId.info = 0;
513
514 if(!enterId.isValid() && !refined)
515 {
516 refined = true;
517 MSG_INFO3((*this->thesolver->spxout),
518 (*this->thesolver->spxout) << "WDEVEX02 trying refinement step..\n";)
519 enterId = selectEnterX(this->theeps / DEVEX_REFINETOL);
520
521 if(enterId.isSPxColId() && this->thesolver->isBasic(SPxColId(enterId)))
522 enterId.info = 0;
523
524 if(enterId.isSPxRowId() && this->thesolver->isBasic(SPxRowId(enterId)))
525 enterId.info = 0;
526 }
527
528 return enterId;
529 }
530
531 // choose the best entering index among columns and rows but prefer sparsity
532 template <class R>
533 SPxId SPxDevexPR<R>::selectEnterX(R tol)
534 {
535 SPxId enterId;
536 SPxId enterCoId;
537 R best;
538 R bestCo;
539
540 best = 0;
541 bestCo = 0;
542 last = 1.0;
543
544 // avoid uninitialized value later on in entered4X()
545 last = 1.0;
546
547 if(this->thesolver->hyperPricingEnter && !refined)
548 {
549 if(bestPrices.size() < 2 || this->thesolver->basis().lastUpdate() == 0)
550 enterCoId = (this->thesolver->sparsePricingEnter) ? buildBestPriceVectorEnterDim(best,
551 tol) : selectEnterDenseDim(best, tol);
552 else
553 enterCoId = (this->thesolver->sparsePricingEnter) ? selectEnterHyperDim(best,
554 tol) : selectEnterDenseDim(best, tol);
555
556 if(bestPricesCo.size() < 2 || this->thesolver->basis().lastUpdate() == 0)
557 enterId = (this->thesolver->sparsePricingEnterCo) ? buildBestPriceVectorEnterCoDim(bestCo,
558 tol) : selectEnterDenseCoDim(bestCo, tol);
559 else
560 enterId = (this->thesolver->sparsePricingEnterCo) ? selectEnterHyperCoDim(bestCo,
561 tol) : selectEnterDenseCoDim(bestCo, tol);
562 }
563 else
564 {
565 enterCoId = (this->thesolver->sparsePricingEnter
566 && !refined) ? selectEnterSparseDim(best, tol) : selectEnterDenseDim(best, tol);
567 enterId = (this->thesolver->sparsePricingEnterCo
568 && !refined) ? selectEnterSparseCoDim(bestCo, tol) : selectEnterDenseCoDim(bestCo, tol);
569 }
570
571 // prefer coIds to increase the number of unit vectors in the basis matrix, i.e., rows in colrep and cols in rowrep
572 if(enterCoId.isValid() && (best > SPARSITY_TRADEOFF * bestCo || !enterId.isValid()))
573 return enterCoId;
574 else
575 return enterId;
576 }
577
578 template <class R>
579 SPxId SPxDevexPR<R>::selectEnterHyperDim(R& best, R feastol)
580 {
581 const R* cTest = this->thesolver->coTest().get_const_ptr();
582 const R* cpen = this->thesolver->coWeights.get_const_ptr();
583 R leastBest = -1;
584 R x;
585 int enterIdx = -1;
586 int idx;
587
588 // find the best price from short candidate list
589 for(int i = bestPrices.size() - 1; i >= 0; --i)
590 {
591 idx = bestPrices.index(i);
592 x = cTest[idx];
593
594 if(x < -feastol)
595 {
596 x = devexpr::computePrice(x, cpen[idx], feastol);
597
598 assert(x >= 0);
599
600 // update the best price of candidate list
601 if(x > best)
602 {
603 best = x;
604 enterIdx = idx;
605 last = cpen[idx];
606 }
607
608 // update the smallest price of candidate list
609 if(x < leastBest || leastBest < 0)
610 leastBest = x;
611 }
612 else
613 {
614 bestPrices.remove(i);
615 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
616 }
617 }
618
619 // scan the updated indeces for a better price
620 for(int i = this->thesolver->updateViols.size() - 1; i >= 0; --i)
621 {
622 idx = this->thesolver->updateViols.index(i);
623
624 // only look at indeces that were not checked already
625 if(this->thesolver->isInfeasible[idx] == this->VIOLATED)
626 {
627 x = cTest[idx];
628
629 if(x < -feastol)
630 {
631 x = devexpr::computePrice(x, cpen[idx], feastol);
632
633 if(x > leastBest)
634 {
635 if(x > best)
636 {
637 best = x;
638 enterIdx = idx;
639 last = cpen[idx];
640 }
641
642 // put index into candidate list
643 this->thesolver->isInfeasible[idx] = this->VIOLATED_AND_CHECKED;
644 bestPrices.addIdx(idx);
645 }
646 }
647 else
648 {
649 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
650 }
651 }
652 }
653
654 if(enterIdx >= 0)
655 return this->thesolver->coId(enterIdx);
656 else
657 return SPxId();
658 }
659
660 template <class R>
661 SPxId SPxDevexPR<R>::selectEnterHyperCoDim(R& best, R feastol)
662 {
663 const R* test = this->thesolver->test().get_const_ptr();
664 const R* pen = this->thesolver->weights.get_const_ptr();
665 R leastBest = -1;
666 R x;
667 int enterIdx = -1;
668 int idx;
669
670 // find the best price from short candidate list
671 for(int i = bestPricesCo.size() - 1; i >= 0; --i)
672 {
673 idx = bestPricesCo.index(i);
674 x = test[idx];
675
676 if(x < -feastol)
677 {
678 x = devexpr::computePrice(x, pen[idx], feastol);
679
680 assert(x >= 0);
681
682 // update the best price of candidate list
683 if(x > best)
684 {
685 best = x;
686 enterIdx = idx;
687 last = pen[idx];
688 }
689
690 // update the smallest price of candidate list
691 if(x < leastBest || leastBest < 0)
692 leastBest = x;
693 }
694 else
695 {
696 bestPricesCo.remove(i);
697 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED;
698 }
699 }
700
701 //scan the updated indeces for a better price
702 for(int i = this->thesolver->updateViolsCo.size() - 1; i >= 0; --i)
703 {
704 idx = this->thesolver->updateViolsCo.index(i);
705
706 // only look at indeces that were not checked already
707 if(this->thesolver->isInfeasibleCo[idx] == this->VIOLATED)
708 {
709 x = test[idx];
710
711 if(x < -feastol)
712 {
713 x = devexpr::computePrice(x, pen[idx], feastol);
714
715 if(x > leastBest)
716 {
717 if(x > best)
718 {
719 best = x;
720 enterIdx = idx;
721 last = pen[idx];
722 }
723
724 // put index into candidate list
725 this->thesolver->isInfeasibleCo[idx] = this->VIOLATED_AND_CHECKED;
726 bestPricesCo.addIdx(idx);
727 }
728 }
729 else
730 {
731 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED;
732 }
733 }
734 }
735
736 if(enterIdx >= 0)
737 return this->thesolver->id(enterIdx);
738 else
739 return SPxId();
740 }
741
742
743 template <class R>
744 SPxId SPxDevexPR<R>::selectEnterSparseDim(R& best, R feastol)
745 {
746 const R* cTest = this->thesolver->coTest().get_const_ptr();
747 const R* cpen = this->thesolver->coWeights.get_const_ptr();
748 int enterIdx = -1;
749 int idx;
750 R x;
751
752 assert(this->thesolver->coWeights.dim() == this->thesolver->coTest().dim());
753
754 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i)
755 {
756 idx = this->thesolver->infeasibilities.index(i);
757 x = cTest[idx];
758
759 if(x < -feastol)
760 {
761 x = devexpr::computePrice(x, cpen[idx], feastol);
762
763 if(x > best)
764 {
765 best = x;
766 enterIdx = idx;
767 last = cpen[idx];
768 }
769 }
770 else
771 {
772 this->thesolver->infeasibilities.remove(i);
773 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
774 }
775 }
776
777 if(enterIdx >= 0)
778 return this->thesolver->coId(enterIdx);
779
780 return SPxId();
781 }
782
783
784 template <class R>
785 SPxId SPxDevexPR<R>::selectEnterSparseCoDim(R& best, R feastol)
786 {
787 const R* test = this->thesolver->test().get_const_ptr();
788 const R* pen = this->thesolver->weights.get_const_ptr();
789 int enterIdx = -1;
790 int idx;
791 R x;
792
793 assert(this->thesolver->weights.dim() == this->thesolver->test().dim());
794
795 for(int i = this->thesolver->infeasibilitiesCo.size() - 1; i >= 0; --i)
796 {
797 idx = this->thesolver->infeasibilitiesCo.index(i);
798 x = test[idx];
799
800 if(x < -feastol)
801 {
802 x = devexpr::computePrice(x, pen[idx], feastol);
803
804 if(x > best)
805 {
806 best = x;
807 enterIdx = idx;
808 last = pen[idx];
809 }
810 }
811 else
812 {
813 this->thesolver->infeasibilitiesCo.remove(i);
814 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED;
815 }
816 }
817
818 if(enterIdx >= 0)
819 return this->thesolver->id(enterIdx);
820
821 return SPxId();
822 }
823
824
825 template <class R>
826 SPxId SPxDevexPR<R>::selectEnterDenseDim(R& best, R feastol, int start, int incr)
827 {
828 const R* cTest = this->thesolver->coTest().get_const_ptr();
829 const R* cpen = this->thesolver->coWeights.get_const_ptr();
830 int end = this->thesolver->coWeights.dim();
831 int enterIdx = -1;
832 R x;
833
834 assert(end == this->thesolver->coTest().dim());
835
836 for(; start < end; start += incr)
837 {
838 x = cTest[start];
839
840 if(x < -feastol)
841 {
842 x = devexpr::computePrice(x, cpen[start], feastol);
843
844 if(x > best)
845 {
846 best = x;
847 enterIdx = start;
848 last = cpen[start];
849 }
850 }
851 }
852
853 if(enterIdx >= 0)
854 return this->thesolver->coId(enterIdx);
855
856 return SPxId();
857 }
858
859 template <class R>
860 SPxId SPxDevexPR<R>::selectEnterDenseCoDim(R& best, R feastol, int start, int incr)
861 {
862 const R* test = this->thesolver->test().get_const_ptr();
863 const R* pen = this->thesolver->weights.get_const_ptr();
864 int end = this->thesolver->weights.dim();
865 int enterIdx = -1;
866 R x;
867
868 assert(end == this->thesolver->test().dim());
869
870 for(; start < end; start += incr)
871 {
872 x = test[start];
873
874 if(test[start] < -feastol)
875 {
876 x = devexpr::computePrice(x, pen[start], feastol);
877
878 if(x > best)
879 {
880 best = x;
881 enterIdx = start;
882 last = pen[start];
883 }
884 }
885 }
886
887 if(enterIdx >= 0)
888 return this->thesolver->id(enterIdx);
889
890 return SPxId();
891 }
892
893
894 /**@todo suspicious: the pricer should be informed, that variable id
895 has entered the basis at position n, but the id is not used here
896 (this is true for all pricers)
897 */
898 template <class R>
899 void SPxDevexPR<R>::entered4(SPxId /*id*/, int n)
900 {
901 VectorBase<R>& weights = this->thesolver->weights;
902 VectorBase<R>& coWeights = this->thesolver->coWeights;
903
904 if(n >= 0 && n < this->thesolver->dim())
905 {
906 const R* pVec = this->thesolver->pVec().delta().values();
907 const IdxSet& pIdx = this->thesolver->pVec().idx();
908 const R* coPvec = this->thesolver->coPvec().delta().values();
909 const IdxSet& coPidx = this->thesolver->coPvec().idx();
910 R xi_p = 1 / this->thesolver->fVec().delta()[n];
911 int i, j;
912
913 assert(this->thesolver->fVec().delta()[n] > this->thesolver->epsilon()
914 || this->thesolver->fVec().delta()[n] < -this->thesolver->epsilon());
915
916 xi_p = xi_p * xi_p * last;
917
918 for(j = coPidx.size() - 1; j >= 0; --j)
919 {
920 i = coPidx.index(j);
921 coWeights[i] += xi_p * coPvec[i] * coPvec[i];
922
923 if(coWeights[i] <= 1 || coWeights[i] > 1e+6)
924 {
925 setupWeights(SPxSolverBase<R>::ENTER);
926 return;
927 }
928 }
929
930 for(j = pIdx.size() - 1; j >= 0; --j)
931 {
932 i = pIdx.index(j);
933 weights[i] += xi_p * pVec[i] * pVec[i];
934
935 if(weights[i] <= 1 || weights[i] > 1e+6)
936 {
937 setupWeights(SPxSolverBase<R>::ENTER);
938 return;
939 }
940 }
941 }
942 }
943
944 template <class R>
945 void SPxDevexPR<R>::addedVecs(int n)
946 {
947 int initval = (this->thesolver->type() == SPxSolverBase<R>::ENTER) ? 2 : 1;
948 VectorBase<R>& weights = this->thesolver->weights;
949 n = weights.dim();
950 weights.reDim(this->thesolver->coDim());
951
952 for(int i = weights.dim() - 1; i >= n; --i)
953 weights[i] = initval;
954 }
955
956 template <class R>
957 void SPxDevexPR<R>::addedCoVecs(int n)
958 {
959 int initval = (this->thesolver->type() == SPxSolverBase<R>::ENTER) ? 2 : 1;
960 VectorBase<R>& coWeights = this->thesolver->coWeights;
961 n = coWeights.dim();
962 coWeights.reDim(this->thesolver->dim());
963
964 for(int i = coWeights.dim() - 1; i >= n; --i)
965 coWeights[i] = initval;
966 }
967
968 } // namespace soplex
969