Scippy

SoPlex

Sequential object-oriented simPlex

spxmainsm.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2016 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 #include <iostream>
17 
18 #include "spxmainsm.h"
19 #include "array.h"
20 #include "dataarray.h"
21 #include "sorter.h"
22 #include "spxout.h"
23 #include <sstream>
24 #include <iostream>
25 #include <fstream>
26 
27 
28 //rows
29 #define FREE_LHS_RHS 1
30 #define FREE_CONSTRAINT 1
31 #define EMPTY_CONSTRAINT 1
32 #define ROW_SINGLETON 1
33 #define FORCE_CONSTRAINT 1
34 //cols
35 #define FREE_BOUNDS 1
36 #define EMPTY_COLUMN 1
37 #define FIX_VARIABLE 1
38 #define FREE_ZERO_OBJ_VARIABLE 1
39 #define ZERO_OBJ_COL_SINGLETON 1
40 #define DOUBLETON_EQUATION 1
41 #define FREE_COL_SINGLETON 1
42 //dual
43 #define DOMINATED_COLUMN 1
44 #define WEAKLY_DOMINATED_COLUMN 1
45 
46 
47 #define EXTREMES 1
48 #define ROWS 1
49 #define COLS 1
50 #define DUAL 1
51 ///@todo check: with this simplification step, the unsimplified basis seems to be slightly suboptimal for some instances
52 #define DUPLICATE_ROWS 1
53 #define DUPLICATE_COLS 1
54 
55 
56 #ifndef NDEBUG
57 #define CHECK_BASIC_DIM
58 #endif // NDEBUG
59 
60 namespace soplex
61 {
64 {
65  int numBasis = 0;
66  for(int rs = 0; rs < nRows; ++rs)
67  {
68  if(rows[rs] == SPxSolver::BASIC)
69  numBasis++;
70  }
71 
72  for(int cs = 0; cs < nCols; ++cs)
73  {
74  if(cols[cs] == SPxSolver::BASIC)
75  numBasis++;
76  }
77  assert(numBasis == nRows);
78  return numBasis == nRows;
79 }
80 
83  DataArray<SPxSolver::VarStatus>& rStatus) const
84 {
85  assert(isZero(s[m_i], 1e-9));
86  s[m_i] = -x[m_j];
87 
88  assert(rStatus[m_i] != SPxSolver::UNDEFINED);
89  assert(cStatus[m_j] != SPxSolver::UNDEFINED);
90  assert(rStatus[m_i] != SPxSolver::BASIC || cStatus[m_j] != SPxSolver::BASIC);
91 
92  MSG_DEBUG( std::cout << "RowObjPS: removing slack column " << m_j << " (" << cStatus[m_j] << ") for row " << m_i << " (" << rStatus[m_i] << ").\n" );
93 
94  if( rStatus[m_i] != SPxSolver::BASIC )
95  {
96  switch( cStatus[m_j] )
97  {
99  rStatus[m_i] = SPxSolver::ON_LOWER;
100  break;
101  case SPxSolver::ON_LOWER:
102  rStatus[m_i] = SPxSolver::ON_UPPER;
103  break;
104  default:
105  rStatus[m_i] = cStatus[m_j];
106  }
107 
108  // otherwise checkBasisDim() may fail
109  cStatus[m_j] = SPxSolver::ZERO;
110  }
111 
112 #ifdef CHECK_BASIC_DIM
113  if (!checkBasisDim(rStatus, cStatus))
114  {
115  assert(false);
116  throw SPxInternalCodeException("XMAISM15 Dimension doesn't match after this step.");
117  }
118 #endif
119 }
120 
123  DataArray<SPxSolver::VarStatus>& rStatus) const
124 {
125  // correcting the change of idx by deletion of the row:
126  s[m_old_i] = s[m_i];
127  y[m_old_i] = y[m_i];
128  rStatus[m_old_i] = rStatus[m_i];
129 
130  // primal:
131  Real slack = 0.0;
132 
133  for (int k = 0; k < m_row.size(); ++k)
134  slack += m_row.value(k) * x[m_row.index(k)];
135 
136  s[m_i] = slack;
137 
138  // dual:
139  y[m_i] = m_row_obj;
140 
141  // basis:
142  rStatus[m_i] = SPxSolver::BASIC;
143 
144 #ifdef CHECK_BASIC_DIM
145  if (!checkBasisDim(rStatus, cStatus))
146  {
147  throw SPxInternalCodeException("XMAISM15 Dimension doesn't match after this step.");
148  }
149 #endif
150 }
151 
154  DataArray<SPxSolver::VarStatus>& rStatus) const
155 {
156  // correcting the change of idx by deletion of the row:
157  s[m_old_i] = s[m_i];
158  y[m_old_i] = y[m_i];
159  rStatus[m_old_i] = rStatus[m_i];
160 
161  // primal:
162  s[m_i] = 0.0;
163 
164  // dual:
165  y[m_i] = m_row_obj;
166 
167  // basis:
168  rStatus[m_i] = SPxSolver::BASIC;
169 
170 #ifdef CHECK_BASIC_DIM
171  if (!checkBasisDim(rStatus, cStatus))
172  {
173  throw SPxInternalCodeException("XMAISM16 Dimension doesn't match after this step.");
174  }
175 #endif
176 }
177 
180  DataArray<SPxSolver::VarStatus>& rStatus) const
181 {
182  // correcting the change of idx by deletion of the row:
183  s[m_old_i] = s[m_i];
184  y[m_old_i] = y[m_i];
185  rStatus[m_old_i] = rStatus[m_i];
186 
187  Real aij = m_col[m_i];
188 
189  // primal:
190  s[m_i] = aij * x[m_j];
191 
192  // dual & basis:
193  Real val = m_obj;
194 
195  for(int k = 0; k < m_col.size(); ++k)
196  {
197  if (m_col.index(k) != m_i)
198  val -= m_col.value(k) * y[m_col.index(k)];
199  }
200 
201  Real newLo = (aij > 0) ? m_lhs/aij : m_rhs/aij; // implicit lhs
202  Real newUp = (aij > 0) ? m_rhs/aij : m_lhs/aij; // implicit rhs
203 
204  switch(cStatus[m_j])
205  {
206  case SPxSolver::FIXED:
207  if(newLo <= m_oldLo && newUp >= m_oldUp)
208  {
209  // this row is totally redundant, has not changed bound of xj
210  rStatus[m_i] = SPxSolver::BASIC;
211  y[m_i] = m_row_obj;
212  }
213  else if(EQrel(newLo, newUp, eps()))
214  {
215  // row is in the type aij * xj = b
216  assert(EQrel(newLo, x[m_j], eps()));
217 
218  if(EQrel(m_oldLo, m_oldUp, eps()))
219  {
220  // xj has been fixed in other row
221  rStatus[m_i] = SPxSolver::BASIC;
222  y[m_i] = m_row_obj;
223  }
224  else if((EQrel(m_oldLo, x[m_j], eps()) && r[m_j] <= -eps())
225  || (EQrel(m_oldUp, x[m_j], eps()) && r[m_j] >= eps())
226  || (!EQrel(m_oldLo, x[m_j], eps()) && !(EQrel(m_oldUp, x[m_j], eps()))))
227  {
228  // 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
229  rStatus[m_i] = (EQrel(m_lhs, x[m_j]*aij, eps())) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
230  cStatus[m_j] = SPxSolver::BASIC;
231  y[m_i] = val / aij;
232  r[m_j] = 0.0;
233  }
234  else
235  {
236  // set x_j on one of the bound
237  assert(EQrel(m_oldLo, x[m_j], eps()) || EQrel(m_oldUp, x[m_j], eps()));
238 
239  cStatus[m_j] = EQrel(m_oldLo, x[m_j], eps()) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
240  rStatus[m_i] = SPxSolver::BASIC;
241  y[m_i] = m_row_obj;
242  r[m_j] = val;
243  }
244  }
245  else if(EQrel(newLo, m_oldUp, eps()))
246  {
247  // row is in the type xj >= b/aij, try to set xj on upper
248  if(r[m_j] >= eps())
249  {
250  // the reduced cost is positive, xj should in the basic
251  assert(EQrel(m_rhs, x[m_j]*aij, eps()) || EQrel(m_lhs, x[m_j]*aij, eps()));
252 
253  rStatus[m_i] = (EQrel(m_lhs, x[m_j]*aij, eps())) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
254  cStatus[m_j] = SPxSolver::BASIC;
255  y[m_i] = val / aij;
256  r[m_j] = 0.0;
257  }
258  else
259  {
260  assert(EQrel(m_oldUp, x[m_j], eps()));
261 
262  cStatus[m_j] = SPxSolver::ON_UPPER;
263  rStatus[m_i] = SPxSolver::BASIC;
264  y[m_i] = m_row_obj;
265  r[m_j] = val;
266  }
267  }
268  else if(EQrel(newUp, m_oldLo, eps()))
269  {
270  // row is in the type xj <= b/aij, try to set xj on lower
271  if(r[m_j] <= -eps())
272  {
273  // the reduced cost is negative, xj should in the basic
274  assert(EQrel(m_rhs, x[m_j]*aij, eps()) || EQrel(m_lhs, x[m_j]*aij, eps()));
275 
276  rStatus[m_i] = (EQrel(m_lhs, x[m_j]*aij, eps())) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
277  cStatus[m_j] = SPxSolver::BASIC;
278  y[m_i] = val / aij;
279  r[m_j] = 0.0;
280  }
281  else
282  {
283  assert(EQrel(m_oldLo, x[m_j], eps()));
284 
285  cStatus[m_j] = SPxSolver::ON_LOWER;
286  rStatus[m_i] = SPxSolver::BASIC;
287  y[m_i] = m_row_obj;
288  r[m_j] = val;
289  }
290  }
291  else
292  {
293  // the variable is set to FIXED by other constraints, i.e., this singleton row is redundant
294  rStatus[m_i] = SPxSolver::BASIC;
295  y[m_i] = m_row_obj;
296  }
297  break;
298  case SPxSolver::BASIC:
299  rStatus[m_i] = SPxSolver::BASIC;
300  y[m_i] = m_row_obj;
301  r[m_j] = 0.0;
302  break;
303  case SPxSolver::ON_LOWER:
304  if(EQrel(m_oldLo, x[m_j], eps())) // xj may stay on lower
305  {
306  rStatus[m_i] = SPxSolver::BASIC;
307  y[m_i] = m_row_obj;
308  r[m_j] = val;
309  }
310  else // if reduced costs are negative or old lower bound not equal to xj, we need to change xj into the basis
311  {
312  assert(EQrel(m_rhs, x[m_j]*aij, eps()) || EQrel(m_lhs, x[m_j]*aij, eps()));
313 
314  cStatus[m_j] = SPxSolver::BASIC;
315  rStatus[m_i] = (EQrel(m_lhs, x[m_j]*aij, eps())) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
316  y[m_i] = val / aij;
317  r[m_j] = 0.0;
318  }
319  break;
320  case SPxSolver::ON_UPPER:
321  if(EQrel(m_oldUp, x[m_j], eps())) // xj may stay on upper
322  {
323  rStatus[m_i] = SPxSolver::BASIC;
324  y[m_i] = m_row_obj;
325  r[m_j] = val;
326  }
327  else // if reduced costs are positive or old upper bound not equal to xj, we need to change xj into the basis
328  {
329  assert(EQrel(m_rhs, x[m_j]*aij, eps()) || EQrel(m_lhs, x[m_j]*aij, eps()));
330 
331  cStatus[m_j] = SPxSolver::BASIC;
332  rStatus[m_i] = (EQrel(m_lhs, x[m_j]*aij, eps())) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
333  y[m_i] = val / aij;
334  r[m_j] = 0.0;
335  }
336  break;
337  case SPxSolver::ZERO:
338  rStatus[m_i] = SPxSolver::BASIC;
339  y[m_i] = m_row_obj;
340  r[m_j] = val;
341  break;
342  default:
343  break;
344  }
345 
346 #ifdef CHECK_BASIC_DIM
347  if (!checkBasisDim(rStatus, cStatus))
348  {
349  throw SPxInternalCodeException("XMAISM17 Dimension doesn't match after this step.");
350  }
351 #endif
352 }
353 
356  DataArray<SPxSolver::VarStatus>& rStatus) const
357 {
358  // correcting the change of idx by deletion of the row:
359  s[m_old_i] = s[m_i];
360  y[m_old_i] = y[m_i];
361  rStatus[m_old_i] = rStatus[m_i];
362 
363  // primal:
364  s[m_i] = m_lRhs;
365 
366  // basis:
367  int cBasisCandidate = -1;
368  Real maxViolation = -1.0;
369  int bas_k = -1;
370 
371  for(int k = 0; k < m_row.size(); ++k)
372  {
373  int cIdx = m_row.index(k);
374  Real aij = m_row.value(k);
375  Real oldLo = m_oldLowers[k];
376  Real oldUp = m_oldUppers[k];
377 
378  switch(cStatus[cIdx])
379  {
380  case SPxSolver::FIXED:
381  if(m_fixed[k])
382  {
383  assert(EQrel(oldLo, x[cIdx], eps()) || EQrel(oldUp, x[cIdx], eps()));
384 
385  Real violation = spxAbs(r[cIdx]/aij);
386 
387  cStatus[cIdx] = EQrel(oldLo, x[cIdx], eps()) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
388 
389  if( violation > maxViolation && ( (EQrel(oldLo, x[cIdx], eps()) && r[cIdx] < -eps()) || (EQrel(oldUp, x[cIdx], eps()) && r[cIdx] > eps()) ) )
390  {
391  maxViolation = violation;
392  cBasisCandidate = cIdx;
393  bas_k = k;
394  }
395  } // do nothing, if the old bounds are equal, i.e. variable has been not fixed in this row
396  break;
397  case SPxSolver::ON_LOWER:
398  case SPxSolver::ON_UPPER:
399  case SPxSolver::BASIC:
400  break;
401  default:
402  break;
403  }
404  }
405 
406  // dual and basis :
407  if(cBasisCandidate >= 0) // one of the variable in the row should in the basis
408  {
409  assert(EQrel(m_lRhs, m_rhs, eps()) || EQrel(m_lRhs, m_lhs, eps()));
410  assert(bas_k >= 0);
411  assert(cBasisCandidate == m_row.index(bas_k));
412 
413  cStatus[cBasisCandidate] = SPxSolver::BASIC;
414  rStatus[m_i] = (EQrel(m_lRhs, m_lhs, eps())) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
415 
416  Real aij = m_row.value(bas_k);
417  Real multiplier = r[cBasisCandidate]/aij;
418  r[cBasisCandidate] = 0.0;
419 
420  for(int k = 0; k < m_row.size(); ++k) // update the reduced cost
421  {
422  if(k == bas_k)
423  {
424  continue;
425  }
426  r[m_row.index(k)] -= m_row.value(k) * multiplier;
427  }
428 
429  // compute the value of new dual variable (because we have a new row)
430  Real val = m_objs[bas_k];
431  DSVector basis_col = m_cols[bas_k];
432 
433  for(int k = 0; k < basis_col.size(); ++k)
434  {
435  if (basis_col.index(k) != m_i)
436  val -= basis_col.value(k) * y[basis_col.index(k)];
437  }
438 
439  y[m_i] = val/aij;
440  }
441  else // slack in the basis
442  {
443  rStatus[m_i] = SPxSolver::BASIC;
444  y[m_i] = m_rowobj;
445  }
446 
447 #ifdef CHECK_BASIC_DIM
448  if (!checkBasisDim(rStatus, cStatus))
449  {
450  throw SPxInternalCodeException("XMAISM18 Dimension doesn't match after this step.");
451  }
452 #endif
453 }
454 
457  DataArray<SPxSolver::VarStatus>& rStatus) const
458 {
459  // update the index mapping; if m_correctIdx is false, we assume that this has happened already
460  if(m_correctIdx)
461  {
462  x[m_old_j] = x[m_j];
463  r[m_old_j] = r[m_j];
464  cStatus[m_old_j] = cStatus[m_j];
465  }
466 
467  // primal:
468  x[m_j] = m_val;
469 
470  for(int k = 0; k < m_col.size(); ++k)
471  s[m_col.index(k)] += m_col.value(k) * x[m_j];
472 
473  // dual:
474  Real val = m_obj;
475 
476  for(int k = 0; k < m_col.size(); ++k)
477  val -= m_col.value(k) * y[m_col.index(k)];
478 
479  r[m_j] = val;
480 
481  // basis:
482  if( m_lower == m_upper )
483  {
484  assert(EQrel(m_lower, m_val));
485 
486  cStatus[m_j] = SPxSolver::FIXED;
487  }
488  else
489  {
490  assert(EQrel(m_val, m_lower) || EQrel(m_val, m_upper) || m_val == 0.0);
491 
492  cStatus[m_j] = EQrel(m_val, m_lower) ? SPxSolver::ON_LOWER : (EQrel(m_val, m_upper) ? SPxSolver::ON_UPPER : SPxSolver::ZERO);
493  }
494 
495 #ifdef CHECK_BASIC_DIM
496  if(m_correctIdx)
497  {
498  if (!checkBasisDim(rStatus, cStatus))
499  {
500  throw SPxInternalCodeException("XMAISM19 Dimension doesn't match after this step.");
501  }
502  }
503 #endif
504 }
505 
509 {
510  // basis:
511  cStatus[m_j] = m_status;
512 }
513 
516  DataArray<SPxSolver::VarStatus>& rStatus) const
517 {
518  // correcting the change of idx by deletion of the column and corresponding rows:
519  x[m_old_j] = x[m_j];
520  r[m_old_j] = r[m_j];
521  cStatus[m_old_j] = cStatus[m_j];
522 
523  int rIdx = m_old_i - m_col.size() + 1;
524 
525  for(int k = 0; k < m_col.size(); ++k)
526  {
527  int rIdx_new = m_col.index(k);
528  s[rIdx] = s[rIdx_new];
529  y[rIdx] = y[rIdx_new];
530  rStatus[rIdx] = rStatus[rIdx_new];
531  rIdx++;
532  }
533 
534  // primal:
535  int domIdx = -1;
536  DSVector slack(m_col.size());
537 
538  if (m_loFree)
539  {
540  Real minRowUp = infinity;
541 
542  for(int k = 0; k < m_rows.size(); ++k)
543  {
544  Real val = 0.0;
545  const SVector& row = m_rows[k];
546 
547  for(int l = 0; l < row.size(); ++l)
548  {
549  if (row.index(l) != m_j)
550  val += row.value(l) * x[row.index(l)];
551  }
552 
553  Real scale = maxAbs(m_lRhs[k], val);
554 
555  if (scale < 1.0)
556  scale = 1.0;
557 
558  Real z = (m_lRhs[k] / scale) - (val / scale);
559 
560  if (isZero(z))
561  z = 0.0;
562 
563  Real up = z * scale / row[m_j];
564  slack.add(k, val);
565 
566  if (up < minRowUp)
567  {
568  minRowUp = up;
569  domIdx = k;
570  }
571  }
572 
573  if (m_bnd < minRowUp)
574  {
575  x[m_j] = m_bnd;
576  domIdx = -1;
577  }
578  else
579  x[m_j] = minRowUp;
580  }
581  else
582  {
583  Real maxRowLo = -infinity;
584 
585  for(int k = 0; k < m_rows.size(); ++k)
586  {
587  Real val = 0.0;
588  const SVector& row = m_rows[k];
589 
590  for(int l = 0; l < row.size(); ++l)
591  {
592  if (row.index(l) != m_j)
593  val += row.value(l) * x[row.index(l)];
594  }
595 
596  Real scale = maxAbs(m_lRhs[k], val);
597 
598  if (scale < 1.0)
599  scale = 1.0;
600 
601  Real z = (m_lRhs[k] / scale) - (val / scale);
602 
603  if (isZero(z))
604  z = 0.0;
605 
606  Real lo = z * scale / row[m_j];
607  slack.add(k, val);
608 
609  if (lo > maxRowLo)
610  {
611  maxRowLo = lo;
612  domIdx = k;
613  }
614  }
615 
616  if (m_bnd > maxRowLo)
617  {
618  x[m_j] = m_bnd;
619  domIdx = -1;
620  }
621  else
622  x[m_j] = maxRowLo;
623  }
624 
625  for(int k = 0; k < m_col.size(); ++k)
626  s[m_col.index(k)] = slack[k] + m_col.value(k) * x[m_j];
627 
628  // dual:
629  r[m_j] = 0.0;
630 
631  for(int k = 0; k < m_col.size(); ++k)
632  y[m_col.index(k)] = m_rowObj.value(k);
633 
634  // basis:
635  for(int k = 0; k < m_col.size(); ++k)
636  {
637  if (k != domIdx)
638  rStatus[m_col.index(k)] = SPxSolver::BASIC;
639 
640  else
641  {
642  cStatus[m_j] = SPxSolver::BASIC;
643  if (m_loFree)
644  rStatus[m_col.index(k)] = (m_col.value(k) > 0) ? SPxSolver::ON_UPPER : SPxSolver::ON_LOWER;
645  else
646  rStatus[m_col.index(k)] = (m_col.value(k) > 0) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
647  }
648  }
649  if (domIdx == -1)
650  {
651  if (m_loFree)
652  cStatus[m_j] = SPxSolver::ON_UPPER;
653  else
654  cStatus[m_j] = SPxSolver::ON_LOWER;
655  }
656 
657 #ifdef CHECK_BASIC_DIM
658  if (!checkBasisDim(rStatus, cStatus))
659  {
660  throw SPxInternalCodeException("XMAISM20 Dimension doesn't match after this step.");
661  }
662 #endif
663 }
664 
667  DataArray<SPxSolver::VarStatus>& rStatus) const
668 {
669  // correcting the change of idx by deletion of the column and corresponding rows:
670  x[m_old_j] = x[m_j];
671  r[m_old_j] = r[m_j];
672  cStatus[m_old_j] = cStatus[m_j];
673 
674  // primal & basis:
675  Real aij = m_row[m_j];
676 
677  if (isZero(s[m_i], 1e-6))
678  s[m_i] = 0.0;
679  else if (s[m_i] >= soplex::infinity)
680  // this is a fix for a highly ill conditioned instance that is "solved" in presolving (ilaser0 from MINLP, mittelmann)
681  throw SPxException("Simplifier: scaling factor is too big - aborting unsimplification");
682 
683  Real scale1 = maxAbs(m_lhs, s[m_i]);
684  Real scale2 = maxAbs(m_rhs, s[m_i]);
685 
686  if (scale1 < 1.0)
687  scale1 = 1.0;
688  if (scale2 < 1.0)
689  scale2 = 1.0;
690 
691  Real z1 = (m_lhs / scale1) - (s[m_i] / scale1);
692  Real z2 = (m_rhs / scale2) - (s[m_i] / scale2);
693 
694  if (isZero(z1))
695  z1 = 0.0;
696  if (isZero(z2))
697  z2 = 0.0;
698 
699  Real lo = (aij > 0) ? z1 * scale1 / aij : z2 * scale2 / aij;
700  Real up = (aij > 0) ? z2 * scale2 / aij : z1 * scale1 / aij;
701 
702  if (isZero(lo, eps()))
703  lo = 0.0;
704  if (isZero(up, eps()))
705  up = 0.0;
706 
707  assert(LErel(lo, up));
708  ASSERT_WARN( "WMAISM01", isNotZero(aij, 1.0 / infinity) );
709 
710  if (rStatus[m_i] == SPxSolver::ON_LOWER)
711  {
712  if ( m_lower == m_upper )
713  {
714  x[m_j] = m_lower;
715  cStatus[m_j] = SPxSolver::FIXED;
716  }
717  else if (aij > 0)
718  {
719  x[m_j] = m_upper;
720  cStatus[m_j] = SPxSolver::ON_UPPER;
721  }
722  else if (aij < 0)
723  {
724  x[m_j] = m_lower;
725  cStatus[m_j] = SPxSolver::ON_LOWER;
726  }
727  else
728  throw SPxInternalCodeException("XMAISM01 This should never happen.");
729  }
730  else if (rStatus[m_i] == SPxSolver::ON_UPPER)
731  {
732  if ( m_lower == m_upper )
733  {
734  x[m_j] = m_lower;
735  cStatus[m_j] = SPxSolver::FIXED;
736  }
737  else if (aij > 0)
738  {
739  x[m_j] = m_lower;
740  cStatus[m_j] = SPxSolver::ON_LOWER;
741  }
742  else if (aij < 0)
743  {
744  x[m_j] = m_upper;
745  cStatus[m_j] = SPxSolver::ON_UPPER;
746  }
747  else
748  throw SPxInternalCodeException("XMAISM02 This should never happen.");
749  }
750  else if (rStatus[m_i] == SPxSolver::FIXED)
751  {
752  assert(EQrel(m_lower, m_upper));
753 
754  x[m_j] = m_lower;
755  cStatus[m_j] = SPxSolver::FIXED;
756  }
757  else if (rStatus[m_i] == SPxSolver::BASIC)
758  {
759  if (GErel(m_lower, lo, eps()) && m_lower > -infinity)
760  {
761  x[m_j] = m_lower;
762  cStatus[m_j] = (m_lower == m_upper) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
763  }
764  else if (LErel(m_upper, up, eps()) && m_upper < infinity)
765  {
766  x[m_j] = m_upper;
767  cStatus[m_j] = (m_lower == m_upper) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
768  }
769  else if (lo > -infinity)
770  {
771  // make m_i non-basic and m_j basic
772  x[m_j] = lo;
773  cStatus[m_j] = SPxSolver::BASIC;
774  rStatus[m_i] = (aij > 0 ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER);
775  }
776  else if (up < infinity)
777  {
778  // make m_i non-basic and m_j basic
779  x[m_j] = up;
780  cStatus[m_j] = SPxSolver::BASIC;
781  rStatus[m_i] = (aij > 0 ? SPxSolver::ON_UPPER : SPxSolver::ON_LOWER);
782  }
783  else
784  throw SPxInternalCodeException("XMAISM03 This should never happen.");
785  }
786  else
787  throw SPxInternalCodeException("XMAISM04 This should never happen.");
788 
789  s[m_i] += aij * x[m_j];
790 
791  // dual:
792  r[m_j] = -1.0 * aij * y[m_i];
793 
794  assert(cStatus[m_j] != SPxSolver::BASIC || isZero(r[m_j], eps()));
795 
796 #ifdef CHECK_BASIC_DIM
797  if (!checkBasisDim(rStatus, cStatus))
798  {
799  throw SPxInternalCodeException("XMAISM21 Dimension doesn't match after this step.");
800  }
801 #endif
802 }
803 
806  DataArray<SPxSolver::VarStatus>& rStatus) const
807 {
808 
809  // correcting the change of idx by deletion of the row:
810  s[m_old_i] = s[m_i];
811  y[m_old_i] = y[m_i];
812  rStatus[m_old_i] = rStatus[m_i];
813 
814  // correcting the change of idx by deletion of the column:
815  x[m_old_j] = x[m_j];
816  r[m_old_j] = r[m_j];
817  cStatus[m_old_j] = cStatus[m_j];
818 
819  // primal:
820  Real val = 0.0;
821  Real aij = m_row[m_j];
822 
823  for(int k = 0; k < m_row.size(); ++k)
824  {
825  if (m_row.index(k) != m_j)
826  val += m_row.value(k) * x[m_row.index(k)];
827  }
828 
829  Real scale = maxAbs(m_lRhs, val);
830 
831  if (scale < 1.0)
832  scale = 1.0;
833 
834  Real z = (m_lRhs / scale) - (val / scale);
835 
836  if (isZero(z))
837  z = 0.0;
838 
839  x[m_j] = z * scale / aij;
840  s[m_i] = m_lRhs;
841 
842  // dual:
843  y[m_i] = m_obj / aij;
844  r[m_j] = 0.0;
845 
846  // basis:
847  cStatus[m_j] = SPxSolver::BASIC;
848 
849  if (m_eqCons)
850  rStatus[m_i] = SPxSolver::FIXED;
851  else if (m_onLhs)
852  rStatus[m_i] = SPxSolver::ON_LOWER;
853  else
854  rStatus[m_i] = SPxSolver::ON_UPPER;
855 
856 #ifdef CHECK_BASIC_DIM
857  if (!checkBasisDim(rStatus, cStatus))
858  {
859  throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
860  }
861 #endif
862 }
863 
866  DataArray<SPxSolver::VarStatus>& rStatus) const
867 {
868  // dual:
869  if ((cStatus[m_k] != SPxSolver::BASIC) &&
870  ((cStatus[m_k] == SPxSolver::ON_LOWER && m_strictLo) ||
871  (cStatus[m_k] == SPxSolver::ON_UPPER && m_strictUp) ||
872  (cStatus[m_k] == SPxSolver::FIXED &&
873  (( m_maxSense && ((r[m_j] > 0 && m_strictUp) || (r[m_j] < 0 && m_strictLo))) ||
874  (!m_maxSense && ((r[m_j] > 0 && m_strictLo) || (r[m_j] < 0 && m_strictUp)))))))
875  {
876  Real val = m_kObj;
877  Real aik = m_col[m_i];
878 
879  for(int _k = 0; _k < m_col.size(); ++_k)
880  {
881  if (m_col.index(_k) != m_i)
882  val -= m_col.value(_k) * y[m_col.index(_k)];
883  }
884 
885  y[m_i] = val / aik;
886  r[m_k] = 0.0;
887 
888  r[m_j] = m_jObj - val * m_aij / aik;
889 
890  ASSERT_WARN( "WMAISM73", isNotZero(m_aij * aik) );
891 
892  // basis:
893  if( m_jFixed)
894  cStatus[m_j] = SPxSolver::FIXED;
895  else
896  cStatus[m_j] = (r[m_j] > 0) ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
897 
898  cStatus[m_k] = SPxSolver::BASIC;
899  }
900 
901 #ifdef CHECK_BASIC_DIM
902  if (!checkBasisDim(rStatus, cStatus))
903  {
904  throw SPxInternalCodeException("XMAISM23 Dimension doesn't match after this step.");
905  }
906 #endif
907 }
908 
911  DataArray<SPxSolver::VarStatus>& rStatus) const
912 {
913  // correcting the change of idx by deletion of the duplicated rows:
914  if(m_isLast)
915  {
916  for(int i = m_perm.size() - 1; i >= 0; --i)
917  {
918  if (m_perm[i] >= 0)
919  {
920  int rIdx_new = m_perm[i];
921  int rIdx = i;
922  s[rIdx] = s[rIdx_new];
923  y[rIdx] = y[rIdx_new];
924  rStatus[rIdx] = rStatus[rIdx_new];
925  }
926  }
927  }
928 
929  // primal:
930  for(int k = 0; k < m_scale.size(); ++k)
931  {
932  if (m_scale.index(k) != m_i)
933  s[m_scale.index(k)] = s[m_i] / m_scale.value(k);
934  }
935 
936  // dual & basis:
937  bool haveSetBasis = false;
938 
939  for(int k = 0; k < m_scale.size(); ++k)
940  {
941  int i = m_scale.index(k);
942 
943  if (rStatus[m_i] == SPxSolver::BASIC || (haveSetBasis && i!=m_i))
944  // if the row with tightest lower and upper bound in the basic, every duplicate row should in basic
945  // or basis status of row m_i has been set, this row should be in basis
946  {
947  y[i] = m_rowObj.value(k);
948  rStatus[i] = SPxSolver::BASIC;
949  continue;
950  }
951 
952  ASSERT_WARN( "WMAISM02", isNotZero(m_scale.value(k)) );
953 
954  if (rStatus[m_i] == SPxSolver::FIXED && (i == m_maxLhsIdx || i == m_minRhsIdx))
955  {
956  // this row leads to the tightest lower or upper bound, slack should not be in the basis
957  y[i] = y[m_i] * m_scale.value(k);
958  y[m_i] = m_i_rowObj;
959 
960  if(m_isLhsEqualRhs[k])
961  {
962  rStatus[i] = SPxSolver::FIXED;
963  }
964  else if(i == m_maxLhsIdx)
965  {
966  rStatus[i] = m_scale.value(k)*m_scale.value(0) > 0 ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
967  }
968  else
969  {
970  assert(i == m_minRhsIdx);
971 
972  rStatus[i] = m_scale.value(k)*m_scale.value(0) > 0 ? SPxSolver::ON_UPPER : SPxSolver::ON_LOWER;
973  }
974  if (i != m_i)
975  rStatus[m_i] = SPxSolver::BASIC;
976  haveSetBasis = true;
977  }
978  else if (i == m_maxLhsIdx && rStatus[m_i] == SPxSolver::ON_LOWER)
979  {
980  // this row leads to the tightest lower bound, slack should not be in the basis
981  y[i] = y[m_i] * m_scale.value(k);
982  y[m_i] = m_i_rowObj;
983 
984  rStatus[i] = m_scale.value(k)*m_scale.value(0) > 0 ? SPxSolver::ON_LOWER : SPxSolver::ON_UPPER;
985  if (i != m_i)
986  rStatus[m_i] = SPxSolver::BASIC;
987  haveSetBasis = true;
988  }
989  else if (i == m_minRhsIdx && rStatus[m_i] == SPxSolver::ON_UPPER)
990  {
991  // this row leads to the tightest upper bound, slack should not be in the basis
992  y[i] = y[m_i] * m_scale.value(k);
993  y[m_i] = m_i_rowObj;
994 
995  rStatus[i] = m_scale.value(k)*m_scale.value(0) > 0 ? SPxSolver::ON_UPPER : SPxSolver::ON_LOWER;
996  if (i != m_i)
997  rStatus[m_i] = SPxSolver::BASIC;
998  haveSetBasis = true;
999  }
1000  else if (i != m_i)
1001  {
1002  // this row does not lead to the tightest lower or upper bound, slack should be in the basis
1003  y[i] = m_rowObj.value(k);
1004  rStatus[i] = SPxSolver::BASIC;
1005  }
1006  }
1007 
1008 #ifdef CHECK_BASIC_DIM
1009  if(m_isFirst && !checkBasisDim(rStatus, cStatus))
1010  {
1011  throw SPxInternalCodeException("XMAISM24 Dimension doesn't match after this step.");
1012  }
1013 #endif
1014 
1015  // nothing to do for the reduced cost values
1016 }
1017 
1019  DVector&,
1020  DVector&,
1021  DVector& r,
1023  DataArray<SPxSolver::VarStatus>& rStatus) const
1024 {
1025 
1026  if(m_isFirst)
1027  {
1028 #ifdef CHECK_BASIC_DIM
1029  if (!checkBasisDim(rStatus, cStatus))
1030  {
1031  throw SPxInternalCodeException("XMAISM25 Dimension doesn't match after this step.");
1032  }
1033 #endif
1034  return;
1035  }
1036 
1037 
1038  // correcting the change of idx by deletion of the columns:
1039  if(m_isLast)
1040  {
1041  for(int i = m_perm.size() - 1; i >= 0; --i)
1042  {
1043  if (m_perm[i] >= 0)
1044  {
1045  int cIdx_new = m_perm[i];
1046  int cIdx = i;
1047  x[cIdx] = x[cIdx_new];
1048  r[cIdx] = r[cIdx_new];
1049  cStatus[cIdx] = cStatus[cIdx_new];
1050  }
1051  }
1052  return;
1053  }
1054 
1055  // primal & basis:
1056  ASSERT_WARN( "WMAISM03", isNotZero(m_scale) );
1057 
1058  if (cStatus[m_k] == SPxSolver::ON_LOWER)
1059  {
1060  x[m_k] = m_loK;
1061 
1062  if (m_scale > 0)
1063  {
1064  x[m_j] = m_loJ;
1065  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
1066  }
1067  else
1068  {
1069  x[m_j] = m_upJ;
1070  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
1071  }
1072  }
1073  else if (cStatus[m_k] == SPxSolver::ON_UPPER)
1074  {
1075  x[m_k] = m_upK;
1076 
1077  if (m_scale > 0)
1078  {
1079  x[m_j] = m_upJ;
1080  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
1081  }
1082  else
1083  {
1084  x[m_j] = m_loJ;
1085  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
1086  }
1087  }
1088  else if (cStatus[m_k] == SPxSolver::FIXED)
1089  {
1090  // => x[m_k] and x[m_j] are also fixed before the corresponding preprocessing step
1091  x[m_j] = m_loJ;
1092  cStatus[m_j] = SPxSolver::FIXED;
1093  }
1094  else if (cStatus[m_k] == SPxSolver::ZERO)
1095  {
1096  /* we only aggregate duplicate columns if 0 is contained in their bounds, so we can handle this case properly */
1097  assert(x[m_k] == 0.0);
1098  assert(LErel(m_loJ, 0.0));
1099  assert(GErel(m_upJ, 0.0));
1100  assert(LErel(m_loK, 0.0));
1101  assert(GErel(m_upK, 0.0));
1102 
1103  if (isZero(m_loK) && isZero(m_upK) && m_loK == m_upK)
1104  cStatus[m_k] = SPxSolver::FIXED;
1105  else if (isZero(m_loK))
1106  cStatus[m_k] = SPxSolver::ON_LOWER;
1107  else if (isZero(m_upK))
1108  cStatus[m_k] = SPxSolver::ON_UPPER;
1109  else if (LErel(m_loK, 0.0) && GErel(m_upK, 0.0))
1110  cStatus[m_k] = SPxSolver::ZERO;
1111  else
1112  throw SPxInternalCodeException("XMAISM05 This should never happen.");
1113 
1114  x[m_j] = 0.0;
1115  if (isZero(m_loJ) && isZero(m_upJ) && m_loJ == m_upJ)
1116  cStatus[m_j] = SPxSolver::FIXED;
1117  else if (isZero(m_loJ))
1118  cStatus[m_j] = SPxSolver::ON_LOWER;
1119  else if (isZero(m_upJ))
1120  cStatus[m_j] = SPxSolver::ON_UPPER;
1121  else if (LErel(m_loJ, 0.0) && GErel(m_upJ, 0.0))
1122  cStatus[m_j] = SPxSolver::ZERO;
1123  else
1124  throw SPxInternalCodeException("XMAISM06 This should never happen.");
1125  }
1126  else if (cStatus[m_k] == SPxSolver::BASIC)
1127  {
1128  Real scale1 = maxAbs(x[m_k], m_loK);
1129  Real scale2 = maxAbs(x[m_k], m_upK);
1130 
1131  if (scale1 < 1.0)
1132  scale1 = 1.0;
1133  if (scale2 < 1.0)
1134  scale2 = 1.0;
1135 
1136  Real z1 = (x[m_k] / scale1) - (m_loK / scale1);
1137  Real z2 = (x[m_k] / scale2) - (m_upK / scale2);
1138 
1139  if (isZero(z1))
1140  z1 = 0.0;
1141  if (isZero(z2))
1142  z2 = 0.0;
1143 
1144  if( m_loJ <= -infinity && m_upJ >= infinity && m_loK <= -infinity && m_upK >= infinity )
1145  {
1146  cStatus[m_j] = SPxSolver::ZERO;
1147  x[m_j] = 0.0;
1148  }
1149  else if( m_scale > 0.0 )
1150  {
1151  if( GErel(x[m_k], m_upK + m_scale * m_upJ) )
1152  {
1153  assert(m_upJ < infinity);
1154  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
1155  x[m_j] = m_upJ;
1156  x[m_k] -= m_scale * x[m_j];
1157  }
1158  else if( GErel(x[m_k], m_loK + m_scale * m_upJ) && m_upJ < infinity )
1159  {
1160  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
1161  x[m_j] = m_upJ;
1162  x[m_k] -= m_scale * x[m_j];
1163  }
1164  else if( GErel(x[m_k], m_upK + m_scale * m_loJ) && m_upK < infinity )
1165  {
1166  cStatus[m_k] = (m_loK == m_upK) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
1167  x[m_k] = m_upK;
1168  cStatus[m_j] = SPxSolver::BASIC;
1169  x[m_j] = z2 * scale2 / m_scale;
1170  }
1171  else if( GErel(x[m_k], m_loK + m_scale * m_loJ) && m_loJ > -infinity )
1172  {
1173  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
1174  x[m_j] = m_loJ;
1175  x[m_k] -= m_scale * x[m_j];
1176  }
1177  else if( GErel(x[m_k], m_loK + m_scale * m_loJ) && m_loK > -infinity )
1178  {
1179  cStatus[m_k] = (m_loK == m_upK) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
1180  x[m_k] = m_loK;
1181  cStatus[m_j] = SPxSolver::BASIC;
1182  x[m_j] = z1 * scale1 / m_scale;
1183  }
1184  else if( LTrel(x[m_k], m_loK + m_scale * m_loJ) )
1185  {
1186  assert(m_loJ > -infinity);
1187  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
1188  x[m_j] = m_loJ;
1189  x[m_k] -= m_scale * x[m_j];
1190  }
1191  else
1192  {
1193  throw SPxInternalCodeException("XMAISM08 This should never happen.");
1194  }
1195  }
1196  else
1197  {
1198  assert(m_scale < 0.0);
1199 
1200  if( GErel(x[m_k], m_upK + m_scale * m_loJ) )
1201  {
1202  assert(m_loJ > -infinity);
1203  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
1204  x[m_j] = m_loJ;
1205  x[m_k] -= m_scale * x[m_j];
1206  }
1207  else if( GErel(x[m_k], m_loK + m_scale * m_loJ) && m_loJ > -infinity )
1208  {
1209  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
1210  x[m_j] = m_loJ;
1211  x[m_k] -= m_scale * x[m_j];
1212  }
1213  else if( GErel(x[m_k], m_upK + m_scale * m_upJ) && m_upK < infinity )
1214  {
1215  cStatus[m_k] = (m_loK == m_upK) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
1216  x[m_k] = m_upK;
1217  cStatus[m_j] = SPxSolver::BASIC;
1218  x[m_j] = z2 * scale2 / m_scale;
1219  }
1220  else if( GErel(x[m_k], m_loK + m_scale * m_upJ) && m_upJ < infinity )
1221  {
1222  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
1223  x[m_j] = m_upJ;
1224  x[m_k] -= m_scale * x[m_j];
1225  }
1226  else if( GErel(x[m_k], m_loK + m_scale * m_upJ) && m_loK > -infinity )
1227  {
1228  cStatus[m_k] = (m_loK == m_upK) ? SPxSolver::FIXED : SPxSolver::ON_LOWER;
1229  x[m_k] = m_loK;
1230  cStatus[m_j] = SPxSolver::BASIC;
1231  x[m_j] = z1 * scale1 / m_scale;
1232  }
1233  else if( LTrel(x[m_k], m_loK + m_scale * m_upJ) )
1234  {
1235  assert(m_upJ < infinity);
1236  cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolver::FIXED : SPxSolver::ON_UPPER;
1237  x[m_j] = m_upJ;
1238  x[m_k] -= m_scale * x[m_j];
1239  }
1240  else
1241  {
1242  throw SPxInternalCodeException("XMAISM09 This should never happen.");
1243  }
1244  }
1245  }
1246 
1247  // dual:
1248  r[m_j] = m_scale * r[m_k];
1249 }
1250 
1252 {
1253  for( int i = lp.nRows() - 1; i >= 0; --i )
1254  {
1255  if( lp.maxRowObj(i) != 0.0 )
1256  {
1257  RowObjPS* RowObjPSptr = 0;
1258  spx_alloc(RowObjPSptr);
1259  m_hist.append(new (RowObjPSptr) RowObjPS(lp, i, lp.nCols()));
1260  lp.addCol(lp.rowObj(i), -lp.rhs(i), UnitVector(i), -lp.lhs(i));
1261  lp.changeRange(i, 0.0, 0.0);
1262  lp.changeRowObj(i, 0.0);
1263  m_addedcols++;
1264  }
1265  }
1266 }
1267 
1269 {
1270 
1271  // This method handles extreme value of the given LP by
1272  //
1273  // 1. setting numbers of very small absolute values to zero and
1274  // 2. setting numbers of very large absolute values to -infinity or +infinity, respectively.
1275 
1276  Real maxVal = infinity / 5.0;
1277  Real tol = feastol() * 1e-2;
1278  tol = (tol < epsZero()) ? epsZero() : tol;
1279  int remRows = 0;
1280  int remNzos = 0;
1281  int chgBnds = 0;
1282  int chgLRhs = 0;
1283  int objCnt = 0;
1284 
1285  for(int i = lp.nRows()-1; i >= 0; --i)
1286  {
1287  // lhs
1288  Real lhs = lp.lhs(i);
1289 
1290  if (lhs != 0.0 && isZero(lhs, epsZero()))
1291  {
1292  lp.changeLhs(i, 0.0);
1293  ++chgLRhs;
1294  }
1295  else if (lhs > -infinity && lhs < -maxVal)
1296  {
1297  lp.changeLhs(i, -infinity);
1298  ++chgLRhs;
1299  }
1300  else if (lhs < infinity && lhs > maxVal)
1301  {
1302  lp.changeLhs(i, infinity);
1303  ++chgLRhs;
1304  }
1305 
1306  // rhs
1307  Real rhs = lp.rhs(i);
1308 
1309  if (rhs != 0.0 && isZero(rhs, epsZero()))
1310  {
1311  lp.changeRhs(i, 0.0);
1312  ++chgLRhs;
1313  }
1314  else if (rhs > -infinity && rhs < -maxVal)
1315  {
1316  lp.changeRhs(i, -infinity);
1317  ++chgLRhs;
1318  }
1319  else if (rhs < infinity && rhs > maxVal)
1320  {
1321  lp.changeRhs(i, infinity);
1322  ++chgLRhs;
1323  }
1324 
1325  if (lp.lhs(i) <= -infinity && lp.rhs(i) >= infinity)
1326  {
1327  FreeConstraintPS* FreeConstraintPSptr = 0;
1328  spx_alloc(FreeConstraintPSptr);
1329  m_hist.append(new (FreeConstraintPSptr) FreeConstraintPS(lp, i));
1330 
1331  removeRow(lp, i);
1332  ++remRows;
1333 
1334  ++m_stat[FREE_ROW];
1335  }
1336  }
1337 
1338  for(int j = 0; j < lp.nCols(); ++j)
1339  {
1340  // lower bound
1341  Real lo = lp.lower(j);
1342 
1343  if (lo != 0.0 && isZero(lo, epsZero()))
1344  {
1345  lp.changeLower(j, 0.0);
1346  ++chgBnds;
1347  }
1348  else if (lo > -infinity && lo < -maxVal)
1349  {
1350  lp.changeLower(j, -infinity);
1351  ++chgBnds;
1352  }
1353  else if (lo < infinity && lo > maxVal)
1354  {
1355  lp.changeLower(j, infinity);
1356  ++chgBnds;
1357  }
1358 
1359  // upper bound
1360  Real up = lp.upper(j);
1361 
1362  if (up != 0.0 && isZero(up, epsZero()))
1363  {
1364  lp.changeUpper(j, 0.0);
1365  ++chgBnds;
1366  }
1367  else if (up > -infinity && up < -maxVal)
1368  {
1369  lp.changeUpper(j, -infinity);
1370  ++chgBnds;
1371  }
1372  else if (up < infinity && up > maxVal)
1373  {
1374  lp.changeUpper(j, infinity);
1375  ++chgBnds;
1376  }
1377 
1378  // fixed columns will be eliminated later
1379  if (NE(lo, up))
1380  {
1381  lo = spxAbs(lo);
1382  up = spxAbs(up);
1383 
1384  Real absBnd = (lo > up) ? lo : up;
1385 
1386  if (absBnd < 1.0)
1387  absBnd = 1.0;
1388 
1389  // non-zeros
1390  SVector& col = lp.colVector_w(j);
1391  int i = 0;
1392 
1393  while(i < col.size())
1394  {
1395  Real aij = spxAbs(col.value(i));
1396 
1397  if (isZero(aij * absBnd, tol))
1398  {
1399  SVector& row = lp.rowVector_w(col.index(i));
1400 
1401  // this changes col.size()
1402  row.remove(row.number(j));
1403  col.remove(i);
1404 
1405  MSG_INFO3( (*spxout), (*spxout) << "IMAISM04 aij=" << aij
1406  << " removed, absBnd=" << absBnd
1407  << std::endl; )
1408  ++remNzos;
1409  }
1410  else
1411  {
1412  if (aij > maxVal)
1413  {
1414  MSG_WARNING( (*spxout), (*spxout) << "WMAISM05 Warning! Big matrix coefficient " << aij << std::endl );
1415  }
1416  else if( isZero(aij, tol) )
1417  {
1418  MSG_WARNING( (*spxout), (*spxout) << "WMAISM06 Warning! Tiny matrix coefficient " << aij << std::endl );
1419  }
1420 
1421  ++i;
1422  }
1423  }
1424  }
1425 
1426  // objective
1427  Real obj = lp.obj(j);
1428 
1429  if (obj != 0.0 && isZero(obj, epsZero()))
1430  {
1431  lp.changeObj(j, 0.0);
1432  ++objCnt;
1433  }
1434  else if (obj > -infinity && obj < -maxVal)
1435  {
1436  lp.changeObj(j, -infinity);
1437  ++objCnt;
1438  }
1439  else if (obj < infinity && obj > maxVal)
1440  {
1441  lp.changeObj(j, infinity);
1442  ++objCnt;
1443  }
1444  }
1445 
1446  if (remRows + remNzos + chgLRhs + chgBnds + objCnt > 0)
1447  {
1448  m_remRows += remRows;
1449  m_remNzos += remNzos;
1450  m_chgLRhs += chgLRhs;
1451  m_chgBnds += chgBnds;
1452 
1453  MSG_INFO2( (*spxout), (*spxout) << "Simplifier (extremes) removed "
1454  << remRows << " rows, "
1455  << remNzos << " non-zeros, "
1456  << chgBnds << " col bounds, "
1457  << chgLRhs << " row bounds, "
1458  << objCnt << " objective coefficients" << std::endl; )
1459  }
1460  assert(lp.isConsistent());
1461 }
1462 
1464 {
1465 
1466  // This method removes empty rows and columns from the LP.
1467 
1468  int remRows = 0;
1469  int remCols = 0;
1470 
1471  for(int i = lp.nRows()-1; i >= 0; --i)
1472  {
1473  const SVector& row = lp.rowVector(i);
1474 
1475  if (row.size() == 0)
1476  {
1477  MSG_INFO3( (*spxout), (*spxout) << "IMAISM07 row " << i
1478  << ": empty ->"; )
1479 
1480  if (LT(lp.rhs(i), 0.0, feastol()) || GT(lp.lhs(i), 0.0, feastol()))
1481  {
1482  MSG_INFO3( (*spxout), (*spxout) << " infeasible lhs=" << lp.lhs(i)
1483  << " rhs=" << lp.rhs(i) << std::endl; )
1484  return INFEASIBLE;
1485  }
1486  MSG_INFO3( (*spxout), (*spxout) << " removed" << std::endl; )
1487 
1488  EmptyConstraintPS* EmptyConstraintPSptr = 0;
1489  spx_alloc(EmptyConstraintPSptr);
1490  m_hist.append(new (EmptyConstraintPSptr) EmptyConstraintPS(lp, i));
1491 
1492  ++remRows;
1493  removeRow(lp, i);
1494 
1495  ++m_stat[EMPTY_ROW];
1496  }
1497  }
1498 
1499  for(int j = lp.nCols()-1; j >= 0; --j)
1500  {
1501  const SVector& col = lp.colVector(j);
1502 
1503  if (col.size() == 0)
1504  {
1505  MSG_INFO3( (*spxout), (*spxout) << "IMAISM08 col " << j
1506  << ": empty -> maxObj=" << lp.maxObj(j)
1507  << " lower=" << lp.lower(j)
1508  << " upper=" << lp.upper(j); )
1509 
1510  Real val;
1511 
1512  if (GT(lp.maxObj(j), 0.0, epsZero()))
1513  {
1514  if (lp.upper(j) >= infinity)
1515  {
1516  MSG_INFO3( (*spxout), (*spxout) << " unbounded" << std::endl; )
1517  return UNBOUNDED;
1518  }
1519  val = lp.upper(j);
1520  }
1521  else if (LT(lp.maxObj(j), 0.0, epsZero()))
1522  {
1523  if (lp.lower(j) <= -infinity)
1524  {
1525  MSG_INFO3( (*spxout), (*spxout) << " unbounded" << std::endl; )
1526  return UNBOUNDED;
1527  }
1528  val = lp.lower(j);
1529  }
1530  else
1531  {
1532  ASSERT_WARN( "WMAISM09", isZero(lp.maxObj(j), epsZero()) );
1533  // any value within the bounds is ok
1534  if (lp.lower(j) > -infinity)
1535  val = lp.lower(j);
1536  else if (lp.upper(j) < infinity)
1537  val = lp.upper(j);
1538  else
1539  val = 0.0;
1540  }
1541  MSG_INFO3( (*spxout), (*spxout) << " removed" << std::endl; )
1542 
1543  FixBoundsPS* FixBoundsPSptr = 0;
1544  FixVariablePS* FixVariablePSptr = 0;
1545  spx_alloc(FixBoundsPSptr);
1546  spx_alloc(FixVariablePSptr);
1547  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j, val));
1548  m_hist.append(new (FixVariablePSptr) FixVariablePS(lp, *this, j, val));
1549 
1550  ++remCols;
1551  removeCol(lp, j);
1552 
1553  ++m_stat[EMPTY_COL];
1554  }
1555  }
1556 
1557  if (remRows + remCols > 0)
1558  {
1559  m_remRows += remRows;
1560  m_remCols += remCols;
1561 
1562  MSG_INFO2( (*spxout), (*spxout) << "Simplifier (empty rows/colums) removed "
1563  << remRows << " rows, "
1564  << remCols << " cols"
1565  << std::endl; )
1566 
1567  }
1568  return OKAY;
1569 }
1570 
1572 {
1573  assert(row.size() == 1);
1574 
1575  Real aij = row.value(0);
1576  int j = row.index(0);
1577  Real lo = -infinity;
1578  Real up = infinity;
1579 
1580  MSG_INFO3( (*spxout), (*spxout) << "IMAISM22 row " << i
1581  << ": singleton -> val=" << aij
1582  << " lhs=" << lp.lhs(i)
1583  << " rhs=" << lp.rhs(i); )
1584 
1585  if (GT(aij, 0.0, epsZero())) // aij > 0
1586  {
1587  lo = (lp.lhs(i) <= -infinity) ? -infinity : lp.lhs(i) / aij;
1588  up = (lp.rhs(i) >= infinity) ? infinity : lp.rhs(i) / aij;
1589  }
1590  else if (LT(aij, 0.0, epsZero())) // aij < 0
1591  {
1592  lo = (lp.rhs(i) >= infinity) ? -infinity : lp.rhs(i) / aij;
1593  up = (lp.lhs(i) <= -infinity) ? infinity : lp.lhs(i) / aij;
1594  }
1595  else if (LT(lp.rhs(i), 0.0, feastol()) || GT(lp.lhs(i), 0.0, feastol()))
1596  {
1597  // aij == 0, rhs < 0 or lhs > 0
1598  MSG_INFO3( (*spxout), (*spxout) << " infeasible" << std::endl; )
1599  return INFEASIBLE;
1600  }
1601 
1602  if (isZero(lo, epsZero()))
1603  lo = 0.0;
1604 
1605  if (isZero(up, epsZero()))
1606  up = 0.0;
1607 
1608  MSG_INFO3( (*spxout), (*spxout) << " removed, lower=" << lo
1609  << " (" << lp.lower(j)
1610  << ") upper=" << up
1611  << " (" << lp.upper(j)
1612  << ")" << std::endl; )
1613 
1614  bool stricterUp = false;
1615  bool stricterLo = false;
1616 
1617  Real oldLo = lp.lower(j);
1618  Real oldUp = lp.upper(j);
1619 
1620  if (LTrel(up, lp.upper(j), feastol()))
1621  {
1622  lp.changeUpper(j, up);
1623  stricterUp = true;
1624  }
1625  if (GTrel(lo, lp.lower(j), feastol()))
1626  {
1627  lp.changeLower(j, lo);
1628  stricterLo = true;
1629  }
1630 
1631  RowSingletonPS* RowSingletonPSptr = 0;
1632  spx_alloc(RowSingletonPSptr);
1633  m_hist.append(new (RowSingletonPSptr) RowSingletonPS(lp, i, j, stricterLo, stricterUp, lp.lower(j), lp.upper(j), oldLo, oldUp));
1634 
1635  removeRow(lp, i);
1636 
1637  m_remRows++;
1638  m_remNzos++;
1639  ++m_stat[SINGLETON_ROW];
1640 
1641  return OKAY;
1642 }
1643 
1645 {
1646 
1647  // This method simplifies the rows of the LP.
1648  //
1649  // The following operations are done:
1650  // 1. detect implied free variables
1651  // 2. detect implied free constraints
1652  // 3. detect infeasible constraints
1653  // 4. remove unconstrained constraints
1654  // 5. remove empty constraints
1655  // 6. remove row singletons and tighten the corresponding variable bounds if necessary
1656  // 7. detect forcing rows and fix the corresponding variables
1657 
1658  int remRows = 0;
1659  int remNzos = 0;
1660  int chgLRhs = 0;
1661  int chgBnds = 0;
1662  int keptBnds = 0;
1663  int keptLRhs = 0;
1664 
1665  int oldRows = lp.nRows();
1666 
1667  bool redundantLower;
1668  bool redundantUpper;
1669  bool redundantLhs;
1670  bool redundantRhs;
1671 
1672  for(int i = lp.nRows()-1; i >= 0; --i)
1673  {
1674  const SVector& row = lp.rowVector(i);
1675 
1676 
1677  // compute bounds on constraint value
1678  Real lhsBnd = 0.0; // minimal activity (finite summands)
1679  Real rhsBnd = 0.0; // maximal activity (finite summands)
1680  int lhsCnt = 0; // number of -infinity summands in minimal activity
1681  int rhsCnt = 0; // number of +infinity summands in maximal activity
1682 
1683  for(int k = 0; k < row.size(); ++k)
1684  {
1685  Real aij = row.value(k);
1686  int j = row.index(k);
1687 
1688  if( !isNotZero(aij, 1.0 / infinity) )
1689  {
1690  MSG_WARNING( (*spxout), (*spxout) << "Warning: tiny nonzero coefficient " << aij << " in row " << i << "\n" );
1691  }
1692 
1693  if (aij > 0.0)
1694  {
1695  if (lp.lower(j) <= -infinity)
1696  ++lhsCnt;
1697  else
1698  lhsBnd += aij * lp.lower(j);
1699 
1700  if (lp.upper(j) >= infinity)
1701  ++rhsCnt;
1702  else
1703  rhsBnd += aij * lp.upper(j);
1704  }
1705  else if (aij < 0.0)
1706  {
1707  if (lp.lower(j) <= -infinity)
1708  ++rhsCnt;
1709  else
1710  rhsBnd += aij * lp.lower(j);
1711 
1712  if (lp.upper(j) >= infinity)
1713  ++lhsCnt;
1714  else
1715  lhsBnd += aij * lp.upper(j);
1716  }
1717  }
1718 
1719 #if FREE_BOUNDS
1720  // 1. detect implied free variables
1721  if (rhsCnt <= 1 || lhsCnt <= 1)
1722  {
1723  for(int k = 0; k < row.size(); ++k)
1724  {
1725  Real aij = row.value(k);
1726  int j = row.index(k);
1727 
1728  redundantLower = false;
1729  redundantUpper = false;
1730 
1731  ASSERT_WARN( "WMAISM12", isNotZero(aij, 1.0 / infinity) );
1732 
1733  if (aij > 0.0)
1734  {
1735  if (lp.lhs(i) > -infinity && lp.lower(j) > -infinity && rhsCnt <= 1 && NErel(lp.lhs(i), rhsBnd, feastol())
1736  // do not perform if strongly different orders of magnitude occur
1737  && spxAbs(lp.lhs(i) / rhsBnd) > Param::epsilon())
1738  {
1739  Real lo = -infinity;
1740  Real scale = maxAbs(lp.lhs(i), rhsBnd);
1741 
1742  if (scale < 1.0)
1743  scale = 1.0;
1744 
1745  Real z = (lp.lhs(i) / scale) - (rhsBnd / scale);
1746 
1747  if (isZero(z, epsZero()))
1748  z = 0.0;
1749 
1750  assert(rhsCnt > 0 || lp.upper(j) < infinity);
1751 
1752  if (rhsCnt == 0)
1753  lo = lp.upper(j) + z * scale / aij;
1754  else if (lp.upper(j) >= infinity)
1755  lo = z * scale / aij;
1756 
1757  if (isZero(lo, epsZero()))
1758  lo = 0.0;
1759 
1760  if (GErel(lo, lp.lower(j), feastol()))
1761  {
1762  MSG_INFO3( (*spxout), (*spxout) << "IMAISM13 row " << i
1763  << ": redundant lower bound on x" << j
1764  << " -> lower=" << lo
1765  << " (" << lp.lower(j)
1766  << ")" << std::endl; )
1767 
1768  redundantLower = true;
1769  }
1770 
1771  }
1772  if (lp.rhs(i) < infinity && lp.upper(j) < infinity && lhsCnt <= 1 && NErel(lp.rhs(i), lhsBnd, feastol())
1773  // do not perform if strongly different orders of magnitude occur
1774  && spxAbs(lp.rhs(i) / lhsBnd) > Param::epsilon())
1775  {
1776  Real up = infinity;
1777  Real scale = maxAbs(lp.rhs(i), lhsBnd);
1778 
1779  if (scale < 1.0)
1780  scale = 1.0;
1781 
1782  Real z = (lp.rhs(i) / scale) - (lhsBnd / scale);
1783 
1784  if (isZero(z, epsZero()))
1785  z = 0.0;
1786 
1787  assert(lhsCnt > 0 || lp.lower(j) > -infinity);
1788 
1789  if (lhsCnt == 0)
1790  up = lp.lower(j) + z * scale / aij;
1791  else if (lp.lower(j) <= -infinity)
1792  up = z * scale / aij;
1793 
1794  if (isZero(up, epsZero()))
1795  up = 0.0;
1796 
1797  if (LErel(up, lp.upper(j), feastol()))
1798  {
1799  MSG_INFO3( (*spxout), (*spxout) << "IMAISM14 row " << i
1800  << ": redundant upper bound on x" << j
1801  << " -> upper=" << up
1802  << " (" << lp.upper(j)
1803  << ")" << std::endl; )
1804 
1805  redundantUpper = true;
1806  }
1807  }
1808  if (redundantLower)
1809  {
1810  // no upper bound on x_j OR redundant upper bound
1811  if ((lp.upper(j) >= infinity) || redundantUpper || (!m_keepbounds))
1812  {
1813  ++lhsCnt;
1814  lhsBnd -= aij * lp.lower(j);
1815 
1816  lp.changeLower(j, -infinity);
1817  ++chgBnds;
1818  }
1819  else
1820  ++keptBnds;
1821  }
1822  if (redundantUpper)
1823  {
1824  // no lower bound on x_j OR redundant lower bound
1825  if ((lp.lower(j) <= -infinity) || redundantLower || (!m_keepbounds))
1826  {
1827  ++rhsCnt;
1828  rhsBnd -= aij * lp.upper(j);
1829 
1830  lp.changeUpper(j, infinity);
1831  ++chgBnds;
1832  }
1833  else
1834  ++keptBnds;
1835  }
1836  }
1837  else if (aij < 0.0)
1838  {
1839  if (lp.lhs(i) > -infinity && lp.upper(j) < infinity && rhsCnt <= 1 && NErel(lp.lhs(i), rhsBnd, feastol())
1840  // do not perform if strongly different orders of magnitude occur
1841  && spxAbs(lp.lhs(i) / rhsBnd) > Param::epsilon())
1842  {
1843  Real up = infinity;
1844  Real scale = maxAbs(lp.lhs(i), rhsBnd);
1845 
1846  if (scale < 1.0)
1847  scale = 1.0;
1848 
1849  Real z = (lp.lhs(i) / scale) - (rhsBnd / scale);
1850 
1851  if (isZero(z, epsZero()))
1852  z = 0.0;
1853 
1854  assert(rhsCnt > 0 || lp.lower(j) > -infinity);
1855 
1856  if (rhsCnt == 0)
1857  up = lp.lower(j) + z * scale / aij;
1858  else if (lp.lower(j) <= -infinity)
1859  up = z * scale / aij;
1860 
1861  if (isZero(up, epsZero()))
1862  up = 0.0;
1863 
1864  if (LErel(up, lp.upper(j), feastol()))
1865  {
1866  MSG_INFO3( (*spxout), (*spxout) << "IMAISM15 row " << i
1867  << ": redundant upper bound on x" << j
1868  << " -> upper=" << up
1869  << " (" << lp.upper(j)
1870  << ")" << std::endl; )
1871 
1872  redundantUpper = true;
1873  }
1874  }
1875  if (lp.rhs(i) < infinity && lp.lower(j) > -infinity && lhsCnt <= 1 && NErel(lp.rhs(i), lhsBnd, feastol())
1876  // do not perform if strongly different orders of magnitude occur
1877  && spxAbs(lp.rhs(i) / lhsBnd) > Param::epsilon())
1878  {
1879  Real lo = -infinity;
1880  Real scale = maxAbs(lp.rhs(i), lhsBnd);
1881 
1882  if (scale < 1.0)
1883  scale = 1.0;
1884 
1885  Real z = (lp.rhs(i) / scale) - (lhsBnd / scale);
1886 
1887  if (isZero(z, epsZero()))
1888  z = 0.0;
1889 
1890  assert(lhsCnt > 0 || lp.upper(j) < infinity);
1891 
1892  if (lhsCnt == 0)
1893  lo = lp.upper(j) + z * scale / aij;
1894  else if (lp.upper(j) >= infinity)
1895  lo = z * scale / aij;
1896 
1897  if (isZero(lo, epsZero()))
1898  lo = 0.0;
1899 
1900  if (GErel(lo, lp.lower(j)))
1901  {
1902  MSG_INFO3( (*spxout), (*spxout) << "IMAISM16 row " << i
1903  << ": redundant lower bound on x" << j
1904  << " -> lower=" << lo
1905  << " (" << lp.lower(j)
1906  << ")" << std::endl; )
1907 
1908  redundantLower = true;
1909  }
1910  }
1911  if (redundantUpper)
1912  {
1913  // no lower bound on x_j OR redundant lower bound
1914  if ((lp.lower(j) <= -infinity) || redundantLower || (!m_keepbounds))
1915  {
1916  ++lhsCnt;
1917  lhsBnd -= aij * lp.upper(j);
1918 
1919  lp.changeUpper(j, infinity);
1920  ++chgBnds;
1921  }
1922  else
1923  ++keptBnds;
1924  }
1925  if (redundantLower)
1926  {
1927  // no upper bound on x_j OR redundant upper bound
1928  if ((lp.upper(j) >= infinity) || redundantUpper || (!m_keepbounds))
1929  {
1930  ++rhsCnt;
1931  rhsBnd -= aij * lp.lower(j);
1932 
1933  lp.changeLower(j, -infinity);
1934  ++chgBnds;
1935  }
1936  else
1937  ++keptBnds;
1938  }
1939  }
1940  }
1941  }
1942 #endif
1943 
1944 #if FREE_LHS_RHS
1945 
1946  redundantLhs = false;
1947  redundantRhs = false;
1948 
1949  // 2. detect implied free constraints
1950  if (lp.lhs(i) > -infinity && lhsCnt == 0 && GErel(lhsBnd, lp.lhs(i), feastol()))
1951  {
1952  MSG_INFO3( (*spxout), (*spxout) << "IMAISM17 row " << i
1953  << ": redundant lhs -> lhsBnd=" << lhsBnd
1954  << " lhs=" << lp.lhs(i)
1955  << std::endl; )
1956 
1957  redundantLhs = true;
1958  }
1959  if (lp.rhs(i) < infinity && rhsCnt == 0 && LErel(rhsBnd, lp.rhs(i), feastol()))
1960  {
1961  MSG_INFO3( (*spxout), (*spxout) << "IMAISM18 row " << i
1962  << ": redundant rhs -> rhsBnd=" << rhsBnd
1963  << " rhs=" << lp.rhs(i)
1964  << std::endl; )
1965 
1966  redundantRhs = true;
1967  }
1968  if (redundantLhs)
1969  {
1970  // no rhs for constraint i OR redundant rhs
1971  if ((lp.rhs(i) >= infinity) || redundantRhs || (!m_keepbounds))
1972  {
1973  lp.changeLhs(i, -infinity);
1974  ++chgLRhs;
1975  }
1976  else
1977  ++keptLRhs;
1978  }
1979  if (redundantRhs)
1980  {
1981  // no lhs for constraint i OR redundant lhs
1982  if ((lp.lhs(i) <= -infinity) || redundantLhs || (!m_keepbounds))
1983  {
1984  lp.changeRhs(i, infinity);
1985  ++chgLRhs;
1986  }
1987  else
1988  ++keptLRhs;
1989  }
1990 #endif
1991 
1992  // 3. infeasible constraint
1993  if (LTrel(lp.rhs(i), lp.lhs(i), feastol()) ||
1994  (LTrel(rhsBnd, lp.lhs(i), feastol()) && rhsCnt == 0) ||
1995  (GTrel(lhsBnd, lp.rhs(i), feastol()) && lhsCnt == 0))
1996  {
1997  MSG_INFO3( (*spxout), (*spxout) << "IMAISM19 row " << std::setprecision(20) << i
1998  << ": infeasible -> lhs=" << lp.lhs(i)
1999  << " rhs=" << lp.rhs(i)
2000  << " lhsBnd=" << lhsBnd
2001  << " rhsBnd=" << rhsBnd
2002  << std::endl; )
2003  return INFEASIBLE;
2004  }
2005 
2006 #if FREE_CONSTRAINT
2007  // 4. unconstrained constraint
2008  if (lp.lhs(i) <= -infinity && lp.rhs(i) >= infinity)
2009  {
2010  MSG_INFO3( (*spxout), (*spxout) << "IMAISM20 row " << i
2011  << ": unconstrained -> removed" << std::endl; )
2012 
2013  FreeConstraintPS* FreeConstraintPSptr = 0;
2014  spx_alloc(FreeConstraintPSptr);
2015  m_hist.append(new (FreeConstraintPSptr) FreeConstraintPS(lp, i));
2016 
2017  ++remRows;
2018  remNzos += row.size();
2019  removeRow(lp, i);
2020 
2021  ++m_stat[FREE_ROW];
2022 
2023  continue;
2024  }
2025 #endif
2026 
2027 #if EMPTY_CONSTRAINT
2028  // 5. empty constraint
2029  if (row.size() == 0)
2030  {
2031  MSG_INFO3( (*spxout), (*spxout) << "IMAISM21 row " << i
2032  << ": empty ->"; )
2033 
2034  if (LT(lp.rhs(i), 0.0, feastol()) || GT(lp.lhs(i), 0.0, feastol()))
2035  {
2036  MSG_INFO3( (*spxout), (*spxout) << " infeasible lhs=" << lp.lhs(i)
2037  << " rhs=" << lp.rhs(i) << std::endl; )
2038  return INFEASIBLE;
2039  }
2040  MSG_INFO3( (*spxout), (*spxout) << " removed" << std::endl; )
2041 
2042  EmptyConstraintPS* EmptyConstraintPSptr = 0;
2043  spx_alloc(EmptyConstraintPSptr);
2044  m_hist.append(new (EmptyConstraintPSptr) EmptyConstraintPS(lp, i));
2045 
2046  ++remRows;
2047  removeRow(lp, i);
2048 
2049  ++m_stat[EMPTY_ROW];
2050 
2051  continue;
2052  }
2053 #endif
2054 
2055 #if ROW_SINGLETON
2056  // 6. row singleton
2057  if (row.size() == 1)
2058  {
2059  removeRowSingleton(lp, row, i);
2060  continue;
2061  }
2062 #endif
2063 
2064 #if FORCE_CONSTRAINT
2065  // 7. forcing constraint (postsolving)
2066  // fix variables to obtain the upper bound on constraint value
2067  if (rhsCnt == 0 && EQrel(rhsBnd, lp.lhs(i), feastol()))
2068  {
2069  MSG_INFO3( (*spxout), (*spxout) << "IMAISM24 row " << i
2070  << ": forcing constraint fix on lhs ->"
2071  << " lhs=" << lp.lhs(i)
2072  << " rhsBnd=" << rhsBnd
2073  << std::endl; )
2074 
2075  DataArray<bool> fixedCol(row.size());
2076  DataArray<Real> lowers(row.size());
2077  DataArray<Real> uppers(row.size());
2078 
2079  for(int k = 0; k < row.size(); ++k)
2080  {
2081  Real aij = row.value(k);
2082  int j = row.index(k);
2083 
2084  fixedCol[k] = !(EQrel(lp.upper(j), lp.lower(j), m_epsilon));
2085 
2086  lowers[k] = lp.lower(j);
2087  uppers[k] = lp.upper(j);
2088 
2089  ASSERT_WARN( "WMAISM25", isNotZero(aij, 1.0 /infinity) );
2090 
2091  if (aij > 0.0)
2092  lp.changeLower(j, lp.upper(j));
2093  else
2094  lp.changeUpper(j, lp.lower(j));
2095  }
2096 
2097  ForceConstraintPS* ForceConstraintPSptr = 0;
2098  spx_alloc(ForceConstraintPSptr);
2099  m_hist.append(new (ForceConstraintPSptr) ForceConstraintPS(lp, i, true, fixedCol, lowers, uppers));
2100 
2101  ++remRows;
2102  remNzos += row.size();
2103  removeRow(lp, i);
2104 
2105  ++m_stat[FORCE_ROW];
2106 
2107  continue;
2108  }
2109  // fix variables to obtain the lower bound on constraint value
2110  if (lhsCnt == 0 && EQrel(lhsBnd, lp.rhs(i), feastol()))
2111  {
2112  MSG_INFO3( (*spxout), (*spxout) << "IMAISM26 row " << i
2113  << ": forcing constraint fix on rhs ->"
2114  << " rhs=" << lp.rhs(i)
2115  << " lhsBnd=" << lhsBnd
2116  << std::endl; )
2117 
2118  DataArray<bool> fixedCol(row.size());
2119  DataArray<Real> lowers(row.size());
2120  DataArray<Real> uppers(row.size());
2121 
2122  for(int k = 0; k < row.size(); ++k)
2123  {
2124  Real aij = row.value(k);
2125  int j = row.index(k);
2126 
2127  fixedCol[k] = !(EQrel(lp.upper(j), lp.lower(j), m_epsilon));
2128 
2129  lowers[k] = lp.lower(j);
2130  uppers[k] = lp.upper(j);
2131 
2132  ASSERT_WARN( "WMAISM27", isNotZero(aij, 1.0 / infinity) );
2133 
2134  if (aij > 0.0)
2135  lp.changeUpper(j, lp.lower(j));
2136  else
2137  lp.changeLower(j, lp.upper(j));
2138  }
2139 
2140  ForceConstraintPS* ForceConstraintPSptr = 0;
2141  spx_alloc(ForceConstraintPSptr);
2142  m_hist.append(new (ForceConstraintPSptr) ForceConstraintPS(lp, i, false, fixedCol, lowers, uppers));
2143 
2144  ++remRows;
2145  remNzos += row.size();
2146  removeRow(lp, i);
2147 
2148  ++m_stat[FORCE_ROW];
2149 
2150  continue;
2151  }
2152 #endif
2153  }
2154 
2155  assert(remRows > 0 || remNzos == 0);
2156 
2157  if (remRows + chgLRhs + chgBnds > 0)
2158  {
2159  m_remRows += remRows;
2160  m_remNzos += remNzos;
2161  m_chgLRhs += chgLRhs;
2162  m_chgBnds += chgBnds;
2163  m_keptBnds += keptBnds;
2164  m_keptLRhs += keptLRhs;
2165 
2166  MSG_INFO2( (*spxout), (*spxout) << "Simplifier (rows) removed "
2167  << remRows << " rows, "
2168  << remNzos << " non-zeros, "
2169  << chgBnds << " col bounds, "
2170  << chgLRhs << " row bounds; kept "
2171  << keptBnds << " column bounds, "
2172  << keptLRhs << " row bounds"
2173  << std::endl; )
2174  if( remRows > m_minReduction * oldRows )
2175  again = true;
2176  }
2177  return OKAY;
2178 }
2179 
2181 {
2182 
2183  // This method simplifies the columns of the LP.
2184  //
2185  // The following operations are done:
2186  // 1. detect empty columns and fix corresponding variables
2187  // 2. detect variables that are unconstrained from below or above
2188  // and fix corresponding variables or remove involved constraints
2189  // 3. fix variables
2190  // 4. use column singleton variables with zero objective to adjust constraint bounds
2191  // 5. free column singleton combined with doubleton equation are
2192  // used to make the column singleton variable free
2193  // 6. substitute (implied) free column singletons
2194 
2195  int remRows = 0;
2196  int remCols = 0;
2197  int remNzos = 0;
2198  int chgBnds = 0;
2199 
2200  int oldCols = lp.nCols();
2201  int oldRows = lp.nRows();
2202 
2203  for(int j = lp.nCols()-1; j >= 0; --j)
2204  {
2205  const SVector& col = lp.colVector(j);
2206 
2207  // infeasible bounds
2208  if (GTrel(lp.lower(j), lp.upper(j), feastol()))
2209  {
2210  MSG_INFO3( (*spxout), (*spxout) << "IMAISM29 col " << j
2211  << ": infeasible bounds on x" << j
2212  << " -> lower=" << lp.lower(j)
2213  << " upper=" << lp.upper(j)
2214  << std::endl; )
2215  return INFEASIBLE;
2216  }
2217 
2218  // 1. empty column
2219  if (col.size() == 0)
2220  {
2221 #if EMPTY_COLUMN
2222  MSG_INFO3( (*spxout), (*spxout) << "IMAISM30 col " << j
2223  << ": empty -> maxObj=" << lp.maxObj(j)
2224  << " lower=" << lp.lower(j)
2225  << " upper=" << lp.upper(j); )
2226 
2227  Real val;
2228 
2229  if (GT(lp.maxObj(j), 0.0, epsZero()))
2230  {
2231  if (lp.upper(j) >= infinity)
2232  {
2233  MSG_INFO3( (*spxout), (*spxout) << " unbounded" << std::endl; )
2234  return UNBOUNDED;
2235  }
2236  val = lp.upper(j);
2237  }
2238  else if (LT(lp.maxObj(j), 0.0, epsZero()))
2239  {
2240  if (lp.lower(j) <= -infinity)
2241  {
2242  MSG_INFO3( (*spxout), (*spxout) << " unbounded" << std::endl; )
2243  return UNBOUNDED;
2244  }
2245  val = lp.lower(j);
2246  }
2247  else
2248  {
2249  assert(isZero(lp.maxObj(j), epsZero()));
2250  // any value within the bounds is ok
2251  if (lp.lower(j) > -infinity)
2252  val = lp.lower(j);
2253  else if (lp.upper(j) < infinity)
2254  val = lp.upper(j);
2255  else
2256  val = 0.0;
2257  }
2258  MSG_INFO3( (*spxout), (*spxout) << " removed" << std::endl; )
2259 
2260  FixBoundsPS* FixBoundsPSptr = 0;
2261  FixVariablePS* FixVariablePSptr = 0;
2262  spx_alloc(FixBoundsPSptr);
2263  spx_alloc(FixVariablePSptr);
2264  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j, val));
2265  m_hist.append(new (FixVariablePSptr) FixVariablePS(lp, *this, j, val));
2266 
2267  ++remCols;
2268  removeCol(lp, j);
2269 
2270  ++m_stat[EMPTY_COL];
2271 
2272  continue;
2273 #endif
2274  }
2275 
2276  if (NErel(lp.lower(j), lp.upper(j), feastol()))
2277  {
2278  // will be set to false if any constraint implies a bound on the variable
2279  bool loFree = true;
2280  bool upFree = true;
2281 
2282  // 1. fix and remove variables
2283  for(int k = 0; k < col.size(); ++k)
2284  {
2285  if (!loFree && !upFree)
2286  break;
2287 
2288  int i = col.index(k);
2289 
2290  // warn since this unhandled case may slip through unnoticed otherwise
2291  ASSERT_WARN( "WMAISM31", isNotZero(col.value(k), 1.0 / infinity) );
2292 
2293  if (col.value(k) > 0.0)
2294  {
2295  if (lp.rhs(i) < infinity)
2296  upFree = false;
2297 
2298  if (lp.lhs(i) > -infinity)
2299  loFree = false;
2300  }
2301  else if (col.value(k) < 0.0)
2302  {
2303  if (lp.rhs(i) < infinity)
2304  loFree = false;
2305 
2306  if (lp.lhs(i) > -infinity)
2307  upFree = false;
2308  }
2309  }
2310 
2311  // 2. detect variables that are unconstrained from below or above
2312  // max 3 x
2313  // s.t. 5 x >= 8
2314  if (GT(lp.maxObj(j), 0.0, epsZero()) && upFree)
2315  {
2316 #if FIX_VARIABLE
2317  MSG_INFO3( (*spxout), (*spxout) << "IMAISM32 col " << j
2318  << ": x" << j
2319  << " unconstrained above ->"; )
2320 
2321  if (lp.upper(j) >= infinity)
2322  {
2323  MSG_INFO3( (*spxout), (*spxout) << " unbounded" << std::endl; )
2324 
2325  return UNBOUNDED;
2326  }
2327  MSG_INFO3( (*spxout), (*spxout) << " fixed at upper=" << lp.upper(j) << std::endl; )
2328 
2329  FixBoundsPS* FixBoundsPSptr = 0;
2330  spx_alloc(FixBoundsPSptr);
2331  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j, lp.upper(j)));
2332  lp.changeLower(j, lp.upper(j));
2333  }
2334  // max -3 x
2335  // s.t. 5 x <= 8
2336  else if (LT(lp.maxObj(j), 0.0, epsZero()) && loFree)
2337  {
2338  MSG_INFO3( (*spxout), (*spxout) << "IMAISM33 col " << j
2339  << ": x" << j
2340  << " unconstrained below ->"; )
2341 
2342  if (lp.lower(j) <= -infinity)
2343  {
2344  MSG_INFO3( (*spxout), (*spxout) << " unbounded" << std::endl; )
2345 
2346  return UNBOUNDED;
2347  }
2348  MSG_INFO3( (*spxout), (*spxout) << " fixed at lower=" << lp.lower(j) << std::endl; )
2349 
2350  FixBoundsPS* FixBoundsPSptr = 0;
2351  spx_alloc(FixBoundsPSptr);
2352  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j, lp.lower(j)));
2353  lp.changeUpper(j, lp.lower(j));
2354 #endif
2355  }
2356  else if (isZero(lp.maxObj(j), epsZero()))
2357  {
2358 #if FREE_ZERO_OBJ_VARIABLE
2359  bool unconstrained_below = loFree && lp.lower(j) <= -infinity;
2360  bool unconstrained_above = upFree && lp.upper(j) >= infinity;
2361 
2362  if (unconstrained_below || unconstrained_above)
2363  {
2364  MSG_INFO3( (*spxout), (*spxout) << "IMAISM34 col " << j
2365  << ": x" << j
2366  << " unconstrained "
2367  << (unconstrained_below ? "below" : "above")
2368  << " with zero objective (" << lp.maxObj(j)
2369  << ")" << std::endl; )
2370 
2371  SVector col_idx_sorted(col);
2372 
2373  // sort col elements by increasing idx
2374  IdxCompare compare;
2375  SPxQuicksort(col_idx_sorted.mem(), col_idx_sorted.size(), compare);
2376 
2377  FreeZeroObjVariablePS* FreeZeroObjVariablePSptr = 0;
2378  spx_alloc(FreeZeroObjVariablePSptr);
2379  m_hist.append(new (FreeZeroObjVariablePSptr) FreeZeroObjVariablePS(lp, j, unconstrained_below, col_idx_sorted));
2380 
2381  // we have to remove the rows with larger idx first, because otherwise the rows are reorder and indices
2382  // are out-of-date
2383  remRows += col.size();
2384  for(int k = col_idx_sorted.size()-1; k >= 0; --k)
2385  removeRow(lp, col_idx_sorted.index(k));
2386 
2387  // remove column
2388  removeCol(lp, j);
2389 
2390  // statistics
2391  for(int k = 0; k < col.size(); ++k)
2392  {
2393  int l = col.index(k);
2394  remNzos += lp.rowVector(l).size();
2395  }
2396 
2397  ++m_stat[FREE_ZOBJ_COL];
2398  ++remCols;
2399 
2400  continue;
2401  }
2402 #endif
2403  }
2404  }
2405 
2406 #if FIX_VARIABLE
2407  // 3. fix variable
2408  if (EQrel(lp.lower(j), lp.upper(j), feastol()))
2409  {
2410  MSG_INFO3( (*spxout), (*spxout) << "IMAISM36 col " << j
2411  << ": x" << j
2412  << " fixed -> lower=" << lp.lower(j)
2413  << " upper=" << lp.upper(j) << std::endl; )
2414 
2415  fixColumn(lp, j);
2416 
2417  ++remCols;
2418  remNzos += col.size();
2419  removeCol(lp, j);
2420 
2421  ++m_stat[FIX_COL];
2422 
2423  continue;
2424  }
2425 #endif
2426 
2427  // handle column singletons
2428  if (col.size() == 1)
2429  {
2430  Real aij = col.value(0);
2431  int i = col.index(0);
2432 
2433  // 4. column singleton with zero objective
2434  if (isZero(lp.maxObj(j), epsZero()))
2435  {
2436 #if ZERO_OBJ_COL_SINGLETON
2437  MSG_INFO3( (*spxout), (*spxout) << "IMAISM37 col " << j
2438  << ": singleton in row " << i
2439  << " with zero objective"; )
2440 
2441  Real lhs = -infinity;
2442  Real rhs = +infinity;
2443 
2444  if (GT(aij, 0.0, epsZero()))
2445  {
2446  if (lp.lhs(i) > -infinity && lp.upper(j) < infinity)
2447  lhs = lp.lhs(i) - aij * lp.upper(j);
2448  if (lp.rhs(i) < infinity && lp.lower(j) > -infinity)
2449  rhs = lp.rhs(i) - aij * lp.lower(j);
2450  }
2451  else if (LT(aij, 0.0, epsZero()))
2452  {
2453  if (lp.lhs(i) > -infinity && lp.lower(j) > -infinity)
2454  lhs = lp.lhs(i) - aij * lp.lower(j);
2455  if (lp.rhs(i) < infinity && lp.upper(j) < infinity)
2456  rhs = lp.rhs(i) - aij * lp.upper(j);
2457  }
2458  else
2459  {
2460  lhs = lp.lhs(i);
2461  rhs = lp.rhs(i);
2462  }
2463 
2464  if (isZero(lhs, epsZero()))
2465  lhs = 0.0;
2466  if (isZero(rhs, epsZero()))
2467  rhs = 0.0;
2468 
2469  MSG_INFO3( (*spxout), (*spxout) << " removed -> lhs=" << lhs
2470  << " (" << lp.lhs(i)
2471  << ") rhs=" << rhs
2472  << " (" << lp.rhs(i)
2473  << ")" << std::endl; )
2474 
2475  ZeroObjColSingletonPS* ZeroObjColSingletonPSptr = 0;
2476  spx_alloc(ZeroObjColSingletonPSptr);
2477  m_hist.append(new (ZeroObjColSingletonPSptr) ZeroObjColSingletonPS(lp, *this, j, i));
2478 
2479  lp.changeRange(i, lhs, rhs);
2480 
2481  ++remCols;
2482  ++remNzos;
2483  removeCol(lp, j);
2484 
2486 
2487  if (lp.lhs(i) <= -infinity && lp.rhs(i) >= infinity)
2488  {
2489  FreeConstraintPS* FreeConstraintPSptr = 0;
2490  spx_alloc(FreeConstraintPSptr);
2491  m_hist.append(new (FreeConstraintPSptr) FreeConstraintPS(lp, i));
2492 
2493  ++remRows;
2494  removeRow(lp, i);
2495 
2496  ++m_stat[FREE_ROW];
2497  }
2498 
2499  continue;
2500 #endif
2501  }
2502 
2503  // 5. not free column singleton combined with doubleton equation
2504  else if (EQrel(lp.lhs(i), lp.rhs(i), feastol()) &&
2505  lp.rowVector(i).size() == 2 &&
2506  (lp.lower(j) > -infinity || lp.upper(j) < infinity))
2507  {
2508 #if DOUBLETON_EQUATION
2509  MSG_INFO3( (*spxout), (*spxout) << "IMAISM38 col " << j
2510  << ": singleton in row " << i
2511  << " with doubleton equation ->"; )
2512 
2513  Real lhs = lp.lhs(i);
2514 
2515  const SVector& row = lp.rowVector(i);
2516 
2517  Real aik;
2518  int k;
2519 
2520  if (row.index(0) == j)
2521  {
2522  aik = row.value(1);
2523  k = row.index(1);
2524  }
2525  else if (row.index(1) == j)
2526  {
2527  aik = row.value(0);
2528  k = row.index(0);
2529  }
2530  else
2531  throw SPxInternalCodeException("XMAISM11 This should never happen.");
2532 
2533  ASSERT_WARN( "WMAISM39", isNotZero(aik, 1.0 / infinity) );
2534 
2535  Real lo, up;
2536  Real oldLower = lp.lower(k);
2537  Real oldUpper = lp.upper(k);
2538 
2539  Real scale1 = maxAbs(lhs, aij * lp.upper(j));
2540  Real scale2 = maxAbs(lhs, aij * lp.lower(j));
2541 
2542  if (scale1 < 1.0)
2543  scale1 = 1.0;
2544  if (scale2 < 1.0)
2545  scale2 = 1.0;
2546 
2547  Real z1 = (lhs / scale1) - (aij * lp.upper(j) / scale1);
2548  Real z2 = (lhs / scale2) - (aij * lp.lower(j) / scale2);
2549 
2550  if (isZero(z1, epsZero()))
2551  z1 = 0.0;
2552  if (isZero(z2, epsZero()))
2553  z2 = 0.0;
2554 
2555  if (GT(aij * aik, 0.0, epsZero()))
2556  {
2557  lo = (lp.upper(j) >= infinity) ? -infinity : z1 * scale1 / aik;
2558  up = (lp.lower(j) <= -infinity) ? infinity : z2 * scale2 / aik;
2559  }
2560  else if (LT(aij * aik, 0.0, epsZero()))
2561  {
2562  lo = (lp.lower(j) <= -infinity) ? -infinity : z2 * scale2 / aik;
2563  up = (lp.upper(j) >= infinity) ? infinity : z1 * scale1 / aik;
2564  }
2565  else
2566  throw SPxInternalCodeException("XMAISM12 This should never happen.");
2567 
2568  if (GTrel(lo, lp.lower(k), epsZero()))
2569  lp.changeLower(k, lo);
2570 
2571  if (LTrel(up, lp.upper(k), epsZero()))
2572  lp.changeUpper(k, up);
2573 
2574  MSG_INFO3( (*spxout), (*spxout) << " made free, bounds on x" << k
2575  << ": lower=" << lp.lower(k)
2576  << " (" << oldLower
2577  << ") upper=" << lp.upper(k)
2578  << " (" << oldUpper
2579  << ")" << std::endl; )
2580 
2581  // infeasible bounds
2582  if (GTrel(lp.lower(k), lp.upper(k), feastol()))
2583  {
2584  MSG_INFO3( (*spxout), (*spxout) << "new bounds are infeasible"
2585  << std::endl; )
2586  return INFEASIBLE;
2587  }
2588 
2589  DoubletonEquationPS* DoubletonEquationPSptr = 0;
2590  spx_alloc(DoubletonEquationPSptr);
2591  m_hist.append(new (DoubletonEquationPSptr) DoubletonEquationPS(lp, j, k, i, oldLower, oldUpper));
2592 
2593  if (lp.lower(j) > -infinity && lp.upper(j) < infinity)
2594  chgBnds += 2;
2595  else
2596  ++chgBnds;
2597 
2598  lp.changeBounds(j, -infinity, infinity);
2599 
2600  ++m_stat[DOUBLETON_ROW];
2601 #endif
2602  }
2603 
2604  // 6. (implied) free column singleton
2605  if (lp.lower(j) <= -infinity && lp.upper(j) >= infinity)
2606  {
2607 #if FREE_COL_SINGLETON
2608  Real slackVal = lp.lhs(i);
2609 
2610  // constraint i is an inequality constraint -> transform into equation type
2611  if (NErel(lp.lhs(i), lp.rhs(i), feastol()))
2612  {
2613  MSG_INFO3( (*spxout), (*spxout) << "IMAISM40 col " << j
2614  << ": free singleton in inequality constraint" << std::endl; )
2615 
2616  // do nothing if constraint i is unconstrained
2617  if (lp.lhs(i) <= -infinity && lp.rhs(i) >= infinity)
2618  continue;
2619 
2620  // introduce slack variable to obtain equality constraint
2621  Real sMaxObj = lp.maxObj(j) / aij; // after substituting variable j in objective
2622  Real sLo = lp.lhs(i);
2623  Real sUp = lp.rhs(i);
2624 
2625  if (GT(sMaxObj, 0.0, epsZero()))
2626  {
2627  if (sUp >= infinity)
2628  {
2629  MSG_INFO3( (*spxout), (*spxout) << " -> problem unbounded" << std::endl; )
2630  return UNBOUNDED;
2631  }
2632  slackVal = sUp;
2633  }
2634  else if (LT(sMaxObj, 0.0, epsZero()))
2635  {
2636  if (sLo <= -infinity)
2637  {
2638  MSG_INFO3( (*spxout), (*spxout) << " -> problem unbounded" << std::endl; )
2639  return UNBOUNDED;
2640  }
2641  slackVal = sLo;
2642  }
2643  else
2644  {
2645  assert(isZero(sMaxObj, epsZero()));
2646  // any value within the bounds is ok
2647  if (sLo > -infinity)
2648  slackVal = sLo;
2649  else if (sUp < infinity)
2650  slackVal = sUp;
2651  else
2652  throw SPxInternalCodeException("XMAISM13 This should never happen.");
2653  }
2654  }
2655 
2656  FreeColSingletonPS* FreeColSingletonPSptr = 0;
2657  spx_alloc(FreeColSingletonPSptr);
2658  m_hist.append(new (FreeColSingletonPSptr) FreeColSingletonPS(lp, *this, j, i, slackVal));
2659 
2660  MSG_INFO3( (*spxout), (*spxout) << "IMAISM41 col " << j
2661  << ": free singleton removed" << std::endl; )
2662 
2663  const SVector& row = lp.rowVector(i);
2664 
2665  for (int h = 0; h < row.size(); ++h)
2666  {
2667  int k = row.index(h);
2668 
2669  if (k != j)
2670  {
2671  Real new_obj = lp.obj(k) - (lp.obj(j) * row.value(h) / aij);
2672  lp.changeObj(k, new_obj);
2673  }
2674  }
2675 
2676  ++remRows;
2677  ++remCols;
2678  remNzos += row.size();
2679  removeRow(lp, i);
2680  removeCol(lp, j);
2681 
2683 #endif
2684  }
2685  }
2686  }
2687 
2688  if (remCols + remRows > 0)
2689  {
2690  m_remRows += remRows;
2691  m_remCols += remCols;
2692  m_remNzos += remNzos;
2693  m_chgBnds += chgBnds;
2694 
2695  MSG_INFO2( (*spxout), (*spxout) << "Simplifier (columns) removed "
2696  << remRows << " rows, "
2697  << remCols << " cols, "
2698  << remNzos << " non-zeros, "
2699  << chgBnds << " col bounds"
2700  << std::endl; )
2701  if( remCols + remRows > m_minReduction * (oldCols + oldRows) )
2702  again = true;
2703  }
2704  return OKAY;
2705 }
2706 
2708 {
2709 
2710  // This method simplifies LP using the following dual structures:
2711  //
2712  // 1. dominated columns
2713  // 2. weakly dominated columns
2714  //
2715  // For constructing the dual variables, it is assumed that the objective sense is max
2716 
2717  int remRows = 0;
2718  int remCols = 0;
2719  int remNzos = 0;
2720 
2721  int oldRows = lp.nRows();
2722  int oldCols = lp.nCols();
2723 
2724  DataArray<bool> colSingleton(lp.nCols());
2725  DVector dualVarLo(lp.nRows());
2726  DVector dualVarUp(lp.nRows());
2727  DVector dualConsLo(lp.nCols());
2728  DVector dualConsUp(lp.nCols());
2729 
2730  // init
2731  for(int i = lp.nRows()-1; i >= 0; --i)
2732  {
2733  // check for unconstrained constraints
2734  if (lp.lhs(i) <= -infinity && lp.rhs(i) >= infinity)
2735  {
2736  MSG_INFO3( (*spxout), (*spxout) << "IMAISM43 row " << i
2737  << ": unconstrained" << std::endl; )
2738 
2739  FreeConstraintPS* FreeConstraintPSptr = 0;
2740  spx_alloc(FreeConstraintPSptr);
2741  m_hist.append(new (FreeConstraintPSptr) FreeConstraintPS(lp, i));
2742 
2743  ++remRows;
2744  remNzos += lp.rowVector(i).size();
2745  removeRow(lp, i);
2746 
2747  ++m_stat[FREE_ROW];
2748 
2749  continue;
2750  }
2751 
2752  // corresponds to maximization sense
2753  dualVarLo[i] = (lp.lhs(i) <= -infinity) ? 0.0 : -infinity;
2754  dualVarUp[i] = (lp.rhs(i) >= infinity) ? 0.0 : infinity;
2755  }
2756 
2757  // compute bounds on the dual variables using column singletons
2758  for(int j = 0; j < lp.nCols(); ++j)
2759  {
2760  if (lp.colVector(j).size() == 1)
2761  {
2762  int i = lp.colVector(j).index(0);
2763  Real aij = lp.colVector(j).value(0);
2764 
2765  ASSERT_WARN( "WMAISM44", isNotZero(aij, 1.0 / infinity) );
2766 
2767  Real bound = lp.maxObj(j) / aij;
2768 
2769  if (aij > 0)
2770  {
2771  if (lp.lower(j) <= -infinity && bound < dualVarUp[i])
2772  dualVarUp[i] = bound;
2773  if (lp.upper(j) >= infinity && bound > dualVarLo[i])
2774  dualVarLo[i] = bound;
2775  }
2776  else if (aij < 0)
2777  {
2778  if (lp.lower(j) <= -infinity && bound > dualVarLo[i])
2779  dualVarLo[i] = bound;
2780  if (lp.upper(j) >= infinity && bound < dualVarUp[i])
2781  dualVarUp[i] = bound;
2782  }
2783  }
2784  }
2785 
2786  // compute bounds on the dual constraints
2787  for(int j = 0; j < lp.nCols(); ++j)
2788  {
2789  dualConsLo[j] = dualConsUp[j] = 0.0;
2790 
2791  const SVector& col = lp.colVector(j);
2792 
2793  for(int k = 0; k < col.size(); ++k)
2794  {
2795  if (dualConsLo[j] <= -infinity && dualConsUp[j] >= infinity)
2796  break;
2797 
2798  Real aij = col.value(k);
2799  int i = col.index(k);
2800 
2801  ASSERT_WARN( "WMAISM45", isNotZero(aij, 1.0 / infinity) );
2802 
2803  if (aij > 0)
2804  {
2805  if (dualVarLo[i] <= -infinity)
2806  dualConsLo[j] = -infinity;
2807  else
2808  dualConsLo[j] += aij * dualVarLo[i];
2809 
2810  if (dualVarUp[i] >= infinity)
2811  dualConsUp[j] = infinity;
2812  else
2813  dualConsUp[j] += aij * dualVarUp[i];
2814  }
2815  else if (aij < 0)
2816  {
2817  if (dualVarLo[i] <= -infinity)
2818  dualConsUp[j] = infinity;
2819  else
2820  dualConsUp[j] += aij * dualVarLo[i];
2821 
2822  if (dualVarUp[i] >= infinity)
2823  dualConsLo[j] = -infinity;
2824  else
2825  dualConsLo[j] += aij * dualVarUp[i];
2826  }
2827  }
2828  }
2829 
2830  for(int j = lp.nCols()-1; j >= 0; --j)
2831  {
2832  if (lp.colVector(j).size() <= 1)
2833  continue;
2834 
2835  // dual infeasibility checks
2836  if (LTrel(dualConsUp[j], dualConsLo[j], opttol()))
2837  {
2838  MSG_INFO3( (*spxout), (*spxout) << "IMAISM46 col " << j
2839  << ": dual infeasible -> dual lhs bound=" << dualConsLo[j]
2840  << " dual rhs bound=" << dualConsUp[j] << std::endl; )
2841  return DUAL_INFEASIBLE;
2842  }
2843 
2844  Real obj = lp.maxObj(j);
2845 
2846  // 1. dominated column
2847  // Is the problem really unbounded in the cases below ??? Or is only dual infeasibility be shown
2848  if (GTrel(obj, dualConsUp[j], opttol()))
2849  {
2850 #if DOMINATED_COLUMN
2851  MSG_INFO3( (*spxout), (*spxout) << "IMAISM47 col " << j
2852  << ": dominated -> maxObj=" << obj
2853  << " dual rhs bound=" << dualConsUp[j] << std::endl; )
2854 
2855  if (lp.upper(j) >= infinity)
2856  {
2857  MSG_INFO2( (*spxout), (*spxout) << " unbounded" << std::endl; )
2858  return UNBOUNDED;
2859  }
2860 
2861  MSG_INFO3( (*spxout), (*spxout) << " fixed at upper=" << lp.upper(j) << std::endl; )
2862 
2863  FixBoundsPS* FixBoundsPSptr = 0;
2864  spx_alloc(FixBoundsPSptr);
2865  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j, lp.upper(j)));
2866  lp.changeLower(j, lp.upper(j));
2867 
2868  ++m_stat[DOMINATED_COL];
2869 #endif
2870  }
2871  else if (LTrel(obj, dualConsLo[j], opttol()))
2872  {
2873 #if DOMINATED_COLUMN
2874  MSG_INFO3( (*spxout), (*spxout) << "IMAISM48 col " << j
2875  << ": dominated -> maxObj=" << obj
2876  << " dual lhs bound=" << dualConsLo[j] << std::endl; )
2877 
2878  if (lp.lower(j) <= -infinity)
2879  {
2880  MSG_INFO2( (*spxout), (*spxout) << " unbounded" << std::endl; )
2881  return UNBOUNDED;
2882  }
2883 
2884  MSG_INFO3( (*spxout), (*spxout) << " fixed at lower=" << lp.lower(j) << std::endl; )
2885 
2886  FixBoundsPS* FixBoundsPSptr = 0;
2887  spx_alloc(FixBoundsPSptr);
2888  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j, lp.lower(j)));
2889  lp.changeUpper(j, lp.lower(j));
2890 
2891  ++m_stat[DOMINATED_COL];
2892 #endif
2893  }
2894 
2895  // 2. weakly dominated column (no postsolving)
2896  else if (lp.upper(j) < infinity && EQrel(obj, dualConsUp[j], opttol()))
2897  {
2898 #if WEAKLY_DOMINATED_COLUMN
2899  MSG_INFO3( (*spxout), (*spxout) << "IMAISM49 col " << j
2900  << ": weakly dominated -> maxObj=" << obj
2901  << " dual rhs bound=" << dualConsUp[j] << std::endl; )
2902 
2903  FixBoundsPS* FixBoundsPSptr = 0;
2904  spx_alloc(FixBoundsPSptr);
2905  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j, lp.upper(j)));
2906  lp.changeLower(j, lp.upper(j));
2907 
2909 #endif
2910  }
2911  else if (lp.lower(j) > -infinity && EQrel(obj, dualConsLo[j], opttol()))
2912  {
2913 #if WEAKLY_DOMINATED_COLUMN
2914  MSG_INFO3( (*spxout), (*spxout) << "IMAISM50 col " << j
2915  << ": weakly dominated -> maxObj=" << obj
2916  << " dual lhs bound=" << dualConsLo[j] << std::endl; )
2917 
2918  FixBoundsPS* FixBoundsPSptr = 0;
2919  spx_alloc(FixBoundsPSptr);
2920  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j, lp.lower(j)));
2921  lp.changeUpper(j, lp.lower(j));
2922 
2924 #endif
2925  }
2926 
2927  // fix column
2928  if (EQrel(lp.lower(j), lp.upper(j), feastol()))
2929  {
2930 #if FIX_VARIABLE
2931  fixColumn(lp, j);
2932 
2933  ++remCols;
2934  remNzos += lp.colVector(j).size();
2935  removeCol(lp, j);
2936 
2937  ++m_stat[FIX_COL];
2938 #endif
2939  }
2940  }
2941 
2942  assert(remRows > 0 || remCols > 0 || remNzos == 0);
2943 
2944  if (remCols + remRows > 0)
2945  {
2946  m_remRows += remRows;
2947  m_remCols += remCols;
2948  m_remNzos += remNzos;
2949 
2950  MSG_INFO2( (*spxout), (*spxout) << "Simplifier (dual) removed "
2951  << remRows << " rows, "
2952  << remCols << " cols, "
2953  << remNzos << " non-zeros"
2954  << std::endl; )
2955  if( remCols + remRows > m_minReduction * (oldCols + oldRows) )
2956  again = true;
2957  }
2958  return OKAY;
2959 }
2960 
2962 {
2963 
2964  // This method simplifies the LP by removing duplicate rows
2965  // Duplicates are detected using the algorithm of Bixby and Wagner [1987]
2966 
2967  // Possible extension: use generalized definition of duplicate rows according to Andersen and Andersen
2968  // However: the resulting sparsification is often very small since the involved rows are usually very sparse
2969 
2970  int remRows = 0;
2971  int remNzos = 0;
2972 
2973  int oldRows = lp.nRows();
2974 
2975  // remove empty rows and columns
2977  if (ret != OKAY)
2978  return ret;
2979 
2980 #if ROW_SINGLETON
2981  int rs_remRows = 0;
2982  for (int i = 0; i < lp.nRows(); ++i)
2983  {
2984  const SVector& row = lp.rowVector(i);
2985 
2986  if (row.size() == 1)
2987  {
2988  removeRowSingleton(lp, row, i);
2989  rs_remRows++;
2990  }
2991  }
2992 
2993  if (rs_remRows > 0)
2994  {
2995  MSG_INFO2( (*spxout), (*spxout) << "Simplifier duplicate rows (row singleton stage) removed "
2996  << rs_remRows << " rows, "
2997  << rs_remRows << " non-zeros"
2998  << std::endl; )
2999  }
3000 #endif
3001 
3002  if (lp.nRows() < 2)
3003  return OKAY;
3004 
3005  DataArray<int> pClass(lp.nRows()); // class of parallel rows
3006  DataArray<int> classSize(lp.nRows()); // size of each class
3007  DataArray<Real> scale(lp.nRows()); // scaling factor for each row
3008  int* idxMem = 0;
3009 
3010  try
3011  {
3012  spx_alloc(idxMem, lp.nRows());
3013  }
3014  catch( const SPxMemoryException& x )
3015  {
3016  spx_free(idxMem);
3017  throw x;
3018  }
3019 
3020  IdxSet idxSet(lp.nRows(), idxMem); // set of feasible indices for new pClass
3021 
3022  // init
3023  pClass[0] = 0;
3024  scale[0] = 0.0;
3025  classSize[0] = lp.nRows();
3026 
3027  for(int i = 1; i < lp.nRows(); ++i)
3028  {
3029  pClass[i] = 0;
3030  scale[i] = 0.0;
3031  classSize[i] = 0;
3032  idxSet.addIdx(i);
3033  }
3034 
3035  Real oldVal;
3036 
3037  // main loop
3038  for(int j = 0; j < lp.nCols(); ++j)
3039  {
3040  const SVector& col = lp.colVector(j);
3041 
3042  for(int k = 0; k < col.size(); ++k)
3043  {
3044  Real aij = col.value(k);
3045  int i = col.index(k);
3046 
3047  if (scale[i] == 0.0)
3048  scale[i] = aij;
3049 
3050  m_classSetRows[pClass[i]].add(i, aij / scale[i]);
3051  if (--classSize[pClass[i]] == 0)
3052  idxSet.addIdx(pClass[i]);
3053  }
3054 
3055  // update each parallel class with non-zero column entry
3056  for(int m = 0; m < col.size(); ++m)
3057  {
3058  int k = pClass[col.index(m)];
3059 
3060  if (m_classSetRows[k].size() > 0)
3061  {
3062  // sort classSet[k] w.r.t. scaled column values
3063  ElementCompare compare;
3064 
3065  if (m_classSetRows[k].size() > 1)
3066  SPxQuicksort(m_classSetRows[k].mem(), m_classSetRows[k].size(), compare);
3067 
3068  // use new index first
3069  int classIdx = idxSet.index(0);
3070  idxSet.remove(0);
3071 
3072  for(int l = 0; l < m_classSetRows[k].size(); ++l)
3073  {
3074  if (l != 0 && NErel(m_classSetRows[k].value(l), oldVal, epsZero()))
3075  {
3076  classIdx = idxSet.index(0);
3077  idxSet.remove(0);
3078  }
3079 
3080  pClass[m_classSetRows[k].index(l)] = classIdx;
3081  ++classSize[classIdx];
3082 
3083  oldVal = m_classSetRows[k].value(l);
3084  }
3085 
3086  m_classSetRows[k].clear();
3087  }
3088  }
3089  }
3090 
3091  spx_free(idxMem);
3092 
3093  DataArray<bool> remRow(lp.nRows());
3094 
3095  for(int k = 0; k < lp.nRows(); ++k )
3096  m_dupRows[k].clear();
3097 
3098  for(int k = 0; k < lp.nRows(); ++k)
3099  {
3100  remRow[k] = false;
3101  m_dupRows[pClass[k]].add(k, 0.0);
3102  }
3103 
3104  const int nRowsOld_tmp = lp.nRows();
3105  int* perm_tmp = 0;
3106  spx_alloc(perm_tmp, nRowsOld_tmp);
3107 
3108  for(int j = 0; j < nRowsOld_tmp; ++j)
3109  {
3110  perm_tmp[j] = 0;
3111  }
3112 
3113  int idxFirstDupRows = -1;
3114  int idxLastDupRows = -1;
3115  int numDelRows = 0;
3116 
3117  for(int k = 0; k < lp.nRows(); ++k)
3118  {
3119  if (m_dupRows[k].size() > 1 && !(lp.rowVector(m_dupRows[k].index(0)).size() == 1))
3120  {
3121  idxLastDupRows = k;
3122 
3123  if(idxFirstDupRows < 0)
3124  {
3125  idxFirstDupRows = k;
3126  }
3127 
3128  for(int l = 1; l < m_dupRows[k].size(); ++l)
3129  {
3130  int i = m_dupRows[k].index(l);
3131  perm_tmp[i] = -1;
3132  }
3133 
3134  numDelRows += (m_dupRows[k].size()-1);
3135  }
3136  }
3137 
3138  {
3139  int k_tmp, j_tmp = -1;
3140  for (k_tmp = j_tmp = 0; k_tmp < nRowsOld_tmp; ++k_tmp)
3141  {
3142  if (perm_tmp[k_tmp] >= 0)
3143  perm_tmp[k_tmp] = j_tmp++;
3144  }
3145  }
3146 
3147  // store rhs and lhs changes for combined update
3148  bool doChangeRanges = false;
3149  DVector newLhsVec(lp.lhs());
3150  DVector newRhsVec(lp.rhs());
3151 
3152  for(int k = 0; k < lp.nRows(); ++k)
3153  {
3154  if (m_dupRows[k].size() > 1 && !(lp.rowVector(m_dupRows[k].index(0)).size() == 1))
3155  {
3156  MSG_INFO3( (*spxout), (*spxout) << "IMAISM53 " << m_dupRows[k].size()
3157  << " duplicate rows found" << std::endl; )
3158 
3159  m_stat[DUPLICATE_ROW] += m_dupRows[k].size()-1;
3160 
3161  // index of one non-column singleton row in dupRows[k]
3162  int rowIdx = -1;
3163  int maxLhsIdx = -1;
3164  int minRhsIdx = -1;
3165  Real maxLhs = -infinity;
3166  Real minRhs = +infinity;
3167 
3168  DataArray<bool> isLhsEqualRhs(m_dupRows[k].size());
3169 
3170  // determine strictest bounds on constraint
3171  for(int l = 0; l < m_dupRows[k].size(); ++l)
3172  {
3173  int i = m_dupRows[k].index(l);
3174  isLhsEqualRhs[l] = (lp.lhs(i) == lp.rhs(i));
3175 
3176  ASSERT_WARN( "WMAISM54", isNotZero(scale[i], 1.0 / infinity) );
3177 
3178  if (rowIdx == -1)
3179  {
3180  rowIdx = i;
3181  maxLhs = lp.lhs(rowIdx);
3182  minRhs = lp.rhs(rowIdx);
3183  }
3184  else
3185  {
3186  Real scaledLhs, scaledRhs;
3187  Real factor = scale[rowIdx] / scale[i];
3188 
3189  if (factor > 0)
3190  {
3191  scaledLhs = (lp.lhs(i) <= -infinity) ? -infinity : lp.lhs(i) * factor;
3192  scaledRhs = (lp.rhs(i) >= infinity) ? infinity : lp.rhs(i) * factor;
3193  }
3194  else
3195  {
3196  scaledLhs = (lp.rhs(i) >= infinity) ? -infinity : lp.rhs(i) * factor;
3197  scaledRhs = (lp.lhs(i) <= -infinity) ? infinity : lp.lhs(i) * factor;
3198  }
3199  if (scaledLhs > maxLhs)
3200  {
3201  maxLhs = scaledLhs;
3202  maxLhsIdx = i;
3203  }
3204  if (scaledRhs < minRhs)
3205  {
3206  minRhs = scaledRhs;
3207  minRhsIdx = i;
3208  }
3209 
3210  remRow[i] = true;
3211  }
3212  }
3213 
3214  if (rowIdx != -1)
3215  {
3216  Real newLhs = (maxLhs > lp.lhs(rowIdx)) ? maxLhs : lp.lhs(rowIdx);
3217  Real newRhs = (minRhs < lp.rhs(rowIdx)) ? minRhs : lp.rhs(rowIdx);
3218 
3219  if(k == idxLastDupRows)
3220  {
3221  DataArray<int> da_perm(nRowsOld_tmp);
3222  for(int j = 0; j < nRowsOld_tmp; ++j)
3223  {
3224  da_perm[j] = perm_tmp[j];
3225  }
3226  DuplicateRowsPS* DuplicateRowsPSptr = 0;
3227  spx_alloc(DuplicateRowsPSptr);
3228  m_hist.append(new (DuplicateRowsPSptr) DuplicateRowsPS(lp, rowIdx, maxLhsIdx, minRhsIdx, m_dupRows[k], scale, da_perm, isLhsEqualRhs, true, EQrel(newLhs, newRhs), k==idxFirstDupRows));
3229  }
3230  else
3231  {
3232  DataArray<int> da_perm_empty(0);
3233  DuplicateRowsPS* DuplicateRowsPSptr = 0;
3234  spx_alloc(DuplicateRowsPSptr);
3235  m_hist.append(new (DuplicateRowsPSptr) DuplicateRowsPS(lp, rowIdx, maxLhsIdx, minRhsIdx, m_dupRows[k], scale, da_perm_empty, isLhsEqualRhs, false, EQrel(newLhs, newRhs), k == idxFirstDupRows));
3236  }
3237 
3238  if (maxLhs > lp.lhs(rowIdx) || minRhs < lp.rhs(rowIdx))
3239  {
3240  // modify lhs and rhs of constraint rowIdx
3241  doChangeRanges = true;
3242 
3243  if (LTrel(newRhs, newLhs, feastol()))
3244  {
3245  MSG_INFO3( (*spxout), (*spxout) << "IMAISM55 duplicate rows yield infeasible bounds:"
3246  << " lhs=" << newLhs
3247  << " rhs=" << newRhs << std::endl; )
3248  spx_free(perm_tmp);
3249  return INFEASIBLE;
3250  }
3251  // if we accept the infeasibility we should clean up the values to avoid problems later
3252  if( newRhs < newLhs )
3253  newRhs = newLhs;
3254 
3255  newLhsVec[rowIdx] = newLhs;
3256  newRhsVec[rowIdx] = newRhs;
3257  }
3258  }
3259  }
3260  }
3261 
3262  // change ranges for all modified constraints by one single call (more efficient)
3263  if (doChangeRanges)
3264  {
3265  lp.changeRange(newLhsVec, newRhsVec);
3266  }
3267 
3268  // remove all rows by one single method call (more efficient)
3269  const int nRowsOld = lp.nRows();
3270  int* perm = 0;
3271  spx_alloc(perm, nRowsOld);
3272 
3273  for(int i = 0; i < nRowsOld; ++i)
3274  {
3275  if (remRow[i])
3276  {
3277  perm[i] = -1;
3278  ++remRows;
3279  remNzos += lp.rowVector(i).size();
3280  }
3281  else
3282  perm[i] = 0;
3283  }
3284  lp.removeRows(perm);
3285 
3286  for(int i = 0; i < nRowsOld; ++i)
3287  {
3288  // assert that the pre-computed permutation was correct
3289  assert(perm[i] == perm_tmp[i]);
3290 
3291  // update the global index mapping
3292  if (perm[i] >= 0)
3293  m_rIdx[perm[i]] = m_rIdx[i];
3294  }
3295 
3296  spx_free(perm);
3297  spx_free(perm_tmp);
3298 
3299  if (remRows + remNzos > 0)
3300  {
3301  m_remRows += remRows;
3302  m_remNzos += remNzos;
3303 
3304  MSG_INFO2( (*spxout), (*spxout) << "Simplifier (duplicate rows) removed "
3305  << remRows << " rows, "
3306  << remNzos << " non-zeros"
3307  << std::endl; )
3308  if( remRows > m_minReduction * oldRows )
3309  again = true;
3310 
3311  }
3312  return OKAY;
3313 }
3314 
3316 {
3317 
3318  // This method simplifies the LP by removing duplicate columns
3319  // Duplicates are detected using the algorithm of Bixby and Wagner [1987]
3320 
3321  int remCols = 0;
3322  int remNzos = 0;
3323 
3324  // remove empty rows and columns
3326  if (ret != OKAY)
3327  return ret;
3328 
3329  if (lp.nCols() < 2)
3330  return OKAY;
3331 
3332  DataArray<int> pClass(lp.nCols()); // class of parallel columns
3333  DataArray<int> classSize(lp.nCols()); // size of each class
3334  DataArray<Real> scale(lp.nCols()); // scaling factor for each column
3335  int* idxMem = 0;
3336 
3337  try
3338  {
3339  spx_alloc(idxMem, lp.nCols());
3340  }
3341  catch( const SPxMemoryException& x )
3342  {
3343  spx_free(idxMem);
3344  throw x;
3345  }
3346 
3347  IdxSet idxSet(lp.nCols(), idxMem); // set of feasible indices for new pClass
3348 
3349  // init
3350  pClass[0] = 0;
3351  scale[0] = 0.0;
3352  classSize[0] = lp.nCols();
3353 
3354  for(int j = 1; j < lp.nCols(); ++j)
3355  {
3356  pClass[j] = 0;
3357  scale[j] = 0.0;
3358  classSize[j] = 0;
3359  idxSet.addIdx(j);
3360  }
3361 
3362  Real oldVal;
3363 
3364  // main loop
3365  for(int i = 0; i < lp.nRows(); ++i)
3366  {
3367  const SVector& row = lp.rowVector(i);
3368 
3369  for(int k = 0; k < row.size(); ++k)
3370  {
3371  Real aij = row.value(k);
3372  int j = row.index(k);
3373 
3374  if (scale[j] == 0.0)
3375  scale[j] = aij;
3376 
3377  m_classSetCols[pClass[j]].add(j, aij / scale[j]);
3378  if (--classSize[pClass[j]] == 0)
3379  idxSet.addIdx(pClass[j]);
3380  }
3381 
3382  // update each parallel class with non-zero row entry
3383  for(int m = 0; m < row.size(); ++m)
3384  {
3385  int k = pClass[row.index(m)];
3386 
3387  if (m_classSetCols[k].size() > 0)
3388  {
3389  // sort classSet[k] w.r.t. scaled row values
3390  ElementCompare compare;
3391 
3392  if (m_classSetCols[k].size() > 1)
3393  SPxQuicksort(m_classSetCols[k].mem(), m_classSetCols[k].size(), compare);
3394 
3395  // use new index first
3396  int classIdx = idxSet.index(0);
3397  idxSet.remove(0);
3398 
3399  for(int l = 0; l < m_classSetCols[k].size(); ++l)
3400  {
3401  if (l != 0 && NErel(m_classSetCols[k].value(l), oldVal, epsZero()))
3402  {
3403  // start new parallel class
3404  classIdx = idxSet.index(0);
3405  idxSet.remove(0);
3406  }
3407 
3408  pClass[m_classSetCols[k].index(l)] = classIdx;
3409  ++classSize[classIdx];
3410 
3411  oldVal = m_classSetCols[k].value(l);
3412  }
3413 
3414  m_classSetCols[k].clear();
3415  }
3416  }
3417  }
3418 
3419  spx_free(idxMem);
3420 
3421  DataArray<bool> remCol(lp.nCols());
3422  DataArray<bool> fixAndRemCol(lp.nCols());
3423 
3424  for( int k = 0; k < lp.nCols(); ++k )
3425  m_dupCols[k].clear();
3426 
3427  for(int k = 0; k < lp.nCols(); ++k)
3428  {
3429  remCol[k] = false;
3430  fixAndRemCol[k] = false;
3431  m_dupCols[pClass[k]].add(k, 0.0);
3432  }
3433 
3434  bool hasDuplicateCol = false;
3435  DataArray<int> m_perm_empty(0);
3436 
3437  for(int k = 0; k < lp.nCols(); ++k)
3438  {
3439  if (m_dupCols[k].size() > 1 && !(lp.colVector(m_dupCols[k].index(0)).size() == 1))
3440  {
3441  MSG_INFO3( (*spxout), (*spxout) << "IMAISM58 " << m_dupCols[k].size()
3442  << " duplicate columns found" << std::endl; )
3443 
3444  if (!hasDuplicateCol)
3445  {
3446  DuplicateColsPS* DuplicateColsPSptr = 0;
3447  spx_alloc(DuplicateColsPSptr);
3448  m_hist.append(new (DuplicateColsPSptr) DuplicateColsPS(lp, 0, 0, 1.0, m_perm_empty, true));
3449  hasDuplicateCol = true;
3450  }
3451 
3452  for(int l = 0; l < m_dupCols[k].size(); ++l)
3453  {
3454  for(int m = 0; m < m_dupCols[k].size(); ++m)
3455  {
3456  int j1 = m_dupCols[k].index(l);
3457  int j2 = m_dupCols[k].index(m);
3458 
3459  if (l != m && !remCol[j1] && !remCol[j2])
3460  {
3461  Real cj1 = lp.maxObj(j1);
3462  Real cj2 = lp.maxObj(j2);
3463 
3464  // A.j1 = factor * A.j2
3465  Real factor = scale[j1] / scale[j2];
3466  Real objDif = cj1 - cj2 * scale[j1] / scale[j2];
3467 
3468  ASSERT_WARN( "WMAISM59", isNotZero(factor, epsZero()) );
3469 
3470  if (isZero(objDif, epsZero()))
3471  {
3472  // case 1: objectives also duplicate
3473 
3474  // if 0 is not within the column bounds, we are not able to postsolve if the aggregated column has
3475  // status ZERO, hence we skip this case
3476  if (LErel(lp.lower(j1), 0.0) && GErel(lp.upper(j1), 0.0)
3477  && LErel(lp.lower(j2), 0.0) && GErel(lp.upper(j2), 0.0))
3478  {
3479  DuplicateColsPS* DuplicateColsPSptr = 0;
3480  spx_alloc(DuplicateColsPSptr);
3481  // variable substitution xj2' := xj2 + factor * xj1 <=> xj2 = -factor * xj1 + xj2'
3482  m_hist.append(new (DuplicateColsPSptr) DuplicateColsPS(lp, j1, j2, factor, m_perm_empty));
3483 
3484  // update bounds of remaining column j2 (new column j2')
3485  if (factor > 0)
3486  {
3487  if (lp.lower(j2) <= -infinity || lp.lower(j1) <= -infinity)
3488  lp.changeLower(j2, -infinity);
3489  else
3490  lp.changeLower(j2, lp.lower(j2) + factor * lp.lower(j1));
3491 
3492  if (lp.upper(j2) >= infinity || lp.upper(j1) >= infinity)
3493  lp.changeUpper(j2, infinity);
3494  else
3495  lp.changeUpper(j2, lp.upper(j2) + factor * lp.upper(j1));
3496  }
3497  else if (factor < 0)
3498  {
3499  if (lp.lower(j2) <= -infinity || lp.upper(j1) >= infinity)
3500  lp.changeLower(j2, -infinity);
3501  else
3502  lp.changeLower(j2, lp.lower(j2) + factor * lp.upper(j1));
3503 
3504  if (lp.upper(j2) >= infinity || lp.lower(j1) <= -infinity)
3505  lp.changeUpper(j2, infinity);
3506  else
3507  lp.changeUpper(j2, lp.upper(j2) + factor * lp.lower(j1));
3508  }
3509 
3510  MSG_INFO3( (*spxout), (*spxout) << "IMAISM60 two duplicate columns " << j1
3511  << ", " << j2
3512  << " replaced by one" << std::endl; )
3513 
3514  remCol[j1] = true;
3515 
3517  }
3518  else
3519  {
3520  MSG_INFO3( (*spxout), (*spxout) << "IMAISM80 not removing two duplicate columns " << j1
3521  << ", " << j2
3522  << " because zero not contained in their bounds" << std::endl; )
3523  }
3524  }
3525  else
3526  {
3527  // case 2: objectives not duplicate
3528  // considered for maximization sense
3529  if (lp.lower(j2) <= -infinity)
3530  {
3531  if (factor > 0 && objDif > 0)
3532  {
3533  if (lp.upper(j1) >= infinity)
3534  {
3535  MSG_INFO3( (*spxout), (*spxout) << "IMAISM75 LP unbounded" << std::endl; )
3536  return UNBOUNDED;
3537  }
3538 
3539  // fix j1 at upper bound
3540  MSG_INFO3( (*spxout), (*spxout) << "IMAISM61 two duplicate columns " << j1
3541  << ", " << j2
3542  << " first one fixed at upper bound=" << lp.upper(j1) << std::endl; )
3543 
3544  FixBoundsPS* FixBoundsPSptr = 0;
3545  spx_alloc(FixBoundsPSptr);
3546  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j1, lp.upper(j1)));
3547  lp.changeLower(j1, lp.upper(j1));
3548  }
3549  else if (factor < 0 && objDif < 0)
3550  {
3551  if (lp.lower(j1) <= -infinity)
3552  {
3553  MSG_INFO3( (*spxout), (*spxout) << "IMAISM76 LP unbounded" << std::endl; )
3554  return UNBOUNDED;
3555  }
3556 
3557  // fix j1 at lower bound
3558  MSG_INFO3( (*spxout), (*spxout) << "IMAISM62 two duplicate columns " << j1
3559  << ", " << j2
3560  << " first one fixed at lower bound=" << lp.lower(j1) << std::endl; )
3561 
3562  FixBoundsPS* FixBoundsPSptr = 0;
3563  spx_alloc(FixBoundsPSptr);
3564  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j1, lp.lower(j1)));
3565  lp.changeUpper(j1, lp.lower(j1));
3566  }
3567  }
3568  else if (lp.upper(j2) >= infinity)
3569  {
3570  // fix j1 at upper bound
3571  if (factor < 0 && objDif > 0)
3572  {
3573  if (lp.upper(j1) >= infinity)
3574  {
3575  MSG_INFO3( (*spxout), (*spxout) << "IMAISM77 LP unbounded" << std::endl; )
3576  return UNBOUNDED;
3577  }
3578 
3579  // fix j1 at upper bound
3580  MSG_INFO3( (*spxout), (*spxout) << "IMAISM63 two duplicate columns " << j1
3581  << ", " << j2
3582  << " first one fixed at upper bound=" << lp.upper(j1) << std::endl; )
3583 
3584  FixBoundsPS* FixBoundsPSptr = 0;
3585  spx_alloc(FixBoundsPSptr);
3586  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j1, lp.upper(j1)));
3587  lp.changeLower(j1, lp.upper(j1));
3588  }
3589 
3590  // fix j1 at lower bound
3591  else if (factor > 0 && objDif < 0)
3592  {
3593  if (lp.lower(j1) <= -infinity)
3594  {
3595  MSG_INFO3( (*spxout), (*spxout) << "IMAISM78 LP unbounded" << std::endl; )
3596  return UNBOUNDED;
3597  }
3598 
3599  // fix j1 at lower bound
3600  MSG_INFO3( (*spxout), (*spxout) << "IMAISM64 two duplicate columns " << j1
3601  << ", " << j2
3602  << " first one fixed at lower bound=" << lp.lower(j1) << std::endl; )
3603 
3604  FixBoundsPS* FixBoundsPSptr = 0;
3605  spx_alloc(FixBoundsPSptr);
3606  m_hist.append(new (FixBoundsPSptr) FixBoundsPS(lp, j1, lp.lower(j1)));
3607  lp.changeUpper(j1, lp.lower(j1));
3608  }
3609  }
3610  if (EQrel(lp.lower(j1), lp.upper(j1), feastol()))
3611  {
3612  remCol[j1] = true;
3613  fixAndRemCol[j1] = true;
3614 
3616  }
3617  }
3618  }
3619  }
3620  }
3621  }
3622  }
3623 
3624  for(int j = 0; j < lp.nCols(); ++j)
3625  {
3626  if(fixAndRemCol[j])
3627  {
3628  assert(remCol[j]);
3629 
3630  // correctIdx == false, because the index mapping will be handled by the postsolving in DuplicateColsPS
3631  fixColumn(lp, j, false);
3632  }
3633  }
3634 
3635  // remove all columns by one single method call (more efficient)
3636  const int nColsOld = lp.nCols();
3637  int* perm = 0;
3638  spx_alloc(perm, nColsOld);
3639 
3640  for(int j = 0; j < nColsOld; ++j)
3641  {
3642  if (remCol[j])
3643  {
3644  perm[j] = -1;
3645  ++remCols;
3646  remNzos += lp.colVector(j).size();
3647  }
3648  else
3649  perm[j] = 0;
3650  }
3651  lp.removeCols(perm);
3652 
3653  for(int j = 0; j < nColsOld; ++j)
3654  {
3655  if (perm[j] >= 0)
3656  m_cIdx[perm[j]] = m_cIdx[j];
3657  }
3658 
3659  DataArray<int> da_perm(nColsOld);
3660  for(int j = 0; j < nColsOld; ++j)
3661  {
3662  da_perm[j] = perm[j];
3663  }
3664 
3665  if (hasDuplicateCol)
3666  {
3667  DuplicateColsPS* DuplicateColsPSptr = 0;
3668  spx_alloc(DuplicateColsPSptr);
3669  m_hist.append(new (DuplicateColsPSptr) DuplicateColsPS(lp, 0, 0, 1.0, da_perm, false, true));
3670  }
3671 
3672  spx_free(perm);
3673 
3674  assert(remCols > 0 || remNzos == 0);
3675 
3676  if (remCols > 0)
3677  {
3678  m_remCols += remCols;
3679  m_remNzos += remNzos;
3680 
3681  MSG_INFO2( (*spxout), (*spxout) << "Simplifier (duplicate columns) removed "
3682  << remCols << " cols, "
3683  << remNzos << " non-zeros"
3684  << std::endl; )
3685  if( remCols > m_minReduction * nColsOld )
3686  again = true;
3687  }
3688  return OKAY;
3689 }
3690 
3691 void SPxMainSM::fixColumn(SPxLP& lp, int j, bool correctIdx)
3692 {
3693 
3694  assert(EQrel(lp.lower(j), lp.upper(j), feastol()));
3695 
3696  Real lo = lp.lower(j);
3697  const SVector& col = lp.colVector(j);
3698 
3699  assert(NE(lo, infinity) && NE(lo, -infinity));
3700 
3701  MSG_INFO3( (*spxout), (*spxout) << "IMAISM66 fix variable x" << j
3702  << ": lower=" << lp.lower(j)
3703  << " upper=" << lp.upper(j)
3704  << std::endl; )
3705 
3706  if (isNotZero(lo, epsZero()))
3707  {
3708  for(int k = 0; k < col.size(); ++k)
3709  {
3710  int i = col.index(k);
3711 
3712  if (lp.rhs(i) < infinity)
3713  {
3714  Real y = lo * col.value(k);
3715  Real scale = maxAbs(lp.rhs(i), y);
3716 
3717  if (scale < 1.0)
3718  scale = 1.0;
3719 
3720  Real rhs = (lp.rhs(i) / scale) - (y / scale);
3721 
3722  if (isZero(rhs, epsZero()))
3723  rhs = 0.0;
3724  else
3725  rhs *= scale;
3726 
3727  MSG_INFO3( (*spxout), (*spxout) << "IMAISM67 row " << i
3728  << ": rhs=" << rhs
3729  << " (" << lp.rhs(i)
3730  << ") aij=" << col.value(k)
3731  << std::endl; )
3732 
3733  lp.changeRhs(i, rhs);
3734  }
3735  if (lp.lhs(i) > -infinity)
3736  {
3737  Real y = lo * col.value(k);
3738  Real scale = maxAbs(lp.lhs(i), y);
3739 
3740  if (scale < 1.0)
3741  scale = 1.0;
3742 
3743  Real lhs = (lp.lhs(i) / scale) - (y / scale);
3744 
3745  if (isZero(lhs, epsZero()))
3746  lhs = 0.0;
3747  else
3748  lhs *= scale;
3749 
3750  MSG_INFO3( (*spxout), (*spxout) << "IMAISM68 row " << i
3751  << ": lhs=" << lhs
3752  << " (" << lp.lhs(i)
3753  << ") aij=" << col.value(k)
3754  << std::endl; )
3755 
3756  lp.changeLhs(i, lhs);
3757  }
3758  assert(lp.lhs(i) <= lp.rhs(i));
3759  }
3760  }
3761 
3762  FixVariablePS* FixVariablePSptr = 0;
3763  spx_alloc(FixVariablePSptr);
3764  m_hist.append(new (FixVariablePSptr) FixVariablePS(lp, *this, j, lp.lower(j), correctIdx));
3765 }
3766 
3768 {
3769  // transfer message handler
3770  spxout = lp.spxout;
3771  assert(spxout != 0);
3772 
3773  m_thesense = lp.spxSense();
3774  m_timeUsed->reset();
3775  m_timeUsed->start();
3776 
3777  m_objoffset = 0.0;
3778 
3779  m_remRows = 0;
3780  m_remCols = 0;
3781  m_remNzos = 0;
3782  m_chgBnds = 0;
3783  m_chgLRhs = 0;
3784  m_keptBnds = 0;
3785  m_keptLRhs = 0;
3786 
3787  m_result = OKAY;
3788  bool again = true;
3789  int numRangedRows = 0;
3790  int numBoxedCols = 0;
3791 
3792  if(m_hist.size() > 0)
3793  {
3794  // delete pointers in old m_hist
3795  for(int k = 0; k < m_hist.size(); ++k)
3796  {
3797  m_hist[k]->~PostStep();
3798  spx_free(m_hist[k]);
3799  }
3800  m_hist.clear();
3801  }
3802 
3803  m_hist.reSize(0);
3804  m_postsolved = false;
3805 
3806  if( eps < 0.0 )
3807  throw SPxInterfaceException("XMAISM30 Cannot use negative epsilon in simplify().");
3808 
3809  if( ftol < 0.0 )
3810  throw SPxInterfaceException("XMAISM31 Cannot use negative feastol in simplify().");
3811 
3812  if( otol < 0.0 )
3813  throw SPxInterfaceException("XMAISM32 Cannot use negative opttol in simplify().");
3814 
3815  m_epsilon = eps;
3816  m_feastol = ftol;
3817  m_opttol = otol;
3818 
3819  for(int i = 0; i < lp.nRows(); ++i)
3820  {
3821  if (lp.lhs(i) > -infinity && lp.rhs(i) < infinity)
3822  ++numRangedRows;
3823  }
3824 
3825  for(int j = 0; j < lp.nCols(); ++j)
3826  {
3827  if (lp.lower(j) > -infinity && lp.upper(j) < infinity)
3828  ++numBoxedCols;
3829  }
3830 
3831  MSG_INFO2( (*spxout), (*spxout) << "LP has "
3832  << numRangedRows << " ranged rows, "
3833  << numBoxedCols << " boxed columns"
3834  << std::endl; )
3835 
3836  m_stat.reSize(15);
3837 
3838  for(int k = 0; k < m_stat.size(); ++k)
3839  m_stat[k] = 0;
3840 
3841  m_addedcols = 0;
3842  handleRowObjectives(lp);
3843 
3844  m_prim.reDim(lp.nCols());
3845  m_slack.reDim(lp.nRows());
3846  m_dual.reDim(lp.nRows());
3847  m_redCost.reDim(lp.nCols());
3848  m_cBasisStat.reSize(lp.nCols());
3849  m_rBasisStat.reSize(lp.nRows());
3850  m_cIdx.reSize(lp.nCols());
3851  m_rIdx.reSize(lp.nRows());
3852 
3853  m_classSetRows.reSize(lp.nRows());
3854  m_classSetCols.reSize(lp.nCols());
3855  m_dupRows.reSize(lp.nRows());
3856  m_dupCols.reSize(lp.nCols());
3857 
3858  m_keepbounds = keepbounds;
3859 
3860  for(int i = 0; i < lp.nRows(); ++i)
3861  m_rIdx[i] = i;
3862 
3863  for(int j = 0; j < lp.nCols(); ++j)
3864  m_cIdx[j] = j;
3865 
3866  // round extreme values (set all values smaller than eps to zero and all values bigger than infinity/5 to infinity)
3867 #if EXTREMES
3868  handleExtremes(lp);
3869 #endif
3870 
3871  // main presolving loop
3872  while(again && m_result == OKAY)
3873  {
3874  again = false;
3875 
3876 #if ROWS
3877  if (m_result == OKAY)
3878  m_result = simplifyRows(lp, again);
3879 #endif
3880 
3881 #if COLS
3882  if (m_result == OKAY)
3883  m_result = simplifyCols(lp, again);
3884 #endif
3885 
3886 #if DUAL
3887  if (m_result == OKAY)
3888  m_result = simplifyDual(lp, again);
3889 #endif
3890 
3891 #if DUPLICATE_ROWS
3892  if (m_result == OKAY)
3893  m_result = duplicateRows(lp, again);
3894 #endif
3895 
3896 #if DUPLICATE_COLS
3897  if (m_result == OKAY)
3898  m_result = duplicateCols(lp, again);
3899 #endif
3900  }
3901 
3902  // preprocessing detected infeasibility or unboundedness
3903  if (m_result != OKAY)
3904  return m_result;
3905 
3908  MSG_INFO1( (*spxout), (*spxout) << "Simplifier removed "
3909  << m_remRows << " rows, "
3910  << m_remCols << " columns, "
3911  << m_remNzos << " nonzeros, "
3912  << m_chgBnds << " col bounds, "
3913  << m_chgLRhs << " row bounds"
3914  << std::endl; )
3915 
3916  if (keepbounds)
3917  MSG_INFO2( (*spxout), (*spxout) << "Simplifier kept "
3918  << m_keptBnds << " column bounds, "
3919  << m_keptLRhs << " row bounds"
3920  << std::endl; )
3921 
3922  MSG_INFO1( (*spxout), (*spxout) << "Reduced LP has "
3923  << lp.nRows() << " rows "
3924  << lp.nCols() << " columns "
3925  << lp.nNzos() << " nonzeros"
3926  << std::endl; )
3927 
3928  // reset counter
3929  numRangedRows = 0;
3930  numBoxedCols = 0;
3931 
3932  for(int i = 0; i < lp.nRows(); ++i)
3933  if (lp.lhs(i) > -infinity && lp.rhs(i) < infinity)
3934  ++numRangedRows;
3935 
3936  for(int j = 0; j < lp.nCols(); ++j)
3937  if (lp.lower(j) > -infinity && lp.upper(j) < infinity)
3938  ++numBoxedCols;
3939 
3940  MSG_INFO2( (*spxout), (*spxout) << "Reduced LP has "
3941  << numRangedRows << " ranged rows, "
3942  << numBoxedCols << " boxed columns"
3943  << std::endl; )
3944 
3945  if (lp.nCols() == 0 && lp.nRows() == 0)
3946  {
3947  MSG_INFO1( (*spxout), (*spxout) << "Simplifier removed all rows and columns" << std::endl; )
3948  m_result = VANISHED;
3949  }
3950 
3951  MSG_INFO2( (*spxout), (*spxout) << "\nSimplifier performed:\n"
3952  << m_stat[EMPTY_ROW] << " empty rows\n"
3953  << m_stat[FREE_ROW] << " free rows\n"
3954  << m_stat[SINGLETON_ROW] << " singleton rows\n"
3955  << m_stat[FORCE_ROW] << " forcing rows\n"
3956  << m_stat[EMPTY_COL] << " empty columns\n"
3957  << m_stat[FIX_COL] << " fixed columns\n"
3958  << m_stat[FREE_ZOBJ_COL] << " free columns with zero objective\n"
3959  << m_stat[ZOBJ_SINGLETON_COL] << " singleton columns with zero objective\n"
3960  << m_stat[DOUBLETON_ROW] << " singleton columns combined with a doubleton equation\n"
3961  << m_stat[FREE_SINGLETON_COL] << " free singleton columns\n"
3962  << m_stat[DOMINATED_COL] << " dominated columns\n"
3963  << m_stat[WEAKLY_DOMINATED_COL] << " weakly dominated columns\n"
3964  << m_stat[DUPLICATE_ROW] << " duplicate rows\n"
3965  << m_stat[FIX_DUPLICATE_COL] << " duplicate columns (fixed)\n"
3966  << m_stat[SUB_DUPLICATE_COL] << " duplicate columns (substituted)\n"
3967  << std::endl; );
3968 
3969  m_timeUsed->stop();
3970 
3971  return m_result;
3972 }
3973 
3974 void SPxMainSM::unsimplify(const Vector& x, const Vector& y, const Vector& s, const Vector& r,
3975  const SPxSolver::VarStatus rows[], const SPxSolver::VarStatus cols[])
3976 {
3977  assert(x.dim() <= m_prim.dim());
3978  assert(y.dim() <= m_dual.dim());
3979  assert(x.dim() == r.dim());
3980  assert(y.dim() == s.dim());
3981 
3982 
3983  // assign values of variables in reduced LP
3984  // NOTE: for maximization problems, we have to switch signs of dual and reduced cost values,
3985  // since simplifier assumes minimization problem
3986  for(int j = 0; j < x.dim(); ++j)
3987  {
3988  m_prim[j] = isZero(x[j], epsZero()) ? 0.0 : x[j];
3989  m_redCost[j] = isZero(r[j], epsZero()) ? 0.0 : (m_thesense == SPxLP::MAXIMIZE ? -r[j] : r[j]);
3990  m_cBasisStat[j] = cols[j];
3991  }
3992  for(int i = 0; i < y.dim(); ++i)
3993  {
3994  m_dual[i] = isZero(y[i], epsZero()) ? 0.0 : (m_thesense == SPxLP::MAXIMIZE ? -y[i] : y[i]);
3995  m_slack[i] = isZero(s[i], epsZero()) ? 0.0 : s[i];
3996  m_rBasisStat[i] = rows[i];
3997  }
3998 
3999  // undo preprocessing
4000  for(int k = m_hist.size()-1; k >= 0; --k)
4001  {
4002  MSG_DEBUG( std::cout << "unsimplifying " << m_hist[k]->getName() << "\n" );
4003 
4004  try
4005  {
4007  }
4008  catch( ... )
4009  {
4010  MSG_INFO1( (*spxout), (*spxout) << "Exception thrown while unsimplifying " << m_hist[k]->getName() << ".\n" );
4011  throw SPxInternalCodeException("XMAISM00 Exception thrown during unsimply().");
4012  }
4013 
4014  m_hist[k]->~PostStep();
4015  spx_free(m_hist[k]);
4016  m_hist.reSize(k);
4017  }
4018 
4019  // for maximization problems, we have to switch signs of dual and reduced cost values back
4021  {
4022  for(int j = 0; j < m_redCost.dim(); ++j)
4023  m_redCost[j] = -m_redCost[j];
4024 
4025  for(int i = 0; i < m_dual.dim(); ++i)
4026  m_dual[i] = -m_dual[i];
4027  }
4028 
4029  if( m_addedcols > 0 )
4030  {
4031  assert(m_prim.dim() >= m_addedcols);
4034  m_cBasisStat.reSize(m_cBasisStat.size() - m_addedcols);
4036  }
4037 
4038 #ifdef CHECK_BASIC_DIM
4039  int numBasis = 0;
4040  for(int rs = 0; rs < m_rBasisStat.size(); ++rs)
4041  {
4042  if(m_rBasisStat[rs] == SPxSolver::BASIC)
4043  numBasis ++;
4044  }
4045 
4046  for(int cs = 0; cs < m_cBasisStat.size(); ++cs)
4047  {
4048  if(m_cBasisStat[cs] == SPxSolver::BASIC)
4049  numBasis ++;
4050  }
4051 
4052  if( numBasis != m_rBasisStat.size())
4053  {
4054  throw SPxInternalCodeException("XMAISM26 Dimension doesn't match after this step.");
4055  }
4056 #endif
4057 
4058  m_hist.clear();
4059  m_postsolved = true;
4060 }
4061 } //namespace soplex
DVector m_prim
unsimplified primal solution vector.
Definition: spxmainsm.h:1076
Rational spxAbs(const Rational &r)
Absolute.
Definition: rational.cpp:3925
bool isNotZero(Real a, Real eps=Param::epsilon())
returns true iff |a| > eps
Definition: spxdefines.h:413
Real m_opttol
dual feasibility tolerance.
Definition: spxmainsm.h:1092
R obj(int i) const
Returns objective value of column i.
Definition: spxlpbase.h:362
DataArray< PostStep * > m_hist
vector of presolve history.
Definition: spxmainsm.h:1084
Exception class for things that should NEVER happen.This class is derived from the SoPlex exception b...
Definition: exceptions.h:109
int cIdx(int j) const
returns for a given column index of the (reduced) LP the corresponding column index in the unsimplifi...
Definition: spxmainsm.h:1152
bool LTrel(Real a, Real b, Real eps=Param::epsilon())
returns true iff relDiff(a,b) <= -eps
Definition: spxdefines.h:431
free variable fixed to zero.
Definition: spxsolver.h:182
Safe arrays of data objects.Class DataArray provides safe arrays of Data Objects. For general C++ obj...
Definition: dataarray.h:63
friend class FreeConstraintPS
Definition: spxmainsm.h:1034
Result removeRowSingleton(SPxLP &lp, const SVector &row, int &i)
remove row singletons.
Definition: spxmainsm.cpp:1571
DataArray< int > m_stat
preprocessing history.
Definition: spxmainsm.h:1093
void reDim(int newdim, const bool setZero=true)
Resets DVectorBase&#39;s dimension to newdim.
Definition: dvectorbase.h:249
int m_keptLRhs
number of kept left- and right-hand sides
Definition: spxsimplifier.h:66
R rowObj(int i) const
Definition: spxlpbase.h:268
UnitVectorBase< Real > UnitVector
Definition: unitvector.h:29
const VectorBase< R > & maxObj() const
Returns objective vector for maximization problem.
Definition: spxlpbase.h:384
virtual void changeRhs(const VectorBase< R > &newRhs)
Changes right hand side vector for constraints to newRhs.
Definition: spxlpbase.h:1371
Real m_epsilon
epsilon zero.
Definition: spxmainsm.h:1090
Postsolves variable bound fixing.
Definition: spxmainsm.h:505
Result
Result of the simplification.
Definition: spxsimplifier.h:81
the problem was so much simplified that it vanished
Definition: spxsimplifier.h:87
virtual void changeObj(const VectorBase< R > &newObj)
Changes objective vector to newObj.
Definition: spxlpbase.h:1182
virtual void unsimplify(const Vector &x, const Vector &y, const Vector &s, const Vector &r, const SPxSolver::VarStatus rows[], const SPxSolver::VarStatus cols[])
reconstructs an optimal solution for the unsimplified LP.
Definition: spxmainsm.cpp:3974
Postsolves duplicate rows.
Definition: spxmainsm.h:875
int m_remNzos
number of removed nonzero coefficients
Definition: spxsimplifier.h:58
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:665
Real m_feastol
primal feasibility tolerance.
Definition: spxmainsm.h:1091
DataArray< SPxSolver::VarStatus > m_rBasisStat
basis status of rows.
Definition: spxmainsm.h:1081
Real maxAbs(Real a, Real b)
returns max(|a|,|b|)
Definition: spxdefines.h:356
virtual void removeCols(int perm[])
Removes multiple columns.
Definition: spxlpbase.h:929
bool LT(Real a, Real b, Real eps=Param::epsilon())
returns true iff a < b + eps
Definition: spxdefines.h:383
SPxLP::SPxSense m_thesense
optimization sense.
Definition: spxmainsm.h:1094
const SVectorBase< R > & colVector(int i) const
Returns column vector of column i.
Definition: spxlpbase.h:341
bool NE(Real a, Real b, Real eps=Param::epsilon())
returns true iff |a-b| > eps
Definition: spxdefines.h:377
virtual void changeLhs(const VectorBase< R > &newLhs)
Changes left hand side vector for constraints to newLhs.
Definition: spxlpbase.h:1342
virtual Result simplify(SPxLP &lp, Real eps, Real delta)
simplify SPxLP lp with identical primal and dual feasibility tolerance.
Definition: spxmainsm.h:1290
int size() const
Number of used indices.
Definition: svectorbase.h:152
friend class DuplicateRowsPS
Definition: spxmainsm.h:1044
virtual const char * getName() const
get name of simplifying step.
Definition: spxmainsm.h:107
const VectorBase< R > & lower() const
Returns lower bound vector.
Definition: spxlpbase.h:420
Real feastol() const
Definition: spxmainsm.h:1162
variable fixed to identical bounds.
Definition: spxsolver.h:181
DVector m_dual
unsimplified dual solution vector.
Definition: spxmainsm.h:1078
R & value(int n)
Reference to value of n &#39;th nonzero.
Definition: svectorbase.h:254
General methods in LP preprocessing.
Postsolves empty constraints.
Definition: spxmainsm.h:222
int m_keptBnds
number of kept bounds
Definition: spxsimplifier.h:64
#define ASSERT_WARN(prefix, expr)
Macro to turn some assertions into warnings.
Definition: spxdefines.h:72
Real m_minReduction
minimal reduction (sum of removed rows/cols) to continue simplification
Definition: spxsimplifier.h:70
void remove(int n, int m)
Remove nonzeros n thru m.
Definition: svectorbase.h:361
void handleExtremes(SPxLP &lp)
handles extreme values by setting them to zero or infinity.
Definition: spxmainsm.cpp:1268
bool m_postsolved
status of postsolving.
Definition: spxmainsm.h:1089
Real epsZero() const
Definition: spxmainsm.h:1157
const VectorBase< R > & upper() const
Returns upper bound vector.
Definition: spxlpbase.h:402
virtual bool checkBasisDim(DataArray< SPxSolver::VarStatus > rows, DataArray< SPxSolver::VarStatus > cols) const
Definition: spxmainsm.cpp:62
SVectorBase< R > & rowVector_w(int i)
Definition: spxlpbase.h:2036
virtual void removeRows(int perm[])
Removes multiple rows.
Definition: spxlpbase.h:832
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:178
void fixColumn(SPxLP &lp, int i, bool correctIdx=true)
handles the fixing of a variable. correctIdx is true iff the index mapping has to be updated...
Definition: spxmainsm.cpp:3691
virtual void changeUpper(const VectorBase< R > &newUpper)
Changes vector of upper bounds to newUpper.
Definition: spxlpbase.h:1278
Wrapper for different output streams and verbosity levels.
virtual void changeRowObj(const VectorBase< R > &newRowObj)
Changes row objective function vector to newRowObj.
Definition: spxlpbase.h:1425
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:354
Exception class for out of memory exceptions.This class is derived from the SoPlex exception base cla...
Definition: exceptions.h:70
virtual void start()=0
start timer, resume accounting user, system and real time.
virtual Real stop()=0
stop timer, return accounted user time.
void spx_alloc(T &p, int n=1)
Allocate memory.
Definition: spxalloc.h:48
simplification could be done
Definition: spxsimplifier.h:83
virtual void addCol(const LPColBase< R > &col)
Definition: spxlpbase.h:650
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:804
Postsolves the case when constraints are removed due to a variable with zero objective that is free i...
Definition: spxmainsm.h:563
virtual void changeLower(const VectorBase< R > &newLower)
Changes vector of lower bounds to newLower.
Definition: spxlpbase.h:1247
int rIdx(int i) const
returns for a given row index of the (reduced) LP the corresponding row index in the unsimplified LP...
Definition: spxmainsm.h:1147
friend class RowSingletonPS
Definition: spxmainsm.h:1036
Timer * m_timeUsed
user time used for simplification
Definition: spxsimplifier.h:51
virtual void changeRange(const VectorBase< R > &newLhs, const VectorBase< R > &newRhs)
Changes left and right hand side vectors.
Definition: spxlpbase.h:1392
nothing known about basis status (possibly due to a singular basis in transformed problem) ...
Definition: spxsolver.h:184
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:455
double Real
SOPLEX_DEBUG.
Definition: spxdefines.h:200
Array< DSVector > m_classSetCols
stores parallel classes with non-zero row entry
Definition: spxmainsm.h:1086
#define MSG_INFO2(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO2.
Definition: spxdefines.h:115
Result duplicateCols(SPxLP &lp, bool &again)
removes duplicate columns
Definition: spxmainsm.cpp:3315
void removeCol(SPxLP &lp, int j)
removes a column in the LP.
Definition: spxmainsm.h:1141
#define MSG_DEBUG(x)
Definition: spxdefines.h:127
int nRows
number of rows
Definition: spxmainsm.h:81
int number(int i) const
Number of index i.
Definition: svectorbase.h:185
bool GT(Real a, Real b, Real eps=Param::epsilon())
returns true iff a > b + eps
Definition: spxdefines.h:395
int & index(int n)
Reference to index of n &#39;th nonzero.
Definition: svectorbase.h:236
Postsolves row singletons.
Definition: spxmainsm.h:269
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:81
Result simplifyRows(SPxLP &lp, bool &again)
performs simplification steps on the rows of the LP.
Definition: spxmainsm.cpp:1644
int nCols() const
Returns number of columns in LP.
Definition: spxlpbase.h:133
friend class DoubletonEquationPS
Definition: spxmainsm.h:1043
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:909
Postsolves doubleton equations combined with a column singleton.
Definition: spxmainsm.h:779
Postsolves free column singletons.
Definition: spxmainsm.h:710
int size() const
return nr. of elements.
Definition: dataarray.h:211
comparator for class SVector::Element: compare nonzeros according to value
Definition: spxmainsm.h:1368
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:506
friend class FixBoundsPS
Definition: spxmainsm.h:1039
DataArray< int > m_cIdx
column index vector in original LP.
Definition: spxmainsm.h:1082
Result m_result
result of the simplification.
Definition: spxmainsm.h:1097
variable set to its upper bound.
Definition: spxsolver.h:179
int nNzos() const
Returns number of nonzeros in LP.
Definition: spxlpbase.h:139
bool GTrel(Real a, Real b, Real eps=Param::epsilon())
returns true iff relDiff(a,b) > eps
Definition: spxdefines.h:443
void SPxQuicksort(T *keys, int end, COMPARATOR &compare, int start=0, bool type=true)
Generic QuickSort implementation.
Definition: sorter.h:73
#define MSG_INFO1(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO1.
Definition: spxdefines.h:113
int dim() const
Dimension of vector.
Definition: vectorbase.h:174
primal infeasibility was detected
Definition: spxsimplifier.h:84
Real opttol() const
Definition: spxmainsm.h:1167
bool GErel(Real a, Real b, Real eps=Param::epsilon())
returns true iff relDiff(a,b) > -eps
Definition: spxdefines.h:449
virtual void changeBounds(const VectorBase< R > &newLower, const VectorBase< R > &newUpper)
Changes variable bounds to newLower and newUpper.
Definition: spxlpbase.h:1309
variable set to its lower bound.
Definition: spxsolver.h:180
Generic QuickSort implementation.
Array< DSVector > m_dupCols
arrange duplicate columns w.r.t. their pClass values
Definition: spxmainsm.h:1088
int m_remRows
number of removed rows
Definition: spxsimplifier.h:54
void add(int i, const R &v)
Append one nonzero (i,v).
Definition: svectorbase.h:272
int nRows() const
Returns number of rows in LP.
Definition: spxlpbase.h:127
Postsolves row objectives.
Definition: spxmainsm.h:128
DataArray< int > m_rIdx
row index vector in original LP.
Definition: spxmainsm.h:1083
Postsolves column singletons with zero objective.
Definition: spxmainsm.h:647
Exception base class.This class implements a base class for our SoPlex exceptions We provide a what()...
Definition: exceptions.h:32
bool m_keepbounds
keep some bounds (for boundflipping)
Definition: spxmainsm.h:1095
const VectorBase< R > & lhs() const
Returns left hand side vector.
Definition: spxlpbase.h:242
Everything should be within this namespace.
bool isConsistent() const
Consistency check.
Definition: spxlpbase.h:1753
bool EQrel(Real a, Real b, Real eps=Param::epsilon())
returns true iff |relDiff(a,b)| <= eps
Definition: spxdefines.h:419
friend class ForceConstraintPS
Definition: spxmainsm.h:1037
void handleRowObjectives(SPxLP &lp)
handle row objectives
Definition: spxmainsm.cpp:1251
primal unboundedness was detected
Definition: spxsimplifier.h:86
#define MSG_WARNING(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::WARNING.
Definition: spxdefines.h:111
Result duplicateRows(SPxLP &lp, bool &again)
removes duplicate rows.
Definition: spxmainsm.cpp:2961
Exception class for incorrect usage of interface methods.
Definition: exceptions.h:126
Save arrays of arbitrary types.
Real m_objoffset
objective offset
Definition: spxsimplifier.h:68
friend class DuplicateColsPS
Definition: spxmainsm.h:1045
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:1018
variable is basic.
Definition: spxsolver.h:183
bool NErel(Real a, Real b, Real eps=Param::epsilon())
returns true iff |relDiff(a,b)| > eps
Definition: spxdefines.h:425
comparator for class SVector::Element: compare nonzeros according to index
Definition: spxmainsm.h:1384
int m_chgBnds
number of changed bounds
Definition: spxsimplifier.h:60
Postsolves duplicate columns.
Definition: spxmainsm.h:970
SPxSense spxSense() const
Returns the optimization sense.
Definition: spxlpbase.h:438
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:121
Result removeEmpty(SPxLP &lp)
removed empty rows and empty columns.
Definition: spxmainsm.cpp:1463
Postsolves forcing constraints.
Definition: spxmainsm.h:352
Result simplifyDual(SPxLP &lp, bool &again)
performs simplification steps on the LP based on dual concepts.
Definition: spxmainsm.cpp:2707
Postsolves unconstraint constraints.
Definition: spxmainsm.h:171
const Real infinity
Definition: spxdefines.cpp:26
SPxOut * spxout
Definition: spxlpbase.h:114
int m_remCols
number of removed columns
Definition: spxsimplifier.h:56
const VectorBase< R > & rhs() const
Returns right hand side vector.
Definition: spxlpbase.h:224
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:864
Postsolves variable fixing.
Definition: spxmainsm.h:442
friend class FreeColSingletonPS
Definition: spxmainsm.h:1042
bool LErel(Real a, Real b, Real eps=Param::epsilon())
returns true iff relDiff(a,b) <= eps
Definition: spxdefines.h:437
const SVectorBase< R > & rowVector(int i) const
Gets row vector of row i.
Definition: spxlpbase.h:212
friend class FreeZeroObjVariablePS
Definition: spxmainsm.h:1040
int m_addedcols
columns added by handleRowObjectives()
Definition: spxmainsm.h:1096
friend class ZeroObjColSingletonPS
Definition: spxmainsm.h:1041
bool isZero(Real a, Real eps=Param::epsilon())
returns true iff |a| <= eps
Definition: spxdefines.h:407
void removeRow(SPxLP &lp, int i)
removes a row in the LP.
Definition: spxmainsm.h:1135
friend class EmptyConstraintPS
Definition: spxmainsm.h:1035
virtual void reset()=0
initialize timer, set timing accounts to zero.
#define MSG_INFO3(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO3.
Definition: spxdefines.h:117
int nCols
number of cols
Definition: spxmainsm.h:79
SPxOut * spxout
message handler
Definition: spxsimplifier.h:72
Save arrays of data objects.
SVectorBase< R > & colVector_w(int i)
Returns the LP as an LPRowSet.
Definition: spxlpbase.h:2030
int m_chgLRhs
number of change right-hand sides
Definition: spxsimplifier.h:62
Result simplifyCols(SPxLP &lp, bool &again)
performs simplification steps on the columns of the LP.
Definition: spxmainsm.cpp:2180
Array< DSVector > m_dupRows
arrange duplicate rows using bucket sort w.r.t. their pClass values
Definition: spxmainsm.h:1087
static Real epsilon()
Definition: spxdefines.h:255
void reSize(int newsize)
reset size to newsize.
Definition: dataarray.h:223
DVector m_slack
unsimplified slack vector.
Definition: spxmainsm.h:1077
Set of indices.Class IdxSet provides a set of indices. At construction it must be given an array of i...
Definition: idxset.h:56
void spx_free(T &p)
Release memory.
Definition: spxalloc.h:109
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:152
const VectorBase< R > & maxRowObj() const
Definition: spxlpbase.h:286
DataArray< SPxSolver::VarStatus > m_cBasisStat
basis status of columns.
Definition: spxmainsm.h:1080
virtual void execute(DVector &x, DVector &y, DVector &s, DVector &r, DataArray< SPxSolver::VarStatus > &cBasis, DataArray< SPxSolver::VarStatus > &rBasis) const
Definition: spxmainsm.cpp:514
dual infeasibility was detected
Definition: spxsimplifier.h:85
Array< DSVector > m_classSetRows
stores parallel classes with non-zero colum entry
Definition: spxmainsm.h:1085
friend class FixVariablePS
Definition: spxmainsm.h:1038
DVector m_redCost
unsimplified reduced cost vector.
Definition: spxmainsm.h:1079