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 <iostream>
26   	
27   	#include "soplex/array.h"
28   	#include "soplex/dataarray.h"
29   	#include "soplex/sorter.h"
30   	#include "soplex/spxout.h"
31   	#include <sstream>
32   	#include <iostream>
33   	#include <fstream>
34   	#include <memory>
35   	
36   	
37   	//rows
38   	#define SOPLEX_FREE_LHS_RHS            1
39   	#define SOPLEX_FREE_CONSTRAINT         1
40   	#define SOPLEX_EMPTY_CONSTRAINT        1
41   	#define SOPLEX_ROW_SINGLETON           1
42   	#define SOPLEX_AGGREGATE_VARS          1
43   	#define SOPLEX_FORCE_CONSTRAINT        1
44   	//cols
45   	#define SOPLEX_FREE_BOUNDS             1
46   	#define SOPLEX_COMPTY_COLUMN           1
47   	#define SOPLEX_FIX_VARIABLE            1
48   	#define SOPLEX_FREE_ZERO_OBJ_VARIABLE  1
49   	#define SOPLEX_ZERO_OBJ_COL_SINGLETON  1
50   	#define SOPLEX_DOUBLETON_EQUATION      1
51   	#define SOPLEX_FREE_COL_SINGLETON      1
52   	//dual
53   	#define SOPLEX_DOMINATED_COLUMN        1
54   	#define SOPLEX_WEAKLY_DOMINATED_COLUMN 1
55   	#define SOPLEX_MULTI_AGGREGATE         1
56   	//other
57   	#define SOPLEX_TRIVIAL_HEURISTICS      1
58   	#define SOPLEX_PSEUDOOBJ               1
59   	
60   	
61   	#define SOPLEX_EXTREMES                1
62   	#define SOPLEX_ROWS_SPXMAINSM          1
63   	#define SOPLEX_COLS_SPXMAINSM          1
64   	#define SOPLEX_DUAL_SPXMAINSM          1
65   	///@todo check: with this simplification step, the unsimplified basis seems to be slightly suboptimal for some instances
66   	#define SOPLEX_DUPLICATE_ROWS          1
67   	#define SOPLEX_DUPLICATE_COLS          1
68   	
69   	
70   	#ifndef NDEBUG
71   	#define SOPLEX_CHECK_BASIS_DIM
72   	#endif  // NDEBUG
73   	
74   	namespace soplex
75   	{
76   	
77   	template <class R>
78   	bool SPxMainSM<R>::PostStep::checkBasisDim(DataArray<typename SPxSolverBase<R>::VarStatus> rows,
79   	      DataArray<typename SPxSolverBase<R>::VarStatus> cols) const
80   	{
81   	   int numBasis = 0;
82   	
83   	   for(int rs = 0; rs < nRows; ++rs)
84   	   {
85   	      if(rows[rs] == SPxSolverBase<R>::BASIC)
86   	         numBasis++;
87   	   }
88   	
89   	   for(int cs = 0; cs < nCols; ++cs)
90   	   {
91   	      if(cols[cs] == SPxSolverBase<R>::BASIC)
92   	         numBasis++;
93   	   }
94   	
95   	   assert(numBasis == nRows);
96   	   return numBasis == nRows;
97   	}
98   	
99   	template <class R>
100  	void SPxMainSM<R>::RowObjPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
101  	                                     VectorBase<R>&,
102  	                                     DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
103  	                                     DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
104  	{
105  	   s[m_i] = s[m_i] - x[m_j];
106  	
107  	   assert(rStatus[m_i] != SPxSolverBase<R>::UNDEFINED);
108  	   assert(cStatus[m_j] != SPxSolverBase<R>::UNDEFINED);
109  	   assert(rStatus[m_i] != SPxSolverBase<R>::BASIC || cStatus[m_j] != SPxSolverBase<R>::BASIC);
110  	
111  	   SPxOut::debug(this, "RowObjPS: removing slack column {} ({}) for row {} ({}).\n", m_j, cStatus[m_j],
112  	                 m_i, rStatus[m_i]);
113  	
114  	   if(rStatus[m_i] != SPxSolverBase<R>::BASIC)
115  	   {
116  	      switch(cStatus[m_j])
117  	      {
118  	      case SPxSolverBase<R>::ON_UPPER:
119  	         rStatus[m_i] = SPxSolverBase<R>::ON_LOWER;
120  	         break;
121  	
122  	      case SPxSolverBase<R>::ON_LOWER:
123  	         rStatus[m_i] = SPxSolverBase<R>::ON_UPPER;
124  	         break;
125  	
126  	      default:
127  	         rStatus[m_i] = cStatus[m_j];
128  	      }
129  	
130  	      // otherwise checkBasisDim() may fail
131  	      cStatus[m_j] = SPxSolverBase<R>::ZERO;
132  	   }
133  	
134  	#ifdef SOPLEX_CHECK_BASIS_DIM
135  	
136  	   if(!this->checkBasisDim(rStatus, cStatus))
137  	   {
138  	      assert(false);
139  	      throw SPxInternalCodeException("XMAISM15 Dimension doesn't match after this step.");
140  	   }
141  	
142  	#endif
143  	}
144  	
145  	template <class R>
146  	void SPxMainSM<R>::FreeConstraintPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
147  	      VectorBase<R>&,
148  	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
149  	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
150  	{
151  	   // correcting the change of idx by deletion of the row:
152  	   if(m_i != m_old_i)
153  	   {
154  	      s[m_old_i] = s[m_i];
155  	      y[m_old_i] = y[m_i];
156  	      rStatus[m_old_i] = rStatus[m_i];
157  	   }
158  	
159  	   // primal:
160  	   R slack = 0.0;
161  	
162  	   for(int k = 0; k < m_row.size(); ++k)
163  	      slack += m_row.value(k) * x[m_row.index(k)];
164  	
165  	   s[m_i] = slack;
166  	
167  	   // dual:
168  	   y[m_i] = m_row_obj;
169  	
170  	   // basis:
171  	   rStatus[m_i] = SPxSolverBase<R>::BASIC;
172  	
173  	#ifdef SOPLEX_CHECK_BASIS_DIM
174  	
175  	   if(!this->checkBasisDim(rStatus, cStatus))
176  	   {
177  	      throw SPxInternalCodeException("XMAISM15 Dimension doesn't match after this step.");
178  	   }
179  	
180  	#endif
181  	}
182  	
183  	template <class R>
184  	void SPxMainSM<R>::EmptyConstraintPS::execute(VectorBase<R>&, VectorBase<R>& y, VectorBase<R>& s,
185  	      VectorBase<R>&,
186  	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
187  	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
188  	{
189  	   // correcting the change of idx by deletion of the row:
190  	   if(m_i != m_old_i)
191  	   {
192  	      s[m_old_i] = s[m_i];
193  	      y[m_old_i] = y[m_i];
194  	      rStatus[m_old_i] = rStatus[m_i];
195  	   }
196  	
197  	   // primal:
198  	   s[m_i] = 0.0;
199  	
200  	   // dual:
201  	   y[m_i] = m_row_obj;
202  	
203  	   // basis:
204  	   rStatus[m_i] = SPxSolverBase<R>::BASIC;
205  	
206  	#ifdef SOPLEX_CHECK_BASIS_DIM
207  	
208  	   if(!this->checkBasisDim(rStatus, cStatus))
209  	   {
210  	      throw SPxInternalCodeException("XMAISM16 Dimension doesn't match after this step.");
211  	   }
212  	
213  	#endif
214  	}
215  	
216  	template <class R>
217  	void SPxMainSM<R>::RowSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
218  	      VectorBase<R>& r,
219  	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
220  	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
221  	{
222  	   // correcting the change of idx by deletion of the row:
223  	   if(m_i != m_old_i)
224  	   {
225  	      y[m_old_i] = y[m_i];
226  	      s[m_old_i] = s[m_i];
227  	      rStatus[m_old_i] = rStatus[m_i];
228  	   }
229  	
230  	   R aij = m_col[m_i];
231  	   assert(aij != 0.0);
232  	
233  	   // primal:
234  	   s[m_i] = aij * x[m_j];
235  	
236  	   // dual & basis:
237  	   R val = m_obj;
238  	
239  	   for(int k = 0; k < m_col.size(); ++k)
240  	   {
241  	      if(m_col.index(k) != m_i)
242  	         val -= m_col.value(k) * y[m_col.index(k)];
243  	   }
244  	
245  	   R newLo = (aij > 0) ? m_lhs / aij : m_rhs / aij; // implicit lhs
246  	   R newUp = (aij > 0) ? m_rhs / aij : m_lhs / aij; // implicit rhs
247  	
248  	   switch(cStatus[m_j])
249  	   {
250  	   case SPxSolverBase<R>::FIXED:
251  	      if(newLo <= m_oldLo && newUp >= m_oldUp)
252  	      {
253  	         // this row is totally redundant, has not changed bound of xj
254  	         rStatus[m_i] = SPxSolverBase<R>::BASIC;
255  	         y[m_i] = m_row_obj;
256  	      }
257  	      else if(EQrel(newLo, newUp, this->feastol()))
258  	      {
259  	         // row is in the type  aij * xj = b
260  	         assert(EQrel(newLo, x[m_j], this->feastol()));
261  	
262  	         if(EQrel(m_oldLo, m_oldUp, this->feastol()))
263  	         {
264  	            // xj has been fixed in other row
265  	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
266  	            y[m_i] = m_row_obj;
267  	         }
268  	         else if((EQrel(m_oldLo, x[m_j], this->feastol()) && r[m_j] <= -this->feastol())
269  	                 || (EQrel(m_oldUp, x[m_j], this->feastol()) && r[m_j] >= this->feastol())
270  	                 || (!EQrel(m_oldLo, x[m_j], this->feastol()) && !(EQrel(m_oldUp, x[m_j], this->feastol()))))
271  	         {
272  	            // if x_j on lower but reduced cost is negative, or x_j on upper but reduced cost is positive, or x_j not on bound: basic
273  	            rStatus[m_i] = (EQrel(m_lhs, x[m_j] * aij,
274  	                                  this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
275  	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
276  	            y[m_i] = val / aij;
277  	            r[m_j] = 0.0;
278  	         }
279  	         else
280  	         {
281  	            // set x_j on one of the bound
282  	            assert(EQrel(m_oldLo, x[m_j], this->feastol()) || EQrel(m_oldUp, x[m_j], this->feastol()));
283  	
284  	            cStatus[m_j] = EQrel(m_oldLo, x[m_j],
285  	                                 this->feastol()) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
286  	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
287  	            y[m_i] = m_row_obj;
288  	            r[m_j] = val;
289  	         }
290  	      }
291  	      else if(EQrel(newLo, m_oldUp, this->feastol()))
292  	      {
293  	         // row is in the type  xj >= b/aij, try to set xj on upper
294  	         if(r[m_j] >= this->feastol())
295  	         {
296  	            // the reduced cost is positive, xj should in the basic
297  	            assert(EQrel(m_rhs / aij, x[m_j], this->feastol()) || EQrel(m_lhs / aij, x[m_j], this->feastol()));
298  	
299  	            rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j],
300  	                                  this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
301  	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
302  	            y[m_i] = val / aij;
303  	            r[m_j] = 0.0;
304  	         }
305  	         else
306  	         {
307  	            assert(EQrel(m_oldUp, x[m_j], this->feastol()));
308  	
309  	            cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
310  	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
311  	            y[m_i] = m_row_obj;
312  	            r[m_j] = val;
313  	         }
314  	      }
315  	      else if(EQrel(newUp, m_oldLo, this->feastol()))
316  	      {
317  	         // row is in the type  xj <= b/aij, try to set xj on lower
318  	         if(r[m_j] <= -this->feastol())
319  	         {
320  	            // the reduced cost is negative, xj should in the basic
321  	            assert(EQrel(m_rhs / aij, x[m_j], this->feastol()) || EQrel(m_lhs / aij, x[m_j], this->feastol()));
322  	
323  	            rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j],
324  	                                  this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
325  	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
326  	            y[m_i] = val / aij;
327  	            r[m_j] = 0.0;
328  	         }
329  	         else
330  	         {
331  	            assert(EQrel(m_oldLo, x[m_j], this->feastol()));
332  	
333  	            cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
334  	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
335  	            y[m_i] = m_row_obj;
336  	            r[m_j] = val;
337  	         }
338  	      }
339  	      else
340  	      {
341  	         // the variable is set to FIXED by other constraints, i.e., this singleton row is redundant
342  	         rStatus[m_i] = SPxSolverBase<R>::BASIC;
343  	         y[m_i] = m_row_obj;
344  	      }
345  	
346  	      break;
347  	
348  	   case SPxSolverBase<R>::BASIC:
349  	      rStatus[m_i] = SPxSolverBase<R>::BASIC;
350  	      y[m_i] = m_row_obj;
351  	      r[m_j] = 0.0;
352  	      break;
353  	
354  	   case SPxSolverBase<R>::ON_LOWER:
355  	      if(EQrel(m_oldLo, x[m_j], this->feastol())) // xj may stay on lower
356  	      {
357  	         rStatus[m_i] = SPxSolverBase<R>::BASIC;
358  	         y[m_i] = m_row_obj;
359  	         r[m_j] = val;
360  	      }
361  	      else // if reduced costs are negative or old lower bound not equal to xj, we need to change xj into the basis
362  	      {
363  	         assert(!isOptimal || EQrel(m_rhs / aij, x[m_j], this->feastol())
364  	                || EQrel(m_lhs / aij, x[m_j], this->feastol()));
365  	
366  	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
367  	         rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j],
368  	                               this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
369  	         y[m_i] = val / aij;
370  	         r[m_j] = 0.0;
371  	      }
372  	
373  	      break;
374  	
375  	   case SPxSolverBase<R>::ON_UPPER:
376  	      if(EQrel(m_oldUp, x[m_j], this->feastol())) // xj may stay on upper
377  	      {
378  	         rStatus[m_i] = SPxSolverBase<R>::BASIC;
379  	         y[m_i] = m_row_obj;
380  	         r[m_j] = val;
381  	      }
382  	      else // if reduced costs are positive or old upper bound not equal to xj, we need to change xj into the basis
383  	      {
384  	         assert(!isOptimal || EQrel(m_rhs / aij, x[m_j], this->feastol())
385  	                || EQrel(m_lhs / aij, x[m_j], this->feastol()));
386  	
387  	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
388  	         rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j],
389  	                               this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
390  	         y[m_i] = val / aij;
391  	         r[m_j] = 0.0;
392  	      }
393  	
394  	      break;
395  	
396  	   case SPxSolverBase<R>::ZERO:
397  	      rStatus[m_i] = SPxSolverBase<R>::BASIC;
398  	      y[m_i] = m_row_obj;
399  	      r[m_j] = val;
400  	      break;
401  	
402  	   default:
403  	      break;
404  	   }
405  	
406  	#ifdef SOPLEX_CHECK_BASIS_DIM
407  	
408  	   if(!this->checkBasisDim(rStatus, cStatus))
409  	   {
410  	      throw SPxInternalCodeException("XMAISM17 Dimension doesn't match after this step.");
411  	   }
412  	
413  	#endif
414  	}
415  	
416  	template <class R>
417  	void SPxMainSM<R>::ForceConstraintPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
418  	      VectorBase<R>& r,
419  	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
420  	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
421  	{
422  	   // correcting the change of idx by deletion of the row:
423  	   if(m_i != m_old_i)
424  	   {
425  	      s[m_old_i] = s[m_i];
426  	      y[m_old_i] = y[m_i];
427  	      rStatus[m_old_i] = rStatus[m_i];
428  	   }
429  	
430  	   // primal:
431  	   s[m_i] = m_lRhs;
432  	
433  	   // basis:
434  	   int cBasisCandidate = -1;
435  	   R maxViolation = -1.0;
436  	   int bas_k = -1;
437  	
438  	   for(int k = 0; k < m_row.size(); ++k)
439  	   {
440  	      int  cIdx  = m_row.index(k);
441  	      R aij   = m_row.value(k);
442  	      R oldLo = m_oldLowers[k];
443  	      R oldUp = m_oldUppers[k];
444  	
445  	      switch(cStatus[cIdx])
446  	      {
447  	      case SPxSolverBase<R>::FIXED:
448  	         if(m_fixed[k])
449  	         {
450  	            assert(EQrel(oldLo, x[cIdx], this->feastol()) || EQrel(oldUp, x[cIdx], this->feastol()));
451  	
452  	            R violation = spxAbs(r[cIdx] / aij);
453  	
454  	            cStatus[cIdx] = EQrel(oldLo, x[cIdx],
455  	                                  this->feastol()) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
456  	
457  	            if(violation > maxViolation && ((EQrel(oldLo, x[cIdx], this->feastol())
458  	                                             && r[cIdx] < -this->feastol())
459  	                                            || (EQrel(oldUp, x[cIdx], this->feastol()) && r[cIdx] > this->feastol())))
460  	            {
461  	               maxViolation = violation;
462  	               cBasisCandidate = cIdx;
463  	               bas_k = k;
464  	            }
465  	         } // do nothing, if the old bounds are equal, i.e. variable has been not fixed in this row
466  	
467  	         break;
468  	
469  	      case SPxSolverBase<R>::ON_LOWER:
470  	      case SPxSolverBase<R>::ON_UPPER:
471  	      case SPxSolverBase<R>::BASIC:
472  	         break;
473  	
474  	      default:
475  	         break;
476  	      }
477  	   }
478  	
479  	   // dual and basis :
480  	   if(cBasisCandidate >= 0)  // one of the variable in the row should in the basis
481  	   {
482  	      assert(EQrel(m_lRhs, m_rhs, this->feastol()) || EQrel(m_lRhs, m_lhs, this->feastol()));
483  	      assert(bas_k >= 0);
484  	      assert(cBasisCandidate == m_row.index(bas_k));
485  	
486  	      cStatus[cBasisCandidate] = SPxSolverBase<R>::BASIC;
487  	      rStatus[m_i] = (EQrel(m_lRhs, m_lhs,
488  	                            this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
489  	
490  	      R aij = m_row.value(bas_k);
491  	      R multiplier = r[cBasisCandidate] / aij;
492  	      r[cBasisCandidate] = 0.0;
493  	
494  	      for(int k = 0; k < m_row.size(); ++k)  // update the reduced cost
495  	      {
496  	         if(k == bas_k)
497  	         {
498  	            continue;
499  	         }
500  	
501  	         r[m_row.index(k)] -= m_row.value(k) * multiplier;
502  	      }
503  	
504  	      // compute the value of new dual variable (because we have a new row)
505  	      R val = m_objs[bas_k];
506  	      DSVectorBase<R> basis_col = m_cols[bas_k];
507  	
508  	      for(int k = 0; k < basis_col.size(); ++k)
509  	      {
510  	         if(basis_col.index(k) != m_i)
511  	            val -= basis_col.value(k) * y[basis_col.index(k)];
512  	      }
513  	
514  	      y[m_i] = val / aij;
515  	   }
516  	   else // slack in the basis
517  	   {
518  	      rStatus[m_i] = SPxSolverBase<R>::BASIC;
519  	      y[m_i] = m_rowobj;
520  	   }
521  	
522  	#ifdef SOPLEX_CHECK_BASIS_DIM
523  	
524  	   if(!this->checkBasisDim(rStatus, cStatus))
525  	   {
526  	      throw SPxInternalCodeException("XMAISM18 Dimension doesn't match after this step.");
527  	   }
528  	
529  	#endif
530  	}
531  	
532  	template <class R>
533  	void SPxMainSM<R>::FixVariablePS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
534  	      VectorBase<R>& r,
535  	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
536  	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
537  	{
538  	   // update the index mapping; if m_correctIdx is false, we assume that this has happened already
539  	   if(m_correctIdx)
540  	   {
541  	      x[m_old_j] = x[m_j];
542  	      r[m_old_j] = r[m_j];
543  	      cStatus[m_old_j] = cStatus[m_j];
544  	   }
545  	
546  	   // primal:
547  	   x[m_j] = m_val;
548  	
549  	   for(int k = 0; k < m_col.size(); ++k)
550  	      s[m_col.index(k)] += m_col.value(k) * x[m_j];
551  	
552  	   // dual:
553  	   R val = m_obj;
554  	
555  	   for(int k = 0; k < m_col.size(); ++k)
556  	      val -= m_col.value(k) * y[m_col.index(k)];
557  	
558  	   r[m_j] = val;
559  	
560  	   // basis:
561  	   if(m_lower == m_upper)
562  	   {
563  	      assert(EQrel(m_lower, m_val, this->epsilon()));
564  	
565  	      cStatus[m_j] = SPxSolverBase<R>::FIXED;
566  	   }
567  	   else
568  	   {
569  	      assert(EQrel(m_val, m_lower, this->epsilon()) || EQrel(m_val, m_upper, this->epsilon())
570  	             || m_val == 0.0);
571  	
572  	      cStatus[m_j] = EQrel(m_val, m_lower, this->epsilon()) ? SPxSolverBase<R>::ON_LOWER : (EQrel(m_val,
573  	                     m_upper, this->epsilon()) ? SPxSolverBase<R>::ON_UPPER : SPxSolverBase<R>::ZERO);
574  	   }
575  	
576  	#ifdef SOPLEX_CHECK_BASIS_DIM
577  	
578  	   if(m_correctIdx)
579  	   {
580  	      if(!this->checkBasisDim(rStatus, cStatus))
581  	      {
582  	         throw SPxInternalCodeException("XMAISM19 Dimension doesn't match after this step.");
583  	      }
584  	   }
585  	
586  	#endif
587  	}
588  	
589  	template <class R>
590  	void SPxMainSM<R>::FixBoundsPS::execute(VectorBase<R>&, VectorBase<R>&, VectorBase<R>&,
591  	                                        VectorBase<R>&,
592  	                                        DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
593  	                                        DataArray<typename SPxSolverBase<R>::VarStatus>&, bool isOptimal) const
594  	{
595  	   // basis:
596  	   cStatus[m_j] = m_status;
597  	}
598  	
599  	template <class R>
600  	void SPxMainSM<R>::FreeZeroObjVariablePS::execute(VectorBase<R>& x, VectorBase<R>& y,
601  	      VectorBase<R>& s, VectorBase<R>& r,
602  	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
603  	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
604  	{
605  	   // correcting the change of idx by deletion of the column and corresponding rows:
606  	   if(m_j != m_old_j)
607  	   {
608  	      x[m_old_j] = x[m_j];
609  	      r[m_old_j] = r[m_j];
610  	      cStatus[m_old_j] = cStatus[m_j];
611  	   }
612  	
613  	   int rIdx = m_old_i - m_col.size() + 1;
614  	
615  	   for(int k = 0; k < m_col.size(); ++k)
616  	   {
617  	      int rIdx_new = m_col.index(k);
618  	      s[rIdx] = s[rIdx_new];
619  	      y[rIdx] = y[rIdx_new];
620  	      rStatus[rIdx] = rStatus[rIdx_new];
621  	      rIdx++;
622  	   }
623  	
624  	   // primal:
625  	   int      domIdx = -1;
626  	   DSVectorBase<R> slack(m_col.size());
627  	
628  	   if(m_loFree)
629  	   {
630  	      R minRowUp = R(infinity);
631  	
632  	      for(int k = 0; k < m_rows.size(); ++k)
633  	      {
634  	         R           val = 0.0;
635  	         const SVectorBase<R>& row = m_rows[k];
636  	
637  	         for(int l = 0; l < row.size(); ++l)
638  	         {
639  	            if(row.index(l) != m_j)
640  	               val += row.value(l) * x[row.index(l)];
641  	         }
642  	
643  	         R scale = maxAbs(m_lRhs[k], val);
644  	
645  	         if(scale < 1.0)
646  	            scale = 1.0;
647  	
648  	         R z = (m_lRhs[k] / scale) - (val / scale);
649  	
650  	         if(isZero(z, this->epsilon()))
651  	            z = 0.0;
652  	
653  	         R up = z * scale / row[m_j];
654  	         slack.add(k, val);
655  	
656  	         if(up < minRowUp)
657  	         {
658  	            minRowUp = up;
659  	            domIdx   = k;
660  	         }
661  	      }
662  	
663  	      if(m_bnd < minRowUp)
664  	      {
665  	         x[m_j] = m_bnd;
666  	         domIdx = -1;
667  	      }
668  	      else
669  	         x[m_j] = minRowUp;
670  	   }
671  	   else
672  	   {
673  	      R maxRowLo = R(-infinity);
674  	
675  	      for(int k = 0; k < m_rows.size(); ++k)
676  	      {
677  	         R val = 0.0;
678  	         const SVectorBase<R>& row = m_rows[k];
679  	
680  	         for(int l = 0; l < row.size(); ++l)
681  	         {
682  	            if(row.index(l) != m_j)
683  	               val += row.value(l) * x[row.index(l)];
684  	         }
685  	
686  	         R scale = maxAbs(m_lRhs[k], val);
687  	
688  	         if(scale < 1.0)
689  	            scale = 1.0;
690  	
691  	         R z = (m_lRhs[k] / scale) - (val / scale);
692  	
693  	         if(isZero(z, this->epsilon()))
694  	            z = 0.0;
695  	
696  	         R lo = z * scale / row[m_j];
697  	         slack.add(k, val);
698  	
699  	         if(lo > maxRowLo)
700  	         {
701  	            maxRowLo = lo;
702  	            domIdx   = k;
703  	         }
704  	      }
705  	
706  	      if(m_bnd > maxRowLo)
707  	      {
708  	         x[m_j] = m_bnd;
709  	         domIdx = -1;
710  	      }
711  	      else
712  	         x[m_j] = maxRowLo;
713  	   }
714  	
715  	   for(int k = 0; k < m_col.size(); ++k)
716  	      s[m_col.index(k)] = slack[k] + m_col.value(k) * x[m_j];
717  	
718  	   // dual:
719  	   r[m_j] = 0.0;
720  	
721  	   for(int k = 0; k < m_col.size(); ++k)
722  	   {
723  	      int idx = m_col.index(k);
724  	      y[idx] = m_rowObj[idx];
725  	   }
726  	
727  	   // basis:
728  	   for(int k = 0; k < m_col.size(); ++k)
729  	   {
730  	      if(k != domIdx)
731  	         rStatus[m_col.index(k)] = SPxSolverBase<R>::BASIC;
732  	
733  	      else
734  	      {
735  	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
736  	
737  	         if(m_loFree)
738  	            rStatus[m_col.index(k)] = (m_col.value(k) > 0) ? SPxSolverBase<R>::ON_UPPER :
739  	                                      SPxSolverBase<R>::ON_LOWER;
740  	         else
741  	            rStatus[m_col.index(k)] = (m_col.value(k) > 0) ? SPxSolverBase<R>::ON_LOWER :
742  	                                      SPxSolverBase<R>::ON_UPPER;
743  	      }
744  	   }
745  	
746  	   if(domIdx == -1)
747  	   {
748  	      if(m_loFree)
749  	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
750  	      else
751  	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
752  	   }
753  	
754  	#ifdef SOPLEX_CHECK_BASIS_DIM
755  	
756  	   if(!this->checkBasisDim(rStatus, cStatus))
757  	   {
758  	      throw SPxInternalCodeException("XMAISM20 Dimension doesn't match after this step.");
759  	   }
760  	
761  	#endif
762  	}
763  	
764  	template <class R>
765  	void SPxMainSM<R>::ZeroObjColSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y,
766  	      VectorBase<R>& s, VectorBase<R>& r,
767  	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
768  	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
769  	{
770  	   // correcting the change of idx by deletion of the column and corresponding rows:
771  	   if(m_j != m_old_j)
772  	   {
773  	      x[m_old_j] = x[m_j];
774  	      r[m_old_j] = r[m_j];
775  	      cStatus[m_old_j] = cStatus[m_j];
776  	   }
777  	
778  	   // primal & basis:
779  	   R aij = m_row[m_j];
780  	
781  	   if(isZero(s[m_i], R(1e-6)))
782  	      s[m_i] = 0.0;
783  	   else if(s[m_i] >= R(infinity))
784  	      // this is a fix for a highly ill conditioned instance that is "solved" in presolving (ilaser0 from MINLP, mittelmann)
785  	      throw SPxException("Simplifier: infinite activities - aborting unsimplification");
786  	
787  	   R scale1 = maxAbs(m_lhs, s[m_i]);
788  	   R scale2 = maxAbs(m_rhs, s[m_i]);
789  	
790  	   if(scale1 < 1.0)
791  	      scale1 = 1.0;
792  	
793  	   if(scale2 < 1.0)
794  	      scale2 = 1.0;
795  	
796  	   R z1 = (m_lhs / scale1) - (s[m_i] / scale1);
797  	   R z2 = (m_rhs / scale2) - (s[m_i] / scale2);
798  	
799  	   if(isZero(z1, this->epsilon()))
800  	      z1 = 0.0;
801  	
802  	   if(isZero(z2, this->epsilon()))
803  	      z2 = 0.0;
804  	
805  	   R lo = (aij > 0) ? z1 * scale1 / aij : z2 * scale2 / aij;
806  	   R up = (aij > 0) ? z2 * scale2 / aij : z1 * scale1 / aij;
807  	
808  	   if(isZero(lo, this->feastol()))
809  	      lo = 0.0;
810  	
811  	   if(isZero(up, this->feastol()))
812  	      up = 0.0;
813  	
814  	   assert(LErel(lo, up, this->epsilon()));
815  	   SOPLEX_ASSERT_WARN("WMAISM01", isNotZero(aij, R(1.0 / R(infinity))));
816  	
817  	   if(rStatus[m_i] == SPxSolverBase<R>::ON_LOWER)
818  	   {
819  	      if(m_lower <= R(-infinity) && m_upper >= R(infinity))
820  	      {
821  	         x[m_j] = 0.0;
822  	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
823  	      }
824  	      else if(m_lower == m_upper)
825  	      {
826  	         x[m_j]       = m_lower;
827  	         cStatus[m_j] = SPxSolverBase<R>::FIXED;
828  	      }
829  	      else if(aij > 0)
830  	      {
831  	         x[m_j]       = m_upper;
832  	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
833  	      }
834  	      else if(aij < 0)
835  	      {
836  	         x[m_j]       = m_lower;
837  	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
838  	      }
839  	      else
840  	         throw SPxInternalCodeException("XMAISM01 This should never happen.");
841  	   }
842  	   else if(rStatus[m_i] == SPxSolverBase<R>::ON_UPPER)
843  	   {
844  	      if(m_lower <= R(-infinity) && m_upper >= R(infinity))
845  	      {
846  	         x[m_j] = 0.0;
847  	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
848  	      }
849  	      else if(m_lower == m_upper)
850  	      {
851  	         x[m_j]       = m_lower;
852  	         cStatus[m_j] = SPxSolverBase<R>::FIXED;
853  	      }
854  	      else if(aij > 0)
855  	      {
856  	         x[m_j]       = m_lower;
857  	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
858  	      }
859  	      else if(aij < 0)
860  	      {
861  	         x[m_j]       = m_upper;
862  	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
863  	      }
864  	      else
865  	         throw SPxInternalCodeException("XMAISM02 This should never happen.");
866  	   }
867  	   else if(rStatus[m_i] == SPxSolverBase<R>::FIXED)
868  	   {
869  	      if(m_lower <= R(-infinity) && m_upper >= R(infinity))
870  	      {
871  	         x[m_j] = 0.0;
872  	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
873  	      }
874  	      else
875  	      {
876  	         assert(EQrel(m_lower, m_upper, this->feastol()));
877  	
878  	         x[m_j]        = (m_lower + m_upper) / 2.0;
879  	         cStatus[m_j]  = SPxSolverBase<R>::FIXED;
880  	      }
881  	   }
882  	   else if(rStatus[m_i] == SPxSolverBase<R>::BASIC)
883  	   {
884  	      if(GErel(m_lower, lo, this->feastol()) && m_lower > R(-infinity))
885  	      {
886  	         x[m_j]       = m_lower;
887  	         cStatus[m_j] = (m_lower == m_upper) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
888  	      }
889  	      else if(LErel(m_upper, up, this->feastol()) && m_upper < R(infinity))
890  	      {
891  	         x[m_j]       = m_upper;
892  	         cStatus[m_j] = (m_lower == m_upper) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
893  	      }
894  	      else if(lo > R(-infinity))
895  	      {
896  	         // make m_i non-basic and m_j basic
897  	         x[m_j]       = lo;
898  	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
899  	         rStatus[m_i] = (aij > 0 ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER);
900  	      }
901  	      else if(up < R(infinity))
902  	      {
903  	         // make m_i non-basic and m_j basic
904  	         x[m_j]       = up;
905  	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
906  	         rStatus[m_i] = (aij > 0 ? SPxSolverBase<R>::ON_UPPER : SPxSolverBase<R>::ON_LOWER);
907  	      }
908  	      else
909  	         throw SPxInternalCodeException("XMAISM03 This should never happen.");
910  	   }
911  	   else
912  	      throw SPxInternalCodeException("XMAISM04 This should never happen.");
913  	
914  	   s[m_i] += aij * x[m_j];
915  	
916  	   // dual:
917  	   r[m_j] = -1.0 * aij * y[m_i];
918  	
919  	   assert(!isOptimal || (cStatus[m_j] != SPxSolverBase<R>::BASIC || isZero(r[m_j], this->feastol())));
920  	
921  	#ifdef SOPLEX_CHECK_BASIS_DIM
922  	
923  	   if(!this->checkBasisDim(rStatus, cStatus))
924  	   {
925  	      throw SPxInternalCodeException("XMAISM21 Dimension doesn't match after this step.");
926  	   }
927  	
928  	#endif
929  	}
930  	
931  	template <class R>
932  	void SPxMainSM<R>::FreeColSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
933  	      VectorBase<R>& r,
934  	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
935  	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
936  	{
937  	
938  	   // correcting the change of idx by deletion of the row:
939  	   if(m_i != m_old_i)
940  	   {
941  	      s[m_old_i] = s[m_i];
942  	      y[m_old_i] = y[m_i];
943  	      rStatus[m_old_i] = rStatus[m_i];
944  	   }
945  	
946  	   // correcting the change of idx by deletion of the column:
947  	   if(m_j != m_old_j)
948  	   {
949  	      x[m_old_j] = x[m_j];
950  	      r[m_old_j] = r[m_j];
951  	      cStatus[m_old_j] = cStatus[m_j];
952  	   }
953  	
954  	   // primal:
955  	   R val = 0.0;
956  	   R aij = m_row[m_j];
957  	
958  	   for(int k = 0; k < m_row.size(); ++k)
959  	   {
960  	      if(m_row.index(k) != m_j)
961  	         val += m_row.value(k) * x[m_row.index(k)];
962  	   }
963  	
964  	   R scale = maxAbs(m_lRhs, val);
965  	
966  	   if(scale < 1.0)
967  	      scale = 1.0;
968  	
969  	   R z = (m_lRhs / scale) - (val / scale);
970  	
971  	   if(isZero(z, this->epsilon()))
972  	      z = 0.0;
973  	
974  	   x[m_j] = z * scale / aij;
975  	   s[m_i] = m_lRhs;
976  	
977  	   // dual:
978  	   y[m_i] = m_obj / aij;
979  	   r[m_j] = 0.0;
980  	
981  	   // basis:
982  	   cStatus[m_j] = SPxSolverBase<R>::BASIC;
983  	
984  	   if(m_eqCons)
985  	      rStatus[m_i] = SPxSolverBase<R>::FIXED;
986  	   else if(m_onLhs)
987  	      rStatus[m_i] = SPxSolverBase<R>::ON_LOWER;
988  	   else
989  	      rStatus[m_i] = SPxSolverBase<R>::ON_UPPER;
990  	
991  	#ifdef SOPLEX_CHECK_BASIS_DIM
992  	
993  	   if(!this->checkBasisDim(rStatus, cStatus))
994  	   {
995  	      throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
996  	   }
997  	
998  	#endif
999  	}
1000 	
1001 	template <class R>
1002 	void SPxMainSM<R>::DoubletonEquationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>&,
1003 	      VectorBase<R>& r,
1004 	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
1005 	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
1006 	{
1007 	   // dual:
1008 	   if((cStatus[m_k]  != SPxSolverBase<R>::BASIC) &&
1009 	         ((cStatus[m_k] == SPxSolverBase<R>::ON_LOWER && m_strictLo) ||
1010 	          (cStatus[m_k] == SPxSolverBase<R>::ON_UPPER && m_strictUp) ||
1011 	          (cStatus[m_k] == SPxSolverBase<R>::FIXED    &&
1012 	           ((m_maxSense && ((r[m_j] > 0 && m_strictUp) || (r[m_j] < 0 && m_strictLo))) ||
1013 	            (!m_maxSense && ((r[m_j] > 0 && m_strictLo) || (r[m_j] < 0 && m_strictUp)))))))
1014 	   {
1015 	      R val  = m_kObj;
1016 	      R aik  = m_col[m_i];
1017 	
1018 	      for(int _k = 0; _k < m_col.size(); ++_k)
1019 	      {
1020 	         if(m_col.index(_k) != m_i)
1021 	            val -= m_col.value(_k) * y[m_col.index(_k)];
1022 	      }
1023 	
1024 	      y[m_i] = val / aik;
1025 	      r[m_k] = 0.0;
1026 	
1027 	      r[m_j] = m_jObj - val * m_aij / aik;
1028 	
1029 	      SOPLEX_ASSERT_WARN("WMAISM73", isNotZero(m_aij * aik, this->epsilon()));
1030 	
1031 	      // basis:
1032 	      if(m_jFixed)
1033 	         cStatus[m_j] = SPxSolverBase<R>::FIXED;
1034 	      else
1035 	      {
1036 	         if(GT(r[m_j], (R) 0, this->epsilon()) || (isZero(r[m_j], this->epsilon())
1037 	               && EQ(x[m_j], m_Lo_j, this->epsilon())))
1038 	            cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
1039 	         else
1040 	            cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
1041 	      }
1042 	
1043 	      cStatus[m_k] = SPxSolverBase<R>::BASIC;
1044 	   }
1045 	
1046 	#ifdef SOPLEX_CHECK_BASIS_DIM
1047 	
1048 	   if(!this->checkBasisDim(rStatus, cStatus))
1049 	   {
1050 	      throw SPxInternalCodeException("XMAISM23 Dimension doesn't match after this step.");
1051 	   }
1052 	
1053 	#endif
1054 	}
1055 	
1056 	template <class R>
1057 	void SPxMainSM<R>::DuplicateRowsPS::execute(VectorBase<R>&, VectorBase<R>& y, VectorBase<R>& s,
1058 	      VectorBase<R>&,
1059 	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
1060 	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
1061 	{
1062 	   // correcting the change of idx by deletion of the duplicated rows:
1063 	   if(m_isLast)
1064 	   {
1065 	      for(int i = m_perm.size() - 1; i >= 0; --i)
1066 	      {
1067 	         if(m_perm[i] >= 0)
1068 	         {
1069 	            int rIdx_new = m_perm[i];
1070 	            int rIdx = i;
1071 	            s[rIdx] = s[rIdx_new];
1072 	            y[rIdx] = y[rIdx_new];
1073 	            rStatus[rIdx] = rStatus[rIdx_new];
1074 	         }
1075 	      }
1076 	   }
1077 	
1078 	   // primal:
1079 	   for(int k = 0; k < m_scale.size(); ++k)
1080 	   {
1081 	      if(m_scale.index(k) != m_i)
1082 	         s[m_scale.index(k)] = s[m_i] / m_scale.value(k);
1083 	   }
1084 	
1085 	   // dual & basis:
1086 	   bool haveSetBasis = false;
1087 	
1088 	   for(int k = 0; k < m_scale.size(); ++k)
1089 	   {
1090 	      int i = m_scale.index(k);
1091 	
1092 	      if(rStatus[m_i] == SPxSolverBase<R>::BASIC || (haveSetBasis && i != m_i))
1093 	         // if the row with tightest lower and upper bound in the basic, every duplicate row should in basic
1094 	         // or basis status of row m_i has been set, this row should be in basis
1095 	      {
1096 	         y[i]       = m_rowObj.value(k);
1097 	         rStatus[i] = SPxSolverBase<R>::BASIC;
1098 	         continue;
1099 	      }
1100 	
1101 	      SOPLEX_ASSERT_WARN("WMAISM02", isNotZero(m_scale.value(k), this->epsilon()));
1102 	
1103 	      if(rStatus[m_i] == SPxSolverBase<R>::FIXED && (i == m_maxLhsIdx || i == m_minRhsIdx))
1104 	      {
1105 	         // this row leads to the tightest lower or upper bound, slack should not be in the basis
1106 	         y[i]   = y[m_i] * m_scale.value(k);
1107 	         y[m_i] = m_i_rowObj;
1108 	
1109 	         if(m_isLhsEqualRhs[k])
1110 	         {
1111 	            rStatus[i] = SPxSolverBase<R>::FIXED;
1112 	         }
1113 	         else if(i == m_maxLhsIdx)
1114 	         {
1115 	            rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_LOWER :
1116 	                         SPxSolverBase<R>::ON_UPPER;
1117 	         }
1118 	         else
1119 	         {
1120 	            assert(i == m_minRhsIdx);
1121 	
1122 	            rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_UPPER :
1123 	                         SPxSolverBase<R>::ON_LOWER;
1124 	         }
1125 	
1126 	         if(i != m_i)
1127 	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
1128 	
1129 	         haveSetBasis = true;
1130 	      }
1131 	      else if(i == m_maxLhsIdx && rStatus[m_i] == SPxSolverBase<R>::ON_LOWER)
1132 	      {
1133 	         // this row leads to the tightest lower bound, slack should not be in the basis
1134 	         y[i]   = y[m_i] * m_scale.value(k);
1135 	         y[m_i] = m_i_rowObj;
1136 	
1137 	         rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_LOWER :
1138 	                      SPxSolverBase<R>::ON_UPPER;
1139 	
1140 	         if(i != m_i)
1141 	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
1142 	
1143 	         haveSetBasis = true;
1144 	      }
1145 	      else if(i == m_minRhsIdx && rStatus[m_i] == SPxSolverBase<R>::ON_UPPER)
1146 	      {
1147 	         // this row leads to the tightest upper bound, slack should not be in the basis
1148 	         y[i]   = y[m_i] * m_scale.value(k);
1149 	         y[m_i] = m_i_rowObj;
1150 	
1151 	         rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_UPPER :
1152 	                      SPxSolverBase<R>::ON_LOWER;
1153 	
1154 	         if(i != m_i)
1155 	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
1156 	
1157 	         haveSetBasis = true;
1158 	      }
1159 	      else if(i != m_i)
1160 	      {
1161 	         // this row does not lead to the tightest lower or upper bound, slack should be in the basis
1162 	         y[i]       = m_rowObj.value(k);
1163 	         rStatus[i] = SPxSolverBase<R>::BASIC;
1164 	      }
1165 	   }
1166 	
1167 	#ifdef SOPLEX_CHECK_BASIS_DIM
1168 	
1169 	   if(m_isFirst && !this->checkBasisDim(rStatus, cStatus))
1170 	   {
1171 	      throw SPxInternalCodeException("XMAISM24 Dimension doesn't match after this step.");
1172 	   }
1173 	
1174 	#endif
1175 	
1176 	   // nothing to do for the reduced cost values
1177 	}
1178 	
1179 	template <class R>
1180 	void SPxMainSM<R>::DuplicateColsPS::execute(VectorBase<R>& x,
1181 	      VectorBase<R>&,
1182 	      VectorBase<R>&,
1183 	      VectorBase<R>& r,
1184 	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
1185 	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
1186 	{
1187 	
1188 	   if(m_isFirst)
1189 	   {
1190 	#ifdef SOPLEX_CHECK_BASIS_DIM
1191 	
1192 	      if(!this->checkBasisDim(rStatus, cStatus))
1193 	      {
1194 	         throw SPxInternalCodeException("XMAISM25 Dimension doesn't match after this step.");
1195 	      }
1196 	
1197 	#endif
1198 	      return;
1199 	   }
1200 	
1201 	
1202 	   // correcting the change of idx by deletion of the columns:
1203 	   if(m_isLast)
1204 	   {
1205 	      for(int i = m_perm.size() - 1; i >= 0; --i)
1206 	      {
1207 	         if(m_perm[i] >= 0)
1208 	         {
1209 	            int cIdx_new = m_perm[i];
1210 	            int cIdx = i;
1211 	            x[cIdx] = x[cIdx_new];
1212 	            r[cIdx] = r[cIdx_new];
1213 	            cStatus[cIdx] = cStatus[cIdx_new];
1214 	         }
1215 	      }
1216 	
1217 	      return;
1218 	   }
1219 	
1220 	   // primal & basis:
1221 	   SOPLEX_ASSERT_WARN("WMAISM03", isNotZero(m_scale, this->epsilon()));
1222 	
1223 	   if(cStatus[m_k] == SPxSolverBase<R>::ON_LOWER)
1224 	   {
1225 	      x[m_k] = m_loK;
1226 	
1227 	      if(m_scale > 0)
1228 	      {
1229 	         x[m_j]       = m_loJ;
1230 	         cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
1231 	      }
1232 	      else
1233 	      {
1234 	         x[m_j]       = m_upJ;
1235 	         cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
1236 	      }
1237 	   }
1238 	   else if(cStatus[m_k] == SPxSolverBase<R>::ON_UPPER)
1239 	   {
1240 	      x[m_k] = m_upK;
1241 	
1242 	      if(m_scale > 0)
1243 	      {
1244 	         x[m_j]       = m_upJ;
1245 	         cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
1246 	      }
1247 	      else
1248 	      {
1249 	         x[m_j]       = m_loJ;
1250 	         cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
1251 	      }
1252 	   }
1253 	   else if(cStatus[m_k] == SPxSolverBase<R>::FIXED)
1254 	   {
1255 	      // => x[m_k] and x[m_j] are also fixed before the corresponding preprocessing step
1256 	      x[m_j]       = m_loJ;
1257 	      cStatus[m_j] = SPxSolverBase<R>::FIXED;
1258 	   }
1259 	   else if(cStatus[m_k] == SPxSolverBase<R>::ZERO)
1260 	   {
1261 	      /* we only aggregate duplicate columns if 0 is contained in their bounds, so we can handle this case properly */
1262 	      assert(isZero(x[m_k], this->epsilon()));
1263 	      assert(LErel(m_loJ, R(0.0), this->epsilon()));
1264 	      assert(GErel(m_upJ, R(0.0), this->epsilon()));
1265 	      assert(LErel(m_loK, R(0.0), this->epsilon()));
1266 	      assert(GErel(m_upK, R(0.0), this->epsilon()));
1267 	
1268 	      if(isZero(m_loK, this->epsilon()) && isZero(m_upK, this->epsilon()) && m_loK == m_upK)
1269 	         cStatus[m_k] = SPxSolverBase<R>::FIXED;
1270 	      else if(isZero(m_loK, this->epsilon()))
1271 	         cStatus[m_k] = SPxSolverBase<R>::ON_LOWER;
1272 	      else if(isZero(m_upK, this->epsilon()))
1273 	         cStatus[m_k] = SPxSolverBase<R>::ON_UPPER;
1274 	      else if(LErel(m_loK, R(0.0), this->epsilon()) && GErel(m_upK, R(0.0), this->epsilon()))
1275 	         cStatus[m_k] = SPxSolverBase<R>::ZERO;
1276 	      else
1277 	         throw SPxInternalCodeException("XMAISM05 This should never happen.");
1278 	
1279 	      x[m_j] = 0.0;
1280 	
1281 	      if(isZero(m_loJ, this->epsilon()) && isZero(m_upJ, this->epsilon()) && m_loJ == m_upJ)
1282 	         cStatus[m_j] = SPxSolverBase<R>::FIXED;
1283 	      else if(isZero(m_loJ, this->epsilon()))
1284 	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
1285 	      else if(isZero(m_upJ, this->epsilon()))
1286 	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
1287 	      else if(LErel(m_loJ, R(0.0), this->epsilon()) && GErel(m_upJ, R(0.0), this->epsilon()))
1288 	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
1289 	      else
1290 	         throw SPxInternalCodeException("XMAISM06 This should never happen.");
1291 	   }
1292 	   else if(cStatus[m_k] == SPxSolverBase<R>::BASIC)
1293 	   {
1294 	      R scale1 = maxAbs(x[m_k], m_loK);
1295 	      R scale2 = maxAbs(x[m_k], m_upK);
1296 	
1297 	      if(scale1 < 1.0)
1298 	         scale1 = 1.0;
1299 	
1300 	      if(scale2 < 1.0)
1301 	         scale2 = 1.0;
1302 	
1303 	      R z1 = (x[m_k] / scale1) - (m_loK / scale1);
1304 	      R z2 = (x[m_k] / scale2) - (m_upK / scale2);
1305 	
1306 	      if(isZero(z1, this->epsilon()))
1307 	         z1 = 0.0;
1308 	
1309 	      if(isZero(z2, this->epsilon()))
1310 	         z2 = 0.0;
1311 	
1312 	      if(m_loJ <= R(-infinity) && m_upJ >= R(infinity) && m_loK <= R(-infinity) && m_upK >= R(infinity))
1313 	      {
1314 	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
1315 	         x[m_j] = 0.0;
1316 	      }
1317 	      else if(m_scale > 0.0)
1318 	      {
1319 	         if(GErel(x[m_k], m_upK + m_scale * m_upJ, this->epsilon()))
1320 	         {
1321 	            assert(m_upJ < R(infinity));
1322 	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
1323 	            x[m_j] = m_upJ;
1324 	            x[m_k] -= m_scale * x[m_j];
1325 	         }
1326 	         else if(GErel(x[m_k], m_loK + m_scale * m_upJ, this->epsilon()) && m_upJ < R(infinity))
1327 	         {
1328 	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
1329 	            x[m_j] = m_upJ;
1330 	            x[m_k] -= m_scale * x[m_j];
1331 	         }
1332 	         else if(GErel(x[m_k], m_upK + m_scale * m_loJ, this->epsilon()) && m_upK < R(infinity))
1333 	         {
1334 	            cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
1335 	            x[m_k] = m_upK;
1336 	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
1337 	            x[m_j] = z2 * scale2 / m_scale;
1338 	         }
1339 	         else if(GErel(x[m_k], m_loK + m_scale * m_loJ, this->epsilon()) && m_loJ > R(-infinity))
1340 	         {
1341 	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
1342 	            x[m_j] = m_loJ;
1343 	            x[m_k] -= m_scale * x[m_j];
1344 	         }
1345 	         else if(GErel(x[m_k], m_loK + m_scale * m_loJ, this->epsilon()) && m_loK > R(-infinity))
1346 	         {
1347 	            cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
1348 	            x[m_k] = m_loK;
1349 	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
1350 	            x[m_j] = z1 * scale1 / m_scale;
1351 	         }
1352 	         else if(LTrel(x[m_k], m_loK + m_scale * m_loJ, this->epsilon()))
1353 	         {
1354 	            assert(m_loJ > R(-infinity));
1355 	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
1356 	            x[m_j] = m_loJ;
1357 	            x[m_k] -= m_scale * x[m_j];
1358 	         }
1359 	         else
1360 	         {
1361 	            throw SPxInternalCodeException("XMAISM08 This should never happen.");
1362 	         }
1363 	      }
1364 	      else
1365 	      {
1366 	         assert(m_scale < 0.0);
1367 	
1368 	         if(GErel(x[m_k], m_upK + m_scale * m_loJ, this->epsilon()))
1369 	         {
1370 	            assert(m_loJ > R(-infinity));
1371 	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
1372 	            x[m_j] = m_loJ;
1373 	            x[m_k] -= m_scale * x[m_j];
1374 	         }
1375 	         else if(GErel(x[m_k], m_loK + m_scale * m_loJ, this->epsilon()) && m_loJ > R(-infinity))
1376 	         {
1377 	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
1378 	            x[m_j] = m_loJ;
1379 	            x[m_k] -= m_scale * x[m_j];
1380 	         }
1381 	         else if(GErel(x[m_k], m_upK + m_scale * m_upJ, this->epsilon()) && m_upK < R(infinity))
1382 	         {
1383 	            cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
1384 	            x[m_k] = m_upK;
1385 	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
1386 	            x[m_j] = z2 * scale2 / m_scale;
1387 	         }
1388 	         else if(GErel(x[m_k], m_loK + m_scale * m_upJ, this->epsilon()) && m_upJ < R(infinity))
1389 	         {
1390 	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
1391 	            x[m_j] = m_upJ;
1392 	            x[m_k] -= m_scale * x[m_j];
1393 	         }
1394 	         else if(GErel(x[m_k], m_loK + m_scale * m_upJ, this->epsilon()) && m_loK > R(-infinity))
1395 	         {
1396 	            cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
1397 	            x[m_k] = m_loK;
1398 	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
1399 	            x[m_j] = z1 * scale1 / m_scale;
1400 	         }
1401 	         else if(LTrel(x[m_k], m_loK + m_scale * m_upJ, this->epsilon()))
1402 	         {
1403 	            assert(m_upJ < R(infinity));
1404 	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
1405 	            x[m_j] = m_upJ;
1406 	            x[m_k] -= m_scale * x[m_j];
1407 	         }
1408 	         else
1409 	         {
1410 	            throw SPxInternalCodeException("XMAISM09 This should never happen.");
1411 	         }
1412 	      }
1413 	   }
1414 	
1415 	   // dual:
1416 	   r[m_j] = m_scale * r[m_k];
1417 	}
1418 	
1419 	template <class R>
1420 	void SPxMainSM<R>::AggregationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
1421 	      VectorBase<R>& r,
1422 	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
1423 	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
1424 	{
1425 	   // correcting the change of idx by deletion of the row:
1426 	   if(m_i != m_old_i)
1427 	   {
1428 	      s[m_old_i] = s[m_i];
1429 	      y[m_old_i] = y[m_i];
1430 	      rStatus[m_old_i] = rStatus[m_i];
1431 	   }
1432 	
1433 	   // correcting the change of idx by deletion of the column:
1434 	   if(m_j != m_old_j)
1435 	   {
1436 	      x[m_old_j] = x[m_j];
1437 	      r[m_old_j] = r[m_j];
1438 	      cStatus[m_old_j] = cStatus[m_j];
1439 	   }
1440 	
1441 	   // primal:
1442 	   R val = 0.0;
1443 	   R aij = m_row[m_j];
1444 	   int active_idx = -1;
1445 	
1446 	   assert(m_row.size() == 2);
1447 	
1448 	   for(int k = 0; k < 2; ++k)
1449 	   {
1450 	      if(m_row.index(k) != m_j)
1451 	      {
1452 	         active_idx = m_row.index(k);
1453 	         val = m_row.value(k) * x[active_idx];
1454 	      }
1455 	   }
1456 	
1457 	   assert(active_idx >= 0);
1458 	
1459 	   R scale = maxAbs(m_rhs, val);
1460 	
1461 	   if(scale < 1.0)
1462 	      scale = 1.0;
1463 	
1464 	   R z = (m_rhs / scale) - (val / scale);
1465 	
1466 	   if(isZero(z, this->epsilon()))
1467 	      z = 0.0;
1468 	
1469 	   x[m_j] = z * scale / aij;
1470 	   s[m_i] = m_rhs;
1471 	
1472 	   if(isOptimal && (LT(x[m_j], m_lower, this->feastol()) || GT(x[m_j], m_upper, this->feastol())))
1473 	   {
1474 	      SPX_MSG_ERROR(std::cerr << "EMAISM: numerical violation after disaggregating variable" << std::endl;
1475 	                   )
1476 	   }
1477 	
1478 	   // dual:
1479 	   R dualVal = 0.0;
1480 	
1481 	   for(int k = 0; k < m_col.size(); ++k)
1482 	   {
1483 	      if(m_col.index(k) != m_i)
1484 	         dualVal += m_col.value(k) * y[m_col.index(k)];
1485 	   }
1486 	
1487 	   z = m_obj - dualVal;
1488 	
1489 	   y[m_i] = z / aij;
1490 	   r[m_j] = 0.0;
1491 	
1492 	   // basis:
1493 	   if(((cStatus[active_idx] == SPxSolverBase<R>::ON_UPPER
1494 	         || cStatus[active_idx] == SPxSolverBase<R>::FIXED)
1495 	         && NE(x[active_idx], m_oldupper, this->feastol())) ||
1496 	         ((cStatus[active_idx] == SPxSolverBase<R>::ON_LOWER
1497 	           || cStatus[active_idx] == SPxSolverBase<R>::FIXED)
1498 	          && NE(x[active_idx], m_oldlower, this->feastol())))
1499 	   {
1500 	      cStatus[active_idx] = SPxSolverBase<R>::BASIC;
1501 	      r[active_idx] = 0.0;
1502 	      assert(NE(m_upper, m_lower, this->epsilon()));
1503 	
1504 	      if(EQ(x[m_j], m_upper, this->feastol()))
1505 	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
1506 	      else if(EQ(x[m_j], m_lower, this->feastol()))
1507 	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
1508 	      else if(m_upper >= R(infinity) && m_lower <= R(-infinity))
1509 	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
1510 	      else
1511 	         throw SPxInternalCodeException("XMAISM unexpected basis status in aggregation unsimplifier.");
1512 	   }
1513 	   else
1514 	   {
1515 	      cStatus[m_j] = SPxSolverBase<R>::BASIC;
1516 	   }
1517 	
1518 	   // sides may not be equal and we always only consider the rhs during aggregation, so set ON_UPPER
1519 	   // (in theory and with exact arithmetic setting it to FIXED would be correct)
1520 	   rStatus[m_i] = SPxSolverBase<R>::ON_UPPER;
1521 	
1522 	#ifdef SOPLEX_CHECK_BASIS_DIM
1523 	
1524 	   if(!this->checkBasisDim(rStatus, cStatus))
1525 	   {
1526 	      throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
1527 	   }
1528 	
1529 	#endif
1530 	}
1531 	
1532 	template <class R>
1533 	void SPxMainSM<R>::MultiAggregationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
1534 	      VectorBase<R>& r,
1535 	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
1536 	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
1537 	{
1538 	
1539 	   // correcting the change of idx by deletion of the row:
1540 	   if(m_i != m_old_i)
1541 	   {
1542 	      s[m_old_i] = s[m_i];
1543 	      y[m_old_i] = y[m_i];
1544 	      rStatus[m_old_i] = rStatus[m_i];
1545 	   }
1546 	
1547 	   // correcting the change of idx by deletion of the column:
1548 	   if(m_j != m_old_j)
1549 	   {
1550 	      x[m_old_j] = x[m_j];
1551 	      r[m_old_j] = r[m_j];
1552 	      cStatus[m_old_j] = cStatus[m_j];
1553 	   }
1554 	
1555 	   // primal:
1556 	   R val = 0.0;
1557 	   R aij = m_row[m_j];
1558 	
1559 	   for(int k = 0; k < m_row.size(); ++k)
1560 	   {
1561 	      if(m_row.index(k) != m_j)
1562 	         val += m_row.value(k) * x[m_row.index(k)];
1563 	   }
1564 	
1565 	   R scale = maxAbs(m_const, val);
1566 	
1567 	   if(scale < 1.0)
1568 	      scale = 1.0;
1569 	
1570 	   R z = (m_const / scale) - (val / scale);
1571 	
1572 	   if(isZero(z, this->epsilon()))
1573 	      z = 0.0;
1574 	
1575 	   x[m_j] = z * scale / aij;
1576 	   s[m_i] = 0.0;
1577 	
1578 	#ifndef NDEBUG
1579 	
1580 	   if(isOptimal && (LT(x[m_j], m_lower, this->feastol()) || GT(x[m_j], m_upper, this->feastol())))
1581 	      SPX_MSG_ERROR(std::cerr << "numerical violation in original space due to MultiAggregation\n";)
1582 	#endif
1583 	
1584 	      // dual:
1585 	      R dualVal = 0.0;
1586 	
1587 	   for(int k = 0; k < m_col.size(); ++k)
1588 	   {
1589 	      if(m_col.index(k) != m_i)
1590 	         dualVal += m_col.value(k) * y[m_col.index(k)];
1591 	   }
1592 	
1593 	   z = m_obj - dualVal;
1594 	
1595 	   y[m_i] = z / aij;
1596 	   r[m_j] = 0.0;
1597 	
1598 	   // basis:
1599 	   cStatus[m_j] = SPxSolverBase<R>::BASIC;
1600 	
1601 	   if(m_eqCons)
1602 	      rStatus[m_i] = SPxSolverBase<R>::FIXED;
1603 	   else if(m_onLhs)
1604 	      rStatus[m_i] = SPxSolverBase<R>::ON_LOWER;
1605 	   else
1606 	      rStatus[m_i] = SPxSolverBase<R>::ON_UPPER;
1607 	
1608 	#ifdef SOPLEX_CHECK_BASIS_DIM
1609 	
1610 	   if(!this->checkBasisDim(rStatus, cStatus))
1611 	   {
1612 	      throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
1613 	   }
1614 	
1615 	#endif
1616 	}
1617 	
1618 	template <class R>
1619 	void SPxMainSM<R>::TightenBoundsPS::execute(VectorBase<R>& x, VectorBase<R>&, VectorBase<R>&,
1620 	      VectorBase<R>&,
1621 	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
1622 	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
1623 	{
1624 	   // basis:
1625 	   switch(cStatus[m_j])
1626 	   {
1627 	   case SPxSolverBase<R>::FIXED:
1628 	      if(LT(x[m_j], m_origupper, this->feastol()) && GT(x[m_j], m_origlower, this->feastol()))
1629 	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
1630 	      else if(LT(x[m_j], m_origupper, this->feastol()))
1631 	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
1632 	      else if(GT(x[m_j], m_origlower, this->feastol()))
1633 	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
1634 	
1635 	      break;
1636 	
1637 	   case SPxSolverBase<R>::ON_LOWER:
1638 	      if(GT(x[m_j], m_origlower, this->feastol()))
1639 	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
1640 	
1641 	      break;
1642 	
1643 	   case SPxSolverBase<R>::ON_UPPER:
1644 	      if(LT(x[m_j], m_origupper, this->feastol()))
1645 	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
1646 	
1647 	      break;
1648 	
1649 	   default:
1650 	      break;
1651 	   }
1652 	
1653 	#ifdef SOPLEX_CHECK_BASIS_DIM
1654 	
1655 	   if(!this->checkBasisDim(rStatus, cStatus))
1656 	   {
1657 	      throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
1658 	   }
1659 	
1660 	#endif
1661 	}
1662 	
1663 	template <class R>
1664 	void SPxMainSM<R>::handleRowObjectives(SPxLPBase<R>& lp)
1665 	{
1666 	   for(int i = lp.nRows() - 1; i >= 0; --i)
1667 	   {
1668 	      if(lp.maxRowObj(i) != 0.0)
1669 	      {
1670 	         std::shared_ptr<PostStep> ptr(new RowObjPS(lp, i, lp.nCols(), this->_tolerances));
1671 	         m_hist.append(ptr);
1672 	         lp.addCol(lp.rowObj(i), -lp.rhs(i), UnitVectorBase<R>(i), -lp.lhs(i));
1673 	         lp.changeRange(i, R(0.0), R(0.0));
1674 	         lp.changeRowObj(i, R(0.0));
1675 	         m_addedcols++;
1676 	      }
1677 	   }
1678 	}
1679 	
1680 	template <class R>
1681 	void SPxMainSM<R>::handleExtremes(SPxLPBase<R>& lp)
1682 	{
1683 	
1684 	   // This method handles extreme value of the given LP by
1685 	   //
1686 	   // 1. setting numbers of very small absolute values to zero and
1687 	   // 2. setting numbers of very large absolute values to R(-infinity) or +R(infinity), respectively.
1688 	
1689 	   R maxVal  = R(infinity) / 5.0;
1690 	   R tol = feastol() * 1e-2;
1691 	   tol = (tol < this->epsZero()) ? this->epsZero() : tol;
1692 	   int  remRows = 0;
1693 	   int  remNzos = 0;
1694 	   int  chgBnds = 0;
1695 	   int  chgLRhs = 0;
1696 	   int  objCnt  = 0;
1697 	
1698 	   for(int i = lp.nRows() - 1; i >= 0; --i)
1699 	   {
1700 	      // lhs
1701 	      R lhs = lp.lhs(i);
1702 	
1703 	      if(lhs != 0.0 && isZero(lhs, this->epsZero()))
1704 	      {
1705 	         lp.changeLhs(i, R(0.0));
1706 	         ++chgLRhs;
1707 	      }
1708 	      else if(lhs > R(-infinity) && lhs < -maxVal)
1709 	      {
1710 	         lp.changeLhs(i, R(-infinity));
1711 	         ++chgLRhs;
1712 	      }
1713 	      else if(lhs <  R(infinity) && lhs >  maxVal)
1714 	      {
1715 	         lp.changeLhs(i,  R(infinity));
1716 	         ++chgLRhs;
1717 	      }
1718 	
1719 	      // rhs
1720 	      R rhs = lp.rhs(i);
1721 	
1722 	      if(rhs != 0.0 && isZero(rhs, this->epsZero()))
1723 	      {
1724 	         lp.changeRhs(i, R(0.0));
1725 	         ++chgLRhs;
1726 	      }
1727 	      else if(rhs > R(-infinity) && rhs < -maxVal)
1728 	      {
1729 	         lp.changeRhs(i, R(-infinity));
1730 	         ++chgLRhs;
1731 	      }
1732 	      else if(rhs <  R(infinity) && rhs >  maxVal)
1733 	      {
1734 	         lp.changeRhs(i,  R(infinity));
1735 	         ++chgLRhs;
1736 	      }
1737 	
1738 	      if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
1739 	      {
1740 	         std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i, this->_tolerances));
1741 	         m_hist.append(ptr);
1742 	
1743 	         removeRow(lp, i);
1744 	         ++remRows;
1745 	
1746 	         ++m_stat[FREE_ROW];
1747 	      }
1748 	   }
1749 	
1750 	   for(int j = 0; j < lp.nCols(); ++j)
1751 	   {
1752 	      // lower bound
1753 	      R lo = lp.lower(j);
1754 	
1755 	      if(lo != 0.0 && isZero(lo, this->epsZero()))
1756 	      {
1757 	         lp.changeLower(j, R(0.0));
1758 	         ++chgBnds;
1759 	      }
1760 	      else if(lo > R(-infinity) && lo < -maxVal)
1761 	      {
1762 	         lp.changeLower(j, R(-infinity));
1763 	         ++chgBnds;
1764 	      }
1765 	      else if(lo <  R(infinity) && lo >  maxVal)
1766 	      {
1767 	         lp.changeLower(j,  R(infinity));
1768 	         ++chgBnds;
1769 	      }
1770 	
1771 	      // upper bound
1772 	      R up = lp.upper(j);
1773 	
1774 	      if(up != 0.0 && isZero(up, this->epsZero()))
1775 	      {
1776 	         lp.changeUpper(j, R(0.0));
1777 	         ++chgBnds;
1778 	      }
1779 	      else if(up > R(-infinity) && up < -maxVal)
1780 	      {
1781 	         lp.changeUpper(j, R(-infinity));
1782 	         ++chgBnds;
1783 	      }
1784 	      else if(up <  R(infinity) && up >  maxVal)
1785 	      {
1786 	         lp.changeUpper(j,  R(infinity));
1787 	         ++chgBnds;
1788 	      }
1789 	
1790 	      // fixed columns will be eliminated later
1791 	      if(NE(lo, up, this->tolerances()->epsilon()))
1792 	      {
1793 	         lo = spxAbs(lo);
1794 	         up = spxAbs(up);
1795 	
1796 	         R absBnd = (lo > up) ? lo : up;
1797 	
1798 	         if(absBnd < 1.0)
1799 	            absBnd = 1.0;
1800 	
1801 	         // non-zeros
1802 	         SVectorBase<R>& col = lp.colVector_w(j);
1803 	         int        i = 0;
1804 	
1805 	         while(i < col.size())
1806 	         {
1807 	            R aij = spxAbs(col.value(i));
1808 	
1809 	            if(isZero(aij * absBnd, tol))
1810 	            {
1811 	               SVectorBase<R>& row = lp.rowVector_w(col.index(i));
1812 	               int row_j = row.pos(j);
1813 	
1814 	               // this changes col.size()
1815 	               if(row_j >= 0)
1816 	                  row.remove(row_j);
1817 	
1818 	               col.remove(i);
1819 	
1820 	               SPxOut::debug(this, "IMAISM04 aij={}  removed, absBnd={} \n", aij, absBnd);
1821 	               ++remNzos;
1822 	            }
1823 	            else
1824 	            {
1825 	               if(aij > maxVal)
1826 	               {
1827 	                  SPX_MSG_WARNING((*this->spxout),
1828 	                                  (*this->spxout) << "WMAISM05 Warning! Big matrix coefficient " << aij
1829 	                                  << std::endl);
1830 	               }
1831 	               else if(isZero(aij, tol))
1832 	               {
1833 	                  SPX_MSG_WARNING((*this->spxout),
1834 	                                  (*this->spxout) << "WMAISM06 Warning! Tiny matrix coefficient " << aij
1835 	                                  << std::endl);
1836 	               }
1837 	
1838 	               ++i;
1839 	            }
1840 	         }
1841 	      }
1842 	
1843 	      // objective
1844 	      R obj = lp.obj(j);
1845 	
1846 	      if(obj != 0.0 && isZero(obj, this->epsZero()))
1847 	      {
1848 	         lp.changeObj(j, R(0.0));
1849 	         ++objCnt;
1850 	      }
1851 	      else if(obj > R(-infinity) && obj < -maxVal)
1852 	      {
1853 	         lp.changeObj(j, R(-infinity));
1854 	         ++objCnt;
1855 	      }
1856 	      else if(obj <  R(infinity) && obj >  maxVal)
1857 	      {
1858 	         lp.changeObj(j,  R(infinity));
1859 	         ++objCnt;
1860 	      }
1861 	   }
1862 	
1863 	   if(remRows + remNzos + chgLRhs + chgBnds + objCnt > 0)
1864 	   {
1865 	      this->m_remRows += remRows;
1866 	      this->m_remNzos += remNzos;
1867 	      this->m_chgLRhs += chgLRhs;
1868 	      this->m_chgBnds += chgBnds;
1869 	
1870 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (extremes) removed "
1871 	                    << remRows << " rows, "
1872 	                    << remNzos << " non-zeros, "
1873 	                    << chgBnds << " col bounds, "
1874 	                    << chgLRhs << " row bounds, "
1875 	                    << objCnt  << " objective coefficients" << std::endl;)
1876 	   }
1877 	
1878 	   assert(lp.isConsistent());
1879 	}
1880 	
1881 	/// computes the minimum and maximum residual activity for a given variable
1882 	template <class R>
1883 	void SPxMainSM<R>::computeMinMaxResidualActivity(SPxLPBase<R>& lp, int rowNumber, int colNumber,
1884 	      R& minAct, R& maxAct)
1885 	{
1886 	   const SVectorBase<R>& row = lp.rowVector(rowNumber);
1887 	   bool minNegInfinite = false;
1888 	   bool maxInfinite = false;
1889 	
1890 	   minAct = 0;   // this is the minimum value that the aggregation can attain
1891 	   maxAct = 0;   // this is the maximum value that the aggregation can attain
1892 	
1893 	   for(int l = 0; l < row.size(); ++l)
1894 	   {
1895 	      if(colNumber < 0 || row.index(l) != colNumber)
1896 	      {
1897 	         // computing the minimum activity of the aggregated variables
1898 	         if(GT(row.value(l), R(0.0), this->tolerances()->epsilon()))
1899 	         {
1900 	            if(lp.lower(row.index(l)) <= R(-infinity))
1901 	               minNegInfinite = true;
1902 	            else
1903 	               minAct += row.value(l) * lp.lower(row.index(l));
1904 	         }
1905 	         else if(LT(row.value(l), R(0.0), this->tolerances()->epsilon()))
1906 	         {
1907 	            if(lp.upper(row.index(l)) >= R(infinity))
1908 	               minNegInfinite = true;
1909 	            else
1910 	               minAct += row.value(l) * lp.upper(row.index(l));
1911 	         }
1912 	
1913 	         // computing the maximum activity of the aggregated variables
1914 	         if(GT(row.value(l), R(0.0), this->tolerances()->epsilon()))
1915 	         {
1916 	            if(lp.upper(row.index(l)) >= R(infinity))
1917 	               maxInfinite = true;
1918 	            else
1919 	               maxAct += row.value(l) * lp.upper(row.index(l));
1920 	         }
1921 	         else if(LT(row.value(l), R(0.0), this->tolerances()->epsilon()))
1922 	         {
1923 	            if(lp.lower(row.index(l)) <= R(-infinity))
1924 	               maxInfinite = true;
1925 	            else
1926 	               maxAct += row.value(l) * lp.lower(row.index(l));
1927 	         }
1928 	      }
1929 	   }
1930 	
1931 	   // if an infinite value exists for the minimum activity, then that it taken
1932 	   if(minNegInfinite)
1933 	      minAct = R(-infinity);
1934 	
1935 	   // if an -infinite value exists for the maximum activity, then that value is taken
1936 	   if(maxInfinite)
1937 	      maxAct = R(infinity);
1938 	}
1939 	
1940 	
1941 	/// calculate min/max value for the multi aggregated variables
1942 	template <class R>
1943 	void SPxMainSM<R>::computeMinMaxValues(SPxLPBase<R>& lp, R side, R val, R minRes, R maxRes,
1944 	                                       R& minVal, R& maxVal)
1945 	{
1946 	   minVal = 0;
1947 	   maxVal = 0;
1948 	
1949 	   if(LT(val, R(0.0), this->tolerances()->epsilon()))
1950 	   {
1951 	      if(minRes <= R(-infinity))
1952 	         minVal = R(-infinity);
1953 	      else
1954 	         minVal = (side - minRes) / val;
1955 	
1956 	      if(maxRes >= R(infinity))
1957 	         maxVal = R(infinity);
1958 	      else
1959 	         maxVal = (side - maxRes) / val;
1960 	   }
1961 	   else if(GT(val, R(0.0), this->tolerances()->epsilon()))
1962 	   {
1963 	      if(maxRes >= R(infinity))
1964 	         minVal = R(-infinity);
1965 	      else
1966 	         minVal = (side - maxRes) / val;
1967 	
1968 	      if(minRes <= R(-infinity))
1969 	         maxVal = R(infinity);
1970 	      else
1971 	         maxVal = (side - minRes) / val;
1972 	   }
1973 	}
1974 	
1975 	
1976 	/// tries to find good lower bound solutions by applying some trivial heuristics
1977 	template <class R>
1978 	void SPxMainSM<R>::trivialHeuristic(SPxLPBase<R>& lp)
1979 	{
1980 	   VectorBase<R>         zerosol(lp.nCols());  // the zero solution VectorBase<R>
1981 	   VectorBase<R>         lowersol(lp.nCols()); // the lower bound solution VectorBase<R>
1982 	   VectorBase<R>         uppersol(lp.nCols()); // the upper bound solution VectorBase<R>
1983 	   VectorBase<R>         locksol(lp.nCols());  // the locks solution VectorBase<R>
1984 	
1985 	   VectorBase<R>         upLocks(lp.nCols());
1986 	   VectorBase<R>         downLocks(lp.nCols());
1987 	
1988 	   R            zeroObj = this->m_objoffset;
1989 	   R            lowerObj = this->m_objoffset;
1990 	   R            upperObj = this->m_objoffset;
1991 	   R            lockObj = this->m_objoffset;
1992 	
1993 	   bool            zerovalid = true;
1994 	
1995 	   R largeValue = R(infinity);
1996 	
1997 	   if(LT(R(1.0 / feastol()), R(infinity), this->tolerances()->epsilon()))
1998 	      largeValue = 1.0 / feastol();
1999 	
2000 	
2001 	
2002 	   for(int j = lp.nCols() - 1; j >= 0; --j)
2003 	   {
2004 	      upLocks[j] = 0;
2005 	      downLocks[j] = 0;
2006 	
2007 	      // computing the locks on the variables
2008 	      const SVectorBase<R>& col = lp.colVector(j);
2009 	
2010 	      for(int k = 0; k < col.size(); ++k)
2011 	      {
2012 	         R aij = col.value(k);
2013 	
2014 	         SOPLEX_ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity))));
2015 	
2016 	         if(GT(lp.lhs(col.index(k)), R(-infinity), this->tolerances()->epsilon())
2017 	               && LT(lp.rhs(col.index(k)), R(infinity), this->tolerances()->epsilon()))
2018 	         {
2019 	            upLocks[j]++;
2020 	            downLocks[j]++;
2021 	         }
2022 	         else if(GT(lp.lhs(col.index(k)), R(-infinity), this->tolerances()->epsilon()))
2023 	         {
2024 	            if(aij > 0)
2025 	               downLocks[j]++;
2026 	            else if(aij < 0)
2027 	               upLocks[j]++;
2028 	         }
2029 	         else if(LT(lp.rhs(col.index(k)), R(infinity), this->tolerances()->epsilon()))
2030 	         {
2031 	            if(aij > 0)
2032 	               upLocks[j]++;
2033 	            else if(aij < 0)
2034 	               downLocks[j]++;
2035 	         }
2036 	      }
2037 	
2038 	      R lower = lp.lower(j);
2039 	      R upper = lp.upper(j);
2040 	
2041 	      if(LE(lower, R(-infinity), this->tolerances()->epsilon()))
2042 	         lower = SOPLEX_MIN(-largeValue, upper);
2043 	
2044 	      if(GE(upper, R(infinity), this->tolerances()->epsilon()))
2045 	         upper = SOPLEX_MAX(lp.lower(j), largeValue);
2046 	
2047 	      if(zerovalid)
2048 	      {
2049 	         if(LE(lower, R(0.0), feastol()) && GE(upper, R(0.0), feastol()))
2050 	            zerosol[j] = 0.0;
2051 	         else
2052 	            zerovalid = false;
2053 	      }
2054 	
2055 	      lowersol[j] = lower;
2056 	      uppersol[j] = upper;
2057 	
2058 	      if(downLocks[j] > upLocks[j])
2059 	         locksol[j] = upper;
2060 	      else if(downLocks[j] < upLocks[j])
2061 	         locksol[j] = lower;
2062 	      else
2063 	         locksol[j] = (lower + upper) / 2.0;
2064 	
2065 	      lowerObj += lp.maxObj(j) * lowersol[j];
2066 	      upperObj += lp.maxObj(j) * uppersol[j];
2067 	      lockObj += lp.maxObj(j) * locksol[j];
2068 	   }
2069 	
2070 	   // trying the lower bound solution
2071 	   if(checkSolution(lp, lowersol))
2072 	   {
2073 	      if(lowerObj > m_cutoffbound)
2074 	         m_cutoffbound = lowerObj;
2075 	   }
2076 	
2077 	   // trying the upper bound solution
2078 	   if(checkSolution(lp, uppersol))
2079 	   {
2080 	      if(upperObj > m_cutoffbound)
2081 	         m_cutoffbound = upperObj;
2082 	   }
2083 	
2084 	   // trying the zero solution
2085 	   if(zerovalid && checkSolution(lp, zerosol))
2086 	   {
2087 	      if(zeroObj > m_cutoffbound)
2088 	         m_cutoffbound = zeroObj;
2089 	   }
2090 	
2091 	   // trying the lock solution
2092 	   if(checkSolution(lp, locksol))
2093 	   {
2094 	      if(lockObj > m_cutoffbound)
2095 	         m_cutoffbound = lockObj;
2096 	   }
2097 	}
2098 	
2099 	
2100 	
2101 	/// checks a solution for feasibility
2102 	template <class R>
2103 	bool SPxMainSM<R>::checkSolution(SPxLPBase<R>& lp, VectorBase<R> sol)
2104 	{
2105 	   for(int i = lp.nRows() - 1; i >= 0; --i)
2106 	   {
2107 	      const SVectorBase<R>& row = lp.rowVector(i);
2108 	      R activity = 0;
2109 	
2110 	      for(int k = 0; k < row.size(); k++)
2111 	         activity += row.value(k) * sol[row.index(k)];
2112 	
2113 	      if(!GE(activity, lp.lhs(i), feastol()) || !LE(activity, lp.rhs(i), feastol()))
2114 	         return false;
2115 	   }
2116 	
2117 	   return true;
2118 	}
2119 	
2120 	
2121 	/// tightens variable bounds by propagating the pseudo objective function value.
2122 	template <class R>
2123 	void SPxMainSM<R>::propagatePseudoobj(SPxLPBase<R>& lp)
2124 	{
2125 	   R pseudoObj = this->m_objoffset;
2126 	
2127 	   for(int j = lp.nCols() - 1; j >= 0; --j)
2128 	   {
2129 	      R val = lp.maxObj(j);
2130 	
2131 	      if(val < 0)
2132 	      {
2133 	         if(lp.lower(j) <= R(-infinity))
2134 	            return;
2135 	
2136 	         pseudoObj += val * lp.lower(j);
2137 	      }
2138 	      else if(val > 0)
2139 	      {
2140 	         if(lp.upper(j) >= R(-infinity))
2141 	            return;
2142 	
2143 	         pseudoObj += val * lp.upper(j);
2144 	      }
2145 	   }
2146 	
2147 	   if(GT(m_cutoffbound, R(-infinity), this->tolerances()->epsilon())
2148 	         && LT(m_cutoffbound, R(infinity), this->tolerances()->epsilon()))
2149 	   {
2150 	      if(pseudoObj > m_pseudoobj)
2151 	         m_pseudoobj = pseudoObj;
2152 	
2153 	      for(int j = lp.nCols() - 1; j >= 0; --j)
2154 	      {
2155 	         R objval = lp.maxObj(j);
2156 	
2157 	         if(EQ(objval, R(0.0), this->tolerances()->epsilon()))
2158 	            continue;
2159 	
2160 	         if(objval < 0.0)
2161 	         {
2162 	            R newbound = lp.lower(j) + (m_cutoffbound - m_pseudoobj) / objval;
2163 	
2164 	            if(LT(newbound, lp.upper(j), this->tolerances()->epsilon()))
2165 	            {
2166 	               std::shared_ptr<PostStep> ptr(new TightenBoundsPS(lp, j, lp.upper(j), lp.lower(j),
2167 	                                             this->_tolerances));
2168 	               m_hist.append(ptr);
2169 	               lp.changeUpper(j, newbound);
2170 	            }
2171 	         }
2172 	         else if(objval > 0.0)
2173 	         {
2174 	            R newbound = lp.upper(j) + (m_cutoffbound - m_pseudoobj) / objval;
2175 	
2176 	            if(GT(newbound, lp.lower(j), this->tolerances()->epsilon()))
2177 	            {
2178 	               std::shared_ptr<PostStep> ptr(new TightenBoundsPS(lp, j, lp.upper(j), lp.lower(j),
2179 	                                             this->_tolerances));
2180 	               m_hist.append(ptr);
2181 	               lp.changeLower(j, newbound);
2182 	            }
2183 	         }
2184 	      }
2185 	   }
2186 	}
2187 	
2188 	
2189 	
2190 	template <class R>
2191 	typename SPxSimplifier<R>::Result SPxMainSM<R>::removeEmpty(SPxLPBase<R>& lp)
2192 	{
2193 	
2194 	   // This method removes empty rows and columns from the LP.
2195 	
2196 	   int remRows = 0;
2197 	   int remCols = 0;
2198 	
2199 	   for(int i = lp.nRows() - 1; i >= 0; --i)
2200 	   {
2201 	      const SVectorBase<R>& row = lp.rowVector(i);
2202 	
2203 	      if(row.size() == 0)
2204 	      {
2205 	         SPxOut::debug(this, "IMAISM07 row {}: empty ->", i);
2206 	
2207 	         if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol()))
2208 	         {
2209 	            SPxOut::debug(this, " infeasible lhs={} rhs={}\n", lp.lhs(i), lp.rhs(i));
2210 	            return this->INFEASIBLE;
2211 	         }
2212 	
2213 	         SPxOut::debug(this, " removed\n");
2214 	
2215 	         std::shared_ptr<PostStep> ptr(new EmptyConstraintPS(lp, i, this->_tolerances));
2216 	         m_hist.append(ptr);
2217 	
2218 	         ++remRows;
2219 	         removeRow(lp, i);
2220 	
2221 	         ++m_stat[EMPTY_ROW];
2222 	      }
2223 	   }
2224 	
2225 	   for(int j = lp.nCols() - 1; j >= 0; --j)
2226 	   {
2227 	      const SVectorBase<R>& col = lp.colVector(j);
2228 	
2229 	      if(col.size() == 0)
2230 	      {
2231 	         SPxOut::debug(this, "IMAISM08 col {}: empty -> maxObj={} lower={} upper={}", j, lp.maxObj(j),
2232 	                       lp.lower(j), lp.upper(j));
2233 	
2234 	         R val;
2235 	
2236 	         if(GT(lp.maxObj(j), R(0.0), this->epsZero()))
2237 	         {
2238 	            if(lp.upper(j) >= R(infinity))
2239 	            {
2240 	               SPxOut::debug(this, " unbounded\n");
2241 	               return this->UNBOUNDED;
2242 	            }
2243 	
2244 	            val = lp.upper(j);
2245 	         }
2246 	         else if(LT(lp.maxObj(j), R(0.0), this->epsZero()))
2247 	         {
2248 	            if(lp.lower(j) <= R(-infinity))
2249 	            {
2250 	               SPxOut::debug(this, " unbounded\n");
2251 	               return this->UNBOUNDED;
2252 	            }
2253 	
2254 	            val = lp.lower(j);
2255 	         }
2256 	         else
2257 	         {
2258 	            SOPLEX_ASSERT_WARN("WMAISM09", isZero(lp.maxObj(j), this->epsZero()));
2259 	
2260 	            // any value within the bounds is ok
2261 	            if(lp.lower(j) > R(-infinity))
2262 	               val = lp.lower(j);
2263 	            else if(lp.upper(j) < R(infinity))
2264 	               val = lp.upper(j);
2265 	            else
2266 	               val = 0.0;
2267 	         }
2268 	
2269 	         SPxOut::debug(this, " removed\n");
2270 	
2271 	         std::shared_ptr<PostStep> ptr1(new FixBoundsPS(lp, j, val, this->_tolerances));
2272 	         std::shared_ptr<PostStep> ptr2(new FixVariablePS(lp, *this, j, val, this->_tolerances));
2273 	         m_hist.append(ptr1);
2274 	         m_hist.append(ptr2);
2275 	
2276 	         ++remCols;
2277 	         removeCol(lp, j);
2278 	
2279 	         ++m_stat[EMPTY_COL];
2280 	      }
2281 	   }
2282 	
2283 	   if(remRows + remCols > 0)
2284 	   {
2285 	      this->m_remRows += remRows;
2286 	      this->m_remCols += remCols;
2287 	
2288 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (empty rows/colums) removed "
2289 	                    << remRows << " rows, "
2290 	                    << remCols << " cols"
2291 	                    << std::endl;)
2292 	
2293 	   }
2294 	
2295 	   return this->OKAY;
2296 	}
2297 	
2298 	template <class R>
2299 	typename SPxSimplifier<R>::Result SPxMainSM<R>::removeRowSingleton(SPxLPBase<R>& lp,
2300 	      const SVectorBase<R>& row, int& i)
2301 	{
2302 	   assert(row.size() == 1);
2303 	
2304 	   R aij = row.value(0);
2305 	   int  j   = row.index(0);
2306 	   R lo  = R(-infinity);
2307 	   R up  =  R(infinity);
2308 	
2309 	   SPxOut::debug(this, "IMAISM22 row {}: singleton -> val={} lhs={} rhs={}", i, aij, lp.lhs(i),
2310 	                 lp.rhs(i));
2311 	
2312 	   if(GT(aij, R(0.0), this->epsZero()))            // aij > 0
2313 	   {
2314 	      lo = (lp.lhs(i) <= R(-infinity)) ? R(-infinity) : lp.lhs(i) / aij;
2315 	      up = (lp.rhs(i) >=  R(infinity)) ?  R(infinity) : lp.rhs(i) / aij;
2316 	   }
2317 	   else if(LT(aij, R(0.0), this->epsZero()))       // aij < 0
2318 	   {
2319 	      lo = (lp.rhs(i) >=  R(infinity)) ? R(-infinity) : lp.rhs(i) / aij;
2320 	      up = (lp.lhs(i) <= R(-infinity)) ?  R(infinity) : lp.lhs(i) / aij;
2321 	   }
2322 	   else if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol()))
2323 	   {
2324 	      // aij == 0, rhs < 0 or lhs > 0
2325 	      SPxOut::debug(this, " infeasible\n");
2326 	      return this->INFEASIBLE;
2327 	   }
2328 	
2329 	   if(isZero(lo, this->epsZero()))
2330 	      lo = 0.0;
2331 	
2332 	   if(isZero(up, this->epsZero()))
2333 	      up = 0.0;
2334 	
2335 	   SPxOut::debug(this, " removed, lower={} ({}) upper={} ({})\n", lo, lp.lower(j), up, lp.upper(j));
2336 	
2337 	   bool stricterUp = false;
2338 	   bool stricterLo = false;
2339 	
2340 	   R oldLo = lp.lower(j);
2341 	   R oldUp = lp.upper(j);
2342 	
2343 	   if(LTrel(up, lp.upper(j), feastol()))
2344 	   {
2345 	      lp.changeUpper(j, up);
2346 	      stricterUp = true;
2347 	   }
2348 	
2349 	   if(GTrel(lo, lp.lower(j), feastol()))
2350 	   {
2351 	      lp.changeLower(j, lo);
2352 	      stricterLo = true;
2353 	   }
2354 	
2355 	   std::shared_ptr<PostStep> ptr(new RowSingletonPS(lp, i, j, stricterLo, stricterUp,
2356 	                                 lp.lower(j),
2357 	                                 lp.upper(j), oldLo, oldUp, this->_tolerances));
2358 	   m_hist.append(ptr);
2359 	
2360 	   removeRow(lp, i);
2361 	
2362 	   this->m_remRows++;
2363 	   this->m_remNzos++;
2364 	   ++m_stat[SINGLETON_ROW];
2365 	
2366 	   return this->OKAY;
2367 	}
2368 	
2369 	/// aggregate variable x_j to x_j = (rhs - aik * x_k) / aij from row i: aij * x_j + aik * x_k = rhs
2370 	template <class R>
2371 	typename SPxSimplifier<R>::Result SPxMainSM<R>::aggregateVars(SPxLPBase<R>& lp,
2372 	      const SVectorBase<R>& row, int& i)
2373 	{
2374 	   assert(row.size() == 2);
2375 	   assert(EQrel(lp.lhs(i), lp.rhs(i), feastol()));
2376 	
2377 	   R rhs = lp.rhs(i);
2378 	   assert(rhs < R(infinity) && rhs > R(-infinity));
2379 	
2380 	   int j = row.index(0);
2381 	   int k = row.index(1);
2382 	   R aij = row.value(0);
2383 	   R aik = row.value(1);
2384 	   R lower_j = lp.lower(j);
2385 	   R upper_j = lp.upper(j);
2386 	   R lower_k = lp.lower(k);
2387 	   R upper_k = lp.upper(k);
2388 	
2389 	   // fixed variables should be removed by simplifyCols()
2390 	   if(EQrel(lower_j, upper_j, feastol()) || EQrel(lower_k, upper_k, feastol()))
2391 	      return this->OKAY;
2392 	
2393 	   assert(isNotZero(aij, this->epsZero()) && isNotZero(aik, this->epsZero()));
2394 	
2395 	   SPxOut::debug(this, "IMAISM23 row {}: doubleton equation -> {} x_{} + {} x_{} = {}", i, aij, j, aik,
2396 	                 k, rhs);
2397 	
2398 	   // determine which variable can be aggregated without requiring bound tightening of the other variable
2399 	   R new_lo_j;
2400 	   R new_up_j;
2401 	   R new_lo_k;
2402 	   R new_up_k;
2403 	
2404 	   if(aij * aik < 0.0)
2405 	   {
2406 	      // orientation persists
2407 	      new_lo_j = (upper_k >=  R(infinity)) ? R(-infinity) : (rhs - aik * upper_k) / aij;
2408 	      new_up_j = (lower_k <= R(-infinity)) ?  R(infinity) : (rhs - aik * lower_k) / aij;
2409 	      new_lo_k = (upper_j >=  R(infinity)) ? R(-infinity) : (rhs - aij * upper_j) / aik;
2410 	      new_up_k = (lower_j <= R(-infinity)) ?  R(infinity) : (rhs - aij * lower_j) / aik;
2411 	   }
2412 	   else if(aij * aik > 0.0)
2413 	   {
2414 	      // orientation is reversed
2415 	      new_lo_j = (lower_k <= R(-infinity)) ? R(-infinity) : (rhs - aik * lower_k) / aij;
2416 	      new_up_j = (upper_k >=  R(infinity)) ?  R(infinity) : (rhs - aik * upper_k) / aij;
2417 	      new_lo_k = (lower_j <= R(-infinity)) ? R(-infinity) : (rhs - aij * lower_j) / aik;
2418 	      new_up_k = (upper_j >=  R(infinity)) ?  R(infinity) : (rhs - aij * upper_j) / aik;
2419 	   }
2420 	   else
2421 	      throw SPxInternalCodeException("XMAISM12 This should never happen.");
2422 	
2423 	   bool flip_jk = false;
2424 	
2425 	   if(new_lo_j <= R(-infinity) && new_up_j >= R(infinity))
2426 	   {
2427 	      // no bound tightening on x_j when x_k is aggregated
2428 	      flip_jk = true;
2429 	   }
2430 	   else if(new_lo_k <= R(-infinity) && new_up_k >= R(infinity))
2431 	   {
2432 	      // no bound tightening on x_k when x_j is aggregated
2433 	      flip_jk = false;
2434 	   }
2435 	   else if(LE(new_lo_j, lower_j, this->tolerances()->epsilon())
2436 	           && GE(new_up_j, upper_j, this->tolerances()->epsilon()))
2437 	   {
2438 	      if(LE(new_lo_k, lower_k, this->tolerances()->epsilon())
2439 	            && GE(new_up_k, upper_k, this->tolerances()->epsilon()))
2440 	      {
2441 	         // both variables' bounds are not affected by aggregation; choose the better aggregation coeff (aik/aij)
2442 	         if(spxAbs(aij) > spxAbs(aik))
2443 	            flip_jk = false;
2444 	         else
2445 	            flip_jk = true;
2446 	      }
2447 	      else
2448 	         flip_jk = false;
2449 	   }
2450 	   else if(LE(new_lo_k, lower_k, this->tolerances()->epsilon())
2451 	           && GE(new_up_k, upper_k, this->tolerances()->epsilon()))
2452 	   {
2453 	      flip_jk = true;
2454 	   }
2455 	   else
2456 	   {
2457 	      if(spxAbs(aij) > spxAbs(aik))
2458 	         flip_jk = false;
2459 	      else
2460 	         flip_jk = true;
2461 	   }
2462 	
2463 	   if(flip_jk)
2464 	   {
2465 	      int _j = j;
2466 	      R _aij = aij;
2467 	      R _lower_j = lower_j;
2468 	      R _upper_j = upper_j;
2469 	      j = k;
2470 	      k = _j;
2471 	      aij = aik;
2472 	      aik = _aij;
2473 	      lower_j = lower_k;
2474 	      lower_k = _lower_j;
2475 	      upper_j = upper_k;
2476 	      upper_k = _upper_j;
2477 	   }
2478 	
2479 	   const SVectorBase<R>& col_j = lp.colVector(j);
2480 	   const SVectorBase<R>& col_k = lp.colVector(k);
2481 	
2482 	   // aggregation coefficients (x_j = aggr_coef * x_k + aggr_const)
2483 	   R aggr_coef = - (aik / aij);
2484 	   R aggr_const = rhs / aij;
2485 	
2486 	   SPxOut::debug(this, " removed, replacing x_{} with {} + {} * x_{}\n", j, aggr_const, aggr_coef, k);
2487 	
2488 	   // replace all occurrences of x_j
2489 	   for(int r = 0; r < col_j.size(); ++r)
2490 	   {
2491 	      int row_r = col_j.index(r);
2492 	      R arj = col_j.value(r);
2493 	
2494 	      // skip row i
2495 	      if(row_r == i)
2496 	         continue;
2497 	
2498 	      // adapt sides of row r
2499 	      R lhs_r = lp.lhs(row_r);
2500 	      R rhs_r = lp.rhs(row_r);
2501 	
2502 	      if(lhs_r > R(-infinity))
2503 	      {
2504 	         lp.changeLhs(row_r, lhs_r - aggr_const * arj);
2505 	         this->m_chgLRhs++;
2506 	      }
2507 	
2508 	      if(rhs_r < R(infinity))
2509 	      {
2510 	         lp.changeRhs(row_r, rhs_r - aggr_const * arj);
2511 	         this->m_chgLRhs++;
2512 	      }
2513 	
2514 	      R newcoef = aggr_coef * arj;
2515 	      int pos_rk = col_k.pos(row_r);
2516 	
2517 	      // check whether x_k is also present in row r and get its coefficient
2518 	      if(pos_rk >= 0)
2519 	      {
2520 	         R ark = col_k.value(pos_rk);
2521 	         newcoef += ark;
2522 	         this->m_remNzos++;
2523 	      }
2524 	
2525 	      // add new column k to row r or adapt the coefficient a_rk
2526 	      lp.changeElement(row_r, k, newcoef);
2527 	   }
2528 	
2529 	   // adapt objective function
2530 	   R obj_j = lp.obj(j);
2531 	
2532 	   if(isNotZero(obj_j, this->epsZero()))
2533 	   {
2534 	      this->addObjoffset(aggr_const * obj_j);
2535 	      R obj_k = lp.obj(k);
2536 	      lp.changeObj(k, obj_k + aggr_coef * obj_j);
2537 	   }
2538 	
2539 	   // adapt bounds of x_k
2540 	   R scale1 = maxAbs(rhs, aij * upper_j);
2541 	   R scale2 = maxAbs(rhs, aij * lower_j);
2542 	
2543 	   if(scale1 < 1.0)
2544 	      scale1 = 1.0;
2545 	
2546 	   if(scale2 < 1.0)
2547 	      scale2 = 1.0;
2548 	
2549 	   R z1 = (rhs / scale1) - (aij * upper_j / scale1);
2550 	   R z2 = (rhs / scale2) - (aij * lower_j / scale2);
2551 	
2552 	   // just some rounding
2553 	   if(isZero(z1, this->epsZero()))
2554 	      z1 = 0.0;
2555 	
2556 	   if(isZero(z2, this->epsZero()))
2557 	      z2 = 0.0;
2558 	
2559 	   // determine which side has to be used for the bounds comparison below
2560 	   if(aik * aij > 0.0)
2561 	   {
2562 	      new_lo_k = (upper_j >=  R(infinity)) ? R(-infinity) : z1 * scale1 / aik;
2563 	      new_up_k = (lower_j <= R(-infinity)) ?  R(infinity) : z2 * scale2 / aik;
2564 	   }
2565 	   else if(aik * aij < 0.0)
2566 	   {
2567 	      new_lo_k = (lower_j <= R(-infinity)) ? R(-infinity) : z2 * scale2 / aik;
2568 	      new_up_k = (upper_j >=  R(infinity)) ?  R(infinity) : z1 * scale1 / aik;
2569 	   }
2570 	   else
2571 	      throw SPxInternalCodeException("XMAISM12 This should never happen.");
2572 	
2573 	   // change bounds of x_k if the new ones are tighter
2574 	   R oldlower_k = lower_k;
2575 	   R oldupper_k = upper_k;
2576 	
2577 	   if(GT(new_lo_k, lower_k, this->epsZero()))
2578 	   {
2579 	      lp.changeLower(k, new_lo_k);
2580 	      this->m_chgBnds++;
2581 	   }
2582 	
2583 	   if(LT(new_up_k, upper_k, this->epsZero()))
2584 	   {
2585 	      lp.changeUpper(k, new_up_k);
2586 	      this->m_chgBnds++;
2587 	   }
2588 	
2589 	   std::shared_ptr<PostStep> ptr(new AggregationPS(lp, i, j, rhs, oldupper_k, oldlower_k,
2590 	                                 this->_tolerances));
2591 	   m_hist.append(ptr);
2592 	
2593 	   removeRow(lp, i);
2594 	   removeCol(lp, j);
2595 	
2596 	   this->m_remRows++;
2597 	   this->m_remCols++;
2598 	   this->m_remNzos += 2;
2599 	
2600 	   ++m_stat[AGGREGATION];
2601 	
2602 	   return this->OKAY;
2603 	}
2604 	
2605 	template <class R>
2606 	typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyRows(SPxLPBase<R>& lp, bool& again)
2607 	{
2608 	
2609 	   // This method simplifies the rows of the LP.
2610 	   //
2611 	   // The following operations are done:
2612 	   // 1. detect implied free variables
2613 	   // 2. detect implied free constraints
2614 	   // 3. detect infeasible constraints
2615 	   // 4. remove unconstrained constraints
2616 	   // 5. remove empty constraints
2617 	   // 6. remove row singletons and tighten the corresponding variable bounds if necessary
2618 	   // 7. remove doubleton equation, aka aggregation
2619 	   // 8. detect forcing rows and fix the corresponding variables
2620 	
2621 	   int remRows = 0;
2622 	   int remNzos = 0;
2623 	   int chgLRhs = 0;
2624 	   int chgBnds = 0;
2625 	   int keptBnds = 0;
2626 	   int keptLRhs = 0;
2627 	
2628 	   int oldRows = lp.nRows();
2629 	
2630 	   bool redundantLower;
2631 	   bool redundantUpper;
2632 	   bool redundantLhs;
2633 	   bool redundantRhs;
2634 	
2635 	   for(int i = lp.nRows() - 1; i >= 0; --i)
2636 	   {
2637 	      const SVectorBase<R>& row = lp.rowVector(i);
2638 	
2639 	
2640 	      // compute bounds on constraint value
2641 	      R lhsBnd = 0.0; // minimal activity (finite summands)
2642 	      R rhsBnd = 0.0; // maximal activity (finite summands)
2643 	      int  lhsCnt = 0; // number of R(-infinity) summands in minimal activity
2644 	      int  rhsCnt = 0; // number of +R(infinity) summands in maximal activity
2645 	
2646 	      for(int k = 0; k < row.size(); ++k)
2647 	      {
2648 	         R aij = row.value(k);
2649 	         int  j   = row.index(k);
2650 	
2651 	         if(!isNotZero(aij, R(1.0 / R(infinity))))
2652 	         {
2653 	            SPX_MSG_WARNING((*this->spxout), (*this->spxout) << "Warning: tiny nonzero coefficient " << aij <<
2654 	                            " in row " << i << "\n");
2655 	         }
2656 	
2657 	         if(aij > 0.0)
2658 	         {
2659 	            if(lp.lower(j) <= R(-infinity))
2660 	               ++lhsCnt;
2661 	            else
2662 	               lhsBnd += aij * lp.lower(j);
2663 	
2664 	            if(lp.upper(j) >= R(infinity))
2665 	               ++rhsCnt;
2666 	            else
2667 	               rhsBnd += aij * lp.upper(j);
2668 	         }
2669 	         else if(aij < 0.0)
2670 	         {
2671 	            if(lp.lower(j) <= R(-infinity))
2672 	               ++rhsCnt;
2673 	            else
2674 	               rhsBnd += aij * lp.lower(j);
2675 	
2676 	            if(lp.upper(j) >= R(infinity))
2677 	               ++lhsCnt;
2678 	            else
2679 	               lhsBnd += aij * lp.upper(j);
2680 	         }
2681 	      }
2682 	
2683 	#if SOPLEX_FREE_BOUNDS
2684 	
2685 	      // 1. detect implied free variables
2686 	      if(rhsCnt <= 1 || lhsCnt <= 1)
2687 	      {
2688 	         for(int k = 0; k < row.size(); ++k)
2689 	         {
2690 	            R aij = row.value(k);
2691 	            int  j   = row.index(k);
2692 	
2693 	            redundantLower = false;
2694 	            redundantUpper = false;
2695 	
2696 	            SOPLEX_ASSERT_WARN("WMAISM12", isNotZero(aij, R(1.0 / R(infinity))));
2697 	
2698 	            if(aij > 0.0)
2699 	            {
2700 	               if(lp.lhs(i) > R(-infinity) && lp.lower(j) > R(-infinity) && rhsCnt <= 1
2701 	                     && NErel(lp.lhs(i), rhsBnd, feastol())
2702 	                     // do not perform if strongly different orders of magnitude occur
2703 	                     && spxAbs(lp.lhs(i) / maxAbs(rhsBnd, R(1.0))) > this->tolerances()->epsilon())
2704 	               {
2705 	                  R lo    = R(-infinity);
2706 	                  R scale = maxAbs(lp.lhs(i), rhsBnd);
2707 	
2708 	                  if(scale < 1.0)
2709 	                     scale = 1.0;
2710 	
2711 	                  R z = (lp.lhs(i) / scale) - (rhsBnd / scale);
2712 	
2713 	                  if(isZero(z, this->epsZero()))
2714 	                     z = 0.0;
2715 	
2716 	                  assert(rhsCnt > 0 || lp.upper(j) < R(infinity));
2717 	
2718 	                  if(rhsCnt == 0)
2719 	                     lo = lp.upper(j) + z * scale / aij;
2720 	                  else if(lp.upper(j) >= R(infinity))
2721 	                     lo = z * scale / aij;
2722 	
2723 	                  if(isZero(lo, this->epsZero()))
2724 	                     lo = 0.0;
2725 	
2726 	                  if(GErel(lo, lp.lower(j), feastol()))
2727 	                  {
2728 	                     SPxOut::debug(this, "IMAISM13 row {}: redundant lower bound on x{} -> lower={} ({})\n", i, j, lo,
2729 	                                   lp.lower(j));
2730 	
2731 	                     redundantLower = true;
2732 	                  }
2733 	
2734 	               }
2735 	
2736 	               if(lp.rhs(i) < R(infinity) && lp.upper(j) < R(infinity) && lhsCnt <= 1
2737 	                     && NErel(lp.rhs(i), lhsBnd, feastol())
2738 	                     // do not perform if strongly different orders of magnitude occur
2739 	                     && spxAbs(lp.rhs(i) / maxAbs(lhsBnd, R(1.0))) > this->tolerances()->epsilon())
2740 	               {
2741 	                  R up    = R(infinity);
2742 	                  R scale = maxAbs(lp.rhs(i), lhsBnd);
2743 	
2744 	                  if(scale < 1.0)
2745 	                     scale = 1.0;
2746 	
2747 	                  R z = (lp.rhs(i) / scale) - (lhsBnd / scale);
2748 	
2749 	                  if(isZero(z, this->epsZero()))
2750 	                     z = 0.0;
2751 	
2752 	                  assert(lhsCnt > 0 || lp.lower(j) > R(-infinity));
2753 	
2754 	                  if(lhsCnt == 0)
2755 	                     up = lp.lower(j) + z * scale / aij;
2756 	                  else if(lp.lower(j) <= R(-infinity))
2757 	                     up = z * scale / aij;
2758 	
2759 	                  if(isZero(up, this->epsZero()))
2760 	                     up = 0.0;
2761 	
2762 	                  if(LErel(up, lp.upper(j), feastol()))
2763 	                  {
2764 	                     SPxOut::debug(this, "IMAISM14 row {}: redundant upper bound on x{} -> upper={} ({})\n", i, j, up,
2765 	                                   lp.upper(j));
2766 	
2767 	                     redundantUpper = true;
2768 	                  }
2769 	               }
2770 	
2771 	               if(redundantLower)
2772 	               {
2773 	                  // no upper bound on x_j OR redundant upper bound
2774 	                  if((lp.upper(j) >= R(infinity)) || redundantUpper || (!m_keepbounds))
2775 	                  {
2776 	                     ++lhsCnt;
2777 	                     lhsBnd -= aij * lp.lower(j);
2778 	
2779 	                     lp.changeLower(j, R(-infinity));
2780 	                     ++chgBnds;
2781 	                  }
2782 	                  else
2783 	                     ++keptBnds;
2784 	               }
2785 	
2786 	               if(redundantUpper)
2787 	               {
2788 	                  // no lower bound on x_j OR redundant lower bound
2789 	                  if((lp.lower(j) <= R(-infinity)) || redundantLower || (!m_keepbounds))
2790 	                  {
2791 	                     ++rhsCnt;
2792 	                     rhsBnd -= aij * lp.upper(j);
2793 	
2794 	                     lp.changeUpper(j, R(infinity));
2795 	                     ++chgBnds;
2796 	                  }
2797 	                  else
2798 	                     ++keptBnds;
2799 	               }
2800 	            }
2801 	            else if(aij < 0.0)
2802 	            {
2803 	               if(lp.lhs(i) > R(-infinity) && lp.upper(j) < R(infinity) && rhsCnt <= 1
2804 	                     && NErel(lp.lhs(i), rhsBnd, feastol())
2805 	                     // do not perform if strongly different orders of magnitude occur
2806 	                     && spxAbs(lp.lhs(i) / maxAbs(rhsBnd, R(1.0))) > this->tolerances()->epsilon())
2807 	               {
2808 	                  R up    = R(infinity);
2809 	                  R scale = maxAbs(lp.lhs(i), rhsBnd);
2810 	
2811 	                  if(scale < 1.0)
2812 	                     scale = 1.0;
2813 	
2814 	                  R z = (lp.lhs(i) / scale) - (rhsBnd / scale);
2815 	
2816 	                  if(isZero(z, this->epsZero()))
2817 	                     z = 0.0;
2818 	
2819 	                  assert(rhsCnt > 0 || lp.lower(j) > R(-infinity));
2820 	
2821 	                  if(rhsCnt == 0)
2822 	                     up = lp.lower(j) + z * scale / aij;
2823 	                  else if(lp.lower(j) <= R(-infinity))
2824 	                     up = z * scale / aij;
2825 	
2826 	                  if(isZero(up, this->epsZero()))
2827 	                     up = 0.0;
2828 	
2829 	                  if(LErel(up, lp.upper(j), feastol()))
2830 	                  {
2831 	                     SPxOut::debug(this, "IMAISM15 row {}: redundant upper bound on x{} -> upper={} ({})\n", i, j, up,
2832 	                                   lp.upper(j));
2833 	
2834 	                     redundantUpper = true;
2835 	                  }
2836 	               }
2837 	
2838 	               if(lp.rhs(i) < R(infinity) && lp.lower(j) > R(-infinity) && lhsCnt <= 1
2839 	                     && NErel(lp.rhs(i), lhsBnd, feastol())
2840 	                     // do not perform if strongly different orders of magnitude occur
2841 	                     && spxAbs(lp.rhs(i) / maxAbs(lhsBnd, R(1.0))) > this->tolerances()->epsilon())
2842 	               {
2843 	                  R lo    = R(-infinity);
2844 	                  R scale = maxAbs(lp.rhs(i), lhsBnd);
2845 	
2846 	                  if(scale < 1.0)
2847 	                     scale = 1.0;
2848 	
2849 	                  R z = (lp.rhs(i) / scale) - (lhsBnd / scale);
2850 	
2851 	                  if(isZero(z, this->epsZero()))
2852 	                     z = 0.0;
2853 	
2854 	                  assert(lhsCnt > 0 || lp.upper(j) < R(infinity));
2855 	
2856 	                  if(lhsCnt == 0)
2857 	                     lo = lp.upper(j) + z * scale / aij;
2858 	                  else if(lp.upper(j) >= R(infinity))
2859 	                     lo = z * scale / aij;
2860 	
2861 	                  if(isZero(lo, this->epsZero()))
2862 	                     lo = 0.0;
2863 	
2864 	                  if(GErel(lo, lp.lower(j), this->tolerances()->epsilon()))
2865 	                  {
2866 	                     SPxOut::debug(this, "IMAISM16 row {}: redundant lower bound on x{} -> lower={} ({})\n", i, j, lo,
2867 	                                   lp.lower(j));
2868 	
2869 	                     redundantLower = true;
2870 	                  }
2871 	               }
2872 	
2873 	               if(redundantUpper)
2874 	               {
2875 	                  // no lower bound on x_j OR redundant lower bound
2876 	                  if((lp.lower(j) <= R(-infinity)) || redundantLower || (!m_keepbounds))
2877 	                  {
2878 	                     ++lhsCnt;
2879 	                     lhsBnd -= aij * lp.upper(j);
2880 	
2881 	                     lp.changeUpper(j, R(infinity));
2882 	                     ++chgBnds;
2883 	                  }
2884 	                  else
2885 	                     ++keptBnds;
2886 	               }
2887 	
2888 	               if(redundantLower)
2889 	               {
2890 	                  // no upper bound on x_j OR redundant upper bound
2891 	                  if((lp.upper(j) >= R(infinity)) || redundantUpper || (!m_keepbounds))
2892 	                  {
2893 	                     ++rhsCnt;
2894 	                     rhsBnd -= aij * lp.lower(j);
2895 	
2896 	                     lp.changeLower(j, R(-infinity));
2897 	                     ++chgBnds;
2898 	                  }
2899 	                  else
2900 	                     ++keptBnds;
2901 	               }
2902 	            }
2903 	         }
2904 	      }
2905 	
2906 	#endif
2907 	
2908 	#if SOPLEX_FREE_LHS_RHS
2909 	
2910 	      redundantLhs = false;
2911 	      redundantRhs = false;
2912 	
2913 	      // 2. detect implied free constraints
2914 	      if(lp.lhs(i) > R(-infinity) && lhsCnt == 0 && GErel(lhsBnd, lp.lhs(i), feastol()))
2915 	      {
2916 	         SPxOut::debug(this, "IMAISM17 row {}: redundant lhs -> lhsBnd={} lhs={}\n", i, lhsBnd, lp.lhs(i));
2917 	
2918 	         redundantLhs = true;
2919 	      }
2920 	
2921 	      if(lp.rhs(i) <  R(infinity) && rhsCnt == 0 && LErel(rhsBnd, lp.rhs(i), feastol()))
2922 	      {
2923 	         SPxOut::debug(this, "IMAISM18 row {}: redundant rhs -> rhsBnd={} rhs={}\n", i, rhsBnd, lp.rhs(i));
2924 	
2925 	         redundantRhs = true;
2926 	      }
2927 	
2928 	      if(redundantLhs)
2929 	      {
2930 	         // no rhs for constraint i OR redundant rhs
2931 	         if((lp.rhs(i) >= R(infinity)) || redundantRhs || (!m_keepbounds))
2932 	         {
2933 	            lp.changeLhs(i, R(-infinity));
2934 	            ++chgLRhs;
2935 	         }
2936 	         else
2937 	            ++keptLRhs;
2938 	      }
2939 	
2940 	      if(redundantRhs)
2941 	      {
2942 	         // no lhs for constraint i OR redundant lhs
2943 	         if((lp.lhs(i) <= R(-infinity)) || redundantLhs || (!m_keepbounds))
2944 	         {
2945 	            lp.changeRhs(i, R(infinity));
2946 	            ++chgLRhs;
2947 	         }
2948 	         else
2949 	            ++keptLRhs;
2950 	      }
2951 	
2952 	#endif
2953 	
2954 	      // 3. infeasible constraint
2955 	      if(LTrel(lp.rhs(i), lp.lhs(i), feastol())                 ||
2956 	            (LTrel(rhsBnd,   lp.lhs(i), feastol()) && rhsCnt == 0) ||
2957 	            (GTrel(lhsBnd,   lp.rhs(i), feastol()) && lhsCnt == 0))
2958 	      {
2959 	         SPxOut::debug(this, "IMAISM19 row {}: infeasible -> lhs={} rhs={} lhsBnd={} rhsBnd={}\n", i,
2960 	                       lp.lhs(i), lp.rhs(i), lhsBnd, rhsBnd);
2961 	         return this->INFEASIBLE;
2962 	      }
2963 	
2964 	#if SOPLEX_FREE_CONSTRAINT
2965 	
2966 	      // 4. unconstrained constraint
2967 	      if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
2968 	      {
2969 	         SPxOut::debug(this, "IMAISM20 row {}: unconstrained -> removed\n", i);
2970 	
2971 	         std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i, this->_tolerances));
2972 	         m_hist.append(ptr);
2973 	
2974 	         ++remRows;
2975 	         remNzos += row.size();
2976 	         removeRow(lp, i);
2977 	
2978 	         ++m_stat[FREE_ROW];
2979 	
2980 	         continue;
2981 	      }
2982 	
2983 	#endif
2984 	
2985 	#if SOPLEX_EMPTY_CONSTRAINT
2986 	
2987 	      // 5. empty constraint
2988 	      if(row.size() == 0)
2989 	      {
2990 	         SPxOut::debug(this, "IMAISM21 row {}: empty ->", i);
2991 	
2992 	         if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol()))
2993 	         {
2994 	            SPxOut::debug(this, " infeasible lhs={} rhs={}\n", lp.lhs(i), lp.rhs(i));
2995 	            return this->INFEASIBLE;
2996 	         }
2997 	
2998 	         SPxOut::debug(this, " removed\n");
2999 	
3000 	         std::shared_ptr<PostStep> ptr(new EmptyConstraintPS(lp, i, this->_tolerances));
3001 	         m_hist.append(ptr);
3002 	
3003 	         ++remRows;
3004 	         removeRow(lp, i);
3005 	
3006 	         ++m_stat[EMPTY_ROW];
3007 	
3008 	         continue;
3009 	      }
3010 	
3011 	#endif
3012 	
3013 	#if SOPLEX_ROW_SINGLETON
3014 	
3015 	      // 6. row singleton
3016 	      if(row.size() == 1)
3017 	      {
3018 	         removeRowSingleton(lp, row, i);
3019 	         continue;
3020 	      }
3021 	
3022 	#endif
3023 	
3024 	#if SOPLEX_AGGREGATE_VARS
3025 	
3026 	      // 7. row doubleton, aka. simple aggregation of two variables in an equation
3027 	      if(row.size() == 2 && EQrel(lp.lhs(i), lp.rhs(i), feastol()))
3028 	      {
3029 	         aggregateVars(lp, row, i);
3030 	         continue;
3031 	      }
3032 	
3033 	#endif
3034 	
3035 	#if SOPLEX_FORCE_CONSTRAINT
3036 	
3037 	      // 8. forcing constraint (postsolving)
3038 	      // fix variables to obtain the upper bound on constraint value
3039 	      if(rhsCnt == 0 && EQrel(rhsBnd, lp.lhs(i), feastol()))
3040 	      {
3041 	         SPxOut::debug(this, "IMAISM24 row {}: forcing constraint fix on lhs -> lhs={} rhsBnd={}\n", i,
3042 	                       lp.lhs(i), rhsBnd);
3043 	
3044 	         DataArray<bool> fixedCol(row.size());
3045 	         Array<R> lowers(row.size());
3046 	         Array<R> uppers(row.size());
3047 	
3048 	         for(int k = 0; k < row.size(); ++k)
3049 	         {
3050 	            R aij = row.value(k);
3051 	            int  j   = row.index(k);
3052 	
3053 	            fixedCol[k] = !(EQrel(lp.upper(j), lp.lower(j), this->epsZero()));
3054 	
3055 	            lowers[k] = lp.lower(j);
3056 	            uppers[k] = lp.upper(j);
3057 	
3058 	            SOPLEX_ASSERT_WARN("WMAISM25", isNotZero(aij, R(1.0 / R(infinity))));
3059 	
3060 	            if(aij > 0.0)
3061 	               lp.changeLower(j, lp.upper(j));
3062 	            else
3063 	               lp.changeUpper(j, lp.lower(j));
3064 	         }
3065 	
3066 	         std::shared_ptr<PostStep> ptr(new ForceConstraintPS(lp, i, true, fixedCol, lowers,
3067 	                                       uppers, this->_tolerances));
3068 	         m_hist.append(ptr);
3069 	
3070 	         ++remRows;
3071 	         remNzos += row.size();
3072 	         removeRow(lp, i);
3073 	
3074 	         ++m_stat[FORCE_ROW];
3075 	
3076 	         continue;
3077 	      }
3078 	
3079 	      // fix variables to obtain the lower bound on constraint value
3080 	      if(lhsCnt == 0 && EQrel(lhsBnd, lp.rhs(i), feastol()))
3081 	      {
3082 	         SPxOut::debug(this, "IMAISM26 row {}: forcing constraint fix on rhs -> rhs={} lhsBnd={}\n", i,
3083 	                       lp.rhs(i), lhsBnd);
3084 	
3085 	         DataArray<bool> fixedCol(row.size());
3086 	         Array<R> lowers(row.size());
3087 	         Array<R> uppers(row.size());
3088 	
3089 	         for(int k = 0; k < row.size(); ++k)
3090 	         {
3091 	            R aij   = row.value(k);
3092 	            int  j     = row.index(k);
3093 	
3094 	            fixedCol[k] = !(EQrel(lp.upper(j), lp.lower(j), this->epsZero()));
3095 	
3096 	            lowers[k] = lp.lower(j);
3097 	            uppers[k] = lp.upper(j);
3098 	
3099 	            SOPLEX_ASSERT_WARN("WMAISM27", isNotZero(aij, R(1.0 / R(infinity))));
3100 	
3101 	            if(aij > 0.0)
3102 	               lp.changeUpper(j, lp.lower(j));
3103 	            else
3104 	               lp.changeLower(j, lp.upper(j));
3105 	         }
3106 	
3107 	         std::shared_ptr<PostStep> ptr(new ForceConstraintPS(lp, i, false, fixedCol, lowers,
3108 	                                       uppers, this->_tolerances));
3109 	         m_hist.append(ptr);
3110 	
3111 	         ++remRows;
3112 	         remNzos += row.size();
3113 	         removeRow(lp, i);
3114 	
3115 	         ++m_stat[FORCE_ROW];
3116 	
3117 	         continue;
3118 	      }
3119 	
3120 	#endif
3121 	   }
3122 	
3123 	   assert(remRows > 0 || remNzos == 0);
3124 	
3125 	   if(remRows + chgLRhs + chgBnds > 0)
3126 	   {
3127 	      this->m_remRows += remRows;
3128 	      this->m_remNzos += remNzos;
3129 	      this->m_chgLRhs += chgLRhs;
3130 	      this->m_chgBnds += chgBnds;
3131 	      this->m_keptBnds += keptBnds;
3132 	      this->m_keptLRhs += keptLRhs;
3133 	
3134 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (rows) removed "
3135 	                    << remRows << " rows, "
3136 	                    << remNzos << " non-zeros, "
3137 	                    << chgBnds << " col bounds, "
3138 	                    << chgLRhs << " row bounds; kept "
3139 	                    << keptBnds << " column bounds, "
3140 	                    << keptLRhs << " row bounds"
3141 	                    << std::endl;)
3142 	
3143 	      if(remRows > this->m_minReduction * oldRows)
3144 	         again = true;
3145 	   }
3146 	
3147 	   return this->OKAY;
3148 	}
3149 	
3150 	template <class R>
3151 	typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyCols(SPxLPBase<R>& lp, bool& again)
3152 	{
3153 	
3154 	   // This method simplifies the columns of the LP.
3155 	   //
3156 	   // The following operations are done:
3157 	   // 1. detect empty columns and fix corresponding variables
3158 	   // 2. detect variables that are unconstrained from below or above
3159 	   //    and fix corresponding variables or remove involved constraints
3160 	   // 3. fix variables
3161 	   // 4. use column singleton variables with zero objective to adjust constraint bounds
3162 	   // 5. (not free) column singleton combined with doubleton equation are
3163 	   //    used to make the column singleton variable free
3164 	   // 6. substitute (implied) free column singletons
3165 	
3166 	   int remRows = 0;
3167 	   int remCols = 0;
3168 	   int remNzos = 0;
3169 	   int chgBnds = 0;
3170 	
3171 	   int oldCols = lp.nCols();
3172 	   int oldRows = lp.nRows();
3173 	
3174 	   for(int j = lp.nCols() - 1; j >= 0; --j)
3175 	   {
3176 	      const SVectorBase<R>& col = lp.colVector(j);
3177 	
3178 	      // infeasible bounds
3179 	      if(GTrel(lp.lower(j), lp.upper(j), feastol()))
3180 	      {
3181 	         SPxOut::debug(this, "IMAISM29 col {}: infeasible bounds on x{} -> lower={} upper={}\n", j, j,
3182 	                       lp.lower(j), lp.upper(j));
3183 	         return this->INFEASIBLE;
3184 	      }
3185 	
3186 	      // 1. empty column
3187 	      if(col.size() == 0)
3188 	      {
3189 	#if SOPLEX_COMPTY_COLUMN
3190 	         SPxOut::debug(this, "IMAISM30 col {}: empty -> maxObj={} lower={} upper={}\n", j, lp.maxObj(j),
3191 	                       lp.lower(j), lp.upper(j));
3192 	
3193 	         R val;
3194 	
3195 	         if(lp.maxObj(j) > R(0.0))
3196 	         {
3197 	            if(lp.upper(j) >= R(infinity))
3198 	            {
3199 	               SPxOut::debug(this, " unbounded\n");
3200 	               return this->UNBOUNDED;
3201 	            }
3202 	
3203 	            val = lp.upper(j);
3204 	         }
3205 	         else if(lp.maxObj(j) < R(0.0))
3206 	         {
3207 	            if(lp.lower(j) <= R(-infinity))
3208 	            {
3209 	               SPxOut::debug(this, " unbounded\n");
3210 	               return this->UNBOUNDED;
3211 	            }
3212 	
3213 	            val = lp.lower(j);
3214 	         }
3215 	         else
3216 	         {
3217 	            assert(isZero(lp.maxObj(j), this->epsZero()));
3218 	
3219 	            // any value within the bounds is ok
3220 	            if(lp.lower(j) > R(-infinity))
3221 	               val = lp.lower(j);
3222 	            else if(lp.upper(j) < R(infinity))
3223 	               val = lp.upper(j);
3224 	            else
3225 	               val = 0.0;
3226 	         }
3227 	
3228 	         SPxOut::debug(this, " removed\n");
3229 	
3230 	         std::shared_ptr<PostStep> ptr1(new FixBoundsPS(lp, j, val, this->_tolerances));
3231 	         std::shared_ptr<PostStep> ptr2(new FixVariablePS(lp, *this, j, val, this->_tolerances));
3232 	         m_hist.append(ptr1);
3233 	         m_hist.append(ptr2);
3234 	
3235 	         ++remCols;
3236 	         removeCol(lp, j);
3237 	
3238 	         ++m_stat[EMPTY_COL];
3239 	
3240 	         continue;
3241 	#endif
3242 	      }
3243 	
3244 	      if(NErel(lp.lower(j), lp.upper(j), feastol()))
3245 	      {
3246 	         // will be set to false if any constraint implies a bound on the variable
3247 	         bool loFree = true;
3248 	         bool upFree = true;
3249 	
3250 	         // 1. fix and remove variables
3251 	         for(int k = 0; k < col.size(); ++k)
3252 	         {
3253 	            if(!loFree && !upFree)
3254 	               break;
3255 	
3256 	            int i = col.index(k);
3257 	
3258 	            // warn since this unhandled case may slip through unnoticed otherwise
3259 	            SOPLEX_ASSERT_WARN("WMAISM31", isNotZero(col.value(k), R(1.0 / R(infinity))));
3260 	
3261 	            if(col.value(k) > 0.0)
3262 	            {
3263 	               if(lp.rhs(i) <  R(infinity))
3264 	                  upFree = false;
3265 	
3266 	               if(lp.lhs(i) > R(-infinity))
3267 	                  loFree = false;
3268 	            }
3269 	            else if(col.value(k) < 0.0)
3270 	            {
3271 	               if(lp.rhs(i) <  R(infinity))
3272 	                  loFree = false;
3273 	
3274 	               if(lp.lhs(i) > R(-infinity))
3275 	                  upFree = false;
3276 	            }
3277 	         }
3278 	
3279 	         // 2. detect variables that are unconstrained from below or above
3280 	         // max  3 x
3281 	         // s.t. 5 x >= 8
3282 	         if(GT(lp.maxObj(j), R(0.0), this->epsZero()) && upFree)
3283 	         {
3284 	#if SOPLEX_FIX_VARIABLE
3285 	            SPxOut::debug(this, "IMAISM32 col {}: x{} unconstrained above ->", j, j);
3286 	
3287 	            if(lp.upper(j) >= R(infinity))
3288 	            {
3289 	               SPxOut::debug(this, " unbounded\n");
3290 	
3291 	               return this->UNBOUNDED;
3292 	            }
3293 	
3294 	            SPxOut::debug(this, " fixed at upper={}\n", lp.upper(j));
3295 	
3296 	            std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j), this->_tolerances));
3297 	            m_hist.append(ptr);
3298 	            lp.changeLower(j, lp.upper(j));
3299 	         }
3300 	         // max -3 x
3301 	         // s.t. 5 x <= 8
3302 	         else if(LT(lp.maxObj(j), R(0.0), this->epsZero()) && loFree)
3303 	         {
3304 	            SPxOut::debug(this, "IMAISM33 col {}: x{} unconstrained below ->", j, j);
3305 	
3306 	            if(lp.lower(j) <= R(-infinity))
3307 	            {
3308 	               SPxOut::debug(this, " unbounded\n");
3309 	
3310 	               return this->UNBOUNDED;
3311 	            }
3312 	
3313 	            SPxOut::debug(this, " fixed at lower={}\n", lp.lower(j));
3314 	
3315 	            std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j), this->_tolerances));
3316 	            m_hist.append(ptr);
3317 	            lp.changeUpper(j, lp.lower(j));
3318 	#endif
3319 	         }
3320 	         else if(isZero(lp.maxObj(j), this->epsZero()))
3321 	         {
3322 	#if SOPLEX_FREE_ZERO_OBJ_VARIABLE
3323 	            bool unconstrained_below = loFree && lp.lower(j) <= R(-infinity);
3324 	            bool unconstrained_above = upFree && lp.upper(j) >= R(infinity);
3325 	
3326 	            if(unconstrained_below || unconstrained_above)
3327 	            {
3328 	               SPxOut::debug(this, "IMAISM34 col {}: x{} unconstrained {}  with zero objective ({})\n", j, j,
3329 	                             (unconstrained_below ? "below" : "above"), lp.maxObj(j));
3330 	
3331 	               DSVectorBase<R> col_idx_sorted(col);
3332 	
3333 	               // sort col elements by increasing idx
3334 	               IdxCompare compare;
3335 	               SPxQuicksort(col_idx_sorted.mem(), col_idx_sorted.size(), compare);
3336 	
3337 	               std::shared_ptr<PostStep> ptr(new FreeZeroObjVariablePS(lp, j, unconstrained_below,
3338 	                                             col_idx_sorted, this->_tolerances));
3339 	               m_hist.append(ptr);
3340 	
3341 	               // we have to remove the rows with larger idx first, because otherwise the rows are reorder and indices
3342 	               // are out-of-date
3343 	               remRows += col.size();
3344 	
3345 	               for(int k = col_idx_sorted.size() - 1; k >= 0; --k)
3346 	                  removeRow(lp, col_idx_sorted.index(k));
3347 	
3348 	               // remove column
3349 	               removeCol(lp, j);
3350 	
3351 	               // statistics
3352 	               for(int k = 0; k < col.size(); ++k)
3353 	               {
3354 	                  int l   =  col.index(k);
3355 	                  remNzos += lp.rowVector(l).size();
3356 	               }
3357 	
3358 	               ++m_stat[FREE_ZOBJ_COL];
3359 	               ++remCols;
3360 	
3361 	               continue;
3362 	            }
3363 	
3364 	#endif
3365 	         }
3366 	      }
3367 	
3368 	#if SOPLEX_FIX_VARIABLE
3369 	
3370 	      // 3. fix variable
3371 	      if(EQrel(lp.lower(j), lp.upper(j), feastol()))
3372 	      {
3373 	         SPxOut::debug(this, "IMAISM36 col {}: x{} fixed -> lower={} upper={}\n", j, j, lp.lower(j),
3374 	                       lp.upper(j));
3375 	
3376 	         fixColumn(lp, j);
3377 	
3378 	         ++remCols;
3379 	         remNzos += col.size();
3380 	         removeCol(lp, j);
3381 	
3382 	         ++m_stat[FIX_COL];
3383 	
3384 	         continue;
3385 	      }
3386 	
3387 	#endif
3388 	
3389 	      // handle column singletons
3390 	      if(col.size() == 1)
3391 	      {
3392 	         R aij = col.value(0);
3393 	         int  i   = col.index(0);
3394 	
3395 	         // 4. column singleton with zero objective
3396 	         if(isZero(lp.maxObj(j), this->epsZero()))
3397 	         {
3398 	#if SOPLEX_ZERO_OBJ_COL_SINGLETON
3399 	            SPxOut::debug(this, "IMAISM37 col {}: singleton in row {} with zero objective", j, i);
3400 	
3401 	            R lhs = R(-infinity);
3402 	            R rhs = +R(infinity);
3403 	
3404 	            if(GT(aij, R(0.0), this->epsZero()))
3405 	            {
3406 	               if(lp.lhs(i) > R(-infinity) && lp.upper(j) <  R(infinity))
3407 	                  lhs = lp.lhs(i) - aij * lp.upper(j);
3408 	
3409 	               if(lp.rhs(i) <  R(infinity) && lp.lower(j) > R(-infinity))
3410 	                  rhs = lp.rhs(i) - aij * lp.lower(j);
3411 	            }
3412 	            else if(LT(aij, R(0.0), this->epsZero()))
3413 	            {
3414 	               if(lp.lhs(i) > R(-infinity) && lp.lower(j) > R(-infinity))
3415 	                  lhs = lp.lhs(i) - aij * lp.lower(j);
3416 	
3417 	               if(lp.rhs(i) <  R(infinity) && lp.upper(j) <  R(infinity))
3418 	                  rhs = lp.rhs(i) - aij * lp.upper(j);
3419 	            }
3420 	            else
3421 	            {
3422 	               lhs = lp.lhs(i);
3423 	               rhs = lp.rhs(i);
3424 	            }
3425 	
3426 	            if(isZero(lhs, this->epsZero()))
3427 	               lhs = 0.0;
3428 	
3429 	            if(isZero(rhs, this->epsZero()))
3430 	               rhs = 0.0;
3431 	
3432 	            SPxOut::debug(this, " removed -> lhs={} ({}) rhs={} ({})\n", lhs, lp.lhs(i), rhs, lp.rhs(i));
3433 	
3434 	            std::shared_ptr<PostStep> ptr(new ZeroObjColSingletonPS(lp, *this, j, i, this->_tolerances));
3435 	            m_hist.append(ptr);
3436 	
3437 	            lp.changeRange(i, lhs, rhs);
3438 	
3439 	            ++remCols;
3440 	            ++remNzos;
3441 	            removeCol(lp, j);
3442 	
3443 	            ++m_stat[ZOBJ_SINGLETON_COL];
3444 	
3445 	            if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
3446 	            {
3447 	               std::shared_ptr<PostStep> ptr2(new FreeConstraintPS(lp, i, this->_tolerances));
3448 	               m_hist.append(ptr2);
3449 	
3450 	               ++remRows;
3451 	               removeRow(lp, i);
3452 	
3453 	               ++m_stat[FREE_ROW];
3454 	            }
3455 	
3456 	            continue;
3457 	#endif
3458 	         }
3459 	
3460 	         // 5. not free column singleton combined with doubleton equation
3461 	         else if(EQrel(lp.lhs(i), lp.rhs(i), feastol())             &&
3462 	                 lp.rowVector(i).size() == 2                         &&
3463 	                 (lp.lower(j) > R(-infinity) || lp.upper(j) < R(infinity)))
3464 	         {
3465 	#if SOPLEX_DOUBLETON_EQUATION
3466 	            SPxOut::debug(this, "IMAISM38 col {}: singleton in row {} with doubleton equation ->", j, i);
3467 	
3468 	            R lhs = lp.lhs(i);
3469 	
3470 	            const SVectorBase<R>& row = lp.rowVector(i);
3471 	
3472 	            R aik;
3473 	            int  k;
3474 	
3475 	            if(row.index(0) == j)
3476 	            {
3477 	               aik = row.value(1);
3478 	               k   = row.index(1);
3479 	            }
3480 	            else if(row.index(1) == j)
3481 	            {
3482 	               aik = row.value(0);
3483 	               k   = row.index(0);
3484 	            }
3485 	            else
3486 	               throw SPxInternalCodeException("XMAISM11 This should never happen.");
3487 	
3488 	            SOPLEX_ASSERT_WARN("WMAISM39", isNotZero(aik, R(1.0 / R(infinity))));
3489 	
3490 	            R lo, up;
3491 	            R oldLower = lp.lower(k);
3492 	            R oldUpper = lp.upper(k);
3493 	
3494 	            R scale1 = maxAbs(lhs, aij * lp.upper(j));
3495 	            R scale2 = maxAbs(lhs, aij * lp.lower(j));
3496 	
3497 	            if(scale1 < 1.0)
3498 	               scale1 = 1.0;
3499 	
3500 	            if(scale2 < 1.0)
3501 	               scale2 = 1.0;
3502 	
3503 	            R z1 = (lhs / scale1) - (aij * lp.upper(j) / scale1);
3504 	            R z2 = (lhs / scale2) - (aij * lp.lower(j) / scale2);
3505 	
3506 	            if(isZero(z1, this->epsZero()))
3507 	               z1 = 0.0;
3508 	
3509 	            if(isZero(z2, this->epsZero()))
3510 	               z2 = 0.0;
3511 	
3512 	            if(aij * aik > 0.0)
3513 	            {
3514 	               lo = (lp.upper(j) >=  R(infinity)) ? R(-infinity) : z1 * scale1 / aik;
3515 	               up = (lp.lower(j) <= R(-infinity)) ?  R(infinity) : z2 * scale2 / aik;
3516 	            }
3517 	            else if(aij * aik < 0.0)
3518 	            {
3519 	               lo = (lp.lower(j) <= R(-infinity)) ? R(-infinity) : z2 * scale2 / aik;
3520 	               up = (lp.upper(j) >=  R(infinity)) ?  R(infinity) : z1 * scale1 / aik;
3521 	            }
3522 	            else
3523 	               throw SPxInternalCodeException("XMAISM12 This should never happen.");
3524 	
3525 	            if(GTrel(lo, lp.lower(k), this->epsZero()))
3526 	               lp.changeLower(k, lo);
3527 	
3528 	            if(LTrel(up, lp.upper(k), this->epsZero()))
3529 	               lp.changeUpper(k, up);
3530 	
3531 	            SPxOut::debug(this, " made free, bounds on x{}: lower={} ({}) upper={} ({})\n", k, lp.lower(k),
3532 	                          oldLower, lp.upper(k), oldUpper);
3533 	
3534 	            // infeasible bounds
3535 	            if(GTrel(lp.lower(k), lp.upper(k), feastol()))
3536 	            {
3537 	               SPxOut::debug(this, "new bounds are infeasible\n");
3538 	               return this->INFEASIBLE;
3539 	            }
3540 	
3541 	            std::shared_ptr<PostStep> ptr(new DoubletonEquationPS(lp, j, k, i, oldLower, oldUpper,
3542 	                                          this->_tolerances));
3543 	            m_hist.append(ptr);
3544 	
3545 	            if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity))
3546 	               chgBnds += 2;
3547 	            else
3548 	               ++chgBnds;
3549 	
3550 	            lp.changeBounds(j, R(-infinity), R(infinity));
3551 	
3552 	            ++m_stat[DOUBLETON_ROW];
3553 	#endif
3554 	         }
3555 	
3556 	         // 6. (implied) free column singleton
3557 	         if(lp.lower(j) <= R(-infinity) && lp.upper(j) >= R(infinity))
3558 	         {
3559 	#if SOPLEX_FREE_COL_SINGLETON
3560 	            R slackVal = lp.lhs(i);
3561 	
3562 	            // constraint i is an inequality constraint -> transform into equation type
3563 	            if(NErel(lp.lhs(i), lp.rhs(i), feastol()))
3564 	            {
3565 	               SPxOut::debug(this, "IMAISM40 col {}: free singleton in inequality constraint\n", j);
3566 	
3567 	               // do nothing if constraint i is unconstrained
3568 	               if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
3569 	                  continue;
3570 	
3571 	               // introduce slack variable to obtain equality constraint
3572 	               R sMaxObj = lp.maxObj(j) / aij; // after substituting variable j in objective
3573 	               R sLo     = lp.lhs(i);
3574 	               R sUp     = lp.rhs(i);
3575 	
3576 	               if(GT(sMaxObj, R(0.0), this->epsZero()))
3577 	               {
3578 	                  if(sUp >= R(infinity))
3579 	                  {
3580 	                     SPxOut::debug(this, " -> problem unbounded\n");
3581 	                     return this->UNBOUNDED;
3582 	                  }
3583 	
3584 	                  slackVal = sUp;
3585 	               }
3586 	               else if(LT(sMaxObj, R(0.0), this->epsZero()))
3587 	               {
3588 	                  if(sLo <= R(-infinity))
3589 	                  {
3590 	                     SPxOut::debug(this, " -> problem unbounded\n");
3591 	                     return this->UNBOUNDED;
3592 	                  }
3593 	
3594 	                  slackVal = sLo;
3595 	               }
3596 	               else
3597 	               {
3598 	                  assert(isZero(sMaxObj, this->epsZero()));
3599 	
3600 	                  // any value within the bounds is ok
3601 	                  if(sLo > R(-infinity))
3602 	                     slackVal = sLo;
3603 	                  else if(sUp < R(infinity))
3604 	                     slackVal = sUp;
3605 	                  else
3606 	                     throw SPxInternalCodeException("XMAISM13 This should never happen.");
3607 	               }
3608 	            }
3609 	
3610 	            std::shared_ptr<PostStep> ptr(new FreeColSingletonPS(lp, *this, j, i, slackVal, this->_tolerances));
3611 	            m_hist.append(ptr);
3612 	
3613 	            SPxOut::debug(this, "IMAISM41 col {}: free singleton removed\n", j);
3614 	
3615 	            const SVectorBase<R>& row = lp.rowVector(i);
3616 	
3617 	            for(int h = 0; h < row.size(); ++h)
3618 	            {
3619 	               int k = row.index(h);
3620 	
3621 	               if(k != j)
3622 	               {
3623 	                  R new_obj = lp.obj(k) - (lp.obj(j) * row.value(h) / aij);
3624 	                  lp.changeObj(k, new_obj);
3625 	               }
3626 	            }
3627 	
3628 	            ++remRows;
3629 	            ++remCols;
3630 	            remNzos += row.size();
3631 	            removeRow(lp, i);
3632 	            removeCol(lp, j);
3633 	
3634 	            ++m_stat[FREE_SINGLETON_COL];
3635 	
3636 	            continue;
3637 	#endif
3638 	         }
3639 	      }
3640 	   }
3641 	
3642 	   if(remCols + remRows > 0)
3643 	   {
3644 	      this->m_remRows += remRows;
3645 	      this->m_remCols += remCols;
3646 	      this->m_remNzos += remNzos;
3647 	      this->m_chgBnds += chgBnds;
3648 	
3649 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (columns) removed "
3650 	                    << remRows << " rows, "
3651 	                    << remCols << " cols, "
3652 	                    << remNzos << " non-zeros, "
3653 	                    << chgBnds << " col bounds"
3654 	                    << std::endl;)
3655 	
3656 	      if(remCols + remRows > this->m_minReduction * (oldCols + oldRows))
3657 	         again = true;
3658 	   }
3659 	
3660 	   return this->OKAY;
3661 	}
3662 	
3663 	template <class R>
3664 	typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyDual(SPxLPBase<R>& lp, bool& again)
3665 	{
3666 	
3667 	   // This method simplifies LP using the following dual structures:
3668 	   //
3669 	   // 1. dominated columns
3670 	   // 2. weakly dominated columns
3671 	   //
3672 	   // For constructing the dual variables, it is assumed that the objective sense is max
3673 	
3674 	   int remRows = 0;
3675 	   int remCols = 0;
3676 	   int remNzos = 0;
3677 	
3678 	   int oldRows = lp.nRows();
3679 	   int oldCols = lp.nCols();
3680 	
3681 	   DataArray<bool> colSingleton(lp.nCols());
3682 	   VectorBase<R>         dualVarLo(lp.nRows());
3683 	   VectorBase<R>         dualVarUp(lp.nRows());
3684 	   VectorBase<R>         dualConsLo(lp.nCols());
3685 	   VectorBase<R>         dualConsUp(lp.nCols());
3686 	
3687 	   // init
3688 	   for(int i = lp.nRows() - 1; i >= 0; --i)
3689 	   {
3690 	      // check for unconstrained constraints
3691 	      if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
3692 	      {
3693 	         SPxOut::debug(this, "IMAISM43 row {}: unconstrained\n", i);
3694 	
3695 	         std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i, this->_tolerances));
3696 	         m_hist.append(ptr);
3697 	
3698 	         ++remRows;
3699 	         remNzos += lp.rowVector(i).size();
3700 	         removeRow(lp, i);
3701 	
3702 	         ++m_stat[FREE_ROW];
3703 	
3704 	         continue;
3705 	      }
3706 	
3707 	      // corresponds to maximization sense
3708 	      dualVarLo[i] = (lp.lhs(i) <= R(-infinity)) ? 0.0 : R(-infinity);
3709 	      dualVarUp[i] = (lp.rhs(i) >=  R(infinity)) ? 0.0 :  R(infinity);
3710 	   }
3711 	
3712 	   // compute bounds on the dual variables using column singletons
3713 	   for(int j = 0; j < lp.nCols(); ++j)
3714 	   {
3715 	      if(lp.colVector(j).size() == 1)
3716 	      {
3717 	         int  i   = lp.colVector(j).index(0);
3718 	         R aij = lp.colVector(j).value(0);
3719 	
3720 	         SOPLEX_ASSERT_WARN("WMAISM44", isNotZero(aij, R(1.0 / R(infinity))));
3721 	
3722 	         R bound = lp.maxObj(j) / aij;
3723 	
3724 	         if(aij > 0)
3725 	         {
3726 	            if(lp.lower(j) <= R(-infinity) && bound < dualVarUp[i])
3727 	               dualVarUp[i] = bound;
3728 	
3729 	            if(lp.upper(j) >=  R(infinity) && bound > dualVarLo[i])
3730 	               dualVarLo[i] = bound;
3731 	         }
3732 	         else if(aij < 0)
3733 	         {
3734 	            if(lp.lower(j) <= R(-infinity) && bound > dualVarLo[i])
3735 	               dualVarLo[i] = bound;
3736 	
3737 	            if(lp.upper(j) >=  R(infinity) && bound < dualVarUp[i])
3738 	               dualVarUp[i] = bound;
3739 	         }
3740 	      }
3741 	
3742 	   }
3743 	
3744 	   // compute bounds on the dual constraints
3745 	   for(int j = 0; j < lp.nCols(); ++j)
3746 	   {
3747 	      dualConsLo[j] = dualConsUp[j] = 0.0;
3748 	
3749 	      const SVectorBase<R>& col = lp.colVector(j);
3750 	
3751 	      for(int k = 0; k < col.size(); ++k)
3752 	      {
3753 	         if(dualConsLo[j] <= R(-infinity) && dualConsUp[j] >= R(infinity))
3754 	            break;
3755 	
3756 	         R aij = col.value(k);
3757 	         int  i   = col.index(k);
3758 	
3759 	         SOPLEX_ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity))));
3760 	
3761 	         if(aij > 0)
3762 	         {
3763 	            if(dualVarLo[i] <= R(-infinity))
3764 	               dualConsLo[j] = R(-infinity);
3765 	            else
3766 	               dualConsLo[j] += aij * dualVarLo[i];
3767 	
3768 	            if(dualVarUp[i] >= R(infinity))
3769 	               dualConsUp[j] = R(infinity);
3770 	            else
3771 	               dualConsUp[j] += aij * dualVarUp[i];
3772 	         }
3773 	         else if(aij < 0)
3774 	         {
3775 	            if(dualVarLo[i] <= R(-infinity))
3776 	               dualConsUp[j] = R(infinity);
3777 	            else
3778 	               dualConsUp[j] += aij * dualVarLo[i];
3779 	
3780 	            if(dualVarUp[i] >= R(infinity))
3781 	               dualConsLo[j] = R(-infinity);
3782 	            else
3783 	               dualConsLo[j] += aij * dualVarUp[i];
3784 	         }
3785 	      }
3786 	   }
3787 	
3788 	   for(int j = lp.nCols() - 1; j >= 0; --j)
3789 	   {
3790 	      if(lp.colVector(j).size() <= 1)
3791 	         continue;
3792 	
3793 	      // dual infeasibility checks
3794 	      if(LTrel(dualConsUp[j], dualConsLo[j], opttol()))
3795 	      {
3796 	         SPxOut::debug(this, "IMAISM46 col {}: dual infeasible -> dual lhs bound={} dual rhs bound={}\n", j,
3797 	                       dualConsLo[j], dualConsUp[j]);
3798 	         return this->DUAL_INFEASIBLE;
3799 	      }
3800 	
3801 	      R obj = lp.maxObj(j);
3802 	
3803 	      // 1. dominated column
3804 	      // Is the problem really unbounded in the cases below ??? Or is only dual infeasibility be shown
3805 	      if(GTrel(obj, dualConsUp[j], opttol()))
3806 	      {
3807 	#if SOPLEX_DOMINATED_COLUMN
3808 	         SPxOut::debug(this, "IMAISM47 col{}: dominated -> maxObj={} dual rhs bound={}\n", j, obj,
3809 	                       dualConsUp[j]);
3810 	
3811 	         if(lp.upper(j) >= R(infinity))
3812 	         {
3813 	            SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " unbounded" << std::endl;)
3814 	            return this->UNBOUNDED;
3815 	         }
3816 	
3817 	         SPxOut::debug(this, " fixed at upper={}\n", lp.upper(j));
3818 	
3819 	         std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j), this->_tolerances));
3820 	         m_hist.append(ptr);
3821 	         lp.changeLower(j, lp.upper(j));
3822 	
3823 	         ++m_stat[DOMINATED_COL];
3824 	#endif
3825 	      }
3826 	      else if(LTrel(obj, dualConsLo[j], opttol()))
3827 	      {
3828 	#if SOPLEX_DOMINATED_COLUMN
3829 	         SPxOut::debug(this, "IMAISM48 col {}: dominated -> maxObj={} dual lhs bound={}\n", j, obj,
3830 	                       dualConsLo[j]);
3831 	
3832 	         if(lp.lower(j) <= R(-infinity))
3833 	         {
3834 	            SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " unbounded" << std::endl;)
3835 	            return this->UNBOUNDED;
3836 	         }
3837 	
3838 	         SPxOut::debug(this, " fixed at lower={}\n", lp.lower(j));
3839 	
3840 	         std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j), this->_tolerances));
3841 	         m_hist.append(ptr);
3842 	         lp.changeUpper(j, lp.lower(j));
3843 	
3844 	         ++m_stat[DOMINATED_COL];
3845 	#endif
3846 	      }
3847 	
3848 	      // 2. weakly dominated column (no postsolving)
3849 	      else if(lp.upper(j) < R(infinity) && EQrel(obj, dualConsUp[j], opttol()))
3850 	      {
3851 	#if SOPLEX_WEAKLY_DOMINATED_COLUMN
3852 	         SPxOut::debug(this, "IMAISM49 col {}: weakly dominated -> maxObj={} dual rhs bound={}\n", j, obj,
3853 	                       dualConsUp[j]);
3854 	
3855 	         std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j), this->_tolerances));
3856 	         m_hist.append(ptr);
3857 	         lp.changeLower(j, lp.upper(j));
3858 	
3859 	         ++m_stat[WEAKLY_DOMINATED_COL];
3860 	#endif
3861 	      }
3862 	      else if(lp.lower(j) > R(-infinity) && EQrel(obj, dualConsLo[j], opttol()))
3863 	      {
3864 	#if SOPLEX_WEAKLY_DOMINATED_COLUMN
3865 	         SPxOut::debug(this, "IMAISM50 col {}: weakly dominated -> maxObj={} dual lhs bound={}\n", j, obj,
3866 	                       dualConsLo[j]);
3867 	
3868 	         std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j), this->_tolerances));
3869 	         m_hist.append(ptr);
3870 	         lp.changeUpper(j, lp.lower(j));
3871 	
3872 	         ++m_stat[WEAKLY_DOMINATED_COL];
3873 	#endif
3874 	      }
3875 	
3876 	      // fix column
3877 	      if(EQrel(lp.lower(j), lp.upper(j), feastol()))
3878 	      {
3879 	#if SOPLEX_FIX_VARIABLE
3880 	         fixColumn(lp, j);
3881 	
3882 	         ++remCols;
3883 	         remNzos += lp.colVector(j).size();
3884 	         removeCol(lp, j);
3885 	
3886 	         ++m_stat[FIX_COL];
3887 	#endif
3888 	      }
3889 	   }
3890 	
3891 	
3892 	   assert(remRows > 0 || remCols > 0 || remNzos == 0);
3893 	
3894 	   if(remCols + remRows > 0)
3895 	   {
3896 	      this->m_remRows += remRows;
3897 	      this->m_remCols += remCols;
3898 	      this->m_remNzos += remNzos;
3899 	
3900 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (dual) removed "
3901 	                    << remRows << " rows, "
3902 	                    << remCols << " cols, "
3903 	                    << remNzos << " non-zeros"
3904 	                    << std::endl;)
3905 	
3906 	      if(remCols + remRows > this->m_minReduction * (oldCols + oldRows))
3907 	         again = true;
3908 	   }
3909 	
3910 	   return this->OKAY;
3911 	}
3912 	
3913 	
3914 	
3915 	template <class R>
3916 	typename SPxSimplifier<R>::Result SPxMainSM<R>::multiaggregation(SPxLPBase<R>& lp, bool& again)
3917 	{
3918 	   // this simplifier eliminates rows and columns by performing multi aggregations as identified by the constraint
3919 	   // activities.
3920 	   int remRows = 0;
3921 	   int remCols = 0;
3922 	   int remNzos = 0;
3923 	
3924 	   int oldRows = lp.nRows();
3925 	   int oldCols = lp.nCols();
3926 	
3927 	   VectorBase<R> upLocks(lp.nCols());
3928 	   VectorBase<R> downLocks(lp.nCols());
3929 	
3930 	   for(int j = lp.nCols() - 1; j >= 0; --j)
3931 	   {
3932 	      // setting the locks on the variables
3933 	      upLocks[j] = 0;
3934 	      downLocks[j] = 0;
3935 	
3936 	      if(lp.colVector(j).size() <= 1)
3937 	         continue;
3938 	
3939 	      const SVectorBase<R>& col = lp.colVector(j);
3940 	
3941 	      for(int k = 0; k < col.size(); ++k)
3942 	      {
3943 	         R aij = col.value(k);
3944 	
3945 	         SOPLEX_ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity))));
3946 	
3947 	         if(lp.lhs(col.index(k)) > R(-infinity) && lp.rhs(col.index(k)) < R(infinity))
3948 	         {
3949 	            upLocks[j]++;
3950 	            downLocks[j]++;
3951 	         }
3952 	         else if(lp.lhs(col.index(k)) > R(-infinity))
3953 	         {
3954 	            if(aij > 0)
3955 	               downLocks[j]++;
3956 	            else if(aij < 0)
3957 	               upLocks[j]++;
3958 	         }
3959 	         else if(lp.rhs(col.index(k)) < R(infinity))
3960 	         {
3961 	            if(aij > 0)
3962 	               upLocks[j]++;
3963 	            else if(aij < 0)
3964 	               downLocks[j]++;
3965 	         }
3966 	      }
3967 	
3968 	      // multi-aggregate column
3969 	      if(upLocks[j] == 1 || downLocks[j] == 1)
3970 	      {
3971 	         R lower = lp.lower(j);
3972 	         R upper = lp.upper(j);
3973 	         int maxOtherLocks;
3974 	         int bestpos = -1;
3975 	         bool bestislhs = true;
3976 	
3977 	         for(int k = 0; k < col.size(); ++k)
3978 	         {
3979 	            int rowNumber;
3980 	            R lhs;
3981 	            R rhs;
3982 	            bool lhsExists;
3983 	            bool rhsExists;
3984 	            bool aggLhs;
3985 	            bool aggRhs;
3986 	
3987 	            R val = col.value(k);
3988 	
3989 	            rowNumber = col.index(k);
3990 	            lhs = lp.lhs(rowNumber);
3991 	            rhs = lp.rhs(rowNumber);
3992 	
3993 	            if(EQ(lhs, rhs, feastol()))
3994 	               continue;
3995 	
3996 	            lhsExists = lhs > R(-infinity);
3997 	            rhsExists = rhs < R(infinity);
3998 	
3999 	            if(lp.rowVector(rowNumber).size() <= 2)
4000 	               maxOtherLocks = INT_MAX;
4001 	            else if(lp.rowVector(rowNumber).size() == 3)
4002 	               maxOtherLocks = 3;
4003 	            else if(lp.rowVector(rowNumber).size() == 4)
4004 	               maxOtherLocks = 2;
4005 	            else
4006 	               maxOtherLocks = 1;
4007 	
4008 	            aggLhs = lhsExists
4009 	                     && ((col.value(k) > 0.0 && lp.maxObj(j) <= 0.0 && downLocks[j] == 1 && upLocks[j] <= maxOtherLocks)
4010 	                         || (col.value(k) < 0.0 && lp.maxObj(j) >= 0.0 && upLocks[j] == 1 && downLocks[j] <= maxOtherLocks));
4011 	            aggRhs = rhsExists
4012 	                     && ((col.value(k) > 0.0 && lp.maxObj(j) >= 0.0 && upLocks[j] == 1 && downLocks[j] <= maxOtherLocks)
4013 	                         || (col.value(k) < 0.0 && lp.maxObj(j) <= 0.0 && downLocks[j] == 1 && upLocks[j] <= maxOtherLocks));
4014 	
4015 	            if(aggLhs || aggRhs)
4016 	            {
4017 	               R minRes = 0;   // this is the minimum value that the aggregation can attain
4018 	               R maxRes = 0;   // this is the maximum value that the aggregation can attain
4019 	
4020 	               // computing the minimum and maximum residuals if variable j is set to zero.
4021 	               computeMinMaxResidualActivity(lp, rowNumber, j, minRes, maxRes);
4022 	
4023 	               // we will try to aggregate to the lhs
4024 	               if(aggLhs)
4025 	               {
4026 	                  R minVal;
4027 	                  R maxVal;
4028 	
4029 	                  // computing the values of the upper and lower bounds for the aggregated variables
4030 	                  computeMinMaxValues(lp, lhs, val, minRes, maxRes, minVal, maxVal);
4031 	
4032 	                  assert(LE(minVal, maxVal, this->tolerances()->epsilon()));
4033 	
4034 	                  // if the bounds of the aggregation and the original variable are equivalent, then we can reduce
4035 	                  if((minVal > R(-infinity) && GT(minVal, lower, feastol()))
4036 	                        && (maxVal < R(infinity) && LT(maxVal, upper, feastol())))
4037 	                  {
4038 	                     bestpos = col.index(k);
4039 	                     bestislhs = true;
4040 	                     break;
4041 	                  }
4042 	               }
4043 	
4044 	               // we will try to aggregate to the rhs
4045 	               if(aggRhs)
4046 	               {
4047 	                  R minVal;
4048 	                  R maxVal;
4049 	
4050 	                  // computing the values of the upper and lower bounds for the aggregated variables
4051 	                  computeMinMaxValues(lp, rhs, val, minRes, maxRes, minVal, maxVal);
4052 	
4053 	                  assert(LE(minVal, maxVal, this->tolerances()->epsilon()));
4054 	
4055 	                  if((minVal > R(-infinity) && GT(minVal, lower, feastol()))
4056 	                        && (maxVal < R(infinity) && LT(maxVal, upper, feastol())))
4057 	                  {
4058 	                     bestpos = col.index(k);
4059 	                     bestislhs = false;
4060 	                     break;
4061 	                  }
4062 	               }
4063 	            }
4064 	         }
4065 	
4066 	         // it is only possible to aggregate if a best position has been found
4067 	         if(bestpos >= 0)
4068 	         {
4069 	            const SVectorBase<R>& bestRow = lp.rowVector(bestpos);
4070 	            // aggregating the variable and applying the fixings to the all other constraints
4071 	            R aggConstant = (bestislhs ? lp.lhs(bestpos) : lp.rhs(
4072 	                                bestpos));   // this is the lhs or rhs of the aggregated row
4073 	            R aggAij =
4074 	               bestRow[j];                                   // this is the coefficient of the deleted col
4075 	
4076 	            SPxOut::debug(this,
4077 	                          "IMAISM51 col {}: Aggregating row: {} Aggregation Constant={} Coefficient of aggregated col={}\n",
4078 	                          j, bestpos, aggConstant, aggAij);
4079 	
4080 	            std::shared_ptr<PostStep> ptr(new MultiAggregationPS(lp, *this, bestpos, j,
4081 	                                          aggConstant, this->_tolerances));
4082 	            m_hist.append(ptr);
4083 	
4084 	            for(int k = 0; k < col.size(); ++k)
4085 	            {
4086 	               if(col.index(k) != bestpos)
4087 	               {
4088 	                  int rowNumber = col.index(k);
4089 	                  VectorBase<R> updateRow(lp.nCols());
4090 	                  R updateRhs = lp.rhs(col.index(k));
4091 	                  R updateLhs = lp.lhs(col.index(k));
4092 	
4093 	                  updateRow = lp.rowVector(col.index(k));
4094 	
4095 	                  // updating the row with the best row
4096 	                  for(int l = 0; l < bestRow.size(); l++)
4097 	                  {
4098 	                     if(bestRow.index(l) != j)
4099 	                     {
4100 	                        if(lp.rowVector(rowNumber).pos(bestRow.index(l)) >= 0)
4101 	                           lp.changeElement(rowNumber, bestRow.index(l), updateRow[bestRow.index(l)]
4102 	                                            - updateRow[j]*bestRow.value(l) / aggAij);
4103 	                        else
4104 	                           lp.changeElement(rowNumber, bestRow.index(l), -1.0 * updateRow[j]*bestRow.value(l) / aggAij);
4105 	                     }
4106 	                  }
4107 	
4108 	                  // NOTE: I don't know whether we should change the LHS and RHS if they are currently at R(infinity)
4109 	                  if(lp.rhs(rowNumber) < R(infinity))
4110 	                     lp.changeRhs(rowNumber, updateRhs - updateRow[j]*aggConstant / aggAij);
4111 	
4112 	                  if(lp.lhs(rowNumber) > R(-infinity))
4113 	                     lp.changeLhs(rowNumber, updateLhs - updateRow[j]*aggConstant / aggAij);
4114 	
4115 	                  assert(LE(lp.lhs(rowNumber), lp.rhs(rowNumber), this->tolerances()->epsilon()));
4116 	               }
4117 	            }
4118 	
4119 	            for(int l = 0; l < bestRow.size(); l++)
4120 	            {
4121 	               if(bestRow.index(l) != j)
4122 	                  lp.changeMaxObj(bestRow.index(l),
4123 	                                  lp.maxObj(bestRow.index(l)) - lp.maxObj(j)*bestRow.value(l) / aggAij);
4124 	            }
4125 	
4126 	            ++remCols;
4127 	            remNzos += lp.colVector(j).size();
4128 	            removeCol(lp, j);
4129 	            ++remRows;
4130 	            remNzos += lp.rowVector(bestpos).size();
4131 	            removeRow(lp, bestpos);
4132 	
4133 	            ++m_stat[MULTI_AGG];
4134 	         }
4135 	      }
4136 	   }
4137 	
4138 	
4139 	   assert(remRows > 0 || remCols > 0 || remNzos == 0);
4140 	
4141 	   if(remCols + remRows > 0)
4142 	   {
4143 	      this->m_remRows += remRows;
4144 	      this->m_remCols += remCols;
4145 	      this->m_remNzos += remNzos;
4146 	
4147 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (multi-aggregation) removed "
4148 	                    << remRows << " rows, "
4149 	                    << remCols << " cols, "
4150 	                    << remNzos << " non-zeros"
4151 	                    << std::endl;)
4152 	
4153 	      if(remCols + remRows > this->m_minReduction * (oldCols + oldRows))
4154 	         again = true;
4155 	   }
4156 	
4157 	   return this->OKAY;
4158 	}
4159 	
4160 	
4161 	
4162 	template <class R>
4163 	typename SPxSimplifier<R>::Result SPxMainSM<R>::duplicateRows(SPxLPBase<R>& lp, bool& again)
4164 	{
4165 	
4166 	   // This method simplifies the LP by removing duplicate rows
4167 	   // Duplicates are detected using the algorithm of Bixby and Wagner [1987]
4168 	
4169 	   // Possible extension: use generalized definition of duplicate rows according to Andersen and Andersen
4170 	   // However: the resulting sparsification is often very small since the involved rows are usually very sparse
4171 	
4172 	   int remRows = 0;
4173 	   int remNzos = 0;
4174 	
4175 	   int oldRows = lp.nRows();
4176 	
4177 	   // remove empty rows and columns
4178 	   typename SPxSimplifier<R>::Result ret = removeEmpty(lp);
4179 	
4180 	   if(ret != this->OKAY)
4181 	      return ret;
4182 	
4183 	#if SOPLEX_ROW_SINGLETON
4184 	   int rs_remRows = 0;
4185 	
4186 	   for(int i = 0; i < lp.nRows(); ++i)
4187 	   {
4188 	      const SVectorBase<R>& row = lp.rowVector(i);
4189 	
4190 	      if(row.size() == 1)
4191 	      {
4192 	         removeRowSingleton(lp, row, i);
4193 	         rs_remRows++;
4194 	      }
4195 	   }
4196 	
4197 	   if(rs_remRows > 0)
4198 	   {
4199 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) <<
4200 	                    "Simplifier duplicate rows (row singleton stage) removed "
4201 	                    << rs_remRows << " rows, "
4202 	                    << rs_remRows << " non-zeros"
4203 	                    << std::endl;)
4204 	   }
4205 	
4206 	#endif
4207 	
4208 	   if(lp.nRows() < 2)
4209 	      return this->OKAY;
4210 	
4211 	   DataArray<int>    pClass(lp.nRows());           // class of parallel rows
4212 	   DataArray<int>    classSize(lp.nRows());        // size of each class
4213 	   Array<R>   scale(lp.nRows());            // scaling factor for each row
4214 	   int*              idxMem = 0;
4215 	
4216 	   try
4217 	   {
4218 	      spx_alloc(idxMem, lp.nRows());
4219 	   }
4220 	   catch(const SPxMemoryException& x)
4221 	   {
4222 	      spx_free(idxMem);
4223 	      throw x;
4224 	   }
4225 	
4226 	   IdxSet idxSet(lp.nRows(), idxMem);           // set of feasible indices for new pClass
4227 	
4228 	   // init
4229 	   pClass[0]    = 0;
4230 	   scale[0]     = 0.0;
4231 	   classSize[0] = lp.nRows();
4232 	
4233 	   for(int i = 1; i < lp.nRows(); ++i)
4234 	   {
4235 	      pClass[i] = 0;
4236 	      scale[i]  = 0.0;
4237 	      classSize[i] = 0;
4238 	      idxSet.addIdx(i);
4239 	   }
4240 	
4241 	   R oldVal = 0.0;
4242 	
4243 	   // main loop
4244 	   for(int j = 0; j < lp.nCols(); ++j)
4245 	   {
4246 	      const SVectorBase<R>& col = lp.colVector(j);
4247 	
4248 	      for(int k = 0; k < col.size(); ++k)
4249 	      {
4250 	         R aij = col.value(k);
4251 	         int  i   = col.index(k);
4252 	
4253 	         if(scale[i] == 0.0)
4254 	            scale[i] = aij;
4255 	
4256 	         m_classSetRows[pClass[i]].add(i, aij / scale[i]);
4257 	
4258 	         if(--classSize[pClass[i]] == 0)
4259 	            idxSet.addIdx(pClass[i]);
4260 	      }
4261 	
4262 	      // update each parallel class with non-zero column entry
4263 	      for(int m = 0; m < col.size(); ++m)
4264 	      {
4265 	         int k = pClass[col.index(m)];
4266 	
4267 	         if(m_classSetRows[k].size() > 0)
4268 	         {
4269 	            // sort classSet[k] w.r.t. scaled column values
4270 	            ElementCompare compare(this->tolerances()->epsilon());
4271 	
4272 	            if(m_classSetRows[k].size() > 1)
4273 	               SPxQuicksort(m_classSetRows[k].mem(), m_classSetRows[k].size(), compare);
4274 	
4275 	            // use new index first
4276 	            int classIdx = idxSet.index(0);
4277 	            idxSet.remove(0);
4278 	
4279 	            for(int l = 0; l < m_classSetRows[k].size(); ++l)
4280 	            {
4281 	               if(l != 0 && NErel(m_classSetRows[k].value(l), oldVal, this->epsZero()))
4282 	               {
4283 	                  classIdx = idxSet.index(0);
4284 	                  idxSet.remove(0);
4285 	               }
4286 	
4287 	               pClass[m_classSetRows[k].index(l)] = classIdx;
4288 	               ++classSize[classIdx];
4289 	
4290 	               oldVal = m_classSetRows[k].value(l);
4291 	            }
4292 	
4293 	            m_classSetRows[k].clear();
4294 	         }
4295 	      }
4296 	   }
4297 	
4298 	   spx_free(idxMem);
4299 	
4300 	   DataArray<bool> remRow(lp.nRows());
4301 	
4302 	   for(int k = 0; k < lp.nRows(); ++k)
4303 	      m_dupRows[k].clear();
4304 	
4305 	   for(int k = 0; k < lp.nRows(); ++k)
4306 	   {
4307 	      remRow[k] = false;
4308 	      m_dupRows[pClass[k]].add(k, 0.0);
4309 	   }
4310 	
4311 	   const int nRowsOld_tmp = lp.nRows();
4312 	   int* perm_tmp = 0;
4313 	   spx_alloc(perm_tmp, nRowsOld_tmp);
4314 	
4315 	   for(int j = 0; j < nRowsOld_tmp; ++j)
4316 	   {
4317 	      perm_tmp[j] = 0;
4318 	   }
4319 	
4320 	   int idxFirstDupRows = -1;
4321 	   int idxLastDupRows = -1;
4322 	   int numDelRows = 0;
4323 	
4324 	   for(int k = 0; k < lp.nRows(); ++k)
4325 	   {
4326 	      if(m_dupRows[k].size() > 1 && !(lp.rowVector(m_dupRows[k].index(0)).size() == 1))
4327 	      {
4328 	         idxLastDupRows = k;
4329 	
4330 	         if(idxFirstDupRows < 0)
4331 	         {
4332 	            idxFirstDupRows = k;
4333 	         }
4334 	
4335 	         for(int l = 1; l < m_dupRows[k].size(); ++l)
4336 	         {
4337 	            int i = m_dupRows[k].index(l);
4338 	            perm_tmp[i] = -1;
4339 	         }
4340 	
4341 	         numDelRows += (m_dupRows[k].size() - 1);
4342 	      }
4343 	   }
4344 	
4345 	   {
4346 	      int k_tmp, j_tmp = -1;
4347 	
4348 	      for(k_tmp = j_tmp = 0; k_tmp < nRowsOld_tmp; ++k_tmp)
4349 	      {
4350 	         if(perm_tmp[k_tmp] >= 0)
4351 	            perm_tmp[k_tmp] = j_tmp++;
4352 	      }
4353 	   }
4354 	
4355 	   // store rhs and lhs changes for combined update
4356 	   bool doChangeRanges = false;
4357 	   VectorBase<R> newLhsVec(lp.lhs());
4358 	   VectorBase<R> newRhsVec(lp.rhs());
4359 	
4360 	   for(int k = 0; k < lp.nRows(); ++k)
4361 	   {
4362 	      if(m_dupRows[k].size() > 1 && !(lp.rowVector(m_dupRows[k].index(0)).size() == 1))
4363 	      {
4364 	         SPxOut::debug(this, "IMAISM53 {} duplicate rows found\n", m_dupRows[k].size());
4365 	
4366 	         m_stat[DUPLICATE_ROW] += m_dupRows[k].size() - 1;
4367 	
4368 	         // index of one non-column singleton row in dupRows[k]
4369 	         int  rowIdx    = -1;
4370 	         int  maxLhsIdx = -1;
4371 	         int  minRhsIdx = -1;
4372 	         R maxLhs    = R(-infinity);
4373 	         R minRhs    = +R(infinity);
4374 	
4375 	         DataArray<bool> isLhsEqualRhs(m_dupRows[k].size());
4376 	
4377 	         // determine strictest bounds on constraint
4378 	         for(int l = 0; l < m_dupRows[k].size(); ++l)
4379 	         {
4380 	            int i = m_dupRows[k].index(l);
4381 	            isLhsEqualRhs[l] = (lp.lhs(i) == lp.rhs(i));
4382 	
4383 	            SOPLEX_ASSERT_WARN("WMAISM54", isNotZero(scale[i], R(1.0 / R(infinity))));
4384 	
4385 	            if(rowIdx == -1)
4386 	            {
4387 	               rowIdx = i;
4388 	               maxLhs = lp.lhs(rowIdx);
4389 	               minRhs = lp.rhs(rowIdx);
4390 	            }
4391 	            else
4392 	            {
4393 	               R scaledLhs, scaledRhs;
4394 	               R factor = scale[rowIdx] / scale[i];
4395 	
4396 	               if(factor > 0)
4397 	               {
4398 	                  scaledLhs = (lp.lhs(i) <= R(-infinity)) ? R(-infinity) : lp.lhs(i) * factor;
4399 	                  scaledRhs = (lp.rhs(i) >=  R(infinity)) ?  R(infinity) : lp.rhs(i) * factor;
4400 	               }
4401 	               else
4402 	               {
4403 	                  scaledLhs = (lp.rhs(i) >=  R(infinity)) ? R(-infinity) : lp.rhs(i) * factor;
4404 	                  scaledRhs = (lp.lhs(i) <= R(-infinity)) ?  R(infinity) : lp.lhs(i) * factor;
4405 	               }
4406 	
4407 	               if(scaledLhs > maxLhs)
4408 	               {
4409 	                  maxLhs    = scaledLhs;
4410 	                  maxLhsIdx = i;
4411 	               }
4412 	
4413 	               if(scaledRhs < minRhs)
4414 	               {
4415 	                  minRhs    = scaledRhs;
4416 	                  minRhsIdx = i;
4417 	               }
4418 	
4419 	               remRow[i] = true;
4420 	            }
4421 	         }
4422 	
4423 	         if(rowIdx != -1)
4424 	         {
4425 	            R newLhs = (maxLhs > lp.lhs(rowIdx)) ? maxLhs : lp.lhs(rowIdx);
4426 	            R newRhs = (minRhs < lp.rhs(rowIdx)) ? minRhs : lp.rhs(rowIdx);
4427 	
4428 	            if(k == idxLastDupRows)
4429 	            {
4430 	               DataArray<int> da_perm(nRowsOld_tmp);
4431 	
4432 	               for(int j = 0; j < nRowsOld_tmp; ++j)
4433 	               {
4434 	                  da_perm[j] = perm_tmp[j];
4435 	               }
4436 	
4437 	               std::shared_ptr<PostStep> ptr(new DuplicateRowsPS(lp, rowIdx, maxLhsIdx, minRhsIdx,
4438 	                                             m_dupRows[k], scale, da_perm, isLhsEqualRhs, true,
4439 	                                             EQrel(newLhs, newRhs, this->tolerances()->epsilon()), this->_tolerances, k == idxFirstDupRows));
4440 	               m_hist.append(ptr);
4441 	            }
4442 	            else
4443 	            {
4444 	               DataArray<int> da_perm_empty(0);
4445 	               std::shared_ptr<PostStep> ptr(new DuplicateRowsPS(lp, rowIdx, maxLhsIdx, minRhsIdx,
4446 	                                             m_dupRows[k], scale, da_perm_empty, isLhsEqualRhs, false, EQrel(newLhs, newRhs,
4447 	                                                   this->tolerances()->epsilon()),
4448 	                                             this->_tolerances, k == idxFirstDupRows));
4449 	               m_hist.append(ptr);
4450 	            }
4451 	
4452 	            if(maxLhs > lp.lhs(rowIdx) || minRhs < lp.rhs(rowIdx))
4453 	            {
4454 	               // modify lhs and rhs of constraint rowIdx
4455 	               doChangeRanges = true;
4456 	
4457 	               if(LTrel(newRhs, newLhs, feastol()))
4458 	               {
4459 	                  SPxOut::debug(this, "IMAISM55 duplicate rows yield infeasible bounds: lhs={} rhs={}\n", newLhs,
4460 	                                newRhs);
4461 	                  spx_free(perm_tmp);
4462 	                  return this->INFEASIBLE;
4463 	               }
4464 	
4465 	               // if we accept the infeasibility we should clean up the values to avoid problems later
4466 	               if(newRhs < newLhs)
4467 	                  newRhs = newLhs;
4468 	
4469 	               newLhsVec[rowIdx] = newLhs;
4470 	               newRhsVec[rowIdx] = newRhs;
4471 	            }
4472 	         }
4473 	      }
4474 	   }
4475 	
4476 	   // change ranges for all modified constraints by one single call (more efficient)
4477 	   if(doChangeRanges)
4478 	   {
4479 	      lp.changeRange(newLhsVec, newRhsVec);
4480 	   }
4481 	
4482 	   // remove all rows by one single method call (more efficient)
4483 	   const int nRowsOld = lp.nRows();
4484 	   int* perm = 0;
4485 	   spx_alloc(perm, nRowsOld);
4486 	
4487 	   for(int i = 0; i < nRowsOld; ++i)
4488 	   {
4489 	      if(remRow[i])
4490 	      {
4491 	         perm[i] = -1;
4492 	         ++remRows;
4493 	         remNzos += lp.rowVector(i).size();
4494 	      }
4495 	      else
4496 	         perm[i] = 0;
4497 	   }
4498 	
4499 	   lp.removeRows(perm);
4500 	
4501 	   for(int i = 0; i < nRowsOld; ++i)
4502 	   {
4503 	      // assert that the pre-computed permutation was correct
4504 	      assert(perm[i] == perm_tmp[i]);
4505 	
4506 	      // update the global index mapping
4507 	      if(perm[i] >= 0)
4508 	         m_rIdx[perm[i]] = m_rIdx[i];
4509 	   }
4510 	
4511 	   spx_free(perm);
4512 	   spx_free(perm_tmp);
4513 	
4514 	   if(remRows + remNzos > 0)
4515 	   {
4516 	      this->m_remRows += remRows;
4517 	      this->m_remNzos += remNzos;
4518 	
4519 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (duplicate rows) removed "
4520 	                    << remRows << " rows, "
4521 	                    << remNzos << " non-zeros"
4522 	                    << std::endl;)
4523 	
4524 	      if(remRows > this->m_minReduction * oldRows)
4525 	         again = true;
4526 	
4527 	   }
4528 	
4529 	   return this->OKAY;
4530 	}
4531 	
4532 	template <class R>
4533 	typename SPxSimplifier<R>::Result SPxMainSM<R>::duplicateCols(SPxLPBase<R>& lp, bool& again)
4534 	{
4535 	
4536 	   // This method simplifies the LP by removing duplicate columns
4537 	   // Duplicates are detected using the algorithm of Bixby and Wagner [1987]
4538 	
4539 	   int remCols = 0;
4540 	   int remNzos = 0;
4541 	
4542 	   // remove empty rows and columns
4543 	   typename SPxSimplifier<R>::Result ret = removeEmpty(lp);
4544 	
4545 	   if(ret != this->OKAY)
4546 	      return ret;
4547 	
4548 	   if(lp.nCols() < 2)
4549 	      return this->OKAY;
4550 	
4551 	   DataArray<int>    pClass(lp.nCols());          // class of parallel columns
4552 	   DataArray<int>    classSize(lp.nCols());       // size of each class
4553 	   Array<R>   scale(lp.nCols());           // scaling factor for each column
4554 	   int*              idxMem = 0;
4555 	
4556 	   try
4557 	   {
4558 	      spx_alloc(idxMem, lp.nCols());
4559 	   }
4560 	   catch(const SPxMemoryException& x)
4561 	   {
4562 	      spx_free(idxMem);
4563 	      throw x;
4564 	   }
4565 	
4566 	   IdxSet idxSet(lp.nCols(), idxMem);  // set of feasible indices for new pClass
4567 	
4568 	   // init
4569 	   pClass[0]    = 0;
4570 	   scale[0]     = 0.0;
4571 	   classSize[0] = lp.nCols();
4572 	
4573 	   for(int j = 1; j < lp.nCols(); ++j)
4574 	   {
4575 	      pClass[j] = 0;
4576 	      scale[j]  = 0.0;
4577 	      classSize[j] = 0;
4578 	      idxSet.addIdx(j);
4579 	   }
4580 	
4581 	   R oldVal = 0.0;
4582 	
4583 	   // main loop
4584 	   for(int i = 0; i < lp.nRows(); ++i)
4585 	   {
4586 	      const SVectorBase<R>& row = lp.rowVector(i);
4587 	
4588 	      for(int k = 0; k < row.size(); ++k)
4589 	      {
4590 	         R aij = row.value(k);
4591 	         int  j   = row.index(k);
4592 	
4593 	         if(scale[j] == 0.0)
4594 	            scale[j] = aij;
4595 	
4596 	         m_classSetCols[pClass[j]].add(j, aij / scale[j]);
4597 	
4598 	         if(--classSize[pClass[j]] == 0)
4599 	            idxSet.addIdx(pClass[j]);
4600 	      }
4601 	
4602 	      // update each parallel class with non-zero row entry
4603 	      for(int m = 0; m < row.size(); ++m)
4604 	      {
4605 	         int k = pClass[row.index(m)];
4606 	
4607 	         if(m_classSetCols[k].size() > 0)
4608 	         {
4609 	            // sort classSet[k] w.r.t. scaled row values
4610 	            ElementCompare compare(this->tolerances()->epsilon());
4611 	
4612 	            if(m_classSetCols[k].size() > 1)
4613 	               SPxQuicksort(m_classSetCols[k].mem(), m_classSetCols[k].size(), compare);
4614 	
4615 	            // use new index first
4616 	            int classIdx = idxSet.index(0);
4617 	            idxSet.remove(0);
4618 	
4619 	            for(int l = 0; l < m_classSetCols[k].size(); ++l)
4620 	            {
4621 	               if(l != 0 && NErel(m_classSetCols[k].value(l), oldVal, this->epsZero()))
4622 	               {
4623 	                  // start new parallel class
4624 	                  classIdx = idxSet.index(0);
4625 	                  idxSet.remove(0);
4626 	               }
4627 	
4628 	               pClass[m_classSetCols[k].index(l)] = classIdx;
4629 	               ++classSize[classIdx];
4630 	
4631 	               oldVal = m_classSetCols[k].value(l);
4632 	            }
4633 	
4634 	            m_classSetCols[k].clear();
4635 	         }
4636 	      }
4637 	   }
4638 	
4639 	   spx_free(idxMem);
4640 	
4641 	   DataArray<bool> remCol(lp.nCols());
4642 	   DataArray<bool> fixAndRemCol(lp.nCols());
4643 	
4644 	   for(int k = 0; k < lp.nCols(); ++k)
4645 	      m_dupCols[k].clear();
4646 	
4647 	   for(int k = 0; k < lp.nCols(); ++k)
4648 	   {
4649 	      remCol[k] = false;
4650 	      fixAndRemCol[k] = false;
4651 	      m_dupCols[pClass[k]].add(k, 0.0);
4652 	   }
4653 	
4654 	   bool hasDuplicateCol = false;
4655 	   DataArray<int>  m_perm_empty(0);
4656 	
4657 	   for(int k = 0; k < lp.nCols(); ++k)
4658 	   {
4659 	      if(m_dupCols[k].size() > 1 && !(lp.colVector(m_dupCols[k].index(0)).size() == 1))
4660 	      {
4661 	         SPxOut::debug(this, "IMAISM58 {} duplicate columns found\n", m_dupCols[k].size());
4662 	
4663 	         if(!hasDuplicateCol)
4664 	         {
4665 	            std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, 0, 0, 1.0, m_perm_empty, this->_tolerances,
4666 	                                          true));
4667 	            m_hist.append(ptr);
4668 	            hasDuplicateCol = true;
4669 	         }
4670 	
4671 	         for(int l = 0; l < m_dupCols[k].size(); ++l)
4672 	         {
4673 	            for(int m = 0; m < m_dupCols[k].size(); ++m)
4674 	            {
4675 	               int j1  = m_dupCols[k].index(l);
4676 	               int j2  = m_dupCols[k].index(m);
4677 	
4678 	               if(l != m && !remCol[j1] && !remCol[j2])
4679 	               {
4680 	                  R cj1 = lp.maxObj(j1);
4681 	                  R cj2 = lp.maxObj(j2);
4682 	
4683 	                  // A.j1 = factor * A.j2
4684 	                  R factor = scale[j1] / scale[j2];
4685 	                  R objDif = cj1 - cj2 * scale[j1] / scale[j2];
4686 	
4687 	                  SOPLEX_ASSERT_WARN("WMAISM59", isNotZero(factor, this->epsZero()));
4688 	
4689 	                  if(isZero(objDif, this->epsZero()))
4690 	                  {
4691 	                     // case 1: objectives also duplicate
4692 	
4693 	                     // if 0 is not within the column bounds, we are not able to postsolve if the aggregated column has
4694 	                     // status ZERO, hence we skip this case
4695 	                     if(LErel(lp.lower(j1), R(0.0), this->tolerances()->epsilon())
4696 	                           && GErel(lp.upper(j1), R(0.0), this->tolerances()->epsilon())
4697 	                           && LErel(lp.lower(j2), R(0.0), this->tolerances()->epsilon())
4698 	                           && GErel(lp.upper(j2), R(0.0), this->tolerances()->epsilon()))
4699 	                     {
4700 	                        std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, j1, j2, factor, m_perm_empty,
4701 	                                                      this->_tolerances));
4702 	                        // variable substitution xj2' := xj2 + factor * xj1 <=> xj2 = -factor * xj1 + xj2'
4703 	                        m_hist.append(ptr);
4704 	
4705 	                        // update bounds of remaining column j2 (new column j2')
4706 	                        if(factor > 0)
4707 	                        {
4708 	                           if(lp.lower(j2) <= R(-infinity) || lp.lower(j1) <= R(-infinity))
4709 	                              lp.changeLower(j2, R(-infinity));
4710 	                           else
4711 	                              lp.changeLower(j2, lp.lower(j2) + factor * lp.lower(j1));
4712 	
4713 	                           if(lp.upper(j2) >= R(infinity) || lp.upper(j1) >= R(infinity))
4714 	                              lp.changeUpper(j2, R(infinity));
4715 	                           else
4716 	                              lp.changeUpper(j2, lp.upper(j2) + factor * lp.upper(j1));
4717 	                        }
4718 	                        else if(factor < 0)
4719 	                        {
4720 	                           if(lp.lower(j2) <= R(-infinity) || lp.upper(j1) >= R(infinity))
4721 	                              lp.changeLower(j2, R(-infinity));
4722 	                           else
4723 	                              lp.changeLower(j2, lp.lower(j2) + factor * lp.upper(j1));
4724 	
4725 	                           if(lp.upper(j2) >= R(infinity) || lp.lower(j1) <= R(-infinity))
4726 	                              lp.changeUpper(j2, R(infinity));
4727 	                           else
4728 	                              lp.changeUpper(j2, lp.upper(j2) + factor * lp.lower(j1));
4729 	                        }
4730 	
4731 	                        SPxOut::debug(this, "IMAISM60 two duplicate columns {}, {} replaced by one\n", j1, j2);
4732 	
4733 	                        remCol[j1] = true;
4734 	
4735 	                        ++m_stat[SUB_DUPLICATE_COL];
4736 	                     }
4737 	                     else
4738 	                     {
4739 	                        SPxOut::debug(this,
4740 	                                      "IMAISM80 not removing two duplicate columns {}, {} because zero not contained in their bounds\n",
4741 	                                      j1, j2);
4742 	                     }
4743 	                  }
4744 	                  else
4745 	                  {
4746 	                     // case 2: objectives not duplicate
4747 	                     // considered for maximization sense
4748 	                     if(lp.lower(j2) <= R(-infinity))
4749 	                     {
4750 	                        if(factor > 0 && objDif > 0)
4751 	                        {
4752 	                           if(lp.upper(j1) >= R(infinity))
4753 	                           {
4754 	                              SPxOut::debug(this, "IMAISM75 LP unbounded\n");
4755 	                              return this->UNBOUNDED;
4756 	                           }
4757 	
4758 	                           // fix j1 at upper bound
4759 	                           SPxOut::debug(this, "IMAISM61 two duplicate columns {}, {} first one fixed at upper bound={}\n", j1,
4760 	                                         j2, lp.upper(j1));
4761 	
4762 	                           std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.upper(j1), this->_tolerances));
4763 	                           m_hist.append(ptr);
4764 	                           lp.changeLower(j1, lp.upper(j1));
4765 	                        }
4766 	                        else if(factor < 0 && objDif < 0)
4767 	                        {
4768 	                           if(lp.lower(j1) <= R(-infinity))
4769 	                           {
4770 	                              SPxOut::debug(this, "IMAISM76 LP unbounded\n");
4771 	                              return this->UNBOUNDED;
4772 	                           }
4773 	
4774 	                           // fix j1 at lower bound
4775 	                           SPxOut::debug(this, "IMAISM62 two duplicate columns {}, {} first one fixed at lower bound={}\n", j1,
4776 	                                         j2, lp.lower(j1));
4777 	
4778 	                           std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.lower(j1), this->_tolerances));
4779 	                           m_hist.append(ptr);
4780 	                           lp.changeUpper(j1, lp.lower(j1));
4781 	                        }
4782 	                     }
4783 	                     else if(lp.upper(j2) >= R(infinity))
4784 	                     {
4785 	                        // fix j1 at upper bound
4786 	                        if(factor < 0 && objDif > 0)
4787 	                        {
4788 	                           if(lp.upper(j1) >= R(infinity))
4789 	                           {
4790 	                              SPxOut::debug(this, "IMAISM77 LP unbounded\n");
4791 	                              return this->UNBOUNDED;
4792 	                           }
4793 	
4794 	                           // fix j1 at upper bound
4795 	                           SPxOut::debug(this, "IMAISM63 two duplicate columns {}, {},  first one fixed at upper bound={}\n",
4796 	                                         j1, j2, lp.upper(j1));
4797 	
4798 	                           std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.upper(j1), this->_tolerances));
4799 	                           m_hist.append(ptr);
4800 	                           lp.changeLower(j1, lp.upper(j1));
4801 	                        }
4802 	
4803 	                        // fix j1 at lower bound
4804 	                        else if(factor > 0 && objDif < 0)
4805 	                        {
4806 	                           if(lp.lower(j1) <= R(-infinity))
4807 	                           {
4808 	                              SPxOut::debug(this, "IMAISM78 LP unbounded\n");
4809 	                              return this->UNBOUNDED;
4810 	                           }
4811 	
4812 	                           // fix j1 at lower bound
4813 	                           SPxOut::debug(this, "IMAISM64 two duplicate columns {}, {} first one fixed at lower bound={}\n", j1,
4814 	                                         j2, lp.lower(j1));
4815 	
4816 	                           std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.lower(j1), this->_tolerances));
4817 	                           m_hist.append(ptr);
4818 	                           lp.changeUpper(j1, lp.lower(j1));
4819 	                        }
4820 	                     }
4821 	
4822 	                     if(EQrel(lp.lower(j1), lp.upper(j1), feastol()))
4823 	                     {
4824 	                        remCol[j1] = true;
4825 	                        fixAndRemCol[j1] = true;
4826 	
4827 	                        ++m_stat[FIX_DUPLICATE_COL];
4828 	                     }
4829 	                  }
4830 	               }
4831 	            }
4832 	         }
4833 	      }
4834 	   }
4835 	
4836 	   for(int j = 0; j < lp.nCols(); ++j)
4837 	   {
4838 	      if(fixAndRemCol[j])
4839 	      {
4840 	         assert(remCol[j]);
4841 	
4842 	         // correctIdx == false, because the index mapping will be handled by the postsolving in DuplicateColsPS
4843 	         fixColumn(lp, j, false);
4844 	      }
4845 	   }
4846 	
4847 	   // remove all columns by one single method call (more efficient)
4848 	   const int nColsOld = lp.nCols();
4849 	   int* perm = 0;
4850 	   spx_alloc(perm, nColsOld);
4851 	
4852 	   for(int j = 0; j < nColsOld; ++j)
4853 	   {
4854 	      if(remCol[j])
4855 	      {
4856 	         perm[j] = -1;
4857 	         ++remCols;
4858 	         remNzos += lp.colVector(j).size();
4859 	      }
4860 	      else
4861 	         perm[j] = 0;
4862 	   }
4863 	
4864 	   lp.removeCols(perm);
4865 	
4866 	   for(int j = 0; j < nColsOld; ++j)
4867 	   {
4868 	      if(perm[j] >= 0)
4869 	         m_cIdx[perm[j]] = m_cIdx[j];
4870 	   }
4871 	
4872 	   DataArray<int> da_perm(nColsOld);
4873 	
4874 	   for(int j = 0; j < nColsOld; ++j)
4875 	   {
4876 	      da_perm[j] = perm[j];
4877 	   }
4878 	
4879 	   if(hasDuplicateCol)
4880 	   {
4881 	      std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, 0, 0, 1.0, da_perm, this->_tolerances, false,
4882 	                                    true));
4883 	      m_hist.append(ptr);
4884 	   }
4885 	
4886 	   spx_free(perm);
4887 	
4888 	   assert(remCols > 0 || remNzos == 0);
4889 	
4890 	   if(remCols > 0)
4891 	   {
4892 	      this->m_remCols += remCols;
4893 	      this->m_remNzos += remNzos;
4894 	
4895 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (duplicate columns) removed "
4896 	                    << remCols << " cols, "
4897 	                    << remNzos << " non-zeros"
4898 	                    << std::endl;)
4899 	
4900 	      if(remCols > this->m_minReduction * nColsOld)
4901 	         again = true;
4902 	   }
4903 	
4904 	   return this->OKAY;
4905 	}
4906 	
4907 	template <class R>
4908 	void SPxMainSM<R>::fixColumn(SPxLPBase<R>& lp, int j, bool correctIdx)
4909 	{
4910 	
4911 	   assert(EQrel(lp.lower(j), lp.upper(j), feastol()));
4912 	
4913 	   R lo            = lp.lower(j);
4914 	   R up            = lp.upper(j);
4915 	   const SVectorBase<R>& col = lp.colVector(j);
4916 	   R mid           = lo;
4917 	
4918 	   // use the center value between slightly different bounds to improve numerics
4919 	   if(NE(lo, up, this->tolerances()->epsilon()))
4920 	      mid = (up + lo) / 2.0;
4921 	
4922 	   assert(lo < R(infinity) && lo > R(-infinity));
4923 	   assert(up < R(infinity) && up > R(-infinity));
4924 	
4925 	   SPxOut::debug(this, "IMAISM66 fix variable x{}: lower={} upper={} to new value: {}\n", j, lo, up,
4926 	                 mid);
4927 	
4928 	   if(isNotZero(lo, this->epsZero()))
4929 	   {
4930 	      for(int k = 0; k < col.size(); ++k)
4931 	      {
4932 	         int i = col.index(k);
4933 	
4934 	         if(lp.rhs(i) < R(infinity))
4935 	         {
4936 	            R y     = mid * col.value(k);
4937 	            R scale = maxAbs(lp.rhs(i), y);
4938 	
4939 	            if(scale < 1.0)
4940 	               scale = 1.0;
4941 	
4942 	            R rhs = (lp.rhs(i) / scale) - (y / scale);
4943 	
4944 	            if(isZero(rhs, this->epsZero()))
4945 	               rhs = 0.0;
4946 	            else
4947 	               rhs *= scale;
4948 	
4949 	            SPxOut::debug(this, "IMAISM67 row {}: rhs={} ({}) aij={}\n", i, rhs, lp.rhs(i), col.value(k));
4950 	
4951 	            lp.changeRhs(i, rhs);
4952 	         }
4953 	
4954 	         if(lp.lhs(i) > R(-infinity))
4955 	         {
4956 	            R y     = mid * col.value(k);
4957 	            R scale = maxAbs(lp.lhs(i), y);
4958 	
4959 	            if(scale < 1.0)
4960 	               scale = 1.0;
4961 	
4962 	            R lhs = (lp.lhs(i) / scale) - (y / scale);
4963 	
4964 	            if(isZero(lhs, this->epsZero()))
4965 	               lhs = 0.0;
4966 	            else
4967 	               lhs *= scale;
4968 	
4969 	            SPxOut::debug(this, "IMAISM68 row {}: lhs={} ({}) aij={}\n", i, lhs, lp.lhs(i), col.value(k));
4970 	
4971 	            lp.changeLhs(i, lhs);
4972 	         }
4973 	
4974 	         assert(lp.lhs(i) <= lp.rhs(i) + feastol());
4975 	      }
4976 	   }
4977 	
4978 	   std::shared_ptr<PostStep> ptr(new FixVariablePS(lp, *this, j, lp.lower(j), this->_tolerances,
4979 	                                 correctIdx));
4980 	   m_hist.append(ptr);
4981 	}
4982 	
4983 	template <class R>
4984 	typename SPxSimplifier<R>::Result SPxMainSM<R>::simplify(SPxLPBase<R>& lp, Real remainingTime,
4985 	      bool keepbounds, uint32_t seed)
4986 	{
4987 	   // transfer message handler
4988 	   this->spxout = lp.spxout;
4989 	   assert(this->spxout != 0);
4990 	
4991 	   m_thesense = lp.spxSense();
4992 	   this->m_timeUsed->reset();
4993 	   this->m_timeUsed->start();
4994 	
4995 	   this->m_objoffset = 0.0;
4996 	   m_cutoffbound = R(-infinity);
4997 	   m_pseudoobj = R(-infinity);
4998 	
4999 	   this->m_remRows = 0;
5000 	   this->m_remCols = 0;
5001 	   this->m_remNzos = 0;
5002 	   this->m_chgBnds = 0;
5003 	   this->m_chgLRhs = 0;
5004 	   this->m_keptBnds = 0;
5005 	   this->m_keptLRhs = 0;
5006 	
5007 	   m_result     = this->OKAY;
5008 	   bool   again = true;
5009 	   int nrounds = 0;
5010 	
5011 	   if(m_hist.size() > 0)
5012 	   {
5013 	      m_hist.clear();
5014 	   }
5015 	
5016 	   m_hist.reSize(0);
5017 	   m_postsolved = false;
5018 	
5019 	   SPX_MSG_INFO2((*this->spxout),
5020 	                 int numRangedRows = 0;
5021 	                 int numBoxedCols = 0;
5022 	                 int numEqualities = 0;
5023 	
5024 	                 for(int i = 0; i < lp.nRows(); ++i)
5025 	{
5026 	   if(lp.lhs(i) > R(-infinity) && lp.rhs(i) < R(infinity))
5027 	      {
5028 	         if(EQ(lp.lhs(i), lp.rhs(i), this->tolerances()->epsilon()))
5029 	            ++numEqualities;
5030 	         else
5031 	            ++numRangedRows;
5032 	      }
5033 	   }
5034 	   for(int j = 0; j < lp.nCols(); ++j)
5035 	   if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity))
5036 	      ++numBoxedCols;
5037 	
5038 	      (*this->spxout) << "LP has "
5039 	                      << numEqualities << " equations, "
5040 	                      << numRangedRows << " ranged rows, "
5041 	                      << numBoxedCols << " boxed columns"
5042 	                      << std::endl;
5043 	                   )
5044 	
5045 	         m_stat.reSize(17);
5046 	
5047 	   for(int k = 0; k < m_stat.size(); ++k)
5048 	      m_stat[k] = 0;
5049 	
5050 	   m_addedcols = 0;
5051 	   handleRowObjectives(lp);
5052 	
5053 	   m_prim.reDim(lp.nCols());
5054 	   m_slack.reDim(lp.nRows());
5055 	   m_dual.reDim(lp.nRows());
5056 	   m_redCost.reDim(lp.nCols());
5057 	   m_cBasisStat.reSize(lp.nCols());
5058 	   m_rBasisStat.reSize(lp.nRows());
5059 	   m_cIdx.reSize(lp.nCols());
5060 	   m_rIdx.reSize(lp.nRows());
5061 	
5062 	   m_classSetRows.reSize(lp.nRows());
5063 	   m_classSetCols.reSize(lp.nCols());
5064 	   m_dupRows.reSize(lp.nRows());
5065 	   m_dupCols.reSize(lp.nCols());
5066 	
5067 	   m_keepbounds = keepbounds;
5068 	
5069 	   for(int i = 0; i < lp.nRows(); ++i)
5070 	      m_rIdx[i] = i;
5071 	
5072 	   for(int j = 0; j < lp.nCols(); ++j)
5073 	      m_cIdx[j] = j;
5074 	
5075 	   // round extreme values (set all values smaller than this->eps to zero and all values bigger than R(infinity)/5 to R(infinity))
5076 	#if SOPLEX_EXTREMES
5077 	   handleExtremes(lp);
5078 	#endif
5079 	
5080 	   // main presolving loop
5081 	   while(again && m_result == this->OKAY)
5082 	   {
5083 	      nrounds++;
5084 	      SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Round " << nrounds << ":" << std::endl;)
5085 	      again = false;
5086 	
5087 	#if SOPLEX_ROWS_SPXMAINSM
5088 	
5089 	      if(m_result == this->OKAY)
5090 	         m_result = simplifyRows(lp, again);
5091 	
5092 	#endif
5093 	
5094 	#if SOPLEX_COLS_SPXMAINSM
5095 	
5096 	      if(m_result == this->OKAY)
5097 	         m_result = simplifyCols(lp, again);
5098 	
5099 	#endif
5100 	
5101 	      if(again)
5102 	         continue;
5103 	
5104 	#if SOPLEX_DUAL_SPXMAINSM
5105 	
5106 	      if(m_result == this->OKAY)
5107 	         m_result = simplifyDual(lp, again);
5108 	
5109 	#endif
5110 	
5111 	      if(again)
5112 	         continue;
5113 	
5114 	#if SOPLEX_DUPLICATE_ROWS
5115 	
5116 	      if(m_result == this->OKAY)
5117 	         m_result = duplicateRows(lp, again);
5118 	
5119 	#endif
5120 	
5121 	#if SOPLEX_DUPLICATE_COLS
5122 	
5123 	      if(m_result == this->OKAY)
5124 	         m_result = duplicateCols(lp, again);
5125 	
5126 	#endif
5127 	
5128 	      if(!again)
5129 	      {
5130 	#if SOPLEX_TRIVIAL_HEURISTICS
5131 	         trivialHeuristic(lp);
5132 	#endif
5133 	
5134 	#if SOPLEX_PSEUDOOBJ
5135 	         propagatePseudoobj(lp);
5136 	#endif
5137 	
5138 	#if SOPLEX_MULTI_AGGREGATE
5139 	
5140 	         if(m_result == this->OKAY)
5141 	            m_result = multiaggregation(lp, again);
5142 	
5143 	#endif
5144 	      }
5145 	
5146 	   }
5147 	
5148 	   // preprocessing detected infeasibility or unboundedness
5149 	   if(m_result != this->OKAY)
5150 	   {
5151 	      SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier result: " << static_cast<int>
5152 	                    (m_result) << std::endl;)
5153 	      return m_result;
5154 	   }
5155 	
5156 	   this->m_remCols -= m_addedcols;
5157 	   this->m_remNzos -= m_addedcols;
5158 	   SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier removed "
5159 	                 << this->m_remRows << " rows, "
5160 	                 << this->m_remCols << " columns, "
5161 	                 << this->m_remNzos << " nonzeros, "
5162 	                 << this->m_chgBnds << " col bounds, "
5163 	                 << this->m_chgLRhs << " row bounds"
5164 	                 << std::endl;)
5165 	
5166 	   if(keepbounds)
5167 	      SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier kept "
5168 	                    << this->m_keptBnds << " column bounds, "
5169 	                    << this->m_keptLRhs << " row bounds"
5170 	                    << std::endl;)
5171 	
5172 	      SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Reduced LP has "
5173 	                    << lp.nRows() << " rows "
5174 	                    << lp.nCols() << " columns "
5175 	                    << lp.nNzos() << " nonzeros"
5176 	                    << std::endl;)
5177 	
5178 	      SPX_MSG_INFO2((*this->spxout),
5179 	                    int numRangedRows = 0;
5180 	                    int numBoxedCols  = 0;
5181 	                    int numEqualities = 0;
5182 	
5183 	                    for(int i = 0; i < lp.nRows(); ++i)
5184 	   {
5185 	      if(lp.lhs(i) > R(-infinity) && lp.rhs(i) < R(infinity))
5186 	         {
5187 	            if(EQ(lp.lhs(i), lp.rhs(i), this->tolerances()->epsilon()))
5188 	               ++numEqualities;
5189 	            else
5190 	               ++numRangedRows;
5191 	         }
5192 	      }
5193 	   for(int j = 0; j < lp.nCols(); ++j)
5194 	   if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity))
5195 	      ++numBoxedCols;
5196 	
5197 	      (*this->spxout) << "Reduced LP has "
5198 	                      << numEqualities << " equations, "
5199 	                      << numRangedRows << " ranged rows, "
5200 	                      << numBoxedCols << " boxed columns"
5201 	                      << std::endl;
5202 	                   )
5203 	
5204 	         if(lp.nCols() == 0 && lp.nRows() == 0)
5205 	         {
5206 	            SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier removed all rows and columns" <<
5207 	                          std::endl;)
5208 	            m_result = this->VANISHED;
5209 	         }
5210 	
5211 	   SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "\nSimplifier performed:\n"
5212 	                 << m_stat[EMPTY_ROW]            << " empty rows\n"
5213 	                 << m_stat[FREE_ROW]             << " free rows\n"
5214 	                 << m_stat[SINGLETON_ROW]        << " singleton rows\n"
5215 	                 << m_stat[FORCE_ROW]            << " forcing rows\n"
5216 	                 << m_stat[EMPTY_COL]            << " empty columns\n"
5217 	                 << m_stat[FIX_COL]              << " fixed columns\n"
5218 	                 << m_stat[FREE_ZOBJ_COL]        << " free columns with zero objective\n"
5219 	                 << m_stat[ZOBJ_SINGLETON_COL]   << " singleton columns with zero objective\n"
5220 	                 << m_stat[DOUBLETON_ROW]        << " singleton columns combined with a doubleton equation\n"
5221 	                 << m_stat[FREE_SINGLETON_COL]   << " free singleton columns\n"
5222 	                 << m_stat[DOMINATED_COL]        << " dominated columns\n"
5223 	                 << m_stat[WEAKLY_DOMINATED_COL] << " weakly dominated columns\n"
5224 	                 << m_stat[DUPLICATE_ROW]        << " duplicate rows\n"
5225 	                 << m_stat[FIX_DUPLICATE_COL]    << " duplicate columns (fixed)\n"
5226 	                 << m_stat[SUB_DUPLICATE_COL]    << " duplicate columns (substituted)\n"
5227 	                 << m_stat[AGGREGATION]          << " variable aggregations\n"
5228 	                 << m_stat[MULTI_AGG]            << " multi aggregations\n"
5229 	                 << std::endl;);
5230 	
5231 	   this->m_timeUsed->stop();
5232 	
5233 	   return m_result;
5234 	}
5235 	
5236 	template <class R>
5237 	void SPxMainSM<R>::unsimplify(const VectorBase<R>& x, const VectorBase<R>& y,
5238 	                              const VectorBase<R>& s, const VectorBase<R>& r,
5239 	                              const typename SPxSolverBase<R>::VarStatus rows[],
5240 	                              const typename SPxSolverBase<R>::VarStatus cols[], bool isOptimal)
5241 	{
5242 	   SPX_MSG_INFO1((*this->spxout),
5243 	                 (*this->spxout) << " --- unsimplifying solution and basis" << std::endl;)
5244 	   assert(x.dim() <= m_prim.dim());
5245 	   assert(y.dim() <= m_dual.dim());
5246 	   assert(x.dim() == r.dim());
5247 	   assert(y.dim() == s.dim());
5248 	
5249 	   // assign values of variables in reduced LP
5250 	   // NOTE: for maximization problems, we have to switch signs of dual and reduced cost values,
5251 	   // since simplifier assumes minimization problem
5252 	   for(int j = 0; j < x.dim(); ++j)
5253 	   {
5254 	      m_prim[j] = isZero(x[j], this->epsZero()) ? 0.0 : x[j];
5255 	      m_redCost[j] = isZero(r[j], this->epsZero()) ? 0.0 : (m_thesense == SPxLPBase<R>::MAXIMIZE ? -r[j] :
5256 	                     r[j]);
5257 	      m_cBasisStat[j] = cols[j];
5258 	   }
5259 	
5260 	   for(int i = 0; i < y.dim(); ++i)
5261 	   {
5262 	      m_dual[i] = isZero(y[i], this->epsZero()) ? 0.0 : (m_thesense == SPxLPBase<R>::MAXIMIZE ? -y[i] :
5263 	                  y[i]);
5264 	      m_slack[i] = isZero(s[i], this->epsZero()) ? 0.0 : s[i];
5265 	      m_rBasisStat[i] = rows[i];
5266 	   }
5267 	
5268 	   // undo preprocessing
5269 	   for(int k = m_hist.size() - 1; k >= 0; --k)
5270 	   {
5271 	      SPxOut::debug(this, "unsimplifying {} \n", m_hist[k]->getName());
5272 	
5273 	      try
5274 	      {
5275 	         m_hist[k]->execute(m_prim, m_dual, m_slack, m_redCost, m_cBasisStat, m_rBasisStat, isOptimal);
5276 	      }
5277 	      catch(const SPxException& ex)
5278 	      {
5279 	         SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Exception thrown while unsimplifying " <<
5280 	                       m_hist[k]->getName() << ":\n" << ex.what() << "\n");
5281 	         throw SPxInternalCodeException("XMAISM00 Exception thrown during unsimply().");
5282 	      }
5283 	
5284 	      m_hist.reSize(k);
5285 	   }
5286 	
5287 	   // for maximization problems, we have to switch signs of dual and reduced cost values back
5288 	   if(m_thesense == SPxLPBase<R>::MAXIMIZE)
5289 	   {
5290 	      for(int j = 0; j < m_redCost.dim(); ++j)
5291 	         m_redCost[j] = -m_redCost[j];
5292 	
5293 	      for(int i = 0; i < m_dual.dim(); ++i)
5294 	         m_dual[i] = -m_dual[i];
5295 	   }
5296 	
5297 	   if(m_addedcols > 0)
5298 	   {
5299 	      assert(m_prim.dim() >= m_addedcols);
5300 	      m_prim.reDim(m_prim.dim() - m_addedcols);
5301 	      m_redCost.reDim(m_redCost.dim() - m_addedcols);
5302 	      m_cBasisStat.reSize(m_cBasisStat.size() - m_addedcols);
5303 	      m_cIdx.reSize(m_cIdx.size() - m_addedcols);
5304 	   }
5305 	
5306 	#ifdef SOPLEX_CHECK_BASIS_DIM
5307 	   int numBasis = 0;
5308 	
5309 	   for(int rs = 0; rs < m_rBasisStat.size(); ++rs)
5310 	   {
5311 	      if(m_rBasisStat[rs] == SPxSolverBase<R>::BASIC)
5312 	         numBasis ++;
5313 	   }
5314 	
5315 	   for(int cs = 0; cs < m_cBasisStat.size(); ++cs)
5316 	   {
5317 	      if(m_cBasisStat[cs] == SPxSolverBase<R>::BASIC)
5318 	         numBasis ++;
5319 	   }
5320 	
5321 	   if(numBasis != m_rBasisStat.size())
5322 	   {
5323 	      throw SPxInternalCodeException("XMAISM26 Dimension doesn't match after this step.");
5324 	   }
5325 	
5326 	#endif
5327 	
5328 	   m_hist.clear();
5329 	   m_postsolved = true;
5330 	}
5331 	
5332 	// Pretty-printing of solver status.
5333 	template <class R>
5334 	std::ostream& operator<<(std::ostream& os, const typename SPxSimplifier<R>::Result& status)
5335 	{
5336 	   switch(status)
5337 	   {
5338 	   case SPxSimplifier<R>::OKAY:
5339 	      os << "SUCCESS";
5340 	      break;
5341 	
5342 	   case SPxSimplifier<R>::INFEASIBLE:
5343 	      os << "INFEASIBLE";
5344 	      break;
5345 	
5346 	   case SPxSimplifier<R>::DUAL_INFEASIBLE:
5347 	      os << "DUAL_INFEASIBLE";
5348 	      break;
5349 	
5350 	   case SPxSimplifier<R>::UNBOUNDED:
5351 	      os << "UNBOUNDED";
5352 	      break;
5353 	
5354 	   case SPxSimplifier<R>::VANISHED:
5355 	      os << "VANISHED";
5356 	      break;
5357 	
5358 	   default:
5359 	      os << "UNKNOWN";
5360 	      break;
5361 	   }
5362 	
5363 	   return os;
5364 	}
5365 	
5366 	} //namespace soplex
5367