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