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   	//TODO may be faster to have a greater zero tolerance for sparse pricing vectors
26   	//     to reduce the number of nonzero entries, e.g. for workVec
27   	
28   	#include <assert.h>
29   	#include <iostream>
30   	
31   	#include "soplex/spxdefines.h"
32   	#include "soplex/spxsteeppr.h"
33   	#include "soplex/random.h"
34   	
35   	#define SOPLEX_STEEP_REFINETOL 2.0
36   	
37   	namespace soplex
38   	{
39   	
40   	template <class R>
41   	void SPxSteepPR<R>::clear()
42   	{
43   	   this->thesolver = 0;
44   	}
45   	
46   	template <class R>
47   	void SPxSteepPR<R>::load(SPxSolverBase<R>* base)
48   	{
49   	   this->thesolver = base;
50   	
51   	   if(base)
52   	   {
53   	      workVec.clear();
54   	      workVec.reDim(base->dim());
55   	      workRhs.clear();
56   	      workRhs.reDim(base->dim());
57   	   }
58   	}
59   	
60   	template <class R>
61   	void SPxSteepPR<R>::setType(typename SPxSolverBase<R>::Type type)
62   	{
63   	   workRhs.setTolerances(this->_tolerances);
64   	
65   	   setupWeights(type);
66   	   workVec.clear();
67   	   workRhs.clear();
68   	   refined = false;
69   	
70   	   bestPrices.clear();
71   	   bestPrices.setMax(this->thesolver->dim());
72   	   prices.reSize(this->thesolver->dim());
73   	
74   	   if(type == SPxSolverBase<R>::ENTER)
75   	   {
76   	      bestPricesCo.clear();
77   	      bestPricesCo.setMax(this->thesolver->coDim());
78   	      pricesCo.reSize(this->thesolver->coDim());
79   	   }
80   	}
81   	
82   	template <class R>
83   	void SPxSteepPR<R>::setupWeights(typename SPxSolverBase<R>::Type type)
84   	{
85   	   int i;
86   	   int endDim = 0;
87   	   int endCoDim = 0;
88   	   VectorBase<R>& weights = this->thesolver->weights;
89   	   VectorBase<R>& coWeights = this->thesolver->coWeights;
90   	
91   	   if(setup == DEFAULT)
92   	   {
93   	      if(type == SPxSolverBase<R>::ENTER)
94   	      {
95   	         if(this->thesolver->weightsAreSetup)
96   	         {
97   	            // check for added/removed rows and adapt norms accordingly
98   	            if(coWeights.dim() < this->thesolver->dim())
99   	               endDim = coWeights.dim();
100  	            else
101  	               endDim = this->thesolver->dim();
102  	
103  	            if(weights.dim() < this->thesolver->coDim())
104  	               endCoDim = weights.dim();
105  	            else
106  	               endCoDim = this->thesolver->coDim();
107  	         }
108  	
109  	         coWeights.reDim(this->thesolver->dim(), false);
110  	
111  	         for(i = this->thesolver->dim() - 1; i >= endDim; --i)
112  	            coWeights[i] = 2.0;
113  	
114  	         weights.reDim(this->thesolver->coDim(), false);
115  	
116  	         for(i = this->thesolver->coDim() - 1; i >= endCoDim; --i)
117  	            weights[i] = 1.0;
118  	      }
119  	      else
120  	      {
121  	         assert(type == SPxSolverBase<R>::LEAVE);
122  	
123  	         if(this->thesolver->weightsAreSetup)
124  	         {
125  	            // check for added/removed rows and adapt norms accordingly
126  	            if(coWeights.dim() < this->thesolver->dim())
127  	               endDim = coWeights.dim();
128  	            else
129  	               endDim = this->thesolver->dim();
130  	         }
131  	
132  	         coWeights.reDim(this->thesolver->dim(), false);
133  	
134  	         for(i = this->thesolver->dim() - 1; i >= endDim; --i)
135  	            coWeights[i]   = 1.0;
136  	      }
137  	   }
138  	   else
139  	   {
140  	      SPX_MSG_INFO1((*this->thesolver->spxout),
141  	                    (*this->thesolver->spxout) << " --- initializing steepest edge multipliers" << std::endl;)
142  	
143  	      if(type == SPxSolverBase<R>::ENTER)
144  	      {
145  	         coWeights.reDim(this->thesolver->dim(), false);
146  	
147  	         for(i = this->thesolver->dim() - 1; i >= endDim; --i)
148  	            coWeights[i] = 1.0;
149  	
150  	         weights.reDim(this->thesolver->coDim(), false);
151  	
152  	         for(i = this->thesolver->coDim() - 1; i >= endCoDim; --i)
153  	            weights[i] = 1.0 + this->thesolver->vector(i).length2();
154  	      }
155  	      else
156  	      {
157  	         assert(type == SPxSolverBase<R>::LEAVE);
158  	         coWeights.reDim(this->thesolver->dim(), false);
159  	         SSVectorBase<R>  tmp(this->thesolver->dim(), this->thesolver->tolerances());
160  	
161  	         for(i = this->thesolver->dim() - 1; i >= endDim && !this->thesolver->isTimeLimitReached(); --i)
162  	         {
163  	            this->thesolver->basis().coSolve(tmp, this->thesolver->unitVector(i));
164  	            coWeights[i] = tmp.length2();
165  	         }
166  	      }
167  	   }
168  	
169  	   this->thesolver->weightsAreSetup = true;
170  	}
171  	
172  	template <class R>
173  	void SPxSteepPR<R>::setRep(typename SPxSolverBase<R>::Representation)
174  	{
175  	   if(workVec.dim() != this->thesolver->dim())
176  	   {
177  	      VectorBase<R> tmp = this->thesolver->weights;
178  	      this->thesolver->weights = this->thesolver->coWeights;
179  	      this->thesolver->coWeights = tmp;
180  	
181  	      workVec.clear();
182  	      workVec.reDim(this->thesolver->dim());
183  	   }
184  	}
185  	
186  	template <class R>
187  	void SPxSteepPR<R>::left4(int n, SPxId id)
188  	{
189  	   assert(this->thesolver->type() == SPxSolverBase<R>::LEAVE);
190  	
191  	   if(id.isValid())
192  	   {
193  	      R        delta         = 0.1 + 1.0 / this->thesolver->basis().iteration();
194  	      R*       coWeights_ptr = this->thesolver->coWeights.get_ptr();
195  	      const R* workVec_ptr   = workVec.get_const_ptr();
196  	      const R* rhoVec        = this->thesolver->fVec().delta().values();
197  	      R        rhov_1        = 1.0 / rhoVec[n];
198  	      R        beta_q        = this->thesolver->coPvec().delta().length2() * rhov_1 * rhov_1;
199  	
200  	#ifndef NDEBUG
201  	
202  	      if(spxAbs(rhoVec[n]) < this->thetolerance * 0.5)
203  	      {
204  	         SPX_MSG_INFO3((*this->thesolver->spxout), (*this->thesolver->spxout) << "WSTEEP04: rhoVec = "
205  	                       << rhoVec[n] << " with smaller absolute value than 0.5*thetolerance = " << 0.5 * this->thetolerance
206  	                       << std::endl;)
207  	      }
208  	
209  	#endif
210  	
211  	      const IdxSet& rhoIdx = this->thesolver->fVec().idx();
212  	      int           len    = this->thesolver->fVec().idx().size();
213  	
214  	      for(int i = 0; i < len; ++i)
215  	      {
216  	         int  j = rhoIdx.index(i);
217  	         coWeights_ptr[j] += rhoVec[j] * (beta_q * rhoVec[j] - 2.0 * rhov_1 * workVec_ptr[j]);
218  	
219  	         if(coWeights_ptr[j] < delta)
220  	            coWeights_ptr[j] = delta; // coWeights_ptr[j] = delta / (1+delta-x);
221  	         else if(coWeights_ptr[j] >= R(infinity))
222  	            coWeights_ptr[j] = 1.0 / this->thetolerance;
223  	      }
224  	
225  	      coWeights_ptr[n] = beta_q;
226  	   }
227  	}
228  	
229  	
230  	namespace steeppr
231  	{
232  	template <class R>
233  	R computePrice(R viol, R weight, R tol)
234  	{
235  	   if(weight < tol)
236  	      return viol * viol / tol;
237  	   else
238  	      return viol * viol / weight;
239  	}
240  	}
241  	
242  	template <class R>
243  	int SPxSteepPR<R>::buildBestPriceVectorLeave(R feastol)
244  	{
245  	   int idx;
246  	   int nsorted;
247  	   R x;
248  	   const R* fTest = this->thesolver->fTest().get_const_ptr();
249  	   const R* cpen = this->thesolver->coWeights.get_const_ptr();
250  	   typename SPxPricer<R>::IdxElement price;
251  	   prices.clear();
252  	   bestPrices.clear();
253  	
254  	   // construct vector of all prices
255  	   for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i)
256  	   {
257  	      idx = this->thesolver->infeasibilities.index(i);
258  	      x = fTest[idx];
259  	
260  	      if(x < -feastol)
261  	      {
262  	         // it might happen that we call the pricer with a tighter tolerance than what was used when computing the violations
263  	         this->thesolver->isInfeasible[idx] = this->VIOLATED;
264  	         price.val = steeppr::computePrice(x, cpen[idx], feastol);
265  	         price.idx = idx;
266  	         prices.append(price);
267  	      }
268  	   }
269  	
270  	   // set up structures for the quicksort implementation
271  	   this->compare.elements = prices.get_const_ptr();
272  	   // do a partial sort to move the best ones to the front
273  	   // TODO this can be done more efficiently, since we only need the indices
274  	   nsorted = SPxQuicksortPart(prices.get_ptr(), this->compare, 0, prices.size(),
275  	                              SOPLEX_HYPERPRICINGSIZE);
276  	
277  	   // copy indices of best values to bestPrices
278  	   for(int i = 0; i < nsorted; ++i)
279  	   {
280  	      bestPrices.addIdx(prices[i].idx);
281  	      this->thesolver->isInfeasible[prices[i].idx] = this->VIOLATED_AND_CHECKED;
282  	   }
283  	
284  	   if(nsorted > 0)
285  	      return prices[0].idx;
286  	   else
287  	      return -1;
288  	}
289  	
290  	
291  	template <class R>
292  	int SPxSteepPR<R>::selectLeave()
293  	{
294  	   assert(isConsistent());
295  	
296  	   int retid;
297  	
298  	   if(this->thesolver->hyperPricingLeave && this->thesolver->sparsePricingLeave)
299  	   {
300  	      if(bestPrices.size() < 2 || this->thesolver->basis().lastUpdate() == 0)
301  	      {
302  	         // call init method to build up price-vector and return index of largest price
303  	         retid = buildBestPriceVectorLeave(this->thetolerance);
304  	      }
305  	      else
306  	         retid = selectLeaveHyper(this->thetolerance);
307  	   }
308  	   else if(this->thesolver->sparsePricingLeave)
309  	      retid = selectLeaveSparse(this->thetolerance);
310  	   else
311  	      retid = selectLeaveX(this->thetolerance);
312  	
313  	   if(retid < 0 && !refined)
314  	   {
315  	      refined = true;
316  	      SPX_MSG_INFO3((*this->thesolver->spxout),
317  	                    (*this->thesolver->spxout) << "WSTEEP03 trying refinement step..\n";)
318  	      retid = selectLeaveX(this->thetolerance / SOPLEX_STEEP_REFINETOL);
319  	   }
320  	
321  	   if(retid >= 0)
322  	   {
323  	      assert(this->thesolver->coPvec().delta().isConsistent());
324  	      // coPvec().delta() might be not setup after the solve when it contains too many nonzeros.
325  	      // This is intended and forcing to keep the sparsity information leads to a slowdown
326  	      // TODO implement a dedicated solve method for unitvectors
327  	      this->thesolver->basis().coSolve(this->thesolver->coPvec().delta(),
328  	                                       this->thesolver->unitVector(retid));
329  	      assert(this->thesolver->coPvec().delta().isConsistent());
330  	      workRhs.setup_and_assign(this->thesolver->coPvec().delta());
331  	      this->thesolver->setup4solve(&workVec, &workRhs);
332  	   }
333  	
334  	   return retid;
335  	}
336  	
337  	template <class R>
338  	int SPxSteepPR<R>::selectLeaveX(R tol)
339  	{
340  	   const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr();
341  	   const R* fTest         = this->thesolver->fTest().get_const_ptr();
342  	   R best = R(-infinity);
343  	   R x;
344  	   int lastIdx = -1;
345  	
346  	   for(int i = this->thesolver->dim() - 1; i >= 0; --i)
347  	   {
348  	      x = fTest[i];
349  	
350  	      if(x < -tol)
351  	      {
352  	         x = steeppr::computePrice(x, coWeights_ptr[i], tol);
353  	
354  	         if(x > best)
355  	         {
356  	            best = x;
357  	            lastIdx = i;
358  	         }
359  	      }
360  	   }
361  	
362  	   return lastIdx;
363  	}
364  	
365  	template <class R>
366  	int SPxSteepPR<R>::selectLeaveSparse(R tol)
367  	{
368  	   const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr();
369  	   const R* fTest         = this->thesolver->fTest().get_const_ptr();
370  	   R best = R(-infinity);
371  	   R x;
372  	   int lastIdx = -1;
373  	   int idx;
374  	
375  	   for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i)
376  	   {
377  	      idx = this->thesolver->infeasibilities.index(i);
378  	      x = fTest[idx];
379  	
380  	      if(x < -tol)
381  	      {
382  	         x = steeppr::computePrice(x, coWeights_ptr[idx], tol);
383  	
384  	         if(x > best)
385  	         {
386  	            best = x;
387  	            lastIdx = idx;
388  	         }
389  	      }
390  	      else
391  	      {
392  	         this->thesolver->infeasibilities.remove(i);
393  	         assert(this->thesolver->isInfeasible[idx] == this->VIOLATED
394  	                || this->thesolver->isInfeasible[idx] == this->VIOLATED_AND_CHECKED);
395  	         this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
396  	      }
397  	   }
398  	
399  	   return lastIdx;
400  	}
401  	
402  	template <class R>
403  	int SPxSteepPR<R>::selectLeaveHyper(R tol)
404  	{
405  	   const R* coPen = this->thesolver->coWeights.get_const_ptr();
406  	   const R* fTest = this->thesolver->fTest().get_const_ptr();
407  	   R leastBest = -1;
408  	   R best = R(-infinity);
409  	   R x;
410  	   int bestIdx = -1;
411  	   int idx = 0;
412  	
413  	   // find the best price from the short candidate list
414  	   for(int i = bestPrices.size() - 1; i >= 0; --i)
415  	   {
416  	      idx = bestPrices.index(i);
417  	      x = fTest[idx];
418  	
419  	      if(x < -tol)
420  	      {
421  	         assert(this->thesolver->isInfeasible[idx] == this->VIOLATED
422  	                || this->thesolver->isInfeasible[idx] == this->VIOLATED_AND_CHECKED);
423  	         x = steeppr::computePrice(x, coPen[idx], tol);
424  	
425  	         assert(x >= 0);
426  	
427  	         // update the best price of candidate list
428  	         if(x > best)
429  	         {
430  	            best = x;
431  	            bestIdx = idx;
432  	         }
433  	
434  	         // update the smallest price of candidate list
435  	         if(x < leastBest || leastBest < 0)
436  	            leastBest = x;
437  	      }
438  	      else
439  	      {
440  	         bestPrices.remove(i);
441  	         this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
442  	      }
443  	   }
444  	
445  	   // scan the updated indices for a better price
446  	   for(int i = this->thesolver->updateViols.size() - 1; i >= 0; --i)
447  	   {
448  	      idx = this->thesolver->updateViols.index(i);
449  	
450  	      // is this index a candidate for bestPrices?
451  	      if(this->thesolver->isInfeasible[idx] == this->VIOLATED)
452  	      {
453  	         x = fTest[idx];
454  	         assert(x < -tol);
455  	         x = steeppr::computePrice(x, coPen[idx], tol);
456  	
457  	         if(x > leastBest)
458  	         {
459  	            if(x > best)
460  	            {
461  	               best = x;
462  	               bestIdx = idx;
463  	            }
464  	
465  	            this->thesolver->isInfeasible[idx] = this->VIOLATED_AND_CHECKED;
466  	            bestPrices.addIdx(idx);
467  	         }
468  	      }
469  	   }
470  	
471  	   return bestIdx;
472  	}
473  	
474  	/* Entering Simplex
475  	 */
476  	template <class R>
477  	void SPxSteepPR<R>::entered4(SPxId /* id */, int n)
478  	{
479  	   assert(this->thesolver->type() == SPxSolverBase<R>::ENTER);
480  	
481  	   if(n >= 0 && n < this->thesolver->dim())
482  	   {
483  	      R delta = 2 + 1.0 / this->thesolver->basis().iteration();
484  	      R* coWeights_ptr = this->thesolver->coWeights.get_ptr();
485  	      R* weights_ptr = this->thesolver->weights.get_ptr();
486  	      const R* workVec_ptr = workVec.get_const_ptr();
487  	      const R* pVec = this->thesolver->pVec().delta().values();
488  	      const IdxSet& pIdx = this->thesolver->pVec().idx();
489  	      const R* coPvec = this->thesolver->coPvec().delta().values();
490  	      const IdxSet& coPidx = this->thesolver->coPvec().idx();
491  	      R xi_p = 1 / this->thesolver->fVec().delta()[n];
492  	      int i, j;
493  	      R xi_ip;
494  	
495  	      assert(this->thesolver->fVec().delta()[n] > this->thesolver->epsilon()
496  	             || this->thesolver->fVec().delta()[n] < -this->thesolver->epsilon());
497  	
498  	      for(j = coPidx.size() - 1; j >= 0; --j)
499  	      {
500  	         i = coPidx.index(j);
501  	         xi_ip = xi_p * coPvec[i];
502  	         coWeights_ptr[i] += xi_ip * (xi_ip * pi_p - 2.0 * workVec_ptr[i]);
503  	
504  	         /*
505  	           if(coWeights_ptr[i] < 1)
506  	           coWeights_ptr[i] = 1 / (2-x);
507  	         */
508  	         if(coWeights_ptr[i] < delta)
509  	            coWeights_ptr[i] = delta;
510  	         // coWeights_ptr[i] = 1;
511  	         else if(coWeights_ptr[i] > R(infinity))
512  	            coWeights_ptr[i] = 1 / this->thesolver->epsilon();
513  	      }
514  	
515  	      for(j = pIdx.size() - 1; j >= 0; --j)
516  	      {
517  	         i = pIdx.index(j);
518  	         xi_ip = xi_p * pVec[i];
519  	         weights_ptr[i] += xi_ip * (xi_ip * pi_p - 2.0 * (this->thesolver->vector(i) * workVec));
520  	
521  	         /*
522  	           if(weights_ptr[i] < 1)
523  	           weights_ptr[i] = 1 / (2-x);
524  	         */
525  	         if(weights_ptr[i] < delta)
526  	            weights_ptr[i] = delta;
527  	         // weights_ptr[i] = 1;
528  	         else if(weights_ptr[i] > R(infinity))
529  	            weights_ptr[i] = 1.0 / this->thesolver->epsilon();
530  	      }
531  	   }
532  	
533  	   /*@
534  	     if(this->thesolver->isId(id))
535  	     weights[   this->thesolver->number(id) ] *= 1.0001;
536  	     else if(this->thesolver->isCoId(id))
537  	     coWeights[ this->thesolver->number(id) ] *= 1.0001;
538  	   */
539  	
540  	}
541  	
542  	
543  	template <class R>
544  	SPxId SPxSteepPR<R>::buildBestPriceVectorEnterDim(R& best, R feastol)
545  	{
546  	   const R* coTest        = this->thesolver->coTest().get_const_ptr();
547  	   const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr();
548  	   int idx;
549  	   int nsorted;
550  	   R x;
551  	   typename SPxPricer<R>::IdxElement price;
552  	
553  	   prices.clear();
554  	   bestPrices.clear();
555  	
556  	   // construct vector of all prices
557  	   for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i)
558  	   {
559  	      idx = this->thesolver->infeasibilities.index(i);
560  	      x = coTest[idx];
561  	
562  	      if(x < -feastol)
563  	      {
564  	         // it might happen that we call the pricer with a tighter tolerance than what was used when computing the violations
565  	         this->thesolver->isInfeasible[idx] = this->VIOLATED;
566  	         price.val = steeppr::computePrice(x, coWeights_ptr[idx], feastol);
567  	         price.idx = idx;
568  	         prices.append(price);
569  	      }
570  	      else
571  	      {
572  	         this->thesolver->infeasibilities.remove(i);
573  	         this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
574  	      }
575  	   }
576  	
577  	   // set up structures for the quicksort implementation
578  	   this->compare.elements = prices.get_const_ptr();
579  	   // do a partial sort to move the best ones to the front
580  	   // TODO this can be done more efficiently, since we only need the indices
581  	   nsorted = SPxQuicksortPart(prices.get_ptr(), this->compare, 0, prices.size(),
582  	                              SOPLEX_HYPERPRICINGSIZE);
583  	
584  	   // copy indices of best values to bestPrices
585  	   for(int i = 0; i < nsorted; ++i)
586  	   {
587  	      bestPrices.addIdx(prices[i].idx);
588  	      this->thesolver->isInfeasible[prices[i].idx] = this->VIOLATED_AND_CHECKED;
589  	   }
590  	
591  	   if(nsorted > 0)
592  	   {
593  	      best = prices[0].val;
594  	      return this->thesolver->coId(prices[0].idx);
595  	   }
596  	   else
597  	      return SPxId();
598  	}
599  	
600  	
601  	template <class R>
602  	SPxId SPxSteepPR<R>::buildBestPriceVectorEnterCoDim(R& best, R feastol)
603  	{
604  	   const R* test        = this->thesolver->test().get_const_ptr();
605  	   const R* weights_ptr = this->thesolver->weights.get_const_ptr();
606  	   int idx;
607  	   int nsorted;
608  	   R x;
609  	   typename SPxPricer<R>::IdxElement price;
610  	
611  	   pricesCo.clear();
612  	   bestPricesCo.clear();
613  	
614  	   // construct vector of all prices
615  	   for(int i = this->thesolver->infeasibilitiesCo.size() - 1; i >= 0; --i)
616  	   {
617  	      idx = this->thesolver->infeasibilitiesCo.index(i);
618  	      x = test[idx];
619  	
620  	      if(x < -feastol)
621  	      {
622  	         // it might happen that we call the pricer with a tighter tolerance than what was used when computing the violations
623  	         this->thesolver->isInfeasibleCo[idx] = this->VIOLATED;
624  	         price.val = steeppr::computePrice(x, weights_ptr[idx], feastol);
625  	         price.idx = idx;
626  	         pricesCo.append(price);
627  	      }
628  	      else
629  	      {
630  	         this->thesolver->infeasibilitiesCo.remove(i);
631  	         this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED;
632  	      }
633  	   }
634  	
635  	   // set up structures for the quicksort implementation
636  	   this->compare.elements = pricesCo.get_const_ptr();
637  	   // do a partial sort to move the best ones to the front
638  	   // TODO this can be done more efficiently, since we only need the indices
639  	   nsorted = SPxQuicksortPart(pricesCo.get_ptr(), this->compare, 0, pricesCo.size(),
640  	                              SOPLEX_HYPERPRICINGSIZE);
641  	
642  	   // copy indices of best values to bestPrices
643  	   for(int i = 0; i < nsorted; ++i)
644  	   {
645  	      bestPricesCo.addIdx(pricesCo[i].idx);
646  	      this->thesolver->isInfeasibleCo[pricesCo[i].idx] = this->VIOLATED_AND_CHECKED;
647  	   }
648  	
649  	   if(nsorted > 0)
650  	   {
651  	      best = pricesCo[0].val;
652  	      return this->thesolver->id(pricesCo[0].idx);
653  	   }
654  	   else
655  	      return SPxId();
656  	}
657  	
658  	
659  	template <class R>
660  	SPxId SPxSteepPR<R>::selectEnter()
661  	{
662  	   assert(this->thesolver != 0);
663  	   SPxId enterId;
664  	
665  	   enterId = selectEnterX(this->thetolerance);
666  	
667  	   if(!enterId.isValid() && !refined)
668  	   {
669  	      refined = true;
670  	      SPX_MSG_INFO3((*this->thesolver->spxout),
671  	                    (*this->thesolver->spxout) << "WSTEEP05 trying refinement step..\n";)
672  	      enterId = selectEnterX(this->thetolerance / SOPLEX_STEEP_REFINETOL);
673  	   }
674  	
675  	   assert(isConsistent());
676  	
677  	   if(enterId.isValid())
678  	   {
679  	      SSVectorBase<R>& delta = this->thesolver->fVec().delta();
680  	
681  	      this->thesolver->basis().solve4update(delta, this->thesolver->vector(enterId));
682  	
683  	      workRhs.setup_and_assign(delta);
684  	      pi_p = 1 + delta.length2();
685  	
686  	      this->thesolver->setup4coSolve(&workVec, &workRhs);
687  	   }
688  	
689  	   return enterId;
690  	}
691  	
692  	template <class R>
693  	SPxId SPxSteepPR<R>::selectEnterX(R tol)
694  	{
695  	   SPxId enterId;
696  	   SPxId enterCoId;
697  	   R best;
698  	   R bestCo;
699  	
700  	   best = R(-infinity);
701  	   bestCo = R(-infinity);
702  	
703  	   if(this->thesolver->hyperPricingEnter && !refined)
704  	   {
705  	      if(bestPrices.size() < 2 || this->thesolver->basis().lastUpdate() == 0)
706  	         enterCoId = (this->thesolver->sparsePricingEnter) ? buildBestPriceVectorEnterDim(best,
707  	                     tol) : selectEnterDenseDim(best, tol);
708  	      else
709  	         enterCoId = (this->thesolver->sparsePricingEnter) ? selectEnterHyperDim(best,
710  	                     tol) : selectEnterDenseDim(best, tol);
711  	
712  	      if(bestPricesCo.size() < 2 || this->thesolver->basis().lastUpdate() == 0)
713  	         enterId = (this->thesolver->sparsePricingEnterCo) ? buildBestPriceVectorEnterCoDim(bestCo,
714  	                   tol) : selectEnterDenseCoDim(bestCo, tol);
715  	      else
716  	         enterId = (this->thesolver->sparsePricingEnterCo) ? selectEnterHyperCoDim(bestCo,
717  	                   tol) : selectEnterDenseCoDim(bestCo, tol);
718  	   }
719  	   else
720  	   {
721  	      enterCoId = (this->thesolver->sparsePricingEnter
722  	                   && !refined) ? selectEnterSparseDim(best, tol) : selectEnterDenseDim(best, tol);
723  	      enterId = (this->thesolver->sparsePricingEnterCo
724  	                 && !refined) ? selectEnterSparseCoDim(bestCo, tol) : selectEnterDenseCoDim(bestCo, tol);
725  	   }
726  	
727  	   // prefer slack indices to reduce nonzeros in basis matrix
728  	   if(enterCoId.isValid() && (best > SOPLEX_SPARSITY_TRADEOFF * bestCo || !enterId.isValid()))
729  	      return enterCoId;
730  	   else
731  	      return enterId;
732  	}
733  	
734  	
735  	template <class R>
736  	SPxId SPxSteepPR<R>::selectEnterHyperDim(R& best, R tol)
737  	{
738  	   const R* coTest        = this->thesolver->coTest().get_const_ptr();
739  	   const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr();
740  	
741  	   R leastBest = -1;
742  	   R x;
743  	   int enterIdx = -1;
744  	   int idx;
745  	
746  	   // find the best price from short candidate list
747  	   for(int i = bestPrices.size() - 1; i >= 0; --i)
748  	   {
749  	      idx = bestPrices.index(i);
750  	      x = coTest[idx];
751  	
752  	      if(x < -tol)
753  	      {
754  	         x = steeppr::computePrice(x, coWeights_ptr[idx], tol);
755  	
756  	         assert(x >= 0);
757  	
758  	         // update the best price of candidate list
759  	         if(x > best)
760  	         {
761  	            best = x;
762  	            enterIdx = idx;
763  	         }
764  	
765  	         // update the smallest price of candidate list
766  	         if(x < leastBest || leastBest < 0)
767  	            leastBest = x;
768  	      }
769  	      else
770  	      {
771  	         bestPrices.remove(i);
772  	         this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
773  	      }
774  	   }
775  	
776  	   // scan the updated indices for a better price
777  	   for(int i = this->thesolver->updateViols.size() - 1; i >= 0; --i)
778  	   {
779  	      idx = this->thesolver->updateViols.index(i);
780  	
781  	      // only look at indices that were not checked already
782  	      if(this->thesolver->isInfeasible[idx] == this->VIOLATED)
783  	      {
784  	         x = coTest[idx];
785  	
786  	         if(x < -tol)
787  	         {
788  	            x = steeppr::computePrice(x, coWeights_ptr[idx], tol);
789  	
790  	            if(x > leastBest)
791  	            {
792  	               if(x > best)
793  	               {
794  	                  best = x;
795  	                  enterIdx = idx;
796  	               }
797  	
798  	               // put index into candidate list
799  	               this->thesolver->isInfeasible[idx] = this->VIOLATED_AND_CHECKED;
800  	               bestPrices.addIdx(idx);
801  	            }
802  	         }
803  	         else
804  	         {
805  	            this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
806  	         }
807  	      }
808  	   }
809  	
810  	   if(enterIdx >= 0)
811  	      return this->thesolver->coId(enterIdx);
812  	   else
813  	      return SPxId();
814  	}
815  	
816  	
817  	template <class R>
818  	SPxId SPxSteepPR<R>::selectEnterHyperCoDim(R& best, R tol)
819  	{
820  	   const R* test        = this->thesolver->test().get_const_ptr();
821  	   const R* weights_ptr = this->thesolver->weights.get_const_ptr();
822  	
823  	   R leastBest = -1;
824  	   R x;
825  	   int enterIdx = -1;
826  	   int idx;
827  	
828  	   // find the best price from short candidate list
829  	   for(int i = bestPricesCo.size() - 1; i >= 0; --i)
830  	   {
831  	      idx = bestPricesCo.index(i);
832  	      x = test[idx];
833  	
834  	      if(x < -tol)
835  	      {
836  	         x = steeppr::computePrice(x, weights_ptr[idx], tol);
837  	
838  	         assert(x >= 0);
839  	
840  	         // update the best price of candidate list
841  	         if(x > best)
842  	         {
843  	            best = x;
844  	            enterIdx = idx;
845  	         }
846  	
847  	         // update the smallest price of candidate list
848  	         if(x < leastBest || leastBest < 0)
849  	            leastBest = x;
850  	      }
851  	      else
852  	      {
853  	         bestPricesCo.remove(i);
854  	         this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED;
855  	      }
856  	   }
857  	
858  	   // scan the updated indices for a better price
859  	   for(int i = this->thesolver->updateViolsCo.size() - 1; i >= 0; --i)
860  	   {
861  	      idx = this->thesolver->updateViolsCo.index(i);
862  	
863  	      // only look at indices that were not checked already
864  	      if(this->thesolver->isInfeasibleCo[idx] == this->VIOLATED)
865  	      {
866  	         x = test[idx];
867  	
868  	         if(x < -tol)
869  	         {
870  	            x = steeppr::computePrice(x, weights_ptr[idx], tol);
871  	
872  	            if(x > leastBest)
873  	            {
874  	               if(x > best)
875  	               {
876  	                  best = x;
877  	                  enterIdx = idx;
878  	               }
879  	
880  	               // put index into candidate list
881  	               this->thesolver->isInfeasibleCo[idx] = this->VIOLATED_AND_CHECKED;
882  	               bestPricesCo.addIdx(idx);
883  	            }
884  	         }
885  	         else
886  	         {
887  	            this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED;
888  	         }
889  	      }
890  	   }
891  	
892  	   if(enterIdx >= 0)
893  	      return this->thesolver->id(enterIdx);
894  	   else
895  	      return SPxId();
896  	}
897  	
898  	
899  	template <class R>
900  	SPxId SPxSteepPR<R>::selectEnterSparseDim(R& best, R tol)
901  	{
902  	   SPxId enterId;
903  	   const R* coTest        = this->thesolver->coTest().get_const_ptr();
904  	   const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr();
905  	
906  	   int idx;
907  	   R x;
908  	
909  	   for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i)
910  	   {
911  	      idx = this->thesolver->infeasibilities.index(i);
912  	      x = coTest[idx];
913  	
914  	      if(x < -tol)
915  	      {
916  	         x = steeppr::computePrice(x, coWeights_ptr[idx], tol);
917  	
918  	         if(x > best)
919  	         {
920  	            best = x;
921  	            enterId = this->thesolver->coId(idx);
922  	         }
923  	      }
924  	      else
925  	      {
926  	         this->thesolver->infeasibilities.remove(i);
927  	         this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED;
928  	      }
929  	   }
930  	
931  	   return enterId;
932  	}
933  	
934  	template <class R>
935  	SPxId SPxSteepPR<R>::selectEnterSparseCoDim(R& best, R tol)
936  	{
937  	   SPxId enterId;
938  	   const R* test          = this->thesolver->test().get_const_ptr();
939  	   const R* weights_ptr   = this->thesolver->weights.get_const_ptr();
940  	
941  	   int idx;
942  	   R x;
943  	
944  	   for(int i = this->thesolver->infeasibilitiesCo.size() - 1; i >= 0; --i)
945  	   {
946  	      idx = this->thesolver->infeasibilitiesCo.index(i);
947  	      x = test[idx];
948  	
949  	      if(x < -tol)
950  	      {
951  	         x = steeppr::computePrice(x, weights_ptr[idx], tol);
952  	
953  	         if(x > best)
954  	         {
955  	            best   = x;
956  	            enterId = this->thesolver->id(idx);
957  	         }
958  	      }
959  	      else
960  	      {
961  	         this->thesolver->infeasibilitiesCo.remove(i);
962  	         this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED;
963  	      }
964  	   }
965  	
966  	   return enterId;
967  	}
968  	
969  	template <class R>
970  	SPxId SPxSteepPR<R>::selectEnterDenseDim(R& best, R tol)
971  	{
972  	   SPxId enterId;
973  	   const R* coTest        = this->thesolver->coTest().get_const_ptr();
974  	   const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr();
975  	
976  	   R x;
977  	
978  	   for(int i = 0, end = this->thesolver->dim(); i < end; ++i)
979  	   {
980  	      x = coTest[i];
981  	
982  	      if(x < -tol)
983  	      {
984  	         x = steeppr::computePrice(x, coWeights_ptr[i], tol);
985  	
986  	         if(x > best)
987  	         {
988  	            best = x;
989  	            enterId = this->thesolver->coId(i);
990  	         }
991  	      }
992  	   }
993  	
994  	   return enterId;
995  	}
996  	
997  	template <class R>
998  	SPxId SPxSteepPR<R>::selectEnterDenseCoDim(R& best, R tol)
999  	{
1000 	   SPxId enterId;
1001 	   const R* test          = this->thesolver->test().get_const_ptr();
1002 	   const R* weights_ptr   = this->thesolver->weights.get_const_ptr();
1003 	
1004 	   R x;
1005 	
1006 	   for(int i = 0, end = this->thesolver->coDim(); i < end; ++i)
1007 	   {
1008 	      x = test[i];
1009 	
1010 	      if(x < -tol)
1011 	      {
1012 	         x = steeppr::computePrice(x, weights_ptr[i], tol);
1013 	
1014 	         if(x > best)
1015 	         {
1016 	            best   = x;
1017 	            enterId = this->thesolver->id(i);
1018 	         }
1019 	      }
1020 	   }
1021 	
1022 	   return enterId;
1023 	}
1024 	
1025 	
1026 	template <class R>
1027 	void SPxSteepPR<R>::addedVecs(int n)
1028 	{
1029 	   VectorBase<R>& weights = this->thesolver->weights;
1030 	   n = weights.dim();
1031 	   weights.reDim(this->thesolver->coDim());
1032 	
1033 	   if(this->thesolver->type() == SPxSolverBase<R>::ENTER)
1034 	   {
1035 	      for(; n < weights.dim(); ++n)
1036 	         weights[n] = 2;
1037 	   }
1038 	}
1039 	
1040 	template <class R>
1041 	void SPxSteepPR<R>::addedCoVecs(int n)
1042 	{
1043 	   VectorBase<R>& coWeights = this->thesolver->coWeights;
1044 	   n = coWeights.dim();
1045 	   workVec.reDim(this->thesolver->dim());
1046 	   coWeights.reDim(this->thesolver->dim());
1047 	
1048 	   for(; n < coWeights.dim(); ++n)
1049 	      coWeights[n] = 1;
1050 	}
1051 	
1052 	template <class R>
1053 	void SPxSteepPR<R>::removedVec(int i)
1054 	{
1055 	   assert(this->thesolver != 0);
1056 	   VectorBase<R>& weights = this->thesolver->weights;
1057 	   weights[i] = weights[weights.dim()];
1058 	   weights.reDim(this->thesolver->coDim());
1059 	}
1060 	
1061 	template <class R>
1062 	void SPxSteepPR<R>::removedVecs(const int perm[])
1063 	{
1064 	   assert(this->thesolver != 0);
1065 	   VectorBase<R>& weights = this->thesolver->weights;
1066 	
1067 	   if(this->thesolver->type() == SPxSolverBase<R>::ENTER)
1068 	   {
1069 	      int i;
1070 	      int j = weights.dim();
1071 	
1072 	      for(i = 0; i < j; ++i)
1073 	      {
1074 	         if(perm[i] >= 0)
1075 	            weights[perm[i]] = weights[i];
1076 	      }
1077 	   }
1078 	
1079 	   weights.reDim(this->thesolver->coDim());
1080 	}
1081 	
1082 	template <class R>
1083 	void SPxSteepPR<R>::removedCoVec(int i)
1084 	{
1085 	   assert(this->thesolver != 0);
1086 	   VectorBase<R>& coWeights = this->thesolver->coWeights;
1087 	   coWeights[i] = coWeights[coWeights.dim()];
1088 	   coWeights.reDim(this->thesolver->dim());
1089 	}
1090 	
1091 	template <class R>
1092 	void SPxSteepPR<R>::removedCoVecs(const int perm[])
1093 	{
1094 	   assert(this->thesolver != 0);
1095 	   VectorBase<R>& coWeights = this->thesolver->coWeights;
1096 	   int i;
1097 	   int j = coWeights.dim();
1098 	
1099 	   for(i = 0; i < j; ++i)
1100 	   {
1101 	      if(perm[i] >= 0)
1102 	         coWeights[perm[i]] = coWeights[i];
1103 	   }
1104 	
1105 	   coWeights.reDim(this->thesolver->dim());
1106 	}
1107 	
1108 	template <class R>
1109 	bool SPxSteepPR<R>::isConsistent() const
1110 	{
1111 	#ifdef ENABLE_CONSISTENCY_CHECKS
1112 	   VectorBase<R>& w = this->thesolver->weights;
1113 	   VectorBase<R>& coW = this->thesolver->coWeights;
1114 	
1115 	   if(this->thesolver != 0 && this->thesolver->type() == SPxSolverBase<R>::LEAVE && setup == EXACT)
1116 	   {
1117 	      int i;
1118 	      SSVectorBase<R>  tmp(this->thesolver->dim(), this->thesolver->tolerances());
1119 	      R x;
1120 	
1121 	      for(i = this->thesolver->dim() - 1; i >= 0; --i)
1122 	      {
1123 	         this->thesolver->basis().coSolve(tmp, this->thesolver->unitVector(i));
1124 	         x = coW[i] - tmp.length2();
1125 	
1126 	         if(x > this->thesolver->leavetol() || -x > this->thesolver->leavetol())
1127 	         {
1128 	            SPX_MSG_ERROR(std::cerr << "ESTEEP03 x[" << i << "] = " << x << std::endl;)
1129 	         }
1130 	      }
1131 	   }
1132 	
1133 	   if(this->thesolver != 0 && this->thesolver->type() == SPxSolverBase<R>::ENTER)
1134 	   {
1135 	      int i;
1136 	
1137 	      for(i = this->thesolver->dim() - 1; i >= 0; --i)
1138 	      {
1139 	         if(coW[i] < this->thesolver->epsilon())
1140 	            return SPX_MSG_INCONSISTENT("SPxSteepPR");
1141 	      }
1142 	
1143 	      for(i = this->thesolver->coDim() - 1; i >= 0; --i)
1144 	      {
1145 	         if(w[i] < this->thesolver->epsilon())
1146 	            return SPX_MSG_INCONSISTENT("SPxSteepPR");
1147 	      }
1148 	   }
1149 	
1150 	#endif
1151 	
1152 	   return true;
1153 	}
1154 	} // namespace soplex
1155