1 /* */
2 /* This file is part of the class library */
3 /* SoPlex --- the Sequential object-oriented simPlex. */
4 /* */
5 /* Copyright (C) 1996-2022 Konrad-Zuse-Zentrum */
6 /* fuer Informationstechnik Berlin */
7 /* */
8 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
9 /* */
10 /* You should have received a copy of the ZIB Academic License */
11 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
12 /* */
13 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
14
15 #include <iostream>
16 #include <assert.h>
17
18 #include "soplex/spxdefines.h"
19 #include "soplex.h"
20 #include "soplex/statistics.h"
21 #include "soplex/sorter.h"
22
23 //#define NO_TOL
24 #define USE_FEASTOL
25 //#define NO_TRANSFORM
26 //#define PERFORM_COMPPROB_CHECK
27
28 #define MAX_DEGENCHECK 20 /**< the maximum number of degen checks that are performed before the DECOMP is abandoned */
29 #define DEGENCHECK_OFFSET 50 /**< the number of iteration before the degeneracy check is reperformed */
30 #define SLACKCOEFF 1.0 /**< the coefficient of the slack variable in the incompatible rows. */
31 #define TIMELIMIT_FRAC 0.5 /**< the fraction of the total time limit given to the setup of the reduced problem */
32
33 /* This file contains the private functions for the Decomposition Based Dual Simplex (DBDS)
34 *
35 * An important note about the implementation of the DBDS is the reliance on the row representation of the basis matrix.
36 * The forming of the reduced and complementary problems is involves identifying rows in the row-form of the basis
37 * matrix that have zero dual multipliers.
38 *
39 * Ideally, this work will be extended such that the DBDS can be executed using the column-form of the basis matrix. */
40
41 namespace soplex
42 {
43
44 /// solves LP using the decomposition dual simplex
45 template <class R>
46 void SoPlexBase<R>::_solveDecompositionDualSimplex()
47 {
48 assert(_solver.rep() == SPxSolverBase<R>::ROW);
49 assert(_solver.type() == SPxSolverBase<R>::LEAVE);
50
51 // flag to indicate that the algorithm must terminate
52 bool stop = false;
53
54 // setting the initial status of the reduced problem
55 _statistics->redProbStatus = SPxSolverBase<R>::NO_PROBLEM;
56
57 // start timing
58 _statistics->solvingTime->start();
59
60 // remember that last solve was in floating-point
61 _lastSolveMode = SOLVEMODE_REAL;
62
63 // setting the current solving mode.
64 _currentProb = DECOMP_ORIG;
65
66 // setting the sense to maximise. This is to make all matrix operations more consistent.
67 // @todo if the objective sense is changed, the output of the current objective value is the negative of the
68 // actual value. This needs to be corrected in future versions.
69 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE)
70 {
71 assert(_solver.spxSense() == SPxLPBase<R>::MINIMIZE);
72
73 _solver.changeObj(_solver.maxObj());
74 _solver.changeSense(SPxLPBase<R>::MAXIMIZE);
75 }
76
77 // it is necessary to solve the initial problem to find a starting basis
78 _solver.setDecompStatus(SPxSolverBase<R>::FINDSTARTBASIS);
79
80 // setting the decomposition iteration limit to the parameter setting
81 _solver.setDecompIterationLimit(intParam(SoPlexBase<R>::DECOMP_ITERLIMIT));
82
83 // variables used in the initialisation phase of the decomposition solve.
84 int numDegenCheck = 0;
85 R degeneracyLevel = 0;
86 _decompFeasVector.reDim(_solver.nCols());
87 bool initSolveFromScratch = true;
88
89
90 /************************/
91 /* Initialisation phase */
92 /************************/
93
94 // arrays to store the basis status for the rows and columns at the interruption of the original problem solve.
95 DataArray< typename SPxSolverBase<R>::VarStatus > basisStatusRows;
96 DataArray< typename SPxSolverBase<R>::VarStatus > basisStatusCols;
97
98 // since the original LP may have been shifted, the dual multiplier will not be correct for the original LP. This
99 // loop will recheck the degeneracy level and compute the proper dual multipliers.
100 do
101 {
102 // solving the instance.
103 _decompSimplifyAndSolve(_solver, _slufactor, initSolveFromScratch, initSolveFromScratch);
104 initSolveFromScratch = false;
105
106 // checking whether the initialisation must terminate and the original problem is solved using the dual simplex.
107 if(_solver.type() == SPxSolverBase<R>::LEAVE || _solver.status() >= SPxSolverBase<R>::OPTIMAL
108 || _solver.status() == SPxSolverBase<R>::ABORT_EXDECOMP || numDegenCheck > MAX_DEGENCHECK)
109 {
110 // decomposition is deemed not useful. Solving the original problem using regular SoPlexBase.
111
112 // returning the sense to minimise
113 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE)
114 {
115 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE);
116
117 _solver.changeObj(-(_solver.maxObj()));
118 _solver.changeSense(SPxLPBase<R>::MINIMIZE);
119 }
120
121 // switching off the starting basis check. This is only required to initialise the decomposition simplex.
122 _solver.setDecompStatus(SPxSolverBase<R>::DONTFINDSTARTBASIS);
123
124 // the basis is not computed correctly is the problem was unfeasible or unbounded.
125 if(_solver.status() == SPxSolverBase<R>::UNBOUNDED
126 || _solver.status() == SPxSolverBase<R>::INFEASIBLE)
127 _hasBasis = false;
128
129
130 // resolving the problem to update the real lp and solve with the correct objective.
131 // TODO: With some infeasible problem (e.g. refinery) the dual is violated. Setting fromScratch to true
132 // avoids this issue. Need to check what the problem is.
133 // @TODO: Should this be _preprocessAndSolveReal similar to the resolve at the end of the algorithm?
134 _decompSimplifyAndSolve(_solver, _slufactor, true, false);
135
136 // retreiving the original problem statistics prior to destroying it.
137 getOriginalProblemStatistics();
138
139 // storing the solution from the reduced problem
140 _storeSolutionReal();
141
142 // stop timing
143 _statistics->solvingTime->stop();
144 return;
145 }
146 else if(_solver.status() == SPxSolverBase<R>::ABORT_TIME
147 || _solver.status() == SPxSolverBase<R>::ABORT_ITER
148 || _solver.status() == SPxSolverBase<R>::ABORT_VALUE)
149 {
150 // This cleans up the problem is an abort is reached.
151
152 // at this point, the _decompSimplifyAndSolve does not store the realLP. It stores the _decompLP. As such, it is
153 // important to reinstall the _realLP to the _solver.
154 if(!_isRealLPLoaded)
155 {
156 _solver.loadLP(*_decompLP);
157 spx_free(_decompLP);
158 _decompLP = &_solver;
159 _isRealLPLoaded = true;
160 }
161
162 // returning the sense to minimise
163 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE)
164 {
165 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE);
166
167 _solver.changeObj(-(_solver.maxObj()));
168 _solver.changeSense(SPxLPBase<R>::MINIMIZE);
169 }
170
171 // resolving the problem to set up all of the solution vectors. This is required because data from the
172 // initialisation solve may remain at termination causing an infeasible solution to be reported.
173 _preprocessAndSolveReal(false);
174
175 // storing the solution from the reduced problem
176 _storeSolutionReal();
177
178 // stop timing
179 _statistics->solvingTime->stop();
180 return;
181 }
182
183 // checking the degeneracy level
184 _solver.basis().solve(_decompFeasVector, _solver.maxObj());
185 degeneracyLevel = _solver.getDegeneracyLevel(_decompFeasVector);
186
187 _solver.setDegenCompOffset(DEGENCHECK_OFFSET);
188
189 numDegenCheck++;
190 }
191 while((degeneracyLevel > 0.9 || degeneracyLevel < 0.1)
192 || !checkBasisDualFeasibility(_decompFeasVector));
193
194 // decomposition simplex will commence, so the start basis does not need to be found
195 _solver.setDecompStatus(SPxSolverBase<R>::DONTFINDSTARTBASIS);
196
197 // if the original problem has a basis, this will be stored
198 if(_hasBasis)
199 {
200 basisStatusRows.reSize(numRows());
201 basisStatusCols.reSize(numCols());
202 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(),
203 basisStatusCols.size());
204 }
205
206 // updating the algorithm iterations statistic
207 _statistics->callsReducedProb++;
208
209 // setting the statistic information
210 numDecompIter = 0;
211 numRedProbIter = _solver.iterations();
212 numCompProbIter = 0;
213
214 // setting the display information
215 _decompDisplayLine = 0;
216
217 /************************/
218 /* Decomposition phase */
219 /************************/
220
221 MSG_INFO1(spxout,
222 spxout << "======== Degeneracy Detected ========" << std::endl
223 << std::endl
224 << "======== Commencing decomposition solve ========" << std::endl
225 );
226
227 //spxout.setVerbosity( SPxOut::DEBUG );
228 MSG_INFO2(spxout, spxout << "Creating the Reduced and Complementary problems." << std::endl);
229
230 // setting the verbosity level
231 const SPxOut::Verbosity orig_verbosity = spxout.getVerbosity();
232 const SPxOut::Verbosity decomp_verbosity = (SPxOut::Verbosity)intParam(
233 SoPlexBase<R>::DECOMP_VERBOSITY);
234
235 if(decomp_verbosity < orig_verbosity)
236 spxout.setVerbosity(decomp_verbosity);
237
238 // creating copies of the original problem that will be manipulated to form the reduced and complementary
239 // problems.
240 _createDecompReducedAndComplementaryProblems();
241
242 // creating the initial reduced problem from the basis information
243 _formDecompReducedProblem(stop);
244
245 // setting flags for the decomposition solve
246 _hasBasis = false;
247 bool hasRedBasis = false;
248 bool redProbError = false;
249 bool noRedprobIter = false;
250 bool explicitviol = boolParam(SoPlexBase<R>::EXPLICITVIOL);
251 int algIterCount = 0;
252
253
254 // a stop will be triggered if the reduced problem exceeded a specified fraction of the time limit
255 if(stop)
256 {
257 MSG_INFO1(spxout, spxout << "==== Error constructing the reduced problem ====" << std::endl);
258 redProbError = true;
259 _statistics->redProbStatus = SPxSolverBase<R>::NOT_INIT;
260 }
261
262 // the main solving loop of the decomposition simplex.
263 // This loop solves the Reduced problem, and if the problem is feasible, the complementary problem is solved.
264 while(!stop)
265 {
266 int previter = _statistics->iterations;
267
268 // setting the current solving mode.
269 _currentProb = DECOMP_RED;
270
271 // solve the reduced problem
272
273 MSG_INFO2(spxout,
274 spxout << std::endl
275 << "=========================" << std::endl
276 << "Solving: Reduced Problem." << std::endl
277 << "=========================" << std::endl
278 << std::endl);
279
280 _hasBasis = hasRedBasis;
281 // solving the reduced problem
282 _decompSimplifyAndSolve(_solver, _slufactor, !algIterCount, !algIterCount);
283
284 stop = decompTerminate(realParam(
285 SoPlexBase<R>::TIMELIMIT)); // checking whether the algorithm should terminate
286
287 // updating the algorithm iterations statistics
288 _statistics->callsReducedProb++;
289 _statistics->redProbStatus = _solver.status();
290
291 assert(_isRealLPLoaded);
292 hasRedBasis = _hasBasis;
293 _hasBasis = false;
294
295 // checking the status of the reduced problem
296 // It is expected that infeasibility and unboundedness will be found in the initialisation solve. If this is
297 // not the case and the reduced problem is infeasible or unbounded the decomposition simplex will terminate.
298 // If the decomposition simplex terminates, then the original problem will be solved using the stored basis.
299 if(_solver.status() != SPxSolverBase<R>::OPTIMAL)
300 {
301 if(_solver.status() == SPxSolverBase<R>::UNBOUNDED)
302 MSG_INFO2(spxout, spxout << "Unbounded reduced problem." << std::endl);
303
304 if(_solver.status() == SPxSolverBase<R>::INFEASIBLE)
305 MSG_INFO2(spxout, spxout << "Infeasible reduced problem." << std::endl);
306
307 MSG_INFO2(spxout, spxout << "Reduced problem status: " << static_cast<int>
308 (_solver.status()) << std::endl);
309
310 redProbError = true;
311 break;
312 }
313
314 // as a final check, if no iterations were performed with the reduced problem, then the algorithm terminates
315 if(_statistics->iterations == previter)
316 {
317 MSG_WARNING(spxout,
318 spxout << "WIMDSM02: reduced problem performed zero iterations. Terminating." << std::endl;);
319
320 noRedprobIter = true;
321 stop = true;
322 break;
323 }
324
325 // printing display line
326 printDecompDisplayLine(_solver, orig_verbosity, !algIterCount, !algIterCount);
327
328 // get the dual solutions from the reduced problem
329 VectorBase<R> reducedLPDualVector(_solver.nRows());
330 VectorBase<R> reducedLPRedcostVector(_solver.nCols());
331 _solver.getDualSol(reducedLPDualVector);
332 _solver.getRedCostSol(reducedLPRedcostVector);
333
334
335 // Using the solution to the reduced problem:
336 // In the first iteration of the algorithm create the complementary problem.
337 // In the subsequent iterations of the algorithm update the complementary problem.
338 if(algIterCount == 0)
339 _formDecompComplementaryProblem();
340 else
341 {
342 if(boolParam(SoPlexBase<R>::USECOMPDUAL))
343 _updateDecompComplementaryDualProblem(false);
344 else
345 _updateDecompComplementaryPrimalProblem(false);
346 }
347
348 // setting the current solving mode.
349 _currentProb = DECOMP_COMP;
350
351 // solve the complementary problem
352 MSG_INFO2(spxout,
353 spxout << std::endl
354 << "=========================" << std::endl
355 << "Solving: Complementary Problem." << std::endl
356 << "=========================" << std::endl
357 << std::endl);
358
359 if(!explicitviol)
360 {
361 //_compSolver.writeFileLPBase("comp.lp");
362
363 _decompSimplifyAndSolve(_compSolver, _compSlufactor, true, true);
364
365 MSG_INFO2(spxout, spxout << "Iteration " << algIterCount
366 << "Objective Value: " << std::setprecision(10) << _compSolver.objValue()
367 << std::endl);
368 }
369
370
371 assert(_isRealLPLoaded);
372
373
374 // Check whether the complementary problem is solved with a non-negative objective function, is infeasible or
375 // unbounded. If true, then stop the algorithm.
376 if(!explicitviol && (GE(_compSolver.objValue(), R(0.0), R(1e-20))
377 || _compSolver.status() == SPxSolverBase<R>::INFEASIBLE
378 || _compSolver.status() == SPxSolverBase<R>::UNBOUNDED))
379 {
380 _statistics->compProbStatus = _compSolver.status();
381 _statistics->finalCompObj = _compSolver.objValue();
382
383 if(_compSolver.status() == SPxSolverBase<R>::UNBOUNDED)
384 MSG_INFO2(spxout, spxout << "Unbounded complementary problem." << std::endl);
385
386 if(_compSolver.status() == SPxSolverBase<R>::INFEASIBLE)
387 MSG_INFO2(spxout, spxout << "Infeasible complementary problem." << std::endl);
388
389 if(_compSolver.status() == SPxSolverBase<R>::INFEASIBLE
390 || _compSolver.status() == SPxSolverBase<R>::UNBOUNDED)
391 explicitviol = true;
392
393 stop = true;
394 }
395
396 if(!stop && !explicitviol)
397 {
398 // get the primal solutions from the complementary problem
399 VectorBase<R> compLPPrimalVector(_compSolver.nCols());
400 _compSolver.getPrimalSol(compLPPrimalVector);
401
402 // get the dual solutions from the complementary problem
403 VectorBase<R> compLPDualVector(_compSolver.nRows());
404 _compSolver.getDualSol(compLPDualVector);
405
406 // updating the reduced problem
407 _updateDecompReducedProblem(_compSolver.objValue(), reducedLPDualVector, reducedLPRedcostVector,
408 compLPPrimalVector, compLPDualVector);
409 }
410 // if the complementary problem is infeasible or unbounded, it is possible that the algorithm can continue.
411 // a check of the original problem is required to determine whether there are any violations.
412 else if(_compSolver.status() == SPxSolverBase<R>::INFEASIBLE
413 || _compSolver.status() == SPxSolverBase<R>::UNBOUNDED
414 || explicitviol)
415 {
416 // getting the primal vector from the reduced problem
417 VectorBase<R> reducedLPPrimalVector(_solver.nCols());
418 _solver.getPrimalSol(reducedLPPrimalVector);
419
420 // checking the optimality of the reduced problem solution with the original problem
421 _checkOriginalProblemOptimality(reducedLPPrimalVector, true);
422
423 // if there are any violated rows or bounds then stop is reset and the algorithm continues.
424 if(_nDecompViolBounds > 0 || _nDecompViolRows > 0)
425 stop = false;
426
427 if(_nDecompViolBounds == 0 && _nDecompViolRows == 0)
428 stop = true;
429
430 // updating the reduced problem with the original problem violated rows
431 if(!stop)
432 _updateDecompReducedProblemViol(false);
433 }
434
435
436 // =============================================================================
437 // Code check completed up to here
438 // =============================================================================
439
440 numDecompIter++;
441 algIterCount++;
442 }
443
444 // if there is an error in solving the reduced problem, i.e. infeasible or unbounded, then we leave the
445 // decomposition solve and resolve the original. Infeasibility should be dealt with in the original problem.
446 if(!redProbError)
447 {
448 #ifndef NDEBUG
449 // computing the solution for the original variables
450 VectorBase<R> reducedLPPrimalVector(_solver.nCols());
451 _solver.getPrimalSol(reducedLPPrimalVector);
452
453 // checking the optimality of the reduced problem solution with the original problem
454 _checkOriginalProblemOptimality(reducedLPPrimalVector, true);
455
456 // resetting the verbosity level
457 spxout.setVerbosity(orig_verbosity);
458 #endif
459
460 // This is an additional check to ensure that the complementary problem is solving correctly.
461 // This is only used for debugging
462 #ifdef PERFORM_COMPPROB_CHECK
463
464 // solving the complementary problem with the original objective function
465 if(boolParam(SoPlexBase<R>::USECOMPDUAL))
466 _setComplementaryDualOriginalObjective();
467 else
468 _setComplementaryPrimalOriginalObjective();
469
470 // for clean up
471 // the real lp has to be set to not loaded.
472 //_isRealLPLoaded = false;
473 // the basis has to be set to false
474 _hasBasis = false;
475
476 if(boolParam(SoPlexBase<R>::USECOMPDUAL))
477 {
478 SPxLPBase<R> compDualLP;
479 _compSolver.buildDualProblem(compDualLP, _decompPrimalRowIDs.get_ptr(),
480 _decompPrimalColIDs.get_ptr(),
481 _decompDualRowIDs.get_ptr(), _decompDualColIDs.get_ptr(), &_nPrimalRows, &_nPrimalCols, &_nDualRows,
482 &_nDualCols);
483
484 _compSolver.loadLP(compDualLP);
485 }
486
487 _decompSimplifyAndSolve(_compSolver, _compSlufactor, true, true);
488
489 _solReal._hasPrimal = true;
490 _hasSolReal = true;
491 // get the primal solutions from the reduced problem
492 VectorBase<R> testPrimalVector(_compSolver.nCols());
493 _compSolver.getPrimalSol(testPrimalVector);
494 _solReal._primal.reDim(_compSolver.nCols());
495 _solReal._primal = testPrimalVector;
496
497 R maxviol = 0;
498 R sumviol = 0;
499
500 // Since the original objective value has now been installed in the complementary problem, the solution will be
501 // a solution to the original problem.
502 // checking the bound violation of the solution from complementary problem
503 if(getDecompBoundViolation(maxviol, sumviol))
504 MSG_INFO1(spxout, spxout << "Bound violation - "
505 << "Max: " << std::setprecision(20) << maxviol << " "
506 << "Sum: " << sumviol << std::endl);
507
508 // checking the row violation of the solution from the complementary problem
509 if(getDecompRowViolation(maxviol, sumviol))
510 MSG_INFO1(spxout, spxout << "Row violation - "
511 << "Max: " << std::setprecision(21) << maxviol << " "
512 << "Sum: " << sumviol << std::endl);
513
514 MSG_INFO1(spxout, spxout << "Objective Value: " << _compSolver.objValue() << std::endl);
515 #endif
516 }
517
518 // resetting the verbosity level
519 spxout.setVerbosity(orig_verbosity);
520
521 MSG_INFO1(spxout,
522 spxout << "======== Decomposition solve completed ========" << std::endl
523 << std::endl
524 << "======== Resolving original problem ========" << std::endl
525 );
526
527 // if there is a reduced problem error in the first iteration the complementary problme has not been
528 // set up. In this case no memory has been allocated for _decompCompProbColIDsIdx and _fixedOrigVars.
529 if((!redProbError && !noRedprobIter) || algIterCount > 0)
530 {
531 spx_free(_decompCompProbColIDsIdx);
532 spx_free(_fixedOrigVars);
533 }
534
535 spx_free(_decompViolatedRows);
536 spx_free(_decompViolatedBounds);
537 spx_free(_decompReducedProbCols);
538 spx_free(_decompReducedProbRows);
539
540 // retreiving the original problem statistics prior to destroying it.
541 getOriginalProblemStatistics();
542
543 // setting the reduced problem statistics
544 _statistics->numRedProbRows = numIncludedRows;
545 _statistics->numRedProbCols = _solver.nCols();
546
547
548 // printing display line for resolve of problem
549 _solver.printDisplayLine(false, true);
550
551 // if there is an error solving the reduced problem the LP must be solved from scratch
552 if(redProbError)
553 {
554 MSG_INFO1(spxout, spxout << "=== Reduced problem error - Solving original ===" << std::endl);
555
556 // the solver is loaded with the realLP to solve the problem from scratch
557 _solver.loadLP(*_realLP);
558 spx_free(_realLP);
559 _realLP = &_solver;
560 _isRealLPLoaded = true;
561
562 // returning the sense to minimise
563 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE)
564 {
565 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE);
566
567 _solver.changeObj(-(_solver.maxObj()));
568 _solver.changeSense(SPxLPBase<R>::MINIMIZE);
569
570 // Need to add commands to multiply the objective solution values by -1
571 }
572
573 //_solver.setBasis(basisStatusRows.get_const_ptr(), basisStatusCols.get_const_ptr());
574 _preprocessAndSolveReal(true);
575 }
576 else
577 {
578 // resolving the problem to update the real lp and solve with the correct objective.
579 // the realLP is updated with the current solver. This is to resolve the problem to get the correct solution
580 // values
581 _realLP->~SPxLPBase<R>();
582 spx_free(_realLP);
583 _realLP = &_solver;
584 _isRealLPLoaded = true;
585
586 // returning the sense to minimise
587 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE)
588 {
589 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE);
590
591 _solver.changeObj(-(_solver.maxObj()));
592 _solver.changeSense(SPxLPBase<R>::MINIMIZE);
593
594 // Need to add commands to multiply the objective solution values by -1
595 }
596
597 _preprocessAndSolveReal(false);
598 }
599
600 // stop timing
601 _statistics->solvingTime->stop();
602 }
603
604
605
606 /// creating copies of the original problem that will be manipulated to form the reduced and complementary problems
607 template <class R>
608 void SoPlexBase<R>::_createDecompReducedAndComplementaryProblems()
609 {
610 // the reduced problem is formed from the current problem
611 // So, we copy the _solver to the _realLP and work on the _solver
612 // NOTE: there is no need to preprocess because we always have a starting basis.
613 _realLP = nullptr;
614 spx_alloc(_realLP);
615 _realLP = new(_realLP) SPxLPBase<R>(_solver);
616
617 // allocating memory for the reduced problem rows and cols flag array
618 _decompReducedProbRows = nullptr;
619 spx_alloc(_decompReducedProbRows, numRows());
620 _decompReducedProbCols = nullptr;
621 spx_alloc(_decompReducedProbCols, numCols());
622
623 // the complementary problem is formulated with all incompatible rows and those from the reduced problem that have
624 // a positive reduced cost.
625 _compSolver = _solver;
626 _compSolver.setOutstream(spxout);
627 _compSolver.setBasisSolver(&_compSlufactor);
628
629 // allocating memory for the violated bounds and rows arrays
630 _decompViolatedBounds = nullptr;
631 _decompViolatedRows = nullptr;
632 spx_alloc(_decompViolatedBounds, numCols());
633 spx_alloc(_decompViolatedRows, numRows());
634 _nDecompViolBounds = 0;
635 _nDecompViolRows = 0;
636 }
637
638
639
640 /// forms the reduced problem
641 template <class R>
642 void SoPlexBase<R>::_formDecompReducedProblem(bool& stop)
643 {
644 MSG_INFO2(spxout, spxout << "Forming the Reduced problem" << std::endl);
645 int* nonposind = 0;
646 int* compatind = 0;
647 int* rowsforremoval = 0;
648 int* colsforremoval = 0;
649 int nnonposind = 0;
650 int ncompatind = 0;
651
652 assert(_solver.nCols() == numCols());
653 assert(_solver.nRows() == numRows());
654
655 // capturing the basis used for the transformation
656 _decompTransBasis = _solver.basis();
657
658 // setting row counter to zero
659 numIncludedRows = 0;
660
661 // the _decompLP is used as a helper LP object. This is used so that _realLP can still hold the original problem.
662 _decompLP = 0;
663 spx_alloc(_decompLP);
664 _decompLP = new(_decompLP) SPxLPBase<R>(_solver);
665
666 // retreiving the basis information
667 _basisStatusRows.reSize(numRows());
668 _basisStatusCols.reSize(numCols());
669 _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr());
670
671 // get the indices of the rows with positive dual multipliers and columns with positive reduced costs.
672 spx_alloc(nonposind, numCols());
673 spx_alloc(colsforremoval, numCols());
674
675 if(!stop)
676 _getZeroDualMultiplierIndices(_decompFeasVector, nonposind, colsforremoval, &nnonposind, stop);
677
678 // get the compatible columns from the constraint matrix w.r.t the current basis matrix
679 MSG_INFO2(spxout, spxout << "Computing the compatible columns" << std::endl
680 << "Solving time: " << solveTime() << std::endl);
681
682 spx_alloc(compatind, _solver.nRows());
683 spx_alloc(rowsforremoval, _solver.nRows());
684
685 if(!stop)
686 _getCompatibleColumns(_decompFeasVector, nonposind, compatind, rowsforremoval, colsforremoval,
687 nnonposind,
688 &ncompatind, true, stop);
689
690 int* compatboundcons = 0;
691 int ncompatboundcons = 0;
692 spx_alloc(compatboundcons, numCols());
693
694 LPRowSetBase<R> boundcons;
695
696 // identifying the compatible bound constraints
697 if(!stop)
698 _getCompatibleBoundCons(boundcons, compatboundcons, nonposind, &ncompatboundcons, nnonposind, stop);
699
700 // delete rows and columns from the LP to form the reduced problem
701 MSG_INFO2(spxout, spxout << "Deleting rows and columns to form the reduced problem" << std::endl
702 << "Solving time: " << solveTime() << std::endl);
703
704 // allocating memory to add bound constraints
705 SPxRowId* addedrowids = 0;
706 spx_alloc(addedrowids, ncompatboundcons);
707
708 // computing the reduced problem obj coefficient vector and updating the problem
709 if(!stop)
710 {
711 _computeReducedProbObjCoeff(stop);
712
713 _solver.loadLP(*_decompLP);
714
715 // the colsforremoval are the columns with a zero reduced cost.
716 // the rowsforremoval are the rows identified as incompatible.
717 _solver.removeRows(rowsforremoval);
718 // adding the rows for the compatible bound constraints
719 _solver.addRows(addedrowids, boundcons);
720
721 for(int i = 0; i < ncompatboundcons; i++)
722 _decompReducedProbColRowIDs[compatboundcons[i]] = addedrowids[i];
723 }
724
725
726 // freeing allocated memory
727 spx_free(addedrowids);
728 spx_free(compatboundcons);
729 spx_free(rowsforremoval);
730 spx_free(compatind);
731 spx_free(colsforremoval);
732 spx_free(nonposind);
733
734 _decompLP->~SPxLPBase<R>();
735 spx_free(_decompLP);
736 }
737
738
739
740 /// forms the complementary problem
741 template <class R>
742 void SoPlexBase<R>::_formDecompComplementaryProblem()
743 {
744 // delete rows and columns from the LP to form the reduced problem
745 _nElimPrimalRows = 0;
746 _decompElimPrimalRowIDs.reSize(
747 _realLP->nRows()); // the number of eliminated rows is less than the number of rows
748 // in the reduced problem
749 _decompCompProbColIDsIdx = 0;
750 spx_alloc(_decompCompProbColIDsIdx, _realLP->nCols());
751
752 _fixedOrigVars = 0;
753 spx_alloc(_fixedOrigVars, _realLP->nCols());
754
755 for(int i = 0; i < _realLP->nCols(); i++)
756 {
757 _decompCompProbColIDsIdx[i] = -1;
758 _fixedOrigVars[i] = -2; // this must be set to a value differet from -1, 0, 1 to ensure the
759 // _updateComplementaryFixedPrimalVars function updates all variables in the problem.
760 }
761
762 int naddedrows = 0;
763 DataArray < SPxRowId > rangedRowIds(numRows());
764 _deleteAndUpdateRowsComplementaryProblem(rangedRowIds.get_ptr(), naddedrows);
765
766 if(boolParam(
767 SoPlexBase<R>::USECOMPDUAL)) // if we use the dual formulation of the complementary problem, we must
768 // perform the transformation
769 {
770 // initialising the arrays to store the row id's from the primal and the col id's from the dual
771 _decompPrimalRowIDs.reSize(3 * _compSolver.nRows());
772 _decompPrimalColIDs.reSize(3 * _compSolver.nCols());
773 _decompDualRowIDs.reSize(3 * _compSolver.nCols());
774 _decompDualColIDs.reSize(3 * _compSolver.nRows());
775 _nPrimalRows = 0;
776 _nPrimalCols = 0;
777 _nDualRows = 0;
778 _nDualCols = 0;
779
780 // convert complementary problem to dual problem
781 SPxLPBase<R> compDualLP;
782 _compSolver.buildDualProblem(compDualLP, _decompPrimalRowIDs.get_ptr(),
783 _decompPrimalColIDs.get_ptr(),
784 _decompDualRowIDs.get_ptr(), _decompDualColIDs.get_ptr(), &_nPrimalRows, &_nPrimalCols, &_nDualRows,
785 &_nDualCols);
786 compDualLP.setOutstream(spxout);
787
788 // setting the id index array for the complementary problem
789 for(int i = 0; i < _nPrimalRows; i++)
790 {
791 // a primal row may result in two dual columns. So the index array points to the first of the indicies for the
792 // primal row.
793 if(i + 1 < _nPrimalRows &&
794 _realLP->number(SPxRowId(_decompPrimalRowIDs[i])) == _realLP->number(SPxRowId(
795 _decompPrimalRowIDs[i + 1])))
796 i++;
797 }
798
799 for(int i = 0; i < _nPrimalCols; i++)
800 {
801 if(_decompPrimalColIDs[i].getIdx() != _compSlackColId.getIdx())
802 _decompCompProbColIDsIdx[_realLP->number(_decompPrimalColIDs[i])] = i;
803 }
804
805 // retrieving the dual row id for the complementary slack column
806 // there should be a one to one relationship between the number of primal columns and the number of dual rows.
807 // hence, it should be possible to equate the dual row id to the related primal column.
808 assert(_nPrimalCols == _nDualRows);
809 assert(_compSolver.nCols() == compDualLP.nRows());
810 _compSlackDualRowId = compDualLP.rId(_compSolver.number(_compSlackColId));
811
812 _compSolver.loadLP(compDualLP);
813
814 _compSolver.init();
815
816 _updateComplementaryDualSlackColCoeff();
817
818 // initalising the array for the dual columns representing the fixed variables in the complementary problem.
819 _decompFixedVarDualIDs.reSize(_realLP->nCols());
820 _decompVarBoundDualIDs.reSize(2 * _realLP->nCols());
821
822
823 // updating the constraints for the complementary problem
824 _updateDecompComplementaryDualProblem(false);
825 }
826 else // when using the primal of the complementary problem, no transformation of the problem is required.
827 {
828 // initialising the arrays to store the row and column id's from the original problem the col id's from the dual
829 // if a row or column is not included in the complementary problem, the ID is set to INVALID in the
830 // _decompPrimalRowIDs or _decompPrimalColIDs respectively.
831 _decompPrimalRowIDs.reSize(_compSolver.nRows() * 2);
832 _decompPrimalColIDs.reSize(_compSolver.nCols());
833 _decompCompPrimalRowIDs.reSize(_compSolver.nRows() * 2);
834 _decompCompPrimalColIDs.reSize(_compSolver.nCols());
835 _nPrimalRows = 0;
836 _nPrimalCols = 0;
837
838 // at this point in time the complementary problem contains all rows and columns from the original problem
839 int addedrangedrows = 0;
840 _nCompPrimalRows = 0;
841 _nCompPrimalCols = 0;
842
843 for(int i = 0; i < _realLP->nRows(); i++)
844 {
845 assert(_nCompPrimalRows < _compSolver.nRows());
846 _decompPrimalRowIDs[_nPrimalRows] = _realLP->rId(i);
847 _decompCompPrimalRowIDs[_nCompPrimalRows] = _compSolver.rId(i);
848 _nPrimalRows++;
849 _nCompPrimalRows++;
850
851 if(_realLP->rowType(i) == LPRowBase<R>::RANGE || _realLP->rowType(i) == LPRowBase<R>::EQUAL)
852 {
853 assert(EQ(_compSolver.lhs(rangedRowIds[addedrangedrows]), _realLP->lhs(i)));
854 assert(LE(_compSolver.rhs(rangedRowIds[addedrangedrows]), R(infinity)));
855 assert(LE(_compSolver.lhs(i), R(-infinity)));
856 assert(LT(_compSolver.rhs(i), R(infinity)));
857
858 _decompPrimalRowIDs[_nPrimalRows] = _realLP->rId(i);
859 _decompCompPrimalRowIDs[_nCompPrimalRows] = rangedRowIds[addedrangedrows];
860 _nPrimalRows++;
861 _nCompPrimalRows++;
862 addedrangedrows++;
863 }
864 }
865
866 for(int i = 0; i < _compSolver.nCols(); i++)
867 {
868
869 if(_compSolver.cId(i).getIdx() != _compSlackColId.getIdx())
870 {
871 _decompPrimalColIDs[i] = _realLP->cId(i);
872 _decompCompPrimalColIDs[i] = _compSolver.cId(i);
873 _nPrimalCols++;
874 _nCompPrimalCols++;
875
876 _decompCompProbColIDsIdx[_realLP->number(_decompPrimalColIDs[i])] = i;
877 }
878 }
879
880 // updating the constraints for the complementary problem
881 _updateDecompComplementaryPrimalProblem(false);
882 }
883 }
884
885
886
887 /// simplifies the problem and solves
888 template <class R>
889 void SoPlexBase<R>::_decompSimplifyAndSolve(SPxSolverBase<R>& solver, SLUFactor<R>& sluFactor,
890 bool fromScratch, bool applyPreprocessing)
891 {
892 if(realParam(SoPlexBase<R>::TIMELIMIT) < realParam(SoPlexBase<R>::INFTY))
893 solver.setTerminationTime(Real(realParam(SoPlexBase<R>::TIMELIMIT)) -
894 _statistics->solvingTime->time());
895
896 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET));
897 _statistics->preprocessingTime->start();
898
899 typename SPxSimplifier<R>::Result result = SPxSimplifier<R>::OKAY;
900
901 /* delete starting basis if solving from scratch */
902 if(fromScratch)
903 {
904 try
905 {
906 solver.reLoad();
907 }
908 catch(const SPxException& E)
909 {
910 MSG_ERROR(spxout << "Caught exception <" << E.what() << "> during simplify and solve.\n");
911 }
912 }
913
914 //assert(!fromScratch || solver.status() == SPxSolverBase<R>::NO_PROBLEM);
915
916 if(applyPreprocessing)
917 {
918 _enableSimplifierAndScaler();
919 solver.setTerminationValue(realParam(SoPlexBase<R>::INFTY));
920 }
921 else
922 {
923 _disableSimplifierAndScaler();
924 ///@todo implement for both objective senses
925 solver.setTerminationValue(solver.spxSense() == SPxLPBase<R>::MINIMIZE
926 ? realParam(SoPlexBase<R>::OBJLIMIT_UPPER) : realParam(SoPlexBase<R>::INFTY));
927 }
928
929 /* store original lp */
930 applyPreprocessing = (_scaler != nullptr || _simplifier != NULL);
931
932 // @TODO The use of _isRealLPLoaded is not correct. An additional parameter would be useful for this algorithm.
933 // Maybe a parameter _isDecompLPLoaded?
934 if(_isRealLPLoaded)
935 {
936 //assert(_decompLP == &solver); // there should be an assert here, but I don't know what to use.
937
938 // preprocessing is always applied to the LP in the solver; hence we have to create a copy of the original LP
939 // if preprocessing is turned on
940 if(applyPreprocessing)
941 {
942 _decompLP = 0;
943 spx_alloc(_decompLP);
944 _decompLP = new(_decompLP) SPxLPBase<R>(solver);
945 _isRealLPLoaded = false;
946 }
947 else
948 _decompLP = &solver;
949 }
950 else
951 {
952 assert(_decompLP != &solver);
953
954 // ensure that the solver has the original problem
955 solver.loadLP(*_decompLP);
956
957 // load basis if available
958 if(_hasBasis)
959 {
960 assert(_basisStatusRows.size() == solver.nRows());
961 assert(_basisStatusCols.size() == solver.nCols());
962
963 ///@todo this should not fail even if the basis is invalid (wrong dimension or wrong number of basic
964 /// entries); fix either in SPxSolverBase or in SPxBasisBase
965 solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr());
966 }
967
968 // if there is no preprocessing, then the original and the transformed problem are identical and it is more
969 // memory-efficient to keep only the problem in the solver
970 if(!applyPreprocessing)
971 {
972 _decompLP->~SPxLPBase<R>();
973 spx_free(_decompLP);
974 _decompLP = &solver;
975 _isRealLPLoaded = true;
976 }
977 else
978 {
979 _decompLP = 0;
980 spx_alloc(_decompLP);
981 _decompLP = new(_decompLP) SPxLPBase<R>(solver);
982 }
983 }
984
985 // assert that we have two problems if and only if we apply preprocessing
986 assert(_decompLP == &solver || applyPreprocessing);
987 assert(_decompLP != &solver || !applyPreprocessing);
988
989 // apply problem simplification
990 if(_simplifier != 0)
991 {
992 Real remainingTime = _solver.getMaxTime() - _solver.time();
993 result = _simplifier->simplify(solver, realParam(SoPlexBase<R>::EPSILON_ZERO),
994 realParam(SoPlexBase<R>::FEASTOL),
995 realParam(SoPlexBase<R>::OPTTOL),
996 remainingTime);
997 solver.changeObjOffset(_simplifier->getObjoffset() + realParam(SoPlexBase<R>::OBJ_OFFSET));
998 }
999
1000 _statistics->preprocessingTime->stop();
1001
1002 // run the simplex method if problem has not been solved by the simplifier
1003 if(result == SPxSimplifier<R>::OKAY)
1004 {
1005 if(_scaler != nullptr)
1006 _scaler->scale(solver);
1007
1008 bool _hadBasis = _hasBasis;
1009
1010 _statistics->simplexTime->start();
1011
1012 try
1013 {
1014 solver.solve();
1015 }
1016 catch(const SPxException& E)
1017 {
1018 MSG_ERROR(std::cerr << "Caught exception <" << E.what() << "> while solving real LP.\n");
1019 _status = SPxSolverBase<R>::ERROR;
1020 }
1021 catch(...)
1022 {
1023 MSG_ERROR(std::cerr << "Caught unknown exception while solving real LP.\n");
1024 _status = SPxSolverBase<R>::ERROR;
1025 }
1026
1027 _statistics->simplexTime->stop();
1028
1029 // record statistics
1030 // only record the main statistics for the original problem and reduced problem.
1031 // the complementary problem is a pivoting rule, so these will be held in other statistics.
1032 // statistics on the number of iterations for each problem is stored in individual variables.
1033 if(_currentProb == DECOMP_ORIG || _currentProb == DECOMP_RED)
1034 {
1035 _statistics->iterations += solver.iterations();
1036 _statistics->iterationsPrimal += solver.primalIterations();
1037 _statistics->iterationsFromBasis += _hadBasis ? solver.iterations() : 0;
1038 _statistics->boundflips += solver.boundFlips();
1039 _statistics->luFactorizationTimeReal += sluFactor.getFactorTime();
1040 _statistics->luSolveTimeReal += sluFactor.getSolveTime();
1041 _statistics->luFactorizationsReal += sluFactor.getFactorCount();
1042 _statistics->luSolvesReal += sluFactor.getSolveCount();
1043 sluFactor.resetCounters();
1044
1045 _statistics->degenPivotsPrimal += solver.primalDegeneratePivots();
1046 _statistics->degenPivotsDual += solver.dualDegeneratePivots();
1047
1048 _statistics->sumDualDegen += _solver.sumDualDegeneracy();
1049 _statistics->sumPrimalDegen += _solver.sumPrimalDegeneracy();
1050
1051 if(_currentProb == DECOMP_ORIG)
1052 _statistics->iterationsInit += solver.iterations();
1053 else
1054 _statistics->iterationsRedProb += solver.iterations();
1055 }
1056
1057 if(_currentProb == DECOMP_COMP)
1058 _statistics->iterationsCompProb += solver.iterations();
1059
1060 }
1061
1062 // check the result and run again without preprocessing if necessary
1063 _evaluateSolutionDecomp(solver, sluFactor, result);
1064 }
1065
1066
1067
1068 /// loads original problem into solver and solves again after it has been solved to optimality with preprocessing
1069 template <class R>
1070 void SoPlexBase<R>::_decompResolveWithoutPreprocessing(SPxSolverBase<R>& solver,
1071 SLUFactor<R>& sluFactor, typename SPxSimplifier<R>::Result result)
1072 {
1073 assert(_simplifier != 0 || _scaler != nullptr);
1074 assert(result == SPxSimplifier<R>::VANISHED
1075 || (result == SPxSimplifier<R>::OKAY
1076 && (solver.status() == SPxSolverBase<R>::OPTIMAL
1077 || solver.status() == SPxSolverBase<R>::ABORT_DECOMP
1078 || solver.status() == SPxSolverBase<R>::ABORT_EXDECOMP)));
1079
1080 // if simplifier is active and LP is solved in presolving or to optimality, then we unsimplify to get the basis
1081 if(_simplifier != 0)
1082 {
1083 assert(!_simplifier->isUnsimplified());
1084
1085 bool vanished = result == SPxSimplifier<R>::VANISHED;
1086
1087 // get solution vectors for transformed problem
1088 VectorBase<R> primal(vanished ? 0 : solver.nCols());
1089 VectorBase<R> slacks(vanished ? 0 : solver.nRows());
1090 VectorBase<R> dual(vanished ? 0 : solver.nRows());
1091 VectorBase<R> redCost(vanished ? 0 : solver.nCols());
1092
1093 assert(!_isRealLPLoaded);
1094 _basisStatusRows.reSize(_decompLP->nRows());
1095 _basisStatusCols.reSize(_decompLP->nCols());
1096 assert(vanished || _basisStatusRows.size() >= solver.nRows());
1097 assert(vanished || _basisStatusCols.size() >= solver.nCols());
1098
1099 if(!vanished)
1100 {
1101 assert(solver.status() == SPxSolverBase<R>::OPTIMAL
1102 || solver.status() == SPxSolverBase<R>::ABORT_DECOMP
1103 || solver.status() == SPxSolverBase<R>::ABORT_EXDECOMP);
1104
1105 solver.getPrimalSol(primal);
1106 solver.getSlacks(slacks);
1107 solver.getDualSol(dual);
1108 solver.getRedCostSol(redCost);
1109
1110 // unscale vectors
1111 if(_scaler && solver.isScaled())
1112 {
1113 _scaler->unscalePrimal(solver, primal);
1114 _scaler->unscaleSlacks(solver, slacks);
1115 _scaler->unscaleDual(solver, dual);
1116 _scaler->unscaleRedCost(solver, redCost);
1117 }
1118
1119 // get basis of transformed problem
1120 solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr());
1121 }
1122
1123 try
1124 {
1125 _simplifier->unsimplify(primal, dual, slacks, redCost, _basisStatusRows.get_ptr(),
1126 _basisStatusCols.get_ptr(), solver.status() == SPxSolverBase<R>::OPTIMAL);
1127 _simplifier->getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr());
1128 _hasBasis = true;
1129 }
1130 catch(const SPxException& E)
1131 {
1132 MSG_ERROR(spxout << "Caught exception <" << E.what() <<
1133 "> during unsimplification. Resolving without simplifier and scaler.\n");
1134 }
1135 catch(...)
1136 {
1137 MSG_ERROR(spxout <<
1138 "Caught unknown exception during unsimplification. Resolving without simplifier and scaler.\n");
1139 _status = SPxSolverBase<R>::ERROR;
1140 }
1141 }
1142 // if the original problem is not in the solver because of scaling, we also need to store the basis
1143 else if(_scaler != nullptr)
1144 {
1145 _basisStatusRows.reSize(numRows());
1146 _basisStatusCols.reSize(numCols());
1147 assert(_basisStatusRows.size() == solver.nRows());
1148 assert(_basisStatusCols.size() == solver.nCols());
1149
1150 solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr());
1151 _hasBasis = true;
1152 }
1153
1154 // resolve the original problem
1155 _decompSimplifyAndSolve(solver, sluFactor, false, false);
1156 return;
1157 }
1158
1159
1160 /// updates the reduced problem with additional rows using the solution to the complementary problem
1161 template <class R>
1162 void SoPlexBase<R>::_updateDecompReducedProblem(R objValue, VectorBase<R> dualVector,
1163 VectorBase<R> redcostVector,
1164 VectorBase<R> compPrimalVector, VectorBase<R> compDualVector)
1165 {
1166 R feastol = realParam(SoPlexBase<R>::FEASTOL);
1167
1168 R maxDualRatio = R(infinity);
1169
1170 bool usecompdual = boolParam(SoPlexBase<R>::USECOMPDUAL);
1171
1172 for(int i = 0; i < _nPrimalRows; i++)
1173 {
1174 R reducedProbDual = 0;
1175 R compProbPrimal = 0;
1176 R dualRatio = 0;
1177 int rownumber = _realLP->number(SPxRowId(_decompPrimalRowIDs[i]));
1178
1179 if(_decompReducedProbRows[rownumber] && _solver.isBasic(_decompReducedProbRowIDs[rownumber]))
1180 {
1181 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rownumber]);
1182
1183 // retreiving the reduced problem dual solutions and the complementary problem primal solutions
1184 reducedProbDual = dualVector[solverRowNum]; // this is y
1185
1186 if(usecompdual)
1187 compProbPrimal = compPrimalVector[_compSolver.number(SPxColId(_decompDualColIDs[i]))]; // this is u
1188 else
1189 compProbPrimal = compDualVector[_compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i]))];
1190
1191 // the variable in the basis is degenerate.
1192 if(EQ(reducedProbDual, R(0.0), feastol))
1193 {
1194 MSG_WARNING(spxout,
1195 spxout << "WIMDSM01: reduced problem dual value is very close to zero." << std::endl;);
1196 continue;
1197 }
1198
1199 // the translation of the complementary primal problem to the dual some rows resulted in two columns.
1200 if(usecompdual && i < _nPrimalRows - 1 &&
1201 _realLP->number(SPxRowId(_decompPrimalRowIDs[i])) == _realLP->number(SPxRowId(
1202 _decompPrimalRowIDs[i + 1])))
1203 {
1204 i++;
1205 compProbPrimal += compPrimalVector[_compSolver.number(SPxColId(_decompDualColIDs[i]))]; // this is u
1206 }
1207
1208
1209
1210 // updating the ratio
1211 SoPlexBase<R>::DualSign varSign = getExpectedDualVariableSign(solverRowNum);
1212
1213 if(varSign == SoPlexBase<R>::IS_FREE || (varSign == SoPlexBase<R>::IS_POS
1214 && LE(compProbPrimal, (R)0, feastol)) ||
1215 (varSign == SoPlexBase<R>::IS_NEG && GE(compProbPrimal, (R)0, feastol)))
1216 {
1217 dualRatio = R(infinity);
1218 }
1219 else
1220 {
1221 dualRatio = reducedProbDual / compProbPrimal;
1222 }
1223
1224 if(LT(dualRatio, maxDualRatio, feastol))
1225 maxDualRatio = dualRatio;
1226 }
1227 else
1228 {
1229 if(usecompdual)
1230 compProbPrimal = compPrimalVector[_compSolver.number(SPxColId(_decompDualColIDs[i]))]; // this is v
1231 else
1232 compProbPrimal = compDualVector[_compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i]))];
1233 }
1234
1235 }
1236
1237 // setting the maxDualRatio to a maximum of 1
1238 if(maxDualRatio > 1.0)
1239 maxDualRatio = 1.0;
1240
1241 VectorBase<R> compProbRedcost(
1242 _compSolver.nCols()); // the reduced costs of the complementary problem
1243
1244 // Retrieving the reduced costs for each original problem row.
1245 _compSolver.getRedCostSol(compProbRedcost);
1246
1247 LPRowSetBase<R> updaterows;
1248
1249
1250 // Identifying the violated rows
1251 Array<RowViolation> violatedrows;
1252 int nviolatedrows = 0;
1253 int* newrowidx = 0;
1254 int nnewrowidx = 0;
1255 spx_alloc(newrowidx, _nPrimalRows);
1256
1257 violatedrows.reSize(_nPrimalRows);
1258
1259 bool ratioTest = true;
1260
1261 //ratioTest = false;
1262 for(int i = 0; i < _nPrimalRows; i++)
1263 {
1264 LPRowBase<R> origlprow;
1265 DSVectorBase<R> rowtoaddVec(_realLP->nCols());
1266 R compProbPrimal = 0;
1267 R compRowRedcost = 0;
1268 int rownumber = _realLP->number(SPxRowId(_decompPrimalRowIDs[i]));
1269
1270 if(!_decompReducedProbRows[_realLP->number(SPxRowId(_decompPrimalRowIDs[i]))])
1271 {
1272 int compRowNumber;
1273
1274 if(usecompdual)
1275 {
1276 compRowNumber = _compSolver.number(_decompDualColIDs[i]);
1277 // retreiving the complementary problem primal solutions
1278 compProbPrimal = compPrimalVector[compRowNumber]; // this is v
1279 compRowRedcost = compProbRedcost[compRowNumber];
1280 }
1281 else
1282 {
1283 compRowNumber = _compSolver.number(_decompCompPrimalRowIDs[i]);
1284 // retreiving the complementary problem primal solutions
1285 compProbPrimal = compDualVector[compRowNumber]; // this is v
1286 }
1287
1288 // the translation of the complementary primal problem to the dual some rows resulted in two columns.
1289 if(usecompdual && i < _nPrimalRows - 1 &&
1290 _realLP->number(SPxRowId(_decompPrimalRowIDs[i])) == _realLP->number(SPxRowId(
1291 _decompPrimalRowIDs[i + 1])))
1292 {
1293 i++;
1294 compRowNumber = _compSolver.number(SPxColId(_decompDualColIDs[i]));
1295 compProbPrimal += compPrimalVector[compRowNumber]; // this is v
1296 compRowRedcost += compProbRedcost[compRowNumber];
1297 }
1298
1299 SoPlexBase<R>::DualSign varSign = getOrigProbDualVariableSign(rownumber);
1300
1301 // add row to the reduced problem if the computed dual is of the correct sign for a feasible dual solution
1302 if(ratioTest && ((varSign == SoPlexBase<R>::IS_FREE && !isZero(compProbPrimal, feastol)) ||
1303 (varSign == SoPlexBase<R>::IS_POS && GT(compProbPrimal * maxDualRatio, (R) 0, feastol)) ||
1304 (varSign == SoPlexBase<R>::IS_NEG && LT(compProbPrimal * maxDualRatio, (R) 0, feastol))))
1305 {
1306 //this set of statements are required to add rows to the reduced problem. This will add a row for every
1307 //violated constraint. I only want to add a row for the most violated constraint. That is why the row
1308 //adding functionality is put outside the for loop.
1309 if(!_decompReducedProbRows[rownumber])
1310 {
1311 numIncludedRows++;
1312 assert(numIncludedRows <= _realLP->nRows());
1313 }
1314
1315 violatedrows[nviolatedrows].idx = rownumber;
1316 violatedrows[nviolatedrows].violation = spxAbs(compProbPrimal * maxDualRatio);
1317 nviolatedrows++;
1318 }
1319 }
1320 }
1321
1322
1323 // sorting the violated rows by the absolute violation.
1324 // Only a predefined number of rows will be added to the reduced problem
1325 RowViolationCompare compare;
1326 compare.entry = violatedrows.get_const_ptr();
1327
1328 // if no rows are identified by the pricing rule, we add rows based upon the constraint violations
1329 if(!ratioTest || nviolatedrows == 0)
1330 {
1331 _findViolatedRows(objValue, violatedrows, nviolatedrows);
1332 }
1333
1334 int sorted = 0;
1335 int sortsize = MINIMUM(intParam(SoPlexBase<R>::DECOMP_MAXADDEDROWS), nviolatedrows);
1336
1337 // only sorting if the sort size is less than the number of violated rows.
1338 if(sortsize > 0 && sortsize < nviolatedrows)
1339 sorted = SPxQuicksortPart(violatedrows.get_ptr(), compare, sorted + 1, nviolatedrows, sortsize);
1340
1341 // adding the violated rows.
1342 for(int i = 0; i < sortsize; i++)
1343 {
1344 updaterows.add(_transformedRows.lhs(violatedrows[i].idx),
1345 _transformedRows.rowVector(violatedrows[i].idx),
1346 _transformedRows.rhs(violatedrows[i].idx));
1347
1348 _decompReducedProbRows[violatedrows[i].idx] = true;
1349 newrowidx[nnewrowidx] = violatedrows[i].idx;
1350 nnewrowidx++;
1351 }
1352
1353 SPxRowId* addedrowids = 0;
1354 spx_alloc(addedrowids, nnewrowidx);
1355 _solver.addRows(addedrowids, updaterows);
1356
1357 for(int i = 0; i < nnewrowidx; i++)
1358 _decompReducedProbRowIDs[newrowidx[i]] = addedrowids[i];
1359
1360 // freeing allocated memory
1361 spx_free(addedrowids);
1362 spx_free(newrowidx);
1363 }
1364
1365
1366
1367 /// update the reduced problem with additional columns and rows based upon the violated original bounds and rows
1368 // TODO: Allow for the case that no rows are added. This should terminate the algorithm.
1369 // TODO: Check to make sure that only rows added to the problem do not currently exist in the reduced problem.
1370 template <class R>
1371 void SoPlexBase<R>::_updateDecompReducedProblemViol(bool allrows)
1372 {
1373 #ifdef NO_TOL
1374 R feastol = 0.0;
1375 #else
1376 #ifdef USE_FEASTOL
1377 R feastol = realParam(SoPlexBase<R>::FEASTOL);
1378 #else
1379 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO);
1380 #endif
1381 #endif
1382 LPRowSetBase<R> updaterows;
1383
1384 int* newrowidx = 0;
1385 int nnewrowidx = 0;
1386 spx_alloc(newrowidx, _nPrimalRows);
1387
1388 int rowNumber;
1389 int bestrow = -1;
1390 R bestrownorm = R(infinity);
1391 R percenttoadd = 1;
1392
1393 int nrowstoadd = MINIMUM(intParam(SoPlexBase<R>::DECOMP_MAXADDEDROWS), _nDecompViolRows);
1394
1395 if(allrows)
1396 nrowstoadd = _nDecompViolRows; // adding all violated rows
1397
1398 SSVectorBase<R> y(_solver.nCols());
1399 y.unSetup();
1400
1401 // identifying the rows not included in the reduced problem that are violated by the current solution.
1402 for(int i = 0; i < nrowstoadd; i++)
1403 {
1404 rowNumber = _decompViolatedRows[i];
1405
1406 if(!allrows)
1407 {
1408 R norm = 0;
1409
1410 // the rhs of this calculation are the rows of the constraint matrix
1411 // so we are solving y B = A_{i,.}
1412 try
1413 {
1414 _solver.basis().solve(y, _solver.vector(rowNumber));
1415 }
1416 catch(const SPxException& E)
1417 {
1418 MSG_ERROR(spxout << "Caught exception <" << E.what() << "> while computing compatability.\n");
1419 }
1420
1421
1422 // comparing the constraints based upon the row norm
1423 if(y.isSetup())
1424 {
1425 for(int j = 0; j < y.size(); j++)
1426 {
1427 if(isZero(_solver.fVec()[i], feastol))
1428 norm += spxAbs(y.value(j)) * spxAbs(y.value(j));
1429 }
1430 }
1431 else
1432 {
1433 for(int j = 0; j < numCols(); j++)
1434 {
1435 if(isZero(_solver.fVec()[i], feastol))
1436 norm += spxAbs(y[j]) * spxAbs(y[j]);
1437 }
1438 }
1439
1440
1441 // the best row is based upon the row norm
1442 // the best row is added if no violated row is found
1443 norm = soplex::spxSqrt(norm);
1444
1445 if(LT(norm, bestrownorm))
1446 {
1447 bestrow = rowNumber;
1448 bestrownorm = norm;
1449 }
1450
1451
1452 // adding the violated row
1453 if(isZero(norm, feastol) && LT(nnewrowidx / R(numRows()), percenttoadd))
1454 {
1455 updaterows.add(_transformedRows.lhs(rowNumber), _transformedRows.rowVector(rowNumber),
1456 _transformedRows.rhs(rowNumber));
1457
1458 _decompReducedProbRows[rowNumber] = true;
1459 newrowidx[nnewrowidx] = rowNumber;
1460 nnewrowidx++;
1461 }
1462
1463 }
1464 else
1465 {
1466 // adding all violated rows
1467 updaterows.add(_transformedRows.lhs(rowNumber), _transformedRows.rowVector(rowNumber),
1468 _transformedRows.rhs(rowNumber));
1469
1470 _decompReducedProbRows[rowNumber] = true;
1471 newrowidx[nnewrowidx] = rowNumber;
1472 nnewrowidx++;
1473 }
1474 }
1475
1476 // if no violated row is found during the first pass of the available rows, then we add all violated rows.
1477 if(nnewrowidx == 0)
1478 {
1479 for(int i = 0; i < nrowstoadd; i++)
1480 {
1481 rowNumber = _decompViolatedRows[i];
1482
1483 updaterows.add(_transformedRows.lhs(rowNumber), _transformedRows.rowVector(rowNumber),
1484 _transformedRows.rhs(rowNumber));
1485
1486 _decompReducedProbRows[rowNumber] = true;
1487 newrowidx[nnewrowidx] = rowNumber;
1488 nnewrowidx++;
1489 }
1490 }
1491
1492 // we will always add the row that is deemed best based upon the row norm.
1493 // TODO: check whether this should be skipped if a row has already been added.
1494 if(!allrows && bestrow >= 0)
1495 {
1496 updaterows.add(_transformedRows.lhs(bestrow), _transformedRows.rowVector(bestrow),
1497 _transformedRows.rhs(bestrow));
1498
1499 _decompReducedProbRows[bestrow] = true;
1500 newrowidx[nnewrowidx] = bestrow;
1501 nnewrowidx++;
1502 }
1503
1504
1505 SPxRowId* addedrowids = 0;
1506 spx_alloc(addedrowids, nnewrowidx);
1507 _solver.addRows(addedrowids, updaterows);
1508
1509 for(int i = 0; i < nnewrowidx; i++)
1510 _decompReducedProbRowIDs[newrowidx[i]] = addedrowids[i];
1511
1512 // freeing allocated memory
1513 spx_free(addedrowids);
1514 spx_free(newrowidx);
1515 }
1516
1517
1518
1519 /// builds the update rows with those violated in the complmentary problem
1520 // A row is violated in the constraint matrix Ax <= b, if b - A_{i}x < 0
1521 // To aid the computation, all violations are translated to <= constraints
1522 template <class R>
1523 void SoPlexBase<R>::_findViolatedRows(R compObjValue, Array<RowViolation>& violatedrows,
1524 int& nviolatedrows)
1525 {
1526 R feastol = realParam(SoPlexBase<R>::FEASTOL);
1527 VectorBase<R> compProbRedcost(
1528 _compSolver.nCols()); // the reduced costs of the complementary problem
1529 VectorBase<R> compProbPrimal(
1530 _compSolver.nCols()); // the primal vector of the complementary problem
1531 VectorBase<R> compProbActivity(
1532 _compSolver.nRows()); // the activity vector of the complementary problem
1533 R compProbSlackVal = 0;
1534
1535 if(boolParam(SoPlexBase<R>::USECOMPDUAL))
1536 {
1537 // Retrieving the slacks for each row.
1538 _compSolver.getRedCostSol(compProbRedcost);
1539 }
1540 else
1541 {
1542 // Retrieving the primal vector for the complementary problem
1543 _compSolver.getPrimalSol(compProbPrimal);
1544 _compSolver.computePrimalActivity(compProbPrimal, compProbActivity);
1545
1546 compProbSlackVal = compProbPrimal[_compSolver.number(_compSlackColId)];
1547 }
1548
1549 // scanning all rows of the complementary problem for violations.
1550 for(int i = 0; i < _nPrimalRows; i++)
1551 {
1552 LPRowBase<R> origlprow;
1553 DSVectorBase<R> rowtoaddVec(_realLP->nCols());
1554 R compProbViol = 0;
1555 R compSlackCoeff = 0;
1556 int rownumber = _realLP->number(SPxRowId(_decompPrimalRowIDs[i]));
1557 int comprownum = _compSolver.number(SPxRowId(_decompPrimalRowIDs[i]));
1558
1559 if(!_decompReducedProbRows[rownumber])
1560 {
1561 // retreiving the violation of the complementary problem primal constraints
1562 if(boolParam(SoPlexBase<R>::USECOMPDUAL))
1563 {
1564 compSlackCoeff = getCompSlackVarCoeff(i);
1565 compProbViol = compProbRedcost[_compSolver.number(SPxColId(
1566 _decompDualColIDs[i]))]; // this is b - Ax
1567 // subtracting the slack variable value
1568 compProbViol += compObjValue * compSlackCoeff; // must add on the slack variable value.
1569 compProbViol *= compSlackCoeff; // translating the violation to a <= constraint
1570 }
1571 else
1572 {
1573 R viol = _compSolver.rhs(comprownum) - (compProbActivity[comprownum] + compProbSlackVal);
1574
1575 if(viol < 0.0)
1576 compProbViol = viol;
1577
1578 viol = (compProbActivity[comprownum] - compProbSlackVal) - _compSolver.lhs(comprownum);
1579
1580 if(viol < 0.0)
1581 compProbViol = viol;
1582
1583 }
1584
1585 // NOTE: if the row was originally a ranged constraint, we are only interest in one of the inequalities.
1586 // If one inequality of the range violates the bounds, then we will add the row.
1587
1588 // the translation of the complementary primal problem to the dual some rows resulted in two columns.
1589 if(boolParam(SoPlexBase<R>::USECOMPDUAL) && i < _nPrimalRows - 1 &&
1590 _realLP->number(SPxRowId(_decompPrimalRowIDs[i])) == _realLP->number(SPxRowId(
1591 _decompPrimalRowIDs[i + 1])))
1592 {
1593 i++;
1594 compSlackCoeff = getCompSlackVarCoeff(i);
1595 R tempViol = compProbRedcost[_compSolver.number(SPxColId(_decompDualColIDs[i]))]; // this is b - Ax
1596 tempViol += compObjValue * compSlackCoeff;
1597 tempViol *= compSlackCoeff;
1598
1599 // if the other side of the range constraint has a larger violation, then this is used for the
1600 // computation.
1601 if(tempViol < compProbViol)
1602 compProbViol = tempViol;
1603 }
1604
1605
1606 // checking the violation of the row.
1607 if(LT(compProbViol, (R) 0, feastol))
1608 {
1609 if(!_decompReducedProbRows[rownumber])
1610 {
1611 numIncludedRows++;
1612 assert(numIncludedRows <= _realLP->nRows());
1613 }
1614
1615 violatedrows[nviolatedrows].idx = rownumber;
1616 violatedrows[nviolatedrows].violation = spxAbs(compProbViol);
1617 nviolatedrows++;
1618 }
1619 }
1620 }
1621 }
1622
1623
1624
1625 /// identifies the columns of the row-form basis that correspond to rows with zero dual multipliers.
1626 // This function assumes that the basis is in the row form.
1627 // @todo extend this to the case when the basis is in the column form.
1628 template <class R>
1629 void SoPlexBase<R>::_getZeroDualMultiplierIndices(VectorBase<R> feasVector, int* nonposind,
1630 int* colsforremoval,
1631 int* nnonposind, bool& stop)
1632 {
1633 assert(_solver.rep() == SPxSolverBase<R>::ROW);
1634
1635 #ifdef NO_TOL
1636 R feastol = 0.0;
1637 #else
1638 #ifdef USE_FEASTOL
1639 R feastol = realParam(SoPlexBase<R>::FEASTOL);
1640 #else
1641 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO);
1642 #endif
1643 #endif
1644
1645 bool delCol;
1646
1647 _decompReducedProbColIDs.reSize(numCols());
1648 *nnonposind = 0;
1649
1650 // iterating over all columns in the basis matrix
1651 // this identifies the basis indices and the indicies that are positive.
1652 for(int i = 0; i < _solver.nCols(); ++i)
1653 {
1654 _decompReducedProbCols[i] = true;
1655 _decompReducedProbColIDs[i].inValidate();
1656 colsforremoval[i] = i;
1657 delCol = false;
1658
1659 if(_solver.basis().baseId(i).isSPxRowId()) // find the row id's for rows in the basis
1660 {
1661 // record the row if the dual multiple is zero.
1662 if(isZero(feasVector[i], feastol))
1663 {
1664 nonposind[*nnonposind] = i;
1665 (*nnonposind)++;
1666
1667 // NOTE: commenting out the delCol flag at this point. The colsforremoval array should indicate the
1668 // columns that have a zero reduced cost. Hence, the delCol flag should only be set in the isSPxColId
1669 // branch of the if statement.
1670 //delCol = true;
1671 }
1672 }
1673 else if(_solver.basis().baseId(
1674 i).isSPxColId()) // get the column id's for the columns in the basis
1675 {
1676 if(isZero(feasVector[i], feastol))
1677 {
1678 nonposind[*nnonposind] = i;
1679 (*nnonposind)++;
1680
1681 delCol = true;
1682 }
1683 }
1684
1685 // setting an array to identify the columns to be removed from the LP to form the reduced problem
1686 if(delCol)
1687 {
1688 colsforremoval[i] = -1;
1689 _decompReducedProbCols[i] = false;
1690 }
1691 else if(_solver.basis().baseId(i).isSPxColId())
1692 {
1693 if(_solver.basis().baseId(i).isSPxColId())
1694 _decompReducedProbColIDs[_solver.number(_solver.basis().baseId(i))] = SPxColId(
1695 _solver.basis().baseId(i));
1696 else
1697 _decompReducedProbCols[_solver.number(_solver.basis().baseId(i))] = false;
1698 }
1699 }
1700
1701 stop = decompTerminate(realParam(SoPlexBase<R>::TIMELIMIT) * TIMELIMIT_FRAC);
1702 }
1703
1704
1705
1706 /// retrieves the compatible columns from the constraint matrix
1707 // This function also updates the constraint matrix of the reduced problem. It is efficient to perform this in the
1708 // following function because the required linear algebra has been performed.
1709 template <class R>
1710 void SoPlexBase<R>::_getCompatibleColumns(VectorBase<R> feasVector, int* nonposind, int* compatind,
1711 int* rowsforremoval,
1712 int* colsforremoval, int nnonposind, int* ncompatind, bool formRedProb, bool& stop)
1713 {
1714 #ifdef NO_TOL
1715 R feastol = 0.0;
1716 #else
1717 #ifdef USE_FEASTOL
1718 R feastol = realParam(SoPlexBase<R>::FEASTOL);
1719 #else
1720 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO);
1721 #endif
1722 #endif
1723
1724 bool compatible;
1725 SSVectorBase<R> y(_solver.nCols());
1726 y.unSetup();
1727
1728 *ncompatind = 0;
1729
1730 #ifndef NDEBUG
1731 int bind = 0;
1732 bool* activerows = 0;
1733 spx_alloc(activerows, numRows());
1734
1735 for(int i = 0; i < numRows(); ++i)
1736 activerows[i] = false;
1737
1738 for(int i = 0; i < numCols(); ++i)
1739 {
1740 if(_solver.basis().baseId(i).isSPxRowId()) // find the row id's for rows in the basis
1741 {
1742 if(!isZero(feasVector[i], feastol))
1743 {
1744 bind = _realLP->number(SPxRowId(_solver.basis().baseId(i))); // getting the corresponding row
1745 // for the original LP.
1746 assert(bind >= 0 && bind < numRows());
1747 activerows[bind] = true;
1748 }
1749 }
1750 }
1751
1752 #endif
1753
1754 // this function is called from other places where to identify the columns for removal. In these places the
1755 // reduced problem is not formed.
1756 if(formRedProb)
1757 {
1758 _decompReducedProbRowIDs.reSize(_solver.nRows());
1759 _decompReducedProbColRowIDs.reSize(_solver.nRows());
1760 }
1761
1762 for(int i = 0; i < numRows(); ++i)
1763 {
1764 rowsforremoval[i] = i;
1765
1766 if(formRedProb)
1767 _decompReducedProbRows[i] = true;
1768
1769 numIncludedRows++;
1770
1771 // the rhs of this calculation are the rows of the constraint matrix
1772 // so we are solving y B = A_{i,.}
1773 // @NOTE: This may not actually be necessary. The solve process is very time consuming and is a point where the
1774 // approach breaks down. It could be simplier if we use a faster solve. Maybe something like:
1775 // Omer, J.; Towhidi, M. & Soumis, F., "The positive edge pricing rule for the dual simplex",
1776 // Computers & Operations Research , 2015, 61, 135-142
1777 try
1778 {
1779 _solver.basis().solve(y, _solver.vector(i));
1780 }
1781 catch(const SPxException& E)
1782 {
1783 MSG_ERROR(spxout << "Caught exception <" << E.what() << "> while computing compatability.\n");
1784 }
1785
1786 compatible = true;
1787
1788 // a compatible row is given by zeros in all columns related to the nonpositive indices
1789 for(int j = 0; j < nnonposind; ++j)
1790 {
1791 // @TODO: getting a tolerance issue with this check. Don't know how to fix it.
1792 if(!isZero(y[nonposind[j]], feastol))
1793 {
1794 compatible = false;
1795 break;
1796 }
1797 }
1798
1799 // checking that the active rows are compatible
1800 assert(!activerows[i] || compatible);
1801
1802 // changing the matrix coefficients
1803 DSVectorBase<R> newRowVector;
1804 LPRowBase<R> rowtoupdate;
1805
1806 if(y.isSetup())
1807 {
1808 for(int j = 0; j < y.size(); j++)
1809 newRowVector.add(y.index(j), y.value(j));
1810 }
1811 else
1812 {
1813 for(int j = 0; j < numCols(); j++)
1814 {
1815 if(!isZero(y[j], feastol))
1816 newRowVector.add(j, y[j]);
1817 }
1818 }
1819
1820 // transforming the original problem rows
1821 _solver.getRow(i, rowtoupdate);
1822
1823 #ifndef NO_TRANSFORM
1824 rowtoupdate.setRowVector(newRowVector);
1825 #endif
1826
1827 if(formRedProb)
1828 _transformedRows.add(rowtoupdate);
1829
1830
1831 // Making all equality constraints compatible, i.e. they are included in the reduced problem
1832 if(EQ(rowtoupdate.lhs(), rowtoupdate.rhs()))
1833 compatible = true;
1834
1835 if(compatible)
1836 {
1837 compatind[*ncompatind] = i;
1838 (*ncompatind)++;
1839
1840 if(formRedProb)
1841 {
1842 _decompReducedProbRowIDs[i] = _solver.rowId(i);
1843
1844 // updating the compatible row
1845 _decompLP->changeRow(i, rowtoupdate);
1846 }
1847 }
1848 else
1849 {
1850 // setting an array to identify the rows to be removed from the LP to form the reduced problem
1851 rowsforremoval[i] = -1;
1852 numIncludedRows--;
1853
1854 if(formRedProb)
1855 _decompReducedProbRows[i] = false;
1856 }
1857
1858 // determine whether the reduced problem setup should be terminated
1859 stop = decompTerminate(realParam(SoPlexBase<R>::TIMELIMIT) * TIMELIMIT_FRAC);
1860
1861 if(stop)
1862 break;
1863 }
1864
1865 assert(numIncludedRows <= _solver.nRows());
1866
1867 #ifndef NDEBUG
1868 spx_free(activerows);
1869 #endif
1870 }
1871
1872
1873
1874 /// computes the reduced problem objective coefficients
1875 template <class R>
1876 void SoPlexBase<R>::_computeReducedProbObjCoeff(bool& stop)
1877 {
1878 #ifdef NO_TOL
1879 R feastol = 0.0;
1880 #else
1881 #ifdef USE_FEASTOL
1882 R feastol = realParam(SoPlexBase<R>::FEASTOL);
1883 #else
1884 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO);
1885 #endif
1886 #endif
1887
1888 SSVectorBase<R> y(numCols());
1889 y.unSetup();
1890
1891 // the rhs of this calculation is the original objective coefficient vector
1892 // so we are solving y B = c
1893 try
1894 {
1895 _solver.basis().solve(y, _solver.maxObj());
1896 }
1897 catch(const SPxException& E)
1898 {
1899 MSG_ERROR(spxout << "Caught exception <" << E.what() << "> while computing compatability.\n");
1900 }
1901
1902 _transformedObj.reDim(numCols());
1903
1904 if(y.isSetup())
1905 {
1906 int ycount = 0;
1907
1908 for(int i = 0; i < numCols(); i++)
1909 {
1910 if(ycount < y.size() && i == y.index(ycount))
1911 {
1912 _transformedObj[i] = y.value(ycount);
1913 ycount++;
1914 }
1915 else
1916 _transformedObj[i] = 0.0;
1917 }
1918 }
1919 else
1920 {
1921 for(int i = 0; i < numCols(); i++)
1922 {
1923 if(isZero(y[i], feastol))
1924 _transformedObj[i] = 0.0;
1925 else
1926 _transformedObj[i] = y[i];
1927 }
1928 }
1929
1930 // setting the updated objective vector
1931 #ifndef NO_TRANSFORM
1932 _decompLP->changeObj(_transformedObj);
1933 #endif
1934
1935 // determine whether the reduced problem setup should be terminated
1936 stop = decompTerminate(realParam(SoPlexBase<R>::TIMELIMIT) * TIMELIMIT_FRAC);
1937 }
1938
1939
1940
1941 /// computes the compatible bound constraints and adds them to the reduced problem
1942 // NOTE: No columns are getting removed from the reduced problem. Only the bound constraints are being removed.
1943 // So in the reduced problem, all variables are free unless the bound constraints are selected as compatible.
1944 template <class R>
1945 void SoPlexBase<R>::_getCompatibleBoundCons(LPRowSetBase<R>& boundcons, int* compatboundcons,
1946 int* nonposind,
1947 int* ncompatboundcons, int nnonposind, bool& stop)
1948 {
1949 #ifdef NO_TOL
1950 R feastol = 0.0;
1951 #else
1952 #ifdef USE_FEASTOL
1953 R feastol = realParam(SoPlexBase<R>::FEASTOL);
1954 #else
1955 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO);
1956 #endif
1957 #endif
1958
1959 bool compatible;
1960 SSVectorBase<R> y(numCols());
1961 y.unSetup();
1962
1963 _decompReducedProbColRowIDs.reSize(numCols());
1964
1965
1966 // identifying the compatible bound constraints
1967 for(int i = 0; i < numCols(); i++)
1968 {
1969 _decompReducedProbColRowIDs[i].inValidate();
1970
1971 // setting all variables to free variables.
1972 // Bound constraints will only be added to the variables with compatible bound constraints.
1973 _decompLP->changeUpper(i, R(infinity));
1974 _decompLP->changeLower(i, R(-infinity));
1975
1976 // the rhs of this calculation are the unit vectors of the bound constraints
1977 // so we are solving y B = I_{i,.}
1978 // this solve should not be required. We only need the column of the basis inverse.
1979 try
1980 {
1981 _solver.basis().solve(y, _solver.unitVector(i));
1982 }
1983 catch(const SPxException& E)
1984 {
1985 MSG_ERROR(spxout << "Caught exception <" << E.what() << "> while computing compatability.\n");
1986 }
1987
1988 compatible = true;
1989
1990 // a compatible row is given by zeros in all columns related to the nonpositive indices
1991 for(int j = 0; j < nnonposind; ++j)
1992 {
1993 // @todo really need to check this part of the code. Run through this with Ambros or Matthias.
1994 if(!isZero(y[nonposind[j]]))
1995 {
1996 compatible = false;
1997 break;
1998 }
1999 }
2000
2001 // changing the matrix coefficients
2002 DSVectorBase<R> newRowVector;
2003 LPRowBase<R> rowtoupdate;
2004
2005 #ifndef NO_TRANSFORM
2006
2007 if(y.isSetup())
2008 {
2009 for(int j = 0; j < y.size(); j++)
2010 newRowVector.add(y.index(j), y.value(j));
2011 }
2012 else
2013 {
2014 for(int j = 0; j < numCols(); j++)
2015 {
2016 if(!isZero(y[j], feastol))
2017 {
2018 newRowVector.add(j, y[j]);
2019 }
2020 }
2021 }
2022
2023 #else
2024 newRowVector.add(i, 1.0);
2025 #endif
2026
2027 // will probably need to map this set of rows.
2028 _transformedRows.add(_solver.lower(i), newRowVector, _solver.upper(i));
2029
2030 // variable bounds are compatible in the current implementation.
2031 // this is still function is still required to transform the bound constraints with respect to the basis matrix.
2032 compatible = true;
2033
2034 // if the bound constraint is compatible, it remains in the reduced problem.
2035 // Otherwise the bound is removed and the variables are free.
2036 if(compatible)
2037 {
2038 R lhs = R(-infinity);
2039 R rhs = R(infinity);
2040
2041 if(GT(_solver.lower(i), R(-infinity)))
2042 lhs = _solver.lower(i);
2043
2044 if(LT(_solver.upper(i), R(infinity)))
2045 rhs = _solver.upper(i);
2046
2047 if(GT(lhs, R(-infinity)) || LT(rhs, R(infinity)))
2048 {
2049 compatboundcons[(*ncompatboundcons)] = i;
2050 (*ncompatboundcons)++;
2051
2052 boundcons.add(lhs, newRowVector, rhs);
2053 }
2054
2055
2056 }
2057
2058 // determine whether the reduced problem setup should be terminated
2059 stop = decompTerminate(realParam(SoPlexBase<R>::TIMELIMIT) * TIMELIMIT_FRAC);
2060
2061 if(stop)
2062 break;
2063 }
2064
2065 }
2066
2067
2068
2069 /// computes the rows to remove from the complementary problem
2070 template <class R>
2071 void SoPlexBase<R>::_getRowsForRemovalComplementaryProblem(int* nonposind, int* bind,
2072 int* rowsforremoval,
2073 int* nrowsforremoval, int nnonposind)
2074 {
2075 *nrowsforremoval = 0;
2076
2077 for(int i = 0; i < nnonposind; i++)
2078 {
2079 if(bind[nonposind[i]] < 0)
2080 {
2081 rowsforremoval[*nrowsforremoval] = -1 - bind[nonposind[i]];
2082 (*nrowsforremoval)++;
2083 }
2084 }
2085 }
2086
2087
2088
2089 /// removing rows from the complementary problem.
2090 // the rows that are removed from decompCompLP are the rows from the reduced problem that have a non-positive dual
2091 // multiplier in the optimal solution.
2092 template <class R>
2093 void SoPlexBase<R>::_deleteAndUpdateRowsComplementaryProblem(SPxRowId rangedRowIds[],
2094 int& naddedrows)
2095 {
2096 DSVectorBase<R> slackColCoeff;
2097
2098 // setting the objective coefficients of the original variables to zero
2099 VectorBase<R> newObjCoeff(numCols());
2100
2101 for(int i = 0; i < numCols(); i++)
2102 {
2103 _compSolver.changeBounds(_realLP->cId(i), R(-infinity), R(infinity));
2104 newObjCoeff[i] = 0;
2105 }
2106
2107 _compSolver.changeObj(newObjCoeff);
2108
2109 // adding the slack column to the complementary problem
2110 SPxColId* addedcolid = 0;
2111
2112 if(boolParam(SoPlexBase<R>::USECOMPDUAL))
2113 {
2114 spx_alloc(addedcolid, 1);
2115 LPColSetBase<R> compSlackCol;
2116 compSlackCol.add(R(1.0), R(-infinity), slackColCoeff, R(infinity));
2117 _compSolver.addCols(addedcolid, compSlackCol);
2118 _compSlackColId = addedcolid[0];
2119 }
2120 else
2121 {
2122 LPRowSetBase<R>
2123 addrangedrows; // the row set of ranged and equality rows that must be added to the complementary problem.
2124 naddedrows = 0;
2125
2126 // finding all of the ranged and equality rows and creating two <= constraints.
2127 for(int i = 0; i < numRows(); i++)
2128 {
2129 if(_realLP->rowType(i) == LPRowBase<R>::RANGE || _realLP->rowType(i) == LPRowBase<R>::EQUAL)
2130 {
2131 assert(GT(_compSolver.lhs(i), R(-infinity)) && LT(_compSolver.rhs(i), R(infinity)));
2132 assert(_compSolver.rowType(i) == LPRowBase<R>::RANGE
2133 || _compSolver.rowType(i) == LPRowBase<R>::EQUAL);
2134
2135 _compSolver.changeLhs(i, R(-infinity));
2136 addrangedrows.add(_realLP->lhs(i), _realLP->rowVector(i), R(infinity));
2137 naddedrows++;
2138 }
2139 }
2140
2141 // adding the rows for the ranged rows to make <= conatraints
2142 SPxRowId* addedrowid = 0;
2143 spx_alloc(addedrowid, naddedrows);
2144 _compSolver.addRows(addedrowid, addrangedrows);
2145
2146 // collecting the row ids
2147 for(int i = 0; i < naddedrows; i++)
2148 rangedRowIds[i] = addedrowid[i];
2149
2150 spx_free(addedrowid);
2151
2152
2153 // adding the slack column
2154 spx_alloc(addedcolid, 1);
2155 LPColSetBase<R> compSlackCol;
2156 compSlackCol.add(R(-1.0), R(0.0), slackColCoeff, R(infinity));
2157
2158 _compSolver.addCols(addedcolid, compSlackCol);
2159 _compSlackColId = addedcolid[0];
2160 }
2161
2162 // freeing allocated memory
2163 spx_free(addedcolid);
2164 }
2165
2166
2167
2168 /// update the dual complementary problem with additional columns and rows
2169 // Given the solution to the updated reduced problem, the complementary problem will be updated with modifications to
2170 // the constraints and the removal of variables
2171 template <class R>
2172 void SoPlexBase<R>::_updateDecompComplementaryDualProblem(bool origObj)
2173 {
2174 R feastol = realParam(SoPlexBase<R>::FEASTOL);
2175
2176 int prevNumCols =
2177 _compSolver.nCols(); // number of columns in the previous formulation of the complementary prob
2178 int prevNumRows = _compSolver.nRows();
2179 int prevPrimalRowIds = _nPrimalRows;
2180 int prevDualColIds = _nDualCols;
2181
2182 LPColSetBase<R> addElimCols(_nElimPrimalRows); // columns previously eliminated from the
2183 // complementary problem that must be added
2184 int numElimColsAdded = 0;
2185 int currnumrows = prevNumRows;
2186 // looping over all rows from the original LP that were eliminated during the formation of the complementary
2187 // problem. The eliminated rows will be added if they are basic in the reduced problem.
2188
2189 for(int i = 0; i < _nElimPrimalRows; i++)
2190 {
2191 int rowNumber = _realLP->number(_decompElimPrimalRowIDs[i]);
2192
2193 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]);
2194 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows());
2195
2196 // checking for the basic rows in the reduced problem
2197 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER ||
2198 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER ||
2199 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED ||
2200 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE ||
2201 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
2202 LE(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol)) ||
2203 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
2204 LE(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol)))
2205 {
2206 LPRowBase<R> origlprow;
2207 DSVectorBase<R> coltoaddVec(_realLP->nCols());
2208
2209 LPRowSetBase<R> additionalrows;
2210 int nnewrows = 0;
2211
2212 _realLP->getRow(rowNumber, origlprow);
2213
2214 for(int j = 0; j < origlprow.rowVector().size(); j++)
2215 {
2216 // the column of the new row may not exist in the current complementary problem.
2217 // if the column does not exist, then it is necessary to create the column.
2218 int colNumber = origlprow.rowVector().index(j);
2219
2220 if(_decompCompProbColIDsIdx[colNumber] == -1)
2221 {
2222 assert(!_decompReducedProbColIDs[colNumber].isValid());
2223 _decompPrimalColIDs[_nPrimalCols] = _realLP->cId(colNumber);
2224 _decompCompProbColIDsIdx[colNumber] = _nPrimalCols;
2225 _fixedOrigVars[colNumber] = -2;
2226 _nPrimalCols++;
2227
2228 // all columns for the complementary problem are converted to unrestricted.
2229 additionalrows.create(1, _realLP->maxObj(colNumber), _realLP->maxObj(colNumber));
2230 nnewrows++;
2231
2232 coltoaddVec.add(currnumrows, origlprow.rowVector().value(j));
2233 currnumrows++;
2234 }
2235 else
2236 coltoaddVec.add(_compSolver.number(_decompDualRowIDs[_decompCompProbColIDsIdx[colNumber]]),
2237 origlprow.rowVector().value(j));
2238 }
2239
2240 SPxRowId* addedrowids = 0;
2241 spx_alloc(addedrowids, nnewrows);
2242 _compSolver.addRows(addedrowids, additionalrows);
2243
2244 for(int j = 0; j < nnewrows; j++)
2245 {
2246 _decompDualRowIDs[_nDualRows] = addedrowids[j];
2247 _nDualRows++;
2248 }
2249
2250 spx_free(addedrowids);
2251
2252
2253
2254 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER ||
2255 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED ||
2256 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE ||
2257 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
2258 LE(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol)))
2259 {
2260 assert(LT(_realLP->rhs(_decompElimPrimalRowIDs[i]), R(infinity)));
2261 addElimCols.add(_realLP->rhs(_decompElimPrimalRowIDs[i]), R(-infinity), coltoaddVec, R(infinity));
2262
2263 if(_nPrimalRows >= _decompPrimalRowIDs.size())
2264 {
2265 _decompPrimalRowIDs.reSize(_nPrimalRows * 2);
2266 _decompDualColIDs.reSize(_nPrimalRows * 2);
2267 }
2268
2269 _decompPrimalRowIDs[_nPrimalRows] = _decompElimPrimalRowIDs[i];
2270 _nPrimalRows++;
2271
2272 _decompElimPrimalRowIDs.remove(i);
2273 _nElimPrimalRows--;
2274 i--;
2275
2276 numElimColsAdded++;
2277 }
2278 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER ||
2279 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
2280 LE(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol)))
2281 {
2282 // this assert should stay, but there is an issue with the status and the dual vector
2283 //assert(LT(dualVector[_solver.number(_decompReducedProbRowIDs[rowNumber])], 0.0));
2284 assert(GT(_realLP->lhs(_decompElimPrimalRowIDs[i]), R(-infinity)));
2285 addElimCols.add(_realLP->lhs(_decompElimPrimalRowIDs[i]), R(-infinity), coltoaddVec, R(infinity));
2286
2287 _decompPrimalRowIDs[_nPrimalRows] = _decompElimPrimalRowIDs[i];
2288 _nPrimalRows++;
2289
2290 _decompElimPrimalRowIDs.remove(i);
2291 _nElimPrimalRows--;
2292 i--;
2293
2294 numElimColsAdded++;
2295 }
2296 }
2297 }
2298
2299 MSG_INFO2(spxout, spxout << "Number of eliminated columns added to complementary problem: "
2300 << numElimColsAdded << std::endl);
2301
2302 // updating the _decompDualColIDs with the additional columns from the eliminated rows.
2303 _compSolver.addCols(addElimCols);
2304
2305 for(int i = prevNumCols; i < _compSolver.nCols(); i++)
2306 {
2307 _decompDualColIDs[prevDualColIds + i - prevNumCols] = _compSolver.colId(i);
2308 _nDualCols++;
2309 }
2310
2311 assert(_nDualCols == _nPrimalRows);
2312
2313 // looping over all rows from the original problem that were originally contained in the complementary problem.
2314 // The basic rows will be set as free variables, the non-basic rows will be eliminated from the complementary
2315 // problem.
2316 DSVectorBase<R> slackRowCoeff(_compSolver.nCols());
2317
2318 int* colsforremoval = 0;
2319 int ncolsforremoval = 0;
2320 spx_alloc(colsforremoval, prevPrimalRowIds);
2321
2322 for(int i = 0; i < prevPrimalRowIds; i++)
2323 {
2324 int rowNumber = _realLP->number(_decompPrimalRowIDs[i]);
2325
2326 // this loop runs over all rows previously in the complementary problem. If rows are added to the reduced
2327 // problem, they will be transfered from the incompatible set to the compatible set in the following if
2328 // statement.
2329 if(_decompReducedProbRows[rowNumber])
2330 {
2331 // rows added to the reduced problem may have been equality constriants. The equality constraints from the
2332 // original problem are converted into <= and >= constraints. Upon adding these constraints to the reduced
2333 // problem, only a single dual column is needed in the complementary problem. Hence, one of the dual columns
2334 // is removed.
2335 //
2336 // 22.06.2015 Testing keeping all constraints in the complementary problem.
2337 // This requires a dual column to be fixed to zero for the range and equality rows.
2338 #ifdef KEEP_ALL_ROWS_IN_COMP_PROB
2339 bool incrementI = false;
2340 #endif
2341
2342 if(i + 1 < prevPrimalRowIds
2343 && _realLP->number(_decompPrimalRowIDs[i]) == _realLP->number(_decompPrimalRowIDs[i + 1]))
2344 {
2345 assert(_decompPrimalRowIDs[i].idx == _decompPrimalRowIDs[i + 1].idx);
2346
2347 #ifdef KEEP_ALL_ROWS_IN_COMP_PROB // 22.06.2015
2348
2349 if(_realLP->rowType(_decompPrimalRowIDs[i]) == LPRowBase<R>::RANGE)
2350 {
2351 _compSolver.changeObj(_decompDualColIDs[i + 1], 0.0);
2352 _compSolver.changeBounds(_decompDualColIDs[i + 1], 0.0, 0.0);
2353 }
2354
2355 incrementI = true;
2356 #else
2357 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompDualColIDs[i + 1]));
2358 ncolsforremoval++;
2359
2360 _decompPrimalRowIDs.remove(i + 1);
2361 _nPrimalRows--;
2362 _decompDualColIDs.remove(i + 1);
2363 _nDualCols--;
2364
2365 prevPrimalRowIds--;
2366 #endif
2367 }
2368
2369 //assert(i + 1 == prevPrimalRowIds || _decompPrimalRowIDs[i].idx != _decompPrimalRowIDs[i+1].idx);
2370
2371 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]);
2372 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows());
2373
2374 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER ||
2375 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED ||
2376 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE ||
2377 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
2378 LE(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol)))
2379 {
2380 //assert(GT(dualVector[solverRowNum], 0.0));
2381 _compSolver.changeObj(_decompDualColIDs[i], _realLP->rhs(SPxRowId(_decompPrimalRowIDs[i])));
2382 _compSolver.changeBounds(_decompDualColIDs[i], R(-infinity), R(infinity));
2383 }
2384 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER ||
2385 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
2386 LE(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol)))
2387 {
2388 //assert(LT(dualVector[solverRowNum], 0.0));
2389 _compSolver.changeObj(_decompDualColIDs[i], _realLP->lhs(SPxRowId(_decompPrimalRowIDs[i])));
2390 _compSolver.changeBounds(_decompDualColIDs[i], R(-infinity), R(infinity));
2391 }
2392 else //if ( _solver.basis().desc().rowStatus(solverRowNum) != SPxBasisBase<R>::Desc::D_FREE )
2393 {
2394 //assert(isZero(dualVector[solverRowNum], 0.0));
2395
2396 // 22.06.2015 Testing keeping all rows in the complementary problem
2397 #ifdef KEEP_ALL_ROWS_IN_COMP_PROB
2398 switch(_realLP->rowType(_decompPrimalRowIDs[i]))
2399 {
2400 case LPRowBase<R>::RANGE:
2401 assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) ==
2402 _realLP->number(SPxColId(_decompPrimalRowIDs[i + 1])));
2403 assert(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) !=
2404 _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))));
2405
2406 _compSolver.changeObj(_decompDualColIDs[i], _realLP->rhs(_decompPrimalRowIDs[i]));
2407 _compSolver.changeBounds(_decompDualColIDs[i], 0.0, R(infinity));
2408 _compSolver.changeObj(_decompDualColIDs[i + 1], _realLP->lhs(_decompPrimalRowIDs[i]));
2409 _compSolver.changeBounds(_decompDualColIDs[i + 1], R(-infinity), 0.0);
2410
2411 i++;
2412 break;
2413
2414 case LPRowBase<R>::GREATER_EQUAL:
2415 _compSolver.changeObj(_decompDualColIDs[i], _realLP->lhs(_decompPrimalRowIDs[i]));
2416 _compSolver.changeBounds(_decompDualColIDs[i], R(-infinity), 0.0);
2417 break;
2418
2419 case LPRowBase<R>::LESS_EQUAL:
2420 _compSolver.changeObj(_decompDualColIDs[i], _realLP->rhs(_decompPrimalRowIDs[i]));
2421 _compSolver.changeBounds(_decompDualColIDs[i], 0.0, R(infinity));
2422 break;
2423
2424 default:
2425 throw SPxInternalCodeException("XDECOMPSL01 This should never happen.");
2426 }
2427
2428 #else // 22.06.2015 testing keeping all rows in the complementary problem
2429 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompDualColIDs[i]));
2430 ncolsforremoval++;
2431
2432 if(_nElimPrimalRows >= _decompElimPrimalRowIDs.size())
2433 _decompElimPrimalRowIDs.reSize(_realLP->nRows());
2434
2435 _decompElimPrimalRowIDs[_nElimPrimalRows] = _decompPrimalRowIDs[i];
2436 _nElimPrimalRows++;
2437 _decompPrimalRowIDs.remove(i);
2438 _nPrimalRows--;
2439 _decompDualColIDs.remove(i);
2440 _nDualCols--;
2441
2442 i--;
2443 prevPrimalRowIds--;
2444 #endif
2445 }
2446
2447 #ifdef KEEP_ALL_ROWS_IN_COMP_PROB
2448
2449 if(incrementI)
2450 i++;
2451
2452 #endif
2453 }
2454 else
2455 {
2456 switch(_realLP->rowType(_decompPrimalRowIDs[i]))
2457 {
2458 case LPRowBase<R>::RANGE:
2459 assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) ==
2460 _realLP->number(SPxColId(_decompPrimalRowIDs[i + 1])));
2461 assert(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) !=
2462 _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))));
2463
2464 if(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) <
2465 _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))))
2466 {
2467 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), -SLACKCOEFF);
2468 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])), SLACKCOEFF);
2469 }
2470 else
2471 {
2472 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), SLACKCOEFF);
2473 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])), -SLACKCOEFF);
2474 }
2475
2476
2477 i++;
2478 break;
2479
2480 case LPRowBase<R>::EQUAL:
2481 assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) ==
2482 _realLP->number(SPxColId(_decompPrimalRowIDs[i + 1])));
2483
2484 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), SLACKCOEFF);
2485 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])), SLACKCOEFF);
2486
2487 i++;
2488 break;
2489
2490 case LPRowBase<R>::GREATER_EQUAL:
2491 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), -SLACKCOEFF);
2492 break;
2493
2494 case LPRowBase<R>::LESS_EQUAL:
2495 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), SLACKCOEFF);
2496 break;
2497
2498 default:
2499 throw SPxInternalCodeException("XDECOMPSL01 This should never happen.");
2500 }
2501
2502 if(origObj)
2503 {
2504 int numRemove = 1;
2505 int removeCount = 0;
2506
2507 if(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) ==
2508 _realLP->number(SPxColId(_decompPrimalRowIDs[i + 1])))
2509 numRemove++;
2510
2511 do
2512 {
2513 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompDualColIDs[i]));
2514 ncolsforremoval++;
2515
2516 if(_nElimPrimalRows >= _decompElimPrimalRowIDs.size())
2517 _decompElimPrimalRowIDs.reSize(_realLP->nRows());
2518
2519 _decompElimPrimalRowIDs[_nElimPrimalRows] = _decompPrimalRowIDs[i];
2520 _nElimPrimalRows++;
2521 _decompPrimalRowIDs.remove(i);
2522 _nPrimalRows--;
2523 _decompDualColIDs.remove(i);
2524 _nDualCols--;
2525
2526 i--;
2527 prevPrimalRowIds--;
2528
2529 removeCount++;
2530 }
2531 while(removeCount < numRemove);
2532 }
2533 }
2534 }
2535
2536 // updating the slack column in the complementary problem
2537 R lhs = 1.0;
2538 R rhs = 1.0;
2539
2540 // it is possible that all rows are included in the reduced problem. In this case, the slack row will be empty. To
2541 // avoid infeasibility, the lhs and rhs are set to 0.
2542 if(slackRowCoeff.size() == 0)
2543 {
2544 lhs = 0.0;
2545 rhs = 0.0;
2546 }
2547
2548 LPRowBase<R> compSlackRow(lhs, slackRowCoeff, rhs);
2549 _compSolver.changeRow(_compSlackDualRowId, compSlackRow);
2550
2551
2552 // if the original objective is used, then all dual columns related to primal rows not in the reduced problem are
2553 // removed from the complementary problem.
2554 // As a result, the slack row becomes an empty row.
2555 int* perm = 0;
2556 spx_alloc(perm, _compSolver.nCols() + numElimColsAdded);
2557 _compSolver.removeCols(colsforremoval, ncolsforremoval, perm);
2558
2559
2560 // updating the dual columns to represent the fixed primal variables.
2561 int* currFixedVars = 0;
2562 spx_alloc(currFixedVars, _realLP->nCols());
2563 _identifyComplementaryDualFixedPrimalVars(currFixedVars);
2564 _removeComplementaryDualFixedPrimalVars(currFixedVars);
2565 _updateComplementaryDualFixedPrimalVars(currFixedVars);
2566
2567
2568 // freeing allocated memory
2569 spx_free(currFixedVars);
2570 spx_free(perm);
2571 spx_free(colsforremoval);
2572 }
2573
2574
2575
2576 /// update the primal complementary problem with additional columns and rows
2577 // Given the solution to the updated reduced problem, the complementary problem will be updated with modifications to
2578 // the constraints and the removal of variables
2579 template <class R>
2580 void SoPlexBase<R>::_updateDecompComplementaryPrimalProblem(bool origObj)
2581 {
2582 R feastol = realParam(SoPlexBase<R>::FEASTOL);
2583
2584 int prevNumRows = _compSolver.nRows();
2585 int prevPrimalRowIds = _nPrimalRows;
2586
2587 assert(_nPrimalRows == _nCompPrimalRows);
2588
2589 LPRowSetBase<R> addElimRows(_nElimPrimalRows); // rows previously eliminated from the
2590 // complementary problem that must be added
2591 int numElimRowsAdded = 0;
2592 // looping over all rows from the original LP that were eliminated during the formation of the complementary
2593 // problem. The eliminated rows will be added if they are basic in the reduced problem.
2594
2595 for(int i = 0; i < _nElimPrimalRows; i++)
2596 {
2597 int rowNumber = _realLP->number(_decompElimPrimalRowIDs[i]);
2598
2599 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]);
2600 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows());
2601
2602 // checking the rows that are basic in the reduced problem that should be added to the complementary problem
2603 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER
2604 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER
2605 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED
2606 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE
2607 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
2608 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol))
2609 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
2610 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol)))
2611 {
2612 LPRowBase<R> origlprow;
2613 _realLP->getRow(rowNumber, origlprow);
2614
2615 // NOTE: 11.02.2016 I am assuming that all columns from the original problem are contained in the
2616 // complementary problem. Will need to check this. Since nothing is happenning in the
2617 // _deleteAndUpdateRowsComplementaryProblem function, I am feeling confident that all columns remain.
2618
2619
2620 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER
2621 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED
2622 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE
2623 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
2624 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol)))
2625 {
2626 assert(LT(_realLP->rhs(_decompElimPrimalRowIDs[i]), R(infinity)));
2627
2628 if(_nPrimalRows >= _decompPrimalRowIDs.size())
2629 {
2630 _decompPrimalRowIDs.reSize(_nPrimalRows * 2);
2631 _decompCompPrimalRowIDs.reSize(_nPrimalRows * 2);
2632 }
2633
2634 addElimRows.add(_realLP->rhs(_decompElimPrimalRowIDs[i]), origlprow.rowVector(),
2635 _realLP->rhs(_decompElimPrimalRowIDs[i]));
2636
2637 _decompPrimalRowIDs[_nPrimalRows] = _decompElimPrimalRowIDs[i];
2638 _nPrimalRows++;
2639
2640 _decompElimPrimalRowIDs.remove(i);
2641 _nElimPrimalRows--;
2642 i--;
2643
2644 numElimRowsAdded++;
2645 }
2646 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER
2647 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
2648 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol)))
2649 {
2650 assert(GT(_realLP->lhs(_decompElimPrimalRowIDs[i]), R(-infinity)));
2651
2652 if(_nPrimalRows >= _decompPrimalRowIDs.size())
2653 {
2654 _decompPrimalRowIDs.reSize(_nPrimalRows * 2);
2655 _decompCompPrimalRowIDs.reSize(_nPrimalRows * 2);
2656 }
2657
2658 addElimRows.add(_realLP->lhs(_decompElimPrimalRowIDs[i]), origlprow.rowVector(),
2659 _realLP->lhs(_decompElimPrimalRowIDs[i]));
2660
2661 _decompPrimalRowIDs[_nPrimalRows] = _decompElimPrimalRowIDs[i];
2662 _nPrimalRows++;
2663
2664 _decompElimPrimalRowIDs.remove(i);
2665 _nElimPrimalRows--;
2666 i--;
2667
2668 numElimRowsAdded++;
2669 }
2670 }
2671 }
2672
2673 MSG_INFO2(spxout, spxout << "Number of eliminated rows added to the complementary problem: "
2674 << numElimRowsAdded << std::endl);
2675
2676 // adding the eliminated rows to the complementary problem.
2677 _compSolver.addRows(addElimRows);
2678
2679 for(int i = prevNumRows; i < _compSolver.nRows(); i++)
2680 {
2681 _decompCompPrimalRowIDs[prevPrimalRowIds + i - prevNumRows] = _compSolver.rowId(i);
2682 _nCompPrimalRows++;
2683 }
2684
2685 assert(_nPrimalRows == _nCompPrimalRows);
2686
2687
2688 // looping over all rows from the original problem that were originally contained in the complementary problem.
2689 // The basic rows will be set as equalities, the non-basic rows will be eliminated from the complementary
2690 // problem.
2691 DSVectorBase<R> slackColCoeff(_compSolver.nRows());
2692
2693 int* rowsforremoval = 0;
2694 int nrowsforremoval = 0;
2695 spx_alloc(rowsforremoval, prevPrimalRowIds);
2696
2697 for(int i = 0; i < prevPrimalRowIds; i++)
2698 {
2699 int rowNumber = _realLP->number(_decompPrimalRowIDs[i]);
2700
2701 // this loop runs over all rows previously in the complementary problem. If rows are added to the reduced
2702 // problem, they will be transfered from the incompatible set to the compatible set in the following if
2703 // statement.
2704 if(_decompReducedProbRows[rowNumber])
2705 {
2706 // rows added to the reduced problem may have been equality constriants. The equality constraints from the
2707 // original problem are converted into <= and >= constraints. Upon adding these constraints to the reduced
2708 // problem, only a single dual column is needed in the complementary problem. Hence, one of the dual columns
2709 // is removed.
2710 //
2711
2712 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]);
2713 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows());
2714
2715 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER
2716 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED
2717 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE
2718 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
2719 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol)))
2720 {
2721 _compSolver.changeLhs(_decompCompPrimalRowIDs[i], _realLP->rhs(SPxRowId(_decompPrimalRowIDs[i])));
2722 // need to also update the RHS because a ranged row could have previously been fixed to LOWER
2723 _compSolver.changeRhs(_decompCompPrimalRowIDs[i], _realLP->rhs(SPxRowId(_decompPrimalRowIDs[i])));
2724 }
2725 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER
2726 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
2727 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol)))
2728 {
2729 _compSolver.changeRhs(_decompCompPrimalRowIDs[i], _realLP->lhs(SPxRowId(_decompPrimalRowIDs[i])));
2730 // need to also update the LHS because a ranged row could have previously been fixed to UPPER
2731 _compSolver.changeLhs(_decompCompPrimalRowIDs[i], _realLP->lhs(SPxRowId(_decompPrimalRowIDs[i])));
2732 }
2733 else //if ( _solver.basis().desc().rowStatus(solverRowNum) != SPxBasisBase<R>::Desc::D_FREE )
2734 {
2735 rowsforremoval[nrowsforremoval] = _compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i]));
2736 nrowsforremoval++;
2737
2738 if(_nElimPrimalRows >= _decompElimPrimalRowIDs.size())
2739 _decompElimPrimalRowIDs.reSize(_realLP->nRows());
2740
2741 _decompElimPrimalRowIDs[_nElimPrimalRows] = _decompPrimalRowIDs[i];
2742 _nElimPrimalRows++;
2743 _decompPrimalRowIDs.remove(i);
2744 _nPrimalRows--;
2745 _decompCompPrimalRowIDs.remove(i);
2746 _nCompPrimalRows--;
2747
2748 i--;
2749 prevPrimalRowIds--;
2750 }
2751 }
2752 else
2753 {
2754 switch(_compSolver.rowType(_decompCompPrimalRowIDs[i]))
2755 {
2756 case LPRowBase<R>::RANGE:
2757 assert(false);
2758 break;
2759
2760 case LPRowBase<R>::EQUAL:
2761 assert(false);
2762 break;
2763
2764 case LPRowBase<R>::LESS_EQUAL:
2765 slackColCoeff.add(_compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i])), -SLACKCOEFF);
2766 break;
2767
2768 case LPRowBase<R>::GREATER_EQUAL:
2769 slackColCoeff.add(_compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i])), SLACKCOEFF);
2770 break;
2771
2772 default:
2773 throw SPxInternalCodeException("XDECOMPSL01 This should never happen.");
2774 }
2775
2776 // this is used as a check at the end of the algorithm. If the original objective function is used, then we
2777 // need to remove all unfixed variables.
2778 if(origObj)
2779 {
2780 rowsforremoval[nrowsforremoval] = _compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i]));
2781 nrowsforremoval++;
2782
2783 if(_nElimPrimalRows >= _decompElimPrimalRowIDs.size())
2784 _decompElimPrimalRowIDs.reSize(_realLP->nRows());
2785
2786 _decompElimPrimalRowIDs[_nElimPrimalRows] = _decompPrimalRowIDs[i];
2787 _nElimPrimalRows++;
2788 _decompPrimalRowIDs.remove(i);
2789 _nPrimalRows--;
2790 _decompCompPrimalRowIDs.remove(i);
2791 _nCompPrimalRows--;
2792
2793 i--;
2794 prevPrimalRowIds--;
2795 }
2796 }
2797 }
2798
2799 // updating the slack column in the complementary problem
2800 LPColBase<R> compSlackCol(-1, slackColCoeff, R(infinity), 0.0);
2801 _compSolver.changeCol(_compSlackColId, compSlackCol);
2802
2803
2804 // if the original objective is used, then all complementary rows related to primal rows not in the reduced problem are
2805 // removed from the complementary problem.
2806 // As a result, the slack column becomes an empty column.
2807 int* perm = 0;
2808 spx_alloc(perm, _compSolver.nRows() + numElimRowsAdded);
2809 _compSolver.removeRows(rowsforremoval, nrowsforremoval, perm);
2810
2811 // updating the dual columns to represent the fixed primal variables.
2812 int* currFixedVars = 0;
2813 spx_alloc(currFixedVars, _realLP->nCols());
2814 _identifyComplementaryPrimalFixedPrimalVars(currFixedVars);
2815 _updateComplementaryPrimalFixedPrimalVars(currFixedVars);
2816
2817
2818 // freeing allocated memory
2819 spx_free(currFixedVars);
2820 spx_free(perm);
2821 spx_free(rowsforremoval);
2822 }
2823
2824
2825
2826 /// checking the optimality of the original problem.
2827 // this function is called if the complementary problem is solved with a non-negative objective value. This implies
2828 // that the rows currently included in the reduced problem are sufficient to identify the optimal solution to the
2829 // original problem.
2830 template <class R>
2831 void SoPlexBase<R>::_checkOriginalProblemOptimality(VectorBase<R> primalVector, bool printViol)
2832 {
2833 SSVectorBase<R> x(_solver.nCols());
2834 x.unSetup();
2835
2836 // multiplying the solution vector of the reduced problem with the transformed basis to identify the original
2837 // solution vector.
2838 _decompTransBasis.coSolve(x, primalVector);
2839
2840 if(printViol)
2841 {
2842 MSG_INFO1(spxout, spxout << std::endl
2843 << "Checking consistency between the reduced problem and the original problem." << std::endl);
2844 }
2845
2846
2847 // checking the objective function values of the reduced problem and the original problem.
2848 R redObjVal = 0;
2849 R objectiveVal = 0;
2850
2851 for(int i = 0; i < _solver.nCols(); i++)
2852 {
2853 redObjVal += _solver.maxObj(i) * primalVector[i];
2854 objectiveVal += _realLP->maxObj(i) * x[i];
2855 }
2856
2857 if(printViol)
2858 {
2859 MSG_INFO1(spxout, spxout << "Reduced Problem Objective Value: " << redObjVal << std::endl
2860 << "Original Problem Objective Value: " << objectiveVal << std::endl);
2861 }
2862
2863 _solReal._isPrimalFeasible = true;
2864 _hasSolReal = true;
2865 // get the primal solutions from the reduced problem
2866 _solReal._primal.reDim(_solver.nCols());
2867 _solReal._primal = x;
2868
2869 R maxviol = 0;
2870 R sumviol = 0;
2871
2872 // checking the bound violations
2873 if(getDecompBoundViolation(maxviol, sumviol))
2874 {
2875 if(printViol)
2876 MSG_INFO1(spxout, spxout << "Bound violation - "
2877 << "Max violation: " << maxviol << " Sum violation: " << sumviol << std::endl);
2878 }
2879
2880 _statistics->totalBoundViol = sumviol;
2881 _statistics->maxBoundViol = maxviol;
2882
2883 // checking the row violations
2884 if(getDecompRowViolation(maxviol, sumviol))
2885 {
2886 if(printViol)
2887 MSG_INFO1(spxout, spxout << "Row violation - "
2888 << "Max violation: " << maxviol << " Sum violation: " << sumviol << std::endl);
2889 }
2890
2891 _statistics->totalRowViol = sumviol;
2892 _statistics->maxRowViol = maxviol;
2893
2894 if(printViol)
2895 MSG_INFO1(spxout, spxout << std::endl);
2896 }
2897
2898
2899
2900 /// updating the slack column coefficients to adjust for equality constraints
2901 template <class R>
2902 void SoPlexBase<R>::_updateComplementaryDualSlackColCoeff()
2903 {
2904 // the slack column for the equality constraints is not handled correctly in the dual conversion. Hence, it is
2905 // necessary to change the equality coefficients of the dual row related to the slack column.
2906 for(int i = 0; i < _nPrimalRows; i++)
2907 {
2908 int rowNumber = _realLP->number(SPxRowId(_decompPrimalRowIDs[i]));
2909
2910 if(!_decompReducedProbRows[rowNumber])
2911 {
2912 if(_realLP->rowType(_decompPrimalRowIDs[i]) == LPRowBase<R>::EQUAL)
2913 {
2914 assert(_realLP->lhs(_decompPrimalRowIDs[i]) == _realLP->rhs(_decompPrimalRowIDs[i]));
2915 _compSolver.changeLower(_decompDualColIDs[i],
2916 0.0); // setting the lower bound of the dual column to zero.
2917
2918 LPColBase<R> addEqualityCol(-_realLP->rhs(_decompPrimalRowIDs[i]),
2919 R(-1.0)*_compSolver.colVector(_decompDualColIDs[i]), R(infinity),
2920 0.0); // adding a new column to the dual
2921
2922 SPxColId newDualCol;
2923 _compSolver.addCol(newDualCol, addEqualityCol);
2924
2925 // inserting the row and col ids for the added column. This is to be next to the original column that has
2926 // been duplicated.
2927 _decompPrimalRowIDs.insert(i + 1, 1, _decompPrimalRowIDs[i]);
2928 _decompDualColIDs.insert(i + 1, 1, newDualCol);
2929 assert(_realLP->number(_decompPrimalRowIDs[i]) == _realLP->number(_decompPrimalRowIDs[i + 1]));
2930
2931 i++;
2932 _nPrimalRows++;
2933 _nDualCols++;
2934 }
2935 }
2936 }
2937 }
2938
2939
2940
2941 /// identify the dual columns related to the fixed variables
2942 template <class R>
2943 void SoPlexBase<R>::_identifyComplementaryDualFixedPrimalVars(int* currFixedVars)
2944 {
2945 R feastol = realParam(SoPlexBase<R>::FEASTOL);
2946
2947 int numFixedVar = 0;
2948
2949 for(int i = 0; i < _realLP->nCols(); i++)
2950 {
2951 currFixedVars[i] = 0;
2952
2953 if(!_decompReducedProbColRowIDs[i].isValid())
2954 continue;
2955
2956 int rowNumber = _solver.number(_decompReducedProbColRowIDs[i]);
2957
2958 if(_decompReducedProbColRowIDs[i].isValid())
2959 {
2960 if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_UPPER ||
2961 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_LOWER ||
2962 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FIXED ||
2963 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_FREE)
2964 {
2965 // setting the value of the _fixedOrigVars array to indicate which variables are at their bounds.
2966 currFixedVars[i] = getOrigVarFixedDirection(i);
2967
2968 numFixedVar++;
2969 }
2970 else
2971 {
2972 // the dual flags do not imply anything about the primal status of the rows.
2973 if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
2974 EQ(_solver.rhs(rowNumber) - _solver.pVec()[rowNumber], R(0.0), feastol))
2975 currFixedVars[i] = 1;
2976 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
2977 EQ(_solver.pVec()[rowNumber] - _solver.lhs(rowNumber), R(0.0), feastol))
2978 currFixedVars[i] = -1;
2979 }
2980 }
2981 }
2982
2983 MSG_INFO3(spxout, spxout << "Number of fixed primal variables in the complementary (dual) problem: "
2984 << numFixedVar << std::endl);
2985 }
2986
2987
2988
2989 /// removing the dual columns related to the fixed variables
2990 template <class R>
2991 void SoPlexBase<R>::_removeComplementaryDualFixedPrimalVars(int* currFixedVars)
2992 {
2993 SPxColId tempId;
2994 int ncolsforremoval = 0;
2995 int* colsforremoval = 0;
2996 spx_alloc(colsforremoval, _realLP->nCols() * 2);
2997
2998 tempId.inValidate();
2999
3000 for(int i = 0; i < _realLP->nCols(); i++)
3001 {
3002 assert(_decompCompProbColIDsIdx[i] != -1); // this should be true in the current implementation
3003
3004 if(_decompCompProbColIDsIdx[i] != -1
3005 && _fixedOrigVars[i] != -2) //&& _fixedOrigVars[i] != currFixedVars[i] )
3006 {
3007 if(_fixedOrigVars[i] != 0)
3008 {
3009 assert(_compSolver.number(SPxColId(_decompFixedVarDualIDs[i])) >= 0);
3010 assert(_fixedOrigVars[i] == -1 || _fixedOrigVars[i] == 1);
3011 assert(_decompFixedVarDualIDs[i].isValid());
3012
3013 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompFixedVarDualIDs[i]));
3014 ncolsforremoval++;
3015
3016 _decompFixedVarDualIDs[i] = tempId;
3017 }
3018 else //if( false && !_decompReducedProbColRowIDs[i].isValid() ) // we want to remove all valid columns
3019 // in the current implementation, the only columns not included in the reduced problem are free columns.
3020 {
3021 assert((LE(_realLP->lower(i), R(-infinity)) && GE(_realLP->upper(i), R(infinity))) ||
3022 _compSolver.number(SPxColId(_decompVarBoundDualIDs[i * 2])) >= 0);
3023 int varcount = 0;
3024
3025 if(GT(_realLP->lower(i), R(-infinity)))
3026 {
3027 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompVarBoundDualIDs[i * 2 +
3028 varcount]));
3029 ncolsforremoval++;
3030
3031 _decompVarBoundDualIDs[i * 2 + varcount] = tempId;
3032 varcount++;
3033 }
3034
3035 if(LT(_realLP->upper(i), R(infinity)))
3036 {
3037 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompVarBoundDualIDs[i * 2 +
3038 varcount]));
3039 ncolsforremoval++;
3040
3041 _decompVarBoundDualIDs[i * 2 + varcount] = tempId;
3042 }
3043
3044 }
3045 }
3046 }
3047
3048 int* perm = 0;
3049 spx_alloc(perm, _compSolver.nCols());
3050 _compSolver.removeCols(colsforremoval, ncolsforremoval, perm);
3051
3052 // freeing allocated memory
3053 spx_free(perm);
3054 spx_free(colsforremoval);
3055 }
3056
3057
3058
3059 /// updating the dual columns related to the fixed primal variables.
3060 template <class R>
3061 void SoPlexBase<R>::_updateComplementaryDualFixedPrimalVars(int* currFixedVars)
3062 {
3063 DSVectorBase<R> col(1);
3064 LPColSetBase<R> boundConsCols;
3065 LPColSetBase<R> fixedVarsDualCols(_nPrimalCols);
3066 int numFixedVars = 0;
3067 // the solution to the reduced problem results in a number of variables at their bounds. If such variables exist
3068 // it is necessary to include a dual column to the complementary problem related to a variable fixing. This is
3069 // equivalent to the tight constraints being converted to equality constraints.
3070 int numBoundConsCols = 0;
3071 int* boundConsColsAdded = 0;
3072 spx_alloc(boundConsColsAdded, _realLP->nCols());
3073
3074 // NOTE: this loop only goes over the primal columns that are included in the complementary problem, i.e. the
3075 // columns from the original problem.
3076 // 29.04.15 in the current implementation, all bound constraints are included in the reduced problem. So, all
3077 // variables (columns) are included in the reduced problem.
3078 for(int i = 0; i < _realLP->nCols(); i++)
3079 {
3080 boundConsColsAdded[i] = 0;
3081 assert(_decompCompProbColIDsIdx[i] != -1);
3082
3083 if(_decompCompProbColIDsIdx[i] != -1)
3084 {
3085 int idIndex = _decompCompProbColIDsIdx[i];
3086 assert(_compSolver.number(SPxRowId(_decompDualRowIDs[idIndex])) >= 0);
3087 col.add(_compSolver.number(SPxRowId(_decompDualRowIDs[idIndex])), 1.0);
3088
3089 if(currFixedVars[i] != 0)
3090 {
3091 assert(currFixedVars[i] == -1 || currFixedVars[i] == 1);
3092
3093 assert(_realLP->lower(i) == _solver.lhs(_decompReducedProbColRowIDs[i]));
3094 assert(_realLP->upper(i) == _solver.rhs(_decompReducedProbColRowIDs[i]));
3095 R colObjCoeff = 0;
3096
3097 if(currFixedVars[i] == -1)
3098 colObjCoeff = _solver.lhs(_decompReducedProbColRowIDs[i]);
3099 else
3100 colObjCoeff = _solver.rhs(_decompReducedProbColRowIDs[i]);
3101
3102 fixedVarsDualCols.add(colObjCoeff, R(-infinity), col, R(infinity));
3103 numFixedVars++;
3104 }
3105 // 09.02.15 I think that the else should only be entered if the column does not exist in the reduced
3106 // prob. I have tested by just leaving this as an else (without the if), but I think that this is wrong.
3107 //else if( !_decompReducedProbColRowIDs[i].isValid() )
3108
3109 // NOTE: in the current implementation all columns, except the free columns, are included int the reduced
3110 // problem. There in no need to include the variable bounds in the complementary problem.
3111 else //if( false && _fixedOrigVars[i] == -2 )
3112 {
3113 bool isRedProbCol = _decompReducedProbColRowIDs[i].isValid();
3114
3115 // 29.04.15 in the current implementation only free variables are not included in the reduced problem
3116 if(GT(_realLP->lower(i), R(-infinity)))
3117 {
3118 if(!isRedProbCol)
3119 col.add(_compSolver.number(SPxRowId(_compSlackDualRowId)), -SLACKCOEFF);
3120
3121 boundConsCols.add(_realLP->lower(i), R(-infinity), col, 0.0);
3122
3123 if(!isRedProbCol)
3124 col.remove(col.size() - 1);
3125
3126 boundConsColsAdded[i]++;
3127 numBoundConsCols++;
3128 }
3129
3130 if(LT(_realLP->upper(i), R(infinity)))
3131 {
3132 if(!isRedProbCol)
3133 col.add(_compSolver.number(SPxRowId(_compSlackDualRowId)), SLACKCOEFF);
3134
3135 boundConsCols.add(_realLP->upper(i), 0.0, col, R(infinity));
3136
3137 if(!isRedProbCol)
3138 col.remove(col.size() - 1);
3139
3140 boundConsColsAdded[i]++;
3141 numBoundConsCols++;
3142 }
3143 }
3144
3145 col.clear();
3146 _fixedOrigVars[i] = currFixedVars[i];
3147 }
3148 }
3149
3150 // adding the fixed var dual columns to the complementary problem
3151 SPxColId* addedcolids = 0;
3152 spx_alloc(addedcolids, numFixedVars);
3153 _compSolver.addCols(addedcolids, fixedVarsDualCols);
3154
3155 SPxColId tempId;
3156 int addedcolcount = 0;
3157
3158 tempId.inValidate();
3159
3160 for(int i = 0; i < _realLP->nCols(); i++)
3161 {
3162 if(_fixedOrigVars[i] != 0)
3163 {
3164 assert(_fixedOrigVars[i] == -1 || _fixedOrigVars[i] == 1);
3165 _decompFixedVarDualIDs[i] = addedcolids[addedcolcount];
3166 addedcolcount++;
3167 }
3168 else
3169 _decompFixedVarDualIDs[i] = tempId;
3170 }
3171
3172 // adding the bound cons dual columns to the complementary problem
3173 SPxColId* addedbndcolids = 0;
3174 spx_alloc(addedbndcolids, numBoundConsCols);
3175 _compSolver.addCols(addedbndcolids, boundConsCols);
3176
3177 addedcolcount = 0;
3178
3179 for(int i = 0; i < _realLP->nCols(); i++)
3180 {
3181 if(boundConsColsAdded[i] > 0)
3182 {
3183 for(int j = 0; j < boundConsColsAdded[i]; j++)
3184 {
3185 _decompVarBoundDualIDs[i * 2 + j] = addedbndcolids[addedcolcount];
3186 addedcolcount++;
3187 }
3188 }
3189
3190 switch(boundConsColsAdded[i])
3191 {
3192 case 0:
3193 _decompVarBoundDualIDs[i * 2] = tempId;
3194
3195 // FALLTHROUGH
3196 case 1:
3197 _decompVarBoundDualIDs[i * 2 + 1] = tempId;
3198 break;
3199 }
3200 }
3201
3202 // freeing allocated memory
3203 spx_free(addedbndcolids);
3204 spx_free(addedcolids);
3205 spx_free(boundConsColsAdded);
3206 }
3207
3208
3209 /// identify the dual columns related to the fixed variables
3210 template <class R>
3211 void SoPlexBase<R>::_identifyComplementaryPrimalFixedPrimalVars(int* currFixedVars)
3212 {
3213 int numFixedVar = 0;
3214
3215 for(int i = 0; i < _realLP->nCols(); i++)
3216 {
3217 currFixedVars[i] = 0;
3218
3219 if(!_decompReducedProbColRowIDs[i].isValid())
3220 continue;
3221
3222 int rowNumber = _solver.number(_decompReducedProbColRowIDs[i]);
3223
3224 if(_decompReducedProbColRowIDs[i].isValid() &&
3225 (_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_UPPER ||
3226 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_LOWER ||
3227 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FIXED))
3228 {
3229 // setting the value of the _fixedOrigVars array to indicate which variables are at their bounds.
3230 currFixedVars[i] = getOrigVarFixedDirection(i);
3231
3232 numFixedVar++;
3233 }
3234 }
3235
3236 MSG_INFO3(spxout, spxout <<
3237 "Number of fixed primal variables in the complementary (primal) problem: "
3238 << numFixedVar << std::endl);
3239 }
3240
3241
3242
3243 /// updating the dual columns related to the fixed primal variables.
3244 template <class R>
3245 void SoPlexBase<R>::_updateComplementaryPrimalFixedPrimalVars(int* currFixedVars)
3246 {
3247 int numFixedVars = 0;
3248
3249 // NOTE: this loop only goes over the primal columns that are included in the complementary problem, i.e. the
3250 // columns from the original problem.
3251 // 29.04.15 in the current implementation, all bound constraints are included in the reduced problem. So, all
3252 // variables (columns) are included in the reduced problem.
3253 for(int i = 0; i < _nCompPrimalCols; i++)
3254 {
3255 int colNumber = _compSolver.number(SPxColId(_decompCompPrimalColIDs[i]));
3256
3257 if(_fixedOrigVars[colNumber] != currFixedVars[colNumber])
3258 {
3259 if(currFixedVars[colNumber] != 0)
3260 {
3261 assert(currFixedVars[colNumber] == -1 || currFixedVars[colNumber] == 1);
3262
3263 if(currFixedVars[colNumber] == -1)
3264 _compSolver.changeBounds(colNumber, _realLP->lower(SPxColId(_decompPrimalColIDs[i])),
3265 _realLP->lower(SPxColId(_decompPrimalColIDs[i])));
3266 else
3267 _compSolver.changeBounds(colNumber, _realLP->upper(SPxColId(_decompPrimalColIDs[i])),
3268 _realLP->upper(SPxColId(_decompPrimalColIDs[i])));
3269
3270 numFixedVars++;
3271 }
3272 else
3273 {
3274 _compSolver.changeBounds(colNumber, R(-infinity), R(infinity));
3275 }
3276 }
3277
3278 _fixedOrigVars[colNumber] = currFixedVars[colNumber];
3279 }
3280 }
3281
3282
3283
3284 /// updating the complementary dual problem with the original objective function
3285 template <class R>
3286 void SoPlexBase<R>::_setComplementaryDualOriginalObjective()
3287 {
3288 for(int i = 0; i < _realLP->nCols(); i++)
3289 {
3290 assert(_decompCompProbColIDsIdx[i] != -1); // this should be true in the current implementation
3291 int idIndex = _decompCompProbColIDsIdx[i];
3292 int compRowNumber = _compSolver.number(_decompDualRowIDs[idIndex]);
3293
3294 if(_decompCompProbColIDsIdx[i] != -1)
3295 {
3296 // In the dual conversion, when a variable has a non-standard bound it is converted to a free variable.
3297 if(LE(_realLP->lower(i), R(-infinity)) && GE(_realLP->upper(i), R(infinity)))
3298 {
3299 // unrestricted variable
3300 _compSolver.changeLhs(compRowNumber, _realLP->obj(i));
3301 _compSolver.changeRhs(compRowNumber, _realLP->obj(i));
3302 assert(LE(_compSolver.lhs(compRowNumber), _compSolver.rhs(compRowNumber)));
3303 }
3304 else if(LE(_realLP->lower(i), R(-infinity)))
3305 {
3306 // variable with a finite upper bound
3307 _compSolver.changeRhs(compRowNumber, _realLP->obj(i));
3308
3309 if(isZero(_realLP->upper(i)))
3310 _compSolver.changeLhs(compRowNumber, R(-infinity));
3311 else
3312 _compSolver.changeLhs(compRowNumber, _realLP->obj(i));
3313 }
3314 else if(GE(_realLP->upper(i), R(infinity)))
3315 {
3316 // variable with a finite lower bound
3317 _compSolver.changeLhs(compRowNumber, _realLP->obj(i));
3318
3319 if(isZero(_realLP->upper(i)))
3320 _compSolver.changeRhs(compRowNumber, R(infinity));
3321 else
3322 _compSolver.changeRhs(compRowNumber, _realLP->obj(i));
3323 }
3324 else if(NE(_realLP->lower(i), _realLP->upper(i)))
3325 {
3326 // variable with a finite upper and lower bound
3327 if(isZero(_realLP->upper(i)))
3328 {
3329 _compSolver.changeLhs(compRowNumber, _realLP->obj(i));
3330 _compSolver.changeRhs(compRowNumber, R(infinity));
3331 }
3332 else if(isZero(_realLP->upper(i)))
3333 {
3334 _compSolver.changeLhs(compRowNumber, R(-infinity));
3335 _compSolver.changeRhs(compRowNumber, _realLP->obj(i));
3336 }
3337 else
3338 {
3339 _compSolver.changeLhs(compRowNumber, _realLP->obj(i));
3340 _compSolver.changeRhs(compRowNumber, _realLP->obj(i));
3341 }
3342 }
3343 else
3344 {
3345 // fixed variable
3346 _compSolver.changeLhs(compRowNumber, _realLP->obj(i));
3347 _compSolver.changeRhs(compRowNumber, _realLP->obj(i));
3348 }
3349 }
3350 }
3351
3352 // removing the complementary problem slack column dual row
3353 _compSolver.removeRow(_compSlackDualRowId);
3354 }
3355
3356
3357
3358 /// updating the complementary primal problem with the original objective function
3359 template <class R>
3360 void SoPlexBase<R>::_setComplementaryPrimalOriginalObjective()
3361 {
3362 // the comp solver has not removed any columns. Only the slack variables have been added.
3363 assert(_realLP->nCols() == _compSolver.nCols() - 1);
3364
3365 for(int i = 0; i < _realLP->nCols(); i++)
3366 {
3367 int colNumber = _realLP->number(_decompPrimalColIDs[i]);
3368 int compColNumber = _compSolver.number(_decompCompPrimalColIDs[i]);
3369 _compSolver.changeObj(compColNumber, _realLP->maxObj(colNumber));
3370 }
3371
3372 // removing the complementary problem slack column dual row
3373 _compSolver.removeCol(_compSlackColId);
3374 }
3375
3376
3377
3378 /// determining which bound the primal variables will be fixed to.
3379 template <class R>
3380 int SoPlexBase<R>::getOrigVarFixedDirection(int colNum)
3381 {
3382 if(!_decompReducedProbColRowIDs[colNum].isValid())
3383 return 0;
3384
3385 int rowNumber = _solver.number(_decompReducedProbColRowIDs[colNum]);
3386
3387 // setting the value of the _fixedOrigVars array to indicate which variables are at their bounds.
3388 if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_UPPER ||
3389 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FIXED ||
3390 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_FREE)
3391 {
3392 assert(_solver.rhs(rowNumber) < R(infinity));
3393 return 1;
3394 }
3395 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_LOWER)
3396 {
3397 assert(_solver.lhs(rowNumber) > R(-infinity));
3398 return -1;
3399 }
3400
3401 return 0;
3402 }
3403
3404
3405
3406 // @todo update this function and related comments. It has only been hacked together.
3407 /// checks result of the solving process and solves again without preprocessing if necessary
3408 // @todo need to evaluate the solution to ensure that it is solved to optimality and then we are able to perform the
3409 // next steps in the algorithm.
3410 template <class R>
3411 void SoPlexBase<R>::_evaluateSolutionDecomp(SPxSolverBase<R>& solver, SLUFactor<R>& sluFactor,
3412 typename SPxSimplifier<R>::Result result)
3413 {
3414 typename SPxSolverBase<R>::Status solverStat = SPxSolverBase<R>::UNKNOWN;
3415
3416 if(result == SPxSimplifier<R>::INFEASIBLE)
3417 solverStat = SPxSolverBase<R>::INFEASIBLE;
3418 else if(result == SPxSimplifier<R>::DUAL_INFEASIBLE)
3419 solverStat = SPxSolverBase<R>::INForUNBD;
3420 else if(result == SPxSimplifier<R>::UNBOUNDED)
3421 solverStat = SPxSolverBase<R>::UNBOUNDED;
3422 else if(result == SPxSimplifier<R>::VANISHED)
3423 solverStat = SPxSolverBase<R>::OPTIMAL;
3424 else if(result == SPxSimplifier<R>::OKAY)
3425 solverStat = solver.status();
3426
3427 // updating the status of SoPlexBase if the problem solved is the reduced problem.
3428 if(_currentProb == DECOMP_ORIG || _currentProb == DECOMP_RED)
3429 _status = solverStat;
3430
3431 // process result
3432 switch(solverStat)
3433 {
3434 case SPxSolverBase<R>::OPTIMAL:
3435 if(!_isRealLPLoaded)
3436 {
3437 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET));
3438 _decompResolveWithoutPreprocessing(solver, sluFactor, result);
3439 // Need to solve the complementary problem
3440 return;
3441 }
3442 else
3443 _hasBasis = true;
3444
3445 break;
3446
3447 case SPxSolverBase<R>::UNBOUNDED:
3448 case SPxSolverBase<R>::INFEASIBLE:
3449 case SPxSolverBase<R>::INForUNBD:
3450
3451 // in case of infeasibility or unboundedness, we currently can not unsimplify, but have to solve the original LP again
3452 if(!_isRealLPLoaded)
3453 {
3454 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET));
3455 _decompSimplifyAndSolve(solver, sluFactor, false, false);
3456 return;
3457 }
3458 else
3459 _hasBasis = (solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM);
3460
3461 break;
3462
3463 case SPxSolverBase<R>::ABORT_DECOMP:
3464 case SPxSolverBase<R>::ABORT_EXDECOMP:
3465
3466 // in the initialisation of the decomposition simplex, we want to keep the current basis.
3467 if(!_isRealLPLoaded)
3468 {
3469 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET));
3470 _decompResolveWithoutPreprocessing(solver, sluFactor, result);
3471 //_decompSimplifyAndSolve(solver, sluFactor, false, false);
3472 return;
3473 }
3474 else
3475 _hasBasis = (solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM);
3476
3477 break;
3478
3479 case SPxSolverBase<R>::SINGULAR:
3480 case SPxSolverBase<R>::ABORT_CYCLING:
3481
3482 // in case of infeasibility or unboundedness, we currently can not unsimplify, but have to solve the original LP again
3483 if(!_isRealLPLoaded)
3484 {
3485 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET));
3486 _decompSimplifyAndSolve(solver, sluFactor, false, false);
3487 return;
3488 }
3489
3490 break;
3491
3492
3493 // else fallthrough
3494 case SPxSolverBase<R>::ABORT_TIME:
3495 case SPxSolverBase<R>::ABORT_ITER:
3496 case SPxSolverBase<R>::ABORT_VALUE:
3497 case SPxSolverBase<R>::REGULAR:
3498 case SPxSolverBase<R>::RUNNING:
3499
3500 // store regular basis if there is no simplifier and the original problem is not in the solver because of
3501 // scaling; non-optimal bases should currently not be unsimplified
3502 if(_simplifier == 0 && solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM)
3503 {
3504 _basisStatusRows.reSize(_decompLP->nRows());
3505 _basisStatusCols.reSize(_decompLP->nCols());
3506 assert(_basisStatusRows.size() == solver.nRows());
3507 assert(_basisStatusCols.size() == solver.nCols());
3508
3509 solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr());
3510 _hasBasis = true;
3511 }
3512 else
3513 _hasBasis = false;
3514
3515 break;
3516
3517 default:
3518 _hasBasis = false;
3519 break;
3520 }
3521 }
3522
3523
3524
3525
3526 /// checks the dual feasibility of the current basis
3527 template <class R>
3528 bool SoPlexBase<R>::checkBasisDualFeasibility(VectorBase<R> feasVec)
3529 {
3530 assert(_solver.rep() == SPxSolverBase<R>::ROW);
3531 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE);
3532
3533 R feastol = realParam(SoPlexBase<R>::FEASTOL);
3534
3535 // iterating through all elements in the basis to determine whether the dual feasibility is satisfied.
3536 for(int i = 0; i < _solver.nCols(); i++)
3537 {
3538 if(_solver.basis().baseId(i).isSPxRowId()) // find the row id's for rows in the basis
3539 {
3540 int rownumber = _solver.number(_solver.basis().baseId(i));
3541
3542 if(_solver.basis().desc().rowStatus(rownumber) != SPxBasisBase<R>::Desc::P_ON_UPPER &&
3543 _solver.basis().desc().rowStatus(rownumber) != SPxBasisBase<R>::Desc::P_FIXED)
3544 {
3545 if(GT(feasVec[i], (R) 0, feastol))
3546 return false;
3547 }
3548
3549 if(_solver.basis().desc().rowStatus(rownumber) != SPxBasisBase<R>::Desc::P_ON_LOWER &&
3550 _solver.basis().desc().rowStatus(rownumber) != SPxBasisBase<R>::Desc::P_FIXED)
3551 {
3552 if(LT(feasVec[i], (R) 0, feastol))
3553 return false;
3554 }
3555
3556 }
3557 else if(_solver.basis().baseId(
3558 i).isSPxColId()) // get the column id's for the columns in the basis
3559 {
3560 int colnumber = _solver.number(_solver.basis().baseId(i));
3561
3562 if(_solver.basis().desc().colStatus(colnumber) != SPxBasisBase<R>::Desc::P_ON_UPPER &&
3563 _solver.basis().desc().colStatus(colnumber) != SPxBasisBase<R>::Desc::P_FIXED)
3564 {
3565 if(GT(feasVec[i], (R) 0, feastol))
3566 return false;
3567 }
3568
3569 if(_solver.basis().desc().colStatus(colnumber) != SPxBasisBase<R>::Desc::P_ON_LOWER &&
3570 _solver.basis().desc().colStatus(colnumber) != SPxBasisBase<R>::Desc::P_FIXED)
3571 {
3572 if(LT(feasVec[i], (R) 0, feastol))
3573 return false;
3574 }
3575 }
3576 }
3577
3578 return true;
3579 }
3580
3581
3582
3583 /// returns the expected sign of the dual variables for the reduced problem
3584 template <class R>
3585 typename SoPlexBase<R>::DualSign SoPlexBase<R>::getExpectedDualVariableSign(int rowNumber)
3586 {
3587 if(_solver.isRowBasic(rowNumber))
3588 {
3589 if(_solver.basis().desc().rowStatus(rowNumber) != SPxBasisBase<R>::Desc::P_ON_UPPER &&
3590 _solver.basis().desc().rowStatus(rowNumber) != SPxBasisBase<R>::Desc::P_FIXED)
3591 return SoPlexBase<R>::IS_NEG;
3592
3593 if(_solver.basis().desc().rowStatus(rowNumber) != SPxBasisBase<R>::Desc::P_ON_LOWER &&
3594 _solver.basis().desc().rowStatus(rowNumber) != SPxBasisBase<R>::Desc::P_FIXED)
3595 return SoPlexBase<R>::IS_POS;
3596 }
3597
3598 return SoPlexBase<R>::IS_FREE;
3599 }
3600
3601
3602
3603 /// returns the expected sign of the dual variables for the original problem
3604 template <class R>
3605 typename SoPlexBase<R>::DualSign SoPlexBase<R>::getOrigProbDualVariableSign(int rowNumber)
3606 {
3607 if(_realLP->rowType(rowNumber) == LPRowBase<R>::LESS_EQUAL)
3608 return IS_POS;
3609
3610 if(_realLP->rowType(rowNumber) == LPRowBase<R>::GREATER_EQUAL)
3611 return IS_NEG;
3612
3613 if(_realLP->rowType(rowNumber) == LPRowBase<R>::EQUAL)
3614 return IS_FREE;
3615
3616 // this final statement needs to be checked.
3617 if(_realLP->rowType(rowNumber) == LPRowBase<R>::RANGE)
3618 return IS_FREE;
3619
3620 return IS_FREE;
3621 }
3622
3623
3624 /// print display line of flying table
3625 template <class R>
3626 void SoPlexBase<R>::printDecompDisplayLine(SPxSolverBase<R>& solver,
3627 const SPxOut::Verbosity origVerb, bool force, bool forceHead)
3628 {
3629 // setting the verbosity level
3630 const SPxOut::Verbosity currVerb = spxout.getVerbosity();
3631 spxout.setVerbosity(origVerb);
3632
3633 int displayFreq = intParam(SoPlexBase<R>::DECOMP_DISPLAYFREQ);
3634
3635 MSG_INFO1(spxout,
3636
3637 if(forceHead || (_decompDisplayLine % (displayFreq * 30) == 0))
3638 {
3639 spxout << "type | time | iters | red iter | alg iter | rows | cols | shift | value\n";
3640 }
3641 if(force || (_decompDisplayLine % displayFreq == 0))
3642 {
3643 Real currentTime = _statistics->solvingTime->time();
3644 (solver.type() == SPxSolverBase<R>::LEAVE) ? spxout << " L |" : spxout << " E |";
3645 spxout << std::fixed << std::setw(7) << std::setprecision(1) << currentTime << " |";
3646 spxout << std::scientific << std::setprecision(2);
3647 spxout << std::setw(8) << _statistics->iterations << " | ";
3648 spxout << std::scientific << std::setprecision(2);
3649 spxout << std::setw(8) << _statistics->iterationsRedProb << " | ";
3650 spxout << std::scientific << std::setprecision(2);
3651 spxout << std::setw(8) << _statistics->callsReducedProb << " | ";
3652 spxout << std::scientific << std::setprecision(2);
3653 spxout << std::setw(8) << numIncludedRows << " | ";
3654 spxout << std::scientific << std::setprecision(2);
3655 spxout << std::setw(8) << solver.nCols() << " | "
3656 << solver.shift() << " | "
3657 << std::setprecision(8) << solver.value() + solver.objOffset()
3658 << std::endl;
3659
3660 }
3661 _decompDisplayLine++;
3662 );
3663
3664 spxout.setVerbosity(currVerb);
3665 }
3666
3667
3668
3669 /// stores the problem statistics of the original problem
3670 template <class R>
3671 void SoPlexBase<R>::getOriginalProblemStatistics()
3672 {
3673 numProbRows = _realLP->nRows();
3674 numProbCols = _realLP->nCols();
3675 nNonzeros = _realLP->nNzos();
3676 minAbsNonzero = _realLP->minAbsNzo();
3677 maxAbsNonzero = _realLP->maxAbsNzo();
3678
3679 origCountLower = 0;
3680 origCountUpper = 0;
3681 origCountBoxed = 0;
3682 origCountFreeCol = 0;
3683
3684 origCountLhs = 0;
3685 origCountRhs = 0;
3686 origCountRanged = 0;
3687 origCountFreeRow = 0;
3688
3689 for(int i = 0; i < _realLP->nCols(); i++)
3690 {
3691 bool hasLower = false;
3692 bool hasUpper = false;
3693
3694 if(_realLP->lower(i) > R(-infinity))
3695 {
3696 origCountLower++;
3697 hasLower = true;
3698 }
3699
3700 if(_realLP->upper(i) < R(infinity))
3701 {
3702 origCountUpper++;
3703 hasUpper = true;
3704 }
3705
3706 if(hasUpper && hasLower)
3707 {
3708 origCountBoxed++;
3709 origCountUpper--;
3710 origCountLower--;
3711 }
3712
3713 if(!hasUpper && !hasLower)
3714 origCountFreeCol++;
3715 }
3716
3717 for(int i = 0; i < _realLP->nRows(); i++)
3718 {
3719 bool hasRhs = false;
3720 bool hasLhs = false;
3721
3722 if(_realLP->lhs(i) > R(-infinity))
3723 {
3724 origCountLhs++;
3725 hasLhs = true;
3726 }
3727
3728 if(_realLP->rhs(i) < R(infinity))
3729 {
3730 origCountRhs++;
3731 hasRhs = true;
3732 }
3733
3734 if(hasRhs && hasLhs)
3735 {
3736 if(EQ(_realLP->rhs(i), _realLP->lhs(i)))
3737 origCountEqual++;
3738 else
3739 origCountRanged++;
3740
3741 origCountLhs--;
3742 origCountRhs--;
3743 }
3744
3745 if(!hasRhs && !hasLhs)
3746 origCountFreeRow++;
3747 }
3748 }
3749
3750 template <class R>
3751 void SoPlexBase<R>::printOriginalProblemStatistics(std::ostream& os)
3752 {
3753 os << " Columns : " << numProbCols << "\n"
3754 << " boxed : " << origCountBoxed << "\n"
3755 << " lower bound : " << origCountLower << "\n"
3756 << " upper bound : " << origCountUpper << "\n"
3757 << " free : " << origCountFreeCol << "\n"
3758 << " Rows : " << numProbRows << "\n"
3759 << " equal : " << origCountEqual << "\n"
3760 << " ranged : " << origCountRanged << "\n"
3761 << " lhs : " << origCountLhs << "\n"
3762 << " rhs : " << origCountRhs << "\n"
3763 << " free : " << origCountFreeRow << "\n"
3764 << " Nonzeros : " << nNonzeros << "\n"
3765 << " per column : " << R(nNonzeros) / R(numProbCols) << "\n"
3766 << " per row : " << R(nNonzeros) / R(numProbRows) << "\n"
3767 << " sparsity : " << R(nNonzeros) / R(numProbCols) / R(numProbRows) << "\n"
3768 << " min. abs. value : " << R(minAbsNonzero) << "\n"
3769 << " max. abs. value : " << R(maxAbsNonzero) << "\n";
3770 }
3771
3772
3773
3774 /// gets the coefficient of the slack variable in the primal complementary problem
3775 template <class R>
3776 R SoPlexBase<R>::getCompSlackVarCoeff(int primalRowNum)
3777 {
3778 int indDir = 1;
3779
3780 switch(_realLP->rowType(_decompPrimalRowIDs[primalRowNum]))
3781 {
3782 // NOTE: check the sign of the slackCoeff for the Range constraints. This will depend on the method of
3783 // dual conversion.
3784 case LPRowBase<R>::RANGE:
3785 assert((primalRowNum < _nPrimalRows - 1
3786 && _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum])) ==
3787 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum + 1]))) ||
3788 (primalRowNum > 0 && _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum - 1])) ==
3789 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum]))));
3790
3791 // determine with primalRowNum and primalRowNum+1 or primalRowNum-1 and primalRowNum have the same row id.
3792 if(_realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum - 1])) ==
3793 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum])))
3794 indDir = -1;
3795
3796 if(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[primalRowNum]))) <
3797 _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[primalRowNum + indDir]))))
3798 return -SLACKCOEFF;
3799 else
3800 return SLACKCOEFF;
3801
3802 break;
3803
3804 case LPRowBase<R>::GREATER_EQUAL:
3805 return -SLACKCOEFF;
3806 break;
3807
3808 case LPRowBase<R>::EQUAL:
3809 assert((primalRowNum < _nPrimalRows - 1
3810 && _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum])) ==
3811 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum + 1]))) ||
3812 (primalRowNum > 0 && _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum - 1])) ==
3813 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum]))));
3814
3815 // FALLTHROUGH
3816 case LPRowBase<R>::LESS_EQUAL:
3817 return SLACKCOEFF;
3818 break;
3819
3820 default:
3821 throw SPxInternalCodeException("XDECOMPSL01 This should never happen.");
3822 }
3823 }
3824
3825
3826
3827 /// gets violation of bounds; returns true on success
3828 template <class R>
3829 bool SoPlexBase<R>::getDecompBoundViolation(R& maxviol, R& sumviol)
3830 {
3831 R feastol = realParam(SoPlexBase<R>::FEASTOL);
3832
3833 VectorBase<R>& primal = _solReal._primal;
3834 assert(primal.dim() == _realLP->nCols());
3835
3836 _nDecompViolBounds = 0;
3837
3838 maxviol = 0.0;
3839 sumviol = 0.0;
3840
3841 bool isViol = false;
3842 bool isMaxViol = false;
3843
3844 for(int i = _realLP->nCols() - 1; i >= 0; i--)
3845 {
3846 R currviol = 0.0;
3847
3848 R viol = _realLP->lower(i) - primal[i];
3849
3850 isViol = false;
3851 isMaxViol = false;
3852
3853 if(viol > 0.0)
3854 {
3855 sumviol += viol;
3856
3857 if(viol > maxviol)
3858 {
3859 maxviol = viol;
3860 isMaxViol = true;
3861 }
3862
3863 if(currviol < viol)
3864 currviol = viol;
3865 }
3866
3867 if(GT(viol, R(0.0), feastol))
3868 isViol = true;
3869
3870 viol = primal[i] - _realLP->upper(i);
3871
3872 if(viol > 0.0)
3873 {
3874 sumviol += viol;
3875
3876 if(viol > maxviol)
3877 {
3878 maxviol = viol;
3879 isMaxViol = true;
3880 }
3881
3882 if(currviol < viol)
3883 currviol = viol;
3884 }
3885
3886 if(GT(viol, R(0.0), feastol))
3887 isViol = true;
3888
3889 if(isViol)
3890 {
3891 // updating the violated bounds list
3892 if(isMaxViol)
3893 {
3894 _decompViolatedBounds[_nDecompViolBounds] = _decompViolatedBounds[0];
3895 _decompViolatedBounds[0] = i;
3896 }
3897 else
3898 _decompViolatedBounds[_nDecompViolBounds] = i;
3899
3900 _nDecompViolBounds++;
3901 }
3902 }
3903
3904 return true;
3905 }
3906
3907
3908
3909 /// gets violation of constraints; returns true on success
3910 template <class R>
3911 bool SoPlexBase<R>::getDecompRowViolation(R& maxviol, R& sumviol)
3912 {
3913 R feastol = realParam(SoPlexBase<R>::FEASTOL);
3914
3915 VectorBase<R>& primal = _solReal._primal;
3916 assert(primal.dim() == _realLP->nCols());
3917
3918 VectorBase<R> activity(_realLP->nRows());
3919 _realLP->computePrimalActivity(primal, activity);
3920
3921 _nDecompViolRows = 0;
3922
3923 maxviol = 0.0;
3924 sumviol = 0.0;
3925
3926 bool isViol = false;
3927 bool isMaxViol = false;
3928
3929 for(int i = _realLP->nRows() - 1; i >= 0; i--)
3930 {
3931 R currviol = 0.0;
3932
3933 isViol = false;
3934 isMaxViol = false;
3935
3936 R viol = _realLP->lhs(i) - activity[i];
3937
3938 if(viol > 0.0)
3939 {
3940 sumviol += viol;
3941
3942 if(viol > maxviol)
3943 {
3944 maxviol = viol;
3945 isMaxViol = true;
3946 }
3947
3948 if(viol > currviol)
3949 currviol = viol;
3950 }
3951
3952 if(GT(viol, R(0.0), feastol))
3953 isViol = true;
3954
3955 viol = activity[i] - _realLP->rhs(i);
3956
3957 if(viol > 0.0)
3958 {
3959 sumviol += viol;
3960
3961 if(viol > maxviol)
3962 {
3963 maxviol = viol;
3964 isMaxViol = true;
3965 }
3966
3967 if(viol > currviol)
3968 currviol = viol;
3969 }
3970
3971 if(GT(viol, R(0.0), feastol))
3972 isViol = true;
3973
3974 if(isViol)
3975 {
3976 // updating the violated rows list
3977 if(isMaxViol)
3978 {
3979 _decompViolatedRows[_nDecompViolRows] = _decompViolatedRows[0];
3980 _decompViolatedRows[0] = i;
3981 }
3982 else
3983 _decompViolatedRows[_nDecompViolRows] = i;
3984
3985 _nDecompViolRows++;
3986 }
3987 }
3988
3989 return true;
3990 }
3991
3992
3993
3994 /// function call to terminate the decomposition simplex
3995 template <class R>
3996 bool SoPlexBase<R>::decompTerminate(R timeLimit)
3997 {
3998 R maxTime = timeLimit;
3999
4000 // check if a time limit is actually set
4001 if(maxTime < 0 || maxTime >= R(infinity))
4002 return false;
4003
4004 Real currentTime = _statistics->solvingTime->time();
4005
4006 if(currentTime >= maxTime)
4007 {
4008 MSG_INFO2(spxout, spxout << " --- timelimit (" << _solver.getMaxTime()
4009 << ") reached" << std::endl;)
4010 _solver.setSolverStatus(SPxSolverBase<R>::ABORT_TIME);
4011 return true;
4012 }
4013
4014 return false;
4015 }
4016
4017
4018
4019 /// function to retrieve the original problem row basis status from the reduced and complementary problems
4020 template <class R>
4021 void SoPlexBase<R>::getOriginalProblemBasisRowStatus(DataArray< int >& degenerateRowNums,
4022 DataArray< typename SPxSolverBase<R>::VarStatus >& degenerateRowStatus, int& nDegenerateRows,
4023 int& nNonBasicRows)
4024 {
4025 R feastol = realParam(SoPlexBase<R>::FEASTOL);
4026 int basicRow = 0;
4027
4028 // looping over all rows not contained in the complementary problem
4029 for(int i = 0; i < _nElimPrimalRows; i++)
4030 {
4031 int rowNumber = _realLP->number(_decompElimPrimalRowIDs[i]);
4032
4033 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]);
4034 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows());
4035
4036 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER) /*||
4037 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
4038 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], 0.0, feastol)) )*/
4039 {
4040 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER;
4041 }
4042 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER) /*||
4043 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
4044 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), 0.0, feastol)) )*/
4045 {
4046 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER;
4047 }
4048 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED)
4049 {
4050 _basisStatusRows[rowNumber] = SPxSolverBase<R>::FIXED;
4051 }
4052 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FREE)
4053 {
4054 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ZERO;
4055 }
4056 else
4057 {
4058 _basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC;
4059 basicRow++;
4060 }
4061 }
4062
4063 for(int i = 0; i < _nPrimalRows; i++)
4064 {
4065 int rowNumber = _realLP->number(_decompPrimalRowIDs[i]);
4066
4067 // this loop runs over all rows previously in the complementary problem.
4068 if(_decompReducedProbRows[rowNumber])
4069 {
4070 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]);
4071 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows());
4072
4073 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER)
4074 {
4075 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER;
4076 }
4077 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
4078 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], 0.0, feastol))
4079 {
4080 // if a row is non basic, but is at its bound then the row number and status is stored
4081 degenerateRowNums[nDegenerateRows] = rowNumber;
4082 degenerateRowStatus[nDegenerateRows] = SPxSolverBase<R>::ON_UPPER;
4083 nDegenerateRows++;
4084 }
4085 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER)
4086 {
4087 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER;
4088 }
4089 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
4090 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), 0.0, feastol))
4091 {
4092 // if a row is non basic, but is at its bound then the row number and status is stored
4093 degenerateRowNums[nDegenerateRows] = rowNumber;
4094 degenerateRowStatus[nDegenerateRows] = SPxSolverBase<R>::ON_LOWER;
4095 nDegenerateRows++;
4096 }
4097 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED)
4098 {
4099 _basisStatusRows[rowNumber] = SPxSolverBase<R>::FIXED;
4100 }
4101 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FREE)
4102 {
4103 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ZERO;
4104 }
4105 else
4106 {
4107 _basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC;
4108 basicRow++;
4109 }
4110 }
4111 else
4112 {
4113 _basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC;
4114 basicRow++;
4115
4116 switch(_realLP->rowType(_decompPrimalRowIDs[i]))
4117 {
4118 case LPRowBase<R>::RANGE:
4119 //assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) ==
4120 //_realLP->number(SPxColId(_decompPrimalRowIDs[i+1])));
4121 //assert(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) !=
4122 //_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))));
4123
4124 //// NOTE: This is not correct at present. Need to modify the inequalities. Look at the dual conversion
4125 //// function to get this correct.
4126 //if( _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) <
4127 //_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))))
4128 //{
4129 //if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i]))) ==
4130 //SPxBasisBase<R>::Desc::D_ON_LOWER )
4131 //{
4132 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER;
4133 //}
4134 //else if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))) ==
4135 //SPxBasisBase<R>::Desc::D_ON_UPPER )
4136 //{
4137 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER;
4138 //}
4139 //else
4140 //{
4141 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC;
4142 //basicRow++;
4143 //}
4144 //}
4145 //else
4146 //{
4147 //if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i]))) ==
4148 //SPxBasisBase<R>::Desc::D_ON_UPPER )
4149 //{
4150 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER;
4151 //}
4152 //else if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))) ==
4153 //SPxBasisBase<R>::Desc::D_ON_LOWER )
4154 //{
4155 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER;
4156 //}
4157 //else
4158 //{
4159 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC;
4160 //basicRow++;
4161 //}
4162 //}
4163
4164 i++;
4165 break;
4166
4167 case LPRowBase<R>::EQUAL:
4168 //assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) ==
4169 //_realLP->number(SPxColId(_decompPrimalRowIDs[i+1])));
4170
4171 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::FIXED;
4172
4173 i++;
4174 break;
4175
4176 case LPRowBase<R>::GREATER_EQUAL:
4177 //if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i]))) ==
4178 //SPxBasisBase<R>::Desc::D_ON_LOWER )
4179 //{
4180 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER;
4181 //}
4182 //else
4183 //{
4184 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC;
4185 //basicRow++;
4186 //}
4187
4188 break;
4189
4190 case LPRowBase<R>::LESS_EQUAL:
4191 //if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i]))) ==
4192 //SPxBasisBase<R>::Desc::D_ON_UPPER )
4193 //{
4194 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER;
4195 //}
4196 //else
4197 //{
4198 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC;
4199 //basicRow++;
4200 //}
4201
4202 break;
4203
4204 default:
4205 throw SPxInternalCodeException("XDECOMPSL01 This should never happen.");
4206 }
4207 }
4208 }
4209
4210 nNonBasicRows = _realLP->nRows() - basicRow - nDegenerateRows;
4211 MSG_INFO2(spxout, spxout << "Number of non-basic rows: " << nNonBasicRows << " (from "
4212 << _realLP->nRows() << ")" << std::endl);
4213 }
4214
4215
4216 /// function to retrieve the column status for the original problem basis from the reduced and complementary problems
4217 // all columns are currently in the reduced problem, so it is only necessary to retrieve the status from there.
4218 // However, a transformation of the columns has been made, so it is only possible to retrieve the status from the
4219 // variable bound constraints.
4220 template <class R>
4221 void SoPlexBase<R>::getOriginalProblemBasisColStatus(int& nNonBasicCols)
4222 {
4223 R feastol = realParam(SoPlexBase<R>::FEASTOL);
4224 int basicCol = 0;
4225 int numZeroDual = 0;
4226
4227 // looping over all variable bound constraints
4228 for(int i = 0; i < _realLP->nCols(); i++)
4229 {
4230 if(!_decompReducedProbColRowIDs[i].isValid())
4231 continue;
4232
4233 int rowNumber = _solver.number(_decompReducedProbColRowIDs[i]);
4234
4235 if(_decompReducedProbColRowIDs[i].isValid())
4236 {
4237 if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_UPPER) /*||
4238 (_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_ON_LOWER &&
4239 EQ(_solver.rhs(rowNumber) - _solver.pVec()[rowNumber], 0.0, feastol)) )*/
4240 {
4241 _basisStatusCols[i] = SPxSolverBase<R>::ON_UPPER;
4242 }
4243 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_LOWER) /*||
4244 (_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_ON_UPPER &&
4245 EQ(_solver.pVec()[rowNumber] - _solver.lhs(rowNumber), 0.0, feastol)) )*/
4246 {
4247 _basisStatusCols[i] = SPxSolverBase<R>::ON_LOWER;
4248 }
4249 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FIXED)
4250 {
4251 _basisStatusCols[i] = SPxSolverBase<R>::FIXED;
4252 }
4253 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FREE)
4254 {
4255 _basisStatusCols[i] = SPxSolverBase<R>::ZERO;
4256 }
4257 else
4258 {
4259 _basisStatusCols[i] = SPxSolverBase<R>::BASIC;
4260 basicCol++;
4261 }
4262 }
4263
4264 if(EQ(_solver.fVec()[i], 0.0, feastol))
4265 numZeroDual++;
4266 }
4267
4268 nNonBasicCols = _realLP->nCols() - basicCol;
4269 MSG_INFO2(spxout, spxout << "Number of non-basic columns: "
4270 << nNonBasicCols << " (from " << _realLP->nCols() << ")" << std::endl
4271 << "Number of zero dual columns: " << numZeroDual << " (from " << _realLP->nCols() << ")" <<
4272 std::endl);
4273 }
4274
4275
4276
4277 /// function to build a basis for the original problem as given by the solution to the reduced problem
4278 template <class R>
4279 void SoPlexBase<R>::_writeOriginalProblemBasis(const char* filename, NameSet* rowNames,
4280 NameSet* colNames, bool cpxFormat)
4281 {
4282 int numrows = _realLP->nRows();
4283 int numcols = _realLP->nCols();
4284 int nNonBasicRows = 0;
4285 int nNonBasicCols = 0;
4286 int nDegenerateRows = 0;
4287 DataArray< int > degenerateRowNums; // array to store the orig prob row number of degenerate rows
4288 DataArray< typename SPxSolverBase<R>::VarStatus >
4289 degenerateRowStatus; // possible status for the degenerate rows in the orig prob
4290
4291 // setting the basis status rows and columns to size of the original problem
4292 _basisStatusRows.reSize(numrows);
4293 _basisStatusCols.reSize(numcols);
4294
4295 degenerateRowNums.reSize(numrows);
4296 degenerateRowStatus.reSize(numrows);
4297
4298 // setting the row and column status for the original problem
4299 getOriginalProblemBasisRowStatus(degenerateRowNums, degenerateRowStatus, nDegenerateRows,
4300 nNonBasicRows);
4301 getOriginalProblemBasisColStatus(nNonBasicCols);
4302
4303 // checking the consistency of the non-basic rows and columns for printing out the basis.
4304 // it is necessary to have numcols of non-basic rows and columns.
4305 // if there are not enought non-basic rows and columns, then the degenerate rows are set as non-basic.
4306 // all degenerate rows that are not changed to non-basic must be set to basic.
4307 assert(nDegenerateRows + nNonBasicRows + nNonBasicCols >= numcols);
4308 int degenRowsSetNonBasic = 0;
4309
4310 for(int i = 0; i < nDegenerateRows; i++)
4311 {
4312 if(nNonBasicRows + nNonBasicCols + degenRowsSetNonBasic < numcols)
4313 {
4314 _basisStatusRows[degenerateRowNums[i]] = degenerateRowStatus[i];
4315 degenRowsSetNonBasic++;
4316 }
4317 else
4318 _basisStatusRows[degenerateRowNums[i]] = SPxSolverBase<R>::BASIC;
4319 }
4320
4321 // writing the basis file for the original problem
4322 bool wasRealLPLoaded = _isRealLPLoaded;
4323 _isRealLPLoaded = false;
4324 writeBasisFile(filename, rowNames, colNames, cpxFormat);
4325 _isRealLPLoaded = wasRealLPLoaded;
4326 }
4327 } // namespace soplex
4328