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