1 | /* $Id: ClpSimplex.cpp 2393 2019-02-21 12:01:58Z forrest $ */ |
---|
2 | // Copyright (C) 2002, International Business Machines |
---|
3 | // Corporation and others. All Rights Reserved. |
---|
4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
5 | |
---|
6 | //#undef NDEBUG |
---|
7 | |
---|
8 | #include "ClpConfig.h" |
---|
9 | |
---|
10 | #include "CoinPragma.hpp" |
---|
11 | #include <math.h> |
---|
12 | |
---|
13 | #if SLIM_CLP == 2 |
---|
14 | #define SLIM_NOIO |
---|
15 | #endif |
---|
16 | #include "CoinHelperFunctions.hpp" |
---|
17 | #include "CoinFloatEqual.hpp" |
---|
18 | #include "ClpSimplex.hpp" |
---|
19 | #include "ClpFactorization.hpp" |
---|
20 | #include "ClpPackedMatrix.hpp" |
---|
21 | #include "CoinIndexedVector.hpp" |
---|
22 | #include "ClpDualRowDantzig.hpp" |
---|
23 | #include "ClpDualRowSteepest.hpp" |
---|
24 | #include "ClpPrimalColumnDantzig.hpp" |
---|
25 | #include "ClpPrimalColumnSteepest.hpp" |
---|
26 | #include "ClpPEDualRowSteepest.hpp" |
---|
27 | #include "ClpPEPrimalColumnSteepest.hpp" |
---|
28 | #include "ClpNonLinearCost.hpp" |
---|
29 | #include "ClpMessage.hpp" |
---|
30 | #include "ClpEventHandler.hpp" |
---|
31 | #include "ClpLinearObjective.hpp" |
---|
32 | #include "ClpHelperFunctions.hpp" |
---|
33 | #include "CoinModel.hpp" |
---|
34 | #include "CoinLpIO.hpp" |
---|
35 | #include <cfloat> |
---|
36 | #if CLP_HAS_ABC |
---|
37 | #include "CoinAbcCommon.hpp" |
---|
38 | #endif |
---|
39 | |
---|
40 | #include <string> |
---|
41 | #include <stdio.h> |
---|
42 | #include <iostream> |
---|
43 | //############################################################################# |
---|
44 | |
---|
45 | ClpSimplex::ClpSimplex(bool emptyMessages) |
---|
46 | : |
---|
47 | |
---|
48 | ClpModel(emptyMessages) |
---|
49 | , bestPossibleImprovement_(0.0) |
---|
50 | , zeroTolerance_(1.0e-13) |
---|
51 | , columnPrimalSequence_(-2) |
---|
52 | , rowPrimalSequence_(-2) |
---|
53 | , bestObjectiveValue_(-COIN_DBL_MAX) |
---|
54 | , moreSpecialOptions_(2) |
---|
55 | , baseIteration_(0) |
---|
56 | , vectorMode_(0) |
---|
57 | , primalToleranceToGetOptimal_(-1.0) |
---|
58 | , largeValue_(1.0e15) |
---|
59 | , largestPrimalError_(0.0) |
---|
60 | , largestDualError_(0.0) |
---|
61 | , alphaAccuracy_(-1.0) |
---|
62 | , dualBound_(1.0e10) |
---|
63 | , alpha_(0.0) |
---|
64 | , theta_(0.0) |
---|
65 | , lowerIn_(0.0) |
---|
66 | , valueIn_(0.0) |
---|
67 | , upperIn_(-COIN_DBL_MAX) |
---|
68 | , dualIn_(0.0) |
---|
69 | , lowerOut_(-1) |
---|
70 | , valueOut_(-1) |
---|
71 | , upperOut_(-1) |
---|
72 | , dualOut_(-1) |
---|
73 | , dualTolerance_(1.0e-7) |
---|
74 | , primalTolerance_(1.0e-7) |
---|
75 | , sumDualInfeasibilities_(0.0) |
---|
76 | , sumPrimalInfeasibilities_(0.0) |
---|
77 | , infeasibilityCost_(1.0e10) |
---|
78 | , sumOfRelaxedDualInfeasibilities_(0.0) |
---|
79 | , sumOfRelaxedPrimalInfeasibilities_(0.0) |
---|
80 | , acceptablePivot_(1.0e-8) |
---|
81 | , lower_(NULL) |
---|
82 | , rowLowerWork_(NULL) |
---|
83 | , columnLowerWork_(NULL) |
---|
84 | , upper_(NULL) |
---|
85 | , rowUpperWork_(NULL) |
---|
86 | , columnUpperWork_(NULL) |
---|
87 | , cost_(NULL) |
---|
88 | , rowObjectiveWork_(NULL) |
---|
89 | , objectiveWork_(NULL) |
---|
90 | , sequenceIn_(-1) |
---|
91 | , directionIn_(-1) |
---|
92 | , sequenceOut_(-1) |
---|
93 | , directionOut_(-1) |
---|
94 | , pivotRow_(-1) |
---|
95 | , lastGoodIteration_(-100) |
---|
96 | , dj_(NULL) |
---|
97 | , rowReducedCost_(NULL) |
---|
98 | , reducedCostWork_(NULL) |
---|
99 | , solution_(NULL) |
---|
100 | , rowActivityWork_(NULL) |
---|
101 | , columnActivityWork_(NULL) |
---|
102 | , numberDualInfeasibilities_(0) |
---|
103 | , numberDualInfeasibilitiesWithoutFree_(0) |
---|
104 | , numberPrimalInfeasibilities_(100) |
---|
105 | , numberRefinements_(0) |
---|
106 | , pivotVariable_(NULL) |
---|
107 | , factorization_(NULL) |
---|
108 | , savedSolution_(NULL) |
---|
109 | , numberTimesOptimal_(0) |
---|
110 | , disasterArea_(NULL) |
---|
111 | , changeMade_(1) |
---|
112 | , algorithm_(0) |
---|
113 | , forceFactorization_(-1) |
---|
114 | , perturbation_(100) |
---|
115 | , nonLinearCost_(NULL) |
---|
116 | , lastBadIteration_(-999999) |
---|
117 | , lastFlaggedIteration_(-999999) |
---|
118 | , numberFake_(0) |
---|
119 | , numberChanged_(0) |
---|
120 | , progressFlag_(0) |
---|
121 | , firstFree_(-1) |
---|
122 | , numberExtraRows_(0) |
---|
123 | , maximumBasic_(0) |
---|
124 | , dontFactorizePivots_(0) |
---|
125 | , incomingInfeasibility_(1.0) |
---|
126 | , allowedInfeasibility_(10.0) |
---|
127 | , automaticScale_(0) |
---|
128 | , maximumPerturbationSize_(0) |
---|
129 | , perturbationArray_(NULL) |
---|
130 | , baseModel_(NULL) |
---|
131 | #ifdef ABC_INHERIT |
---|
132 | , abcSimplex_(NULL) |
---|
133 | , abcState_(0) |
---|
134 | #endif |
---|
135 | { |
---|
136 | int i; |
---|
137 | for (i = 0; i < 6; i++) { |
---|
138 | rowArray_[i] = NULL; |
---|
139 | columnArray_[i] = NULL; |
---|
140 | } |
---|
141 | for (i = 0; i < 4; i++) { |
---|
142 | spareIntArray_[i] = 0; |
---|
143 | spareDoubleArray_[i] = 0.0; |
---|
144 | } |
---|
145 | saveStatus_ = NULL; |
---|
146 | // get an empty factorization so we can set tolerances etc |
---|
147 | getEmptyFactorization(); |
---|
148 | // Say sparse |
---|
149 | factorization_->sparseThreshold(1); |
---|
150 | // say Steepest pricing |
---|
151 | dualRowPivot_ = new ClpDualRowSteepest(); |
---|
152 | // say Steepest pricing |
---|
153 | primalColumnPivot_ = new ClpPrimalColumnSteepest(); |
---|
154 | solveType_ = 1; // say simplex based life form |
---|
155 | } |
---|
156 | |
---|
157 | // Subproblem constructor |
---|
158 | ClpSimplex::ClpSimplex(const ClpModel *rhs, |
---|
159 | int numberRows, const int *whichRow, |
---|
160 | int numberColumns, const int *whichColumn, |
---|
161 | bool dropNames, bool dropIntegers, bool fixOthers) |
---|
162 | : ClpModel(rhs, numberRows, whichRow, |
---|
163 | numberColumns, whichColumn, dropNames, dropIntegers) |
---|
164 | , bestPossibleImprovement_(0.0) |
---|
165 | , zeroTolerance_(1.0e-13) |
---|
166 | , columnPrimalSequence_(-2) |
---|
167 | , rowPrimalSequence_(-2) |
---|
168 | , bestObjectiveValue_(-COIN_DBL_MAX) |
---|
169 | , moreSpecialOptions_(2) |
---|
170 | , baseIteration_(0) |
---|
171 | , vectorMode_(0) |
---|
172 | , primalToleranceToGetOptimal_(-1.0) |
---|
173 | , largeValue_(1.0e15) |
---|
174 | , largestPrimalError_(0.0) |
---|
175 | , largestDualError_(0.0) |
---|
176 | , alphaAccuracy_(-1.0) |
---|
177 | , dualBound_(1.0e10) |
---|
178 | , alpha_(0.0) |
---|
179 | , theta_(0.0) |
---|
180 | , lowerIn_(0.0) |
---|
181 | , valueIn_(0.0) |
---|
182 | , upperIn_(-COIN_DBL_MAX) |
---|
183 | , dualIn_(0.0) |
---|
184 | , lowerOut_(-1) |
---|
185 | , valueOut_(-1) |
---|
186 | , upperOut_(-1) |
---|
187 | , dualOut_(-1) |
---|
188 | , dualTolerance_(1.0e-7) |
---|
189 | , primalTolerance_(1.0e-7) |
---|
190 | , sumDualInfeasibilities_(0.0) |
---|
191 | , sumPrimalInfeasibilities_(0.0) |
---|
192 | , infeasibilityCost_(1.0e10) |
---|
193 | , sumOfRelaxedDualInfeasibilities_(0.0) |
---|
194 | , sumOfRelaxedPrimalInfeasibilities_(0.0) |
---|
195 | , acceptablePivot_(1.0e-8) |
---|
196 | , lower_(NULL) |
---|
197 | , rowLowerWork_(NULL) |
---|
198 | , columnLowerWork_(NULL) |
---|
199 | , upper_(NULL) |
---|
200 | , rowUpperWork_(NULL) |
---|
201 | , columnUpperWork_(NULL) |
---|
202 | , cost_(NULL) |
---|
203 | , rowObjectiveWork_(NULL) |
---|
204 | , objectiveWork_(NULL) |
---|
205 | , sequenceIn_(-1) |
---|
206 | , directionIn_(-1) |
---|
207 | , sequenceOut_(-1) |
---|
208 | , directionOut_(-1) |
---|
209 | , pivotRow_(-1) |
---|
210 | , lastGoodIteration_(-100) |
---|
211 | , dj_(NULL) |
---|
212 | , rowReducedCost_(NULL) |
---|
213 | , reducedCostWork_(NULL) |
---|
214 | , solution_(NULL) |
---|
215 | , rowActivityWork_(NULL) |
---|
216 | , columnActivityWork_(NULL) |
---|
217 | , numberDualInfeasibilities_(0) |
---|
218 | , numberDualInfeasibilitiesWithoutFree_(0) |
---|
219 | , numberPrimalInfeasibilities_(100) |
---|
220 | , numberRefinements_(0) |
---|
221 | , pivotVariable_(NULL) |
---|
222 | , factorization_(NULL) |
---|
223 | , savedSolution_(NULL) |
---|
224 | , numberTimesOptimal_(0) |
---|
225 | , disasterArea_(NULL) |
---|
226 | , changeMade_(1) |
---|
227 | , algorithm_(0) |
---|
228 | , forceFactorization_(-1) |
---|
229 | , perturbation_(100) |
---|
230 | , nonLinearCost_(NULL) |
---|
231 | , lastBadIteration_(-999999) |
---|
232 | , lastFlaggedIteration_(-999999) |
---|
233 | , numberFake_(0) |
---|
234 | , numberChanged_(0) |
---|
235 | , progressFlag_(0) |
---|
236 | , firstFree_(-1) |
---|
237 | , numberExtraRows_(0) |
---|
238 | , maximumBasic_(0) |
---|
239 | , dontFactorizePivots_(0) |
---|
240 | , incomingInfeasibility_(1.0) |
---|
241 | , allowedInfeasibility_(10.0) |
---|
242 | , automaticScale_(0) |
---|
243 | , maximumPerturbationSize_(0) |
---|
244 | , perturbationArray_(NULL) |
---|
245 | , baseModel_(NULL) |
---|
246 | #ifdef ABC_INHERIT |
---|
247 | , abcSimplex_(NULL) |
---|
248 | , abcState_(0) |
---|
249 | #endif |
---|
250 | { |
---|
251 | int i; |
---|
252 | for (i = 0; i < 6; i++) { |
---|
253 | rowArray_[i] = NULL; |
---|
254 | columnArray_[i] = NULL; |
---|
255 | } |
---|
256 | for (i = 0; i < 4; i++) { |
---|
257 | spareIntArray_[i] = 0; |
---|
258 | spareDoubleArray_[i] = 0.0; |
---|
259 | } |
---|
260 | saveStatus_ = NULL; |
---|
261 | // get an empty factorization so we can set tolerances etc |
---|
262 | getEmptyFactorization(); |
---|
263 | // say Steepest pricing |
---|
264 | dualRowPivot_ = new ClpDualRowSteepest(); |
---|
265 | // say Steepest pricing |
---|
266 | primalColumnPivot_ = new ClpPrimalColumnSteepest(); |
---|
267 | solveType_ = 1; // say simplex based life form |
---|
268 | if (fixOthers) { |
---|
269 | int numberOtherColumns = rhs->numberColumns(); |
---|
270 | int numberOtherRows = rhs->numberRows(); |
---|
271 | double *solution = new double[numberOtherColumns]; |
---|
272 | CoinZeroN(solution, numberOtherColumns); |
---|
273 | int i; |
---|
274 | for (i = 0; i < numberColumns; i++) { |
---|
275 | int iColumn = whichColumn[i]; |
---|
276 | if (solution[iColumn]) |
---|
277 | fixOthers = false; // duplicates |
---|
278 | solution[iColumn] = 1.0; |
---|
279 | } |
---|
280 | if (fixOthers) { |
---|
281 | const double *otherSolution = rhs->primalColumnSolution(); |
---|
282 | const double *objective = rhs->objective(); |
---|
283 | double offset = 0.0; |
---|
284 | for (i = 0; i < numberOtherColumns; i++) { |
---|
285 | if (solution[i]) { |
---|
286 | solution[i] = 0.0; // in |
---|
287 | } else { |
---|
288 | solution[i] = otherSolution[i]; |
---|
289 | offset += objective[i] * otherSolution[i]; |
---|
290 | } |
---|
291 | } |
---|
292 | double *rhsModification = new double[numberOtherRows]; |
---|
293 | CoinZeroN(rhsModification, numberOtherRows); |
---|
294 | rhs->matrix()->times(solution, rhsModification); |
---|
295 | for (i = 0; i < numberRows; i++) { |
---|
296 | int iRow = whichRow[i]; |
---|
297 | if (rowLower_[i] > -1.0e20) |
---|
298 | rowLower_[i] -= rhsModification[iRow]; |
---|
299 | if (rowUpper_[i] < 1.0e20) |
---|
300 | rowUpper_[i] -= rhsModification[iRow]; |
---|
301 | } |
---|
302 | delete[] rhsModification; |
---|
303 | setObjectiveOffset(rhs->objectiveOffset() - offset); |
---|
304 | // And set objective value to match |
---|
305 | setObjectiveValue(rhs->objectiveValue()); |
---|
306 | } |
---|
307 | delete[] solution; |
---|
308 | } |
---|
309 | } |
---|
310 | // Subproblem constructor |
---|
311 | ClpSimplex::ClpSimplex(const ClpSimplex *rhs, |
---|
312 | int numberRows, const int *whichRow, |
---|
313 | int numberColumns, const int *whichColumn, |
---|
314 | bool dropNames, bool dropIntegers, bool fixOthers) |
---|
315 | : ClpModel(rhs, numberRows, whichRow, |
---|
316 | numberColumns, whichColumn, dropNames, dropIntegers) |
---|
317 | , bestPossibleImprovement_(0.0) |
---|
318 | , zeroTolerance_(1.0e-13) |
---|
319 | , columnPrimalSequence_(-2) |
---|
320 | , rowPrimalSequence_(-2) |
---|
321 | , bestObjectiveValue_(-COIN_DBL_MAX) |
---|
322 | , moreSpecialOptions_(2) |
---|
323 | , baseIteration_(0) |
---|
324 | , vectorMode_(0) |
---|
325 | , primalToleranceToGetOptimal_(-1.0) |
---|
326 | , largeValue_(1.0e15) |
---|
327 | , largestPrimalError_(0.0) |
---|
328 | , largestDualError_(0.0) |
---|
329 | , alphaAccuracy_(-1.0) |
---|
330 | , dualBound_(1.0e10) |
---|
331 | , alpha_(0.0) |
---|
332 | , theta_(0.0) |
---|
333 | , lowerIn_(0.0) |
---|
334 | , valueIn_(0.0) |
---|
335 | , upperIn_(-COIN_DBL_MAX) |
---|
336 | , dualIn_(0.0) |
---|
337 | , lowerOut_(-1) |
---|
338 | , valueOut_(-1) |
---|
339 | , upperOut_(-1) |
---|
340 | , dualOut_(-1) |
---|
341 | , dualTolerance_(rhs->dualTolerance_) |
---|
342 | , primalTolerance_(rhs->primalTolerance_) |
---|
343 | , sumDualInfeasibilities_(0.0) |
---|
344 | , sumPrimalInfeasibilities_(0.0) |
---|
345 | , infeasibilityCost_(1.0e10) |
---|
346 | , sumOfRelaxedDualInfeasibilities_(0.0) |
---|
347 | , sumOfRelaxedPrimalInfeasibilities_(0.0) |
---|
348 | , acceptablePivot_(1.0e-8) |
---|
349 | , lower_(NULL) |
---|
350 | , rowLowerWork_(NULL) |
---|
351 | , columnLowerWork_(NULL) |
---|
352 | , upper_(NULL) |
---|
353 | , rowUpperWork_(NULL) |
---|
354 | , columnUpperWork_(NULL) |
---|
355 | , cost_(NULL) |
---|
356 | , rowObjectiveWork_(NULL) |
---|
357 | , objectiveWork_(NULL) |
---|
358 | , sequenceIn_(-1) |
---|
359 | , directionIn_(-1) |
---|
360 | , sequenceOut_(-1) |
---|
361 | , directionOut_(-1) |
---|
362 | , pivotRow_(-1) |
---|
363 | , lastGoodIteration_(-100) |
---|
364 | , dj_(NULL) |
---|
365 | , rowReducedCost_(NULL) |
---|
366 | , reducedCostWork_(NULL) |
---|
367 | , solution_(NULL) |
---|
368 | , rowActivityWork_(NULL) |
---|
369 | , columnActivityWork_(NULL) |
---|
370 | , numberDualInfeasibilities_(0) |
---|
371 | , numberDualInfeasibilitiesWithoutFree_(0) |
---|
372 | , numberPrimalInfeasibilities_(100) |
---|
373 | , numberRefinements_(0) |
---|
374 | , pivotVariable_(NULL) |
---|
375 | , factorization_(NULL) |
---|
376 | , savedSolution_(NULL) |
---|
377 | , numberTimesOptimal_(0) |
---|
378 | , disasterArea_(NULL) |
---|
379 | , changeMade_(1) |
---|
380 | , algorithm_(0) |
---|
381 | , forceFactorization_(-1) |
---|
382 | , perturbation_(100) |
---|
383 | , nonLinearCost_(NULL) |
---|
384 | , lastBadIteration_(-999999) |
---|
385 | , lastFlaggedIteration_(-999999) |
---|
386 | , numberFake_(0) |
---|
387 | , numberChanged_(0) |
---|
388 | , progressFlag_(0) |
---|
389 | , firstFree_(-1) |
---|
390 | , numberExtraRows_(0) |
---|
391 | , maximumBasic_(0) |
---|
392 | , dontFactorizePivots_(0) |
---|
393 | , incomingInfeasibility_(1.0) |
---|
394 | , allowedInfeasibility_(10.0) |
---|
395 | , automaticScale_(0) |
---|
396 | , maximumPerturbationSize_(0) |
---|
397 | , perturbationArray_(NULL) |
---|
398 | , baseModel_(NULL) |
---|
399 | #ifdef ABC_INHERIT |
---|
400 | , abcSimplex_(NULL) |
---|
401 | , abcState_(rhs->abcState_) |
---|
402 | #endif |
---|
403 | { |
---|
404 | int i; |
---|
405 | for (i = 0; i < 6; i++) { |
---|
406 | rowArray_[i] = NULL; |
---|
407 | columnArray_[i] = NULL; |
---|
408 | } |
---|
409 | for (i = 0; i < 4; i++) { |
---|
410 | spareIntArray_[i] = 0; |
---|
411 | spareDoubleArray_[i] = 0.0; |
---|
412 | } |
---|
413 | saveStatus_ = NULL; |
---|
414 | factorization_ = new ClpFactorization(*rhs->factorization_, -numberRows_); |
---|
415 | //factorization_ = new ClpFactorization(*rhs->factorization_, |
---|
416 | // rhs->factorization_->goDenseThreshold()); |
---|
417 | ClpPEDualRowSteepest *pivotDualPE = dynamic_cast< ClpPEDualRowSteepest * >(rhs->dualRowPivot_); |
---|
418 | if (pivotDualPE) { |
---|
419 | dualRowPivot_ = new ClpPEDualRowSteepest(pivotDualPE->psi()); |
---|
420 | } else { |
---|
421 | ClpDualRowDantzig *pivot = dynamic_cast< ClpDualRowDantzig * >(rhs->dualRowPivot_); |
---|
422 | // say Steepest pricing |
---|
423 | if (!pivot) |
---|
424 | dualRowPivot_ = new ClpDualRowSteepest(); |
---|
425 | else |
---|
426 | dualRowPivot_ = new ClpDualRowDantzig(); |
---|
427 | } |
---|
428 | ClpPEPrimalColumnSteepest *pivotPrimalPE = dynamic_cast< ClpPEPrimalColumnSteepest * >(rhs->primalColumnPivot_); |
---|
429 | if (pivotPrimalPE) { |
---|
430 | primalColumnPivot_ = new ClpPEPrimalColumnSteepest(pivotPrimalPE->psi()); |
---|
431 | } else { |
---|
432 | // say Steepest pricing |
---|
433 | primalColumnPivot_ = new ClpPrimalColumnSteepest(); |
---|
434 | } |
---|
435 | solveType_ = 1; // say simplex based life form |
---|
436 | if (fixOthers) { |
---|
437 | int numberOtherColumns = rhs->numberColumns(); |
---|
438 | int numberOtherRows = rhs->numberRows(); |
---|
439 | double *solution = new double[numberOtherColumns]; |
---|
440 | CoinZeroN(solution, numberOtherColumns); |
---|
441 | int i; |
---|
442 | for (i = 0; i < numberColumns; i++) { |
---|
443 | int iColumn = whichColumn[i]; |
---|
444 | if (solution[iColumn]) |
---|
445 | fixOthers = false; // duplicates |
---|
446 | solution[iColumn] = 1.0; |
---|
447 | } |
---|
448 | if (fixOthers) { |
---|
449 | const double *otherSolution = rhs->primalColumnSolution(); |
---|
450 | const double *objective = rhs->objective(); |
---|
451 | double offset = 0.0; |
---|
452 | for (i = 0; i < numberOtherColumns; i++) { |
---|
453 | if (solution[i]) { |
---|
454 | solution[i] = 0.0; // in |
---|
455 | } else { |
---|
456 | solution[i] = otherSolution[i]; |
---|
457 | offset += objective[i] * otherSolution[i]; |
---|
458 | } |
---|
459 | } |
---|
460 | double *rhsModification = new double[numberOtherRows]; |
---|
461 | CoinZeroN(rhsModification, numberOtherRows); |
---|
462 | rhs->matrix()->times(solution, rhsModification); |
---|
463 | for (i = 0; i < numberRows; i++) { |
---|
464 | int iRow = whichRow[i]; |
---|
465 | if (rowLower_[i] > -1.0e20) |
---|
466 | rowLower_[i] -= rhsModification[iRow]; |
---|
467 | if (rowUpper_[i] < 1.0e20) |
---|
468 | rowUpper_[i] -= rhsModification[iRow]; |
---|
469 | } |
---|
470 | delete[] rhsModification; |
---|
471 | setObjectiveOffset(rhs->objectiveOffset() - offset); |
---|
472 | // And set objective value to match |
---|
473 | setObjectiveValue(rhs->objectiveValue()); |
---|
474 | } |
---|
475 | delete[] solution; |
---|
476 | } |
---|
477 | if (rhs->maximumPerturbationSize_) { |
---|
478 | maximumPerturbationSize_ = 2 * numberColumns; |
---|
479 | perturbationArray_ = new double[maximumPerturbationSize_]; |
---|
480 | for (i = 0; i < numberColumns; i++) { |
---|
481 | int iColumn = whichColumn[i]; |
---|
482 | perturbationArray_[2 * i] = rhs->perturbationArray_[2 * iColumn]; |
---|
483 | perturbationArray_[2 * i + 1] = rhs->perturbationArray_[2 * iColumn + 1]; |
---|
484 | } |
---|
485 | } |
---|
486 | } |
---|
487 | // Puts solution back small model |
---|
488 | void ClpSimplex::getbackSolution(const ClpSimplex &smallModel, const int *whichRow, const int *whichColumn) |
---|
489 | { |
---|
490 | setSumDualInfeasibilities(smallModel.sumDualInfeasibilities()); |
---|
491 | setNumberDualInfeasibilities(smallModel.numberDualInfeasibilities()); |
---|
492 | setSumPrimalInfeasibilities(smallModel.sumPrimalInfeasibilities()); |
---|
493 | setNumberPrimalInfeasibilities(smallModel.numberPrimalInfeasibilities()); |
---|
494 | setNumberIterations(smallModel.numberIterations()); |
---|
495 | setProblemStatus(smallModel.status()); |
---|
496 | setObjectiveValue(smallModel.objectiveValue()); |
---|
497 | const double *solution2 = smallModel.primalColumnSolution(); |
---|
498 | int i; |
---|
499 | int numberRows2 = smallModel.numberRows(); |
---|
500 | int numberColumns2 = smallModel.numberColumns(); |
---|
501 | const double *dj2 = smallModel.dualColumnSolution(); |
---|
502 | for (i = 0; i < numberColumns2; i++) { |
---|
503 | int iColumn = whichColumn[i]; |
---|
504 | columnActivity_[iColumn] = solution2[i]; |
---|
505 | reducedCost_[iColumn] = dj2[i]; |
---|
506 | setStatus(iColumn, smallModel.getStatus(i)); |
---|
507 | } |
---|
508 | const double *dual2 = smallModel.dualRowSolution(); |
---|
509 | memset(dual_, 0, numberRows_ * sizeof(double)); |
---|
510 | for (i = 0; i < numberRows2; i++) { |
---|
511 | int iRow = whichRow[i]; |
---|
512 | setRowStatus(iRow, smallModel.getRowStatus(i)); |
---|
513 | dual_[iRow] = dual2[i]; |
---|
514 | } |
---|
515 | CoinZeroN(rowActivity_, numberRows_); |
---|
516 | #if 0 |
---|
517 | if (!problemStatus_) { |
---|
518 | ClpDisjointCopyN(smallModel.objective(), smallModel.numberColumns_, smallModel.reducedCost_); |
---|
519 | smallModel.matrix_->transposeTimes(-1.0, smallModel.dual_, smallModel.reducedCost_); |
---|
520 | for (int i = 0; i < smallModel.numberColumns_; i++) { |
---|
521 | if (smallModel.getColumnStatus(i) == basic) |
---|
522 | assert (fabs(smallModel.reducedCost_[i]) < 1.0e-5); |
---|
523 | } |
---|
524 | ClpDisjointCopyN(objective(), numberColumns_, reducedCost_); |
---|
525 | matrix_->transposeTimes(-1.0, dual_, reducedCost_); |
---|
526 | for (int i = 0; i < numberColumns_; i++) { |
---|
527 | if (getColumnStatus(i) == basic) |
---|
528 | assert (fabs(reducedCost_[i]) < 1.0e-5); |
---|
529 | } |
---|
530 | } |
---|
531 | #endif |
---|
532 | matrix()->times(columnActivity_, rowActivity_); |
---|
533 | } |
---|
534 | |
---|
535 | //----------------------------------------------------------------------------- |
---|
536 | |
---|
537 | ClpSimplex::~ClpSimplex() |
---|
538 | { |
---|
539 | setPersistenceFlag(0); |
---|
540 | gutsOfDelete(0); |
---|
541 | delete nonLinearCost_; |
---|
542 | } |
---|
543 | //############################################################################# |
---|
544 | void ClpSimplex::setLargeValue(double value) |
---|
545 | { |
---|
546 | if (value > 0.0 && value < COIN_DBL_MAX) |
---|
547 | largeValue_ = value; |
---|
548 | } |
---|
549 | int ClpSimplex::gutsOfSolution(double *givenDuals, |
---|
550 | const double *givenPrimals, |
---|
551 | bool valuesPass) |
---|
552 | { |
---|
553 | |
---|
554 | // if values pass, save values of basic variables |
---|
555 | double *save = NULL; |
---|
556 | double oldValue = 0.0; |
---|
557 | double oldLargestPrimalError = largestPrimalError_; |
---|
558 | double oldLargestDualError = largestDualError_; |
---|
559 | if (valuesPass) { |
---|
560 | assert(algorithm_ > 0); // only primal at present |
---|
561 | assert(nonLinearCost_); |
---|
562 | int iRow; |
---|
563 | checkPrimalSolution(rowActivityWork_, columnActivityWork_); |
---|
564 | // get correct bounds on all variables |
---|
565 | nonLinearCost_->checkInfeasibilities(primalTolerance_); |
---|
566 | oldValue = nonLinearCost_->largestInfeasibility(); |
---|
567 | save = new double[numberRows_]; |
---|
568 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
569 | int iPivot = pivotVariable_[iRow]; |
---|
570 | save[iRow] = solution_[iPivot]; |
---|
571 | } |
---|
572 | } |
---|
573 | // do work |
---|
574 | computePrimals(rowActivityWork_, columnActivityWork_); |
---|
575 | // If necessary - override results |
---|
576 | if (givenPrimals) { |
---|
577 | CoinMemcpyN(givenPrimals, numberColumns_, columnActivityWork_); |
---|
578 | memset(rowActivityWork_, 0, numberRows_ * sizeof(double)); |
---|
579 | times(-1.0, columnActivityWork_, rowActivityWork_); |
---|
580 | } |
---|
581 | double objectiveModification = 0.0; |
---|
582 | if (algorithm_ > 0 && nonLinearCost_ != NULL) { |
---|
583 | #ifdef CLP_USER_DRIVEN |
---|
584 | eventHandler_->eventWithInfo(ClpEventHandler::goodFactorization, NULL); |
---|
585 | #endif |
---|
586 | // primal algorithm |
---|
587 | // get correct bounds on all variables |
---|
588 | // If 4 bit set - Force outgoing variables to exact bound (primal) |
---|
589 | if ((specialOptions_ & 4) == 0) |
---|
590 | nonLinearCost_->checkInfeasibilities(primalTolerance_); |
---|
591 | else |
---|
592 | nonLinearCost_->checkInfeasibilities(0.0); |
---|
593 | objectiveModification += nonLinearCost_->changeInCost(); |
---|
594 | if (nonLinearCost_->numberInfeasibilities()) |
---|
595 | if (handler_->detail(CLP_SIMPLEX_NONLINEAR, messages_) < 100) { |
---|
596 | handler_->message(CLP_SIMPLEX_NONLINEAR, messages_) |
---|
597 | << nonLinearCost_->changeInCost() |
---|
598 | << nonLinearCost_->numberInfeasibilities() |
---|
599 | << CoinMessageEol; |
---|
600 | } |
---|
601 | } |
---|
602 | if (valuesPass) { |
---|
603 | double badInfeasibility = nonLinearCost_->largestInfeasibility(); |
---|
604 | #ifdef CLP_DEBUG |
---|
605 | std::cout << "Largest given infeasibility " << oldValue |
---|
606 | << " now " << nonLinearCost_->largestInfeasibility() << std::endl; |
---|
607 | #endif |
---|
608 | int numberOut = 0; |
---|
609 | // But may be very large rhs etc |
---|
610 | double useError = CoinMin(largestPrimalError_, |
---|
611 | 1.0e5 / maximumAbsElement(solution_, numberRows_ + numberColumns_)); |
---|
612 | if ((oldValue < incomingInfeasibility_ || badInfeasibility > (CoinMax(10.0 * allowedInfeasibility_, 100.0 * oldValue))) |
---|
613 | && (badInfeasibility > CoinMax(incomingInfeasibility_, allowedInfeasibility_) || useError > 1.0e-3)) { |
---|
614 | if (algorithm_ > 1) { |
---|
615 | // nonlinear |
---|
616 | //printf("Original largest infeas %g, now %g, primalError %g\n", |
---|
617 | // oldValue,nonLinearCost_->largestInfeasibility(), |
---|
618 | // largestPrimalError_); |
---|
619 | //printf("going to all slack\n"); |
---|
620 | allSlackBasis(true); |
---|
621 | CoinIotaN(pivotVariable_, numberRows_, numberColumns_); |
---|
622 | return 1; |
---|
623 | } |
---|
624 | //printf("Original largest infeas %g, now %g, primalError %g\n", |
---|
625 | // oldValue,nonLinearCost_->largestInfeasibility(), |
---|
626 | // largestPrimalError_); |
---|
627 | // throw out up to 1000 structurals |
---|
628 | int maxOut = (allowedInfeasibility_ == 10.0) ? 1000 : 100; |
---|
629 | int iRow; |
---|
630 | int *sort = new int[numberRows_]; |
---|
631 | // first put back solution and store difference |
---|
632 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
633 | int iPivot = pivotVariable_[iRow]; |
---|
634 | double difference = fabs(solution_[iPivot] - save[iRow]); |
---|
635 | solution_[iPivot] = save[iRow]; |
---|
636 | save[iRow] = difference; |
---|
637 | } |
---|
638 | int numberBasic = 0; |
---|
639 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
640 | int iPivot = pivotVariable_[iRow]; |
---|
641 | |
---|
642 | if (iPivot < numberColumns_) { |
---|
643 | // column |
---|
644 | double difference = save[iRow]; |
---|
645 | if (difference > 1.0e-4) { |
---|
646 | sort[numberOut] = iRow; |
---|
647 | save[numberOut++] = -difference; |
---|
648 | if (getStatus(iPivot) == basic) |
---|
649 | numberBasic++; |
---|
650 | } |
---|
651 | } |
---|
652 | } |
---|
653 | if (!numberBasic) { |
---|
654 | //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut); |
---|
655 | #if 0 |
---|
656 | allSlackBasis(true); |
---|
657 | CoinIotaN(pivotVariable_, numberRows_, numberColumns_); |
---|
658 | #else |
---|
659 | // allow |
---|
660 | numberOut = 0; |
---|
661 | #endif |
---|
662 | } |
---|
663 | CoinSort_2(save, save + numberOut, sort); |
---|
664 | numberOut = CoinMin(maxOut, numberOut); |
---|
665 | for (iRow = 0; iRow < numberOut; iRow++) { |
---|
666 | int jRow = sort[iRow]; |
---|
667 | int iColumn = pivotVariable_[jRow]; |
---|
668 | setColumnStatus(iColumn, superBasic); |
---|
669 | setRowStatus(jRow, basic); |
---|
670 | pivotVariable_[jRow] = jRow + numberColumns_; |
---|
671 | if (fabs(solution_[iColumn]) > 1.0e10) { |
---|
672 | if (upper_[iColumn] < 0.0) { |
---|
673 | solution_[iColumn] = upper_[iColumn]; |
---|
674 | } else if (lower_[iColumn] > 0.0) { |
---|
675 | solution_[iColumn] = lower_[iColumn]; |
---|
676 | } else { |
---|
677 | solution_[iColumn] = 0.0; |
---|
678 | } |
---|
679 | } |
---|
680 | } |
---|
681 | delete[] sort; |
---|
682 | } |
---|
683 | delete[] save; |
---|
684 | if (numberOut) |
---|
685 | return numberOut; |
---|
686 | } |
---|
687 | if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) { |
---|
688 | //printf("trying feas pump\n"); |
---|
689 | const char *integerType = integerInformation(); |
---|
690 | assert(integerType); |
---|
691 | assert(perturbationArray_); |
---|
692 | CoinZeroN(cost_, numberRows_ + numberColumns_); |
---|
693 | for (int i = 0; i < numberRows_ - numberRows_; i++) { |
---|
694 | int iSequence = pivotVariable_[i]; |
---|
695 | if (iSequence < numberColumns_ && integerType[iSequence]) { |
---|
696 | double lower = lower_[iSequence]; |
---|
697 | double upper = upper_[iSequence]; |
---|
698 | double value = solution_[iSequence]; |
---|
699 | if (value >= lower - primalTolerance_ && value <= upper + primalTolerance_) { |
---|
700 | double sign; |
---|
701 | if (value - lower < upper - value) |
---|
702 | sign = 1.0; |
---|
703 | else |
---|
704 | sign = -1.0; |
---|
705 | cost_[iSequence] = sign * perturbationArray_[iSequence]; |
---|
706 | } |
---|
707 | } |
---|
708 | } |
---|
709 | } |
---|
710 | #if CAN_HAVE_ZERO_OBJ > 1 |
---|
711 | if ((specialOptions_ & 16777216) == 0) { |
---|
712 | #endif |
---|
713 | computeDuals(givenDuals); |
---|
714 | if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) { |
---|
715 | const char *integerType = integerInformation(); |
---|
716 | // Need to do columns and rows to stay dual feasible |
---|
717 | for (int iSequence = 0; iSequence < numberColumns_; iSequence++) { |
---|
718 | if (integerType[iSequence] && getStatus(iSequence) != basic) { |
---|
719 | double djValue = dj_[iSequence]; |
---|
720 | double change = 0.0; |
---|
721 | if (getStatus(iSequence) == atLowerBound) |
---|
722 | change = CoinMax(-djValue, 10.0 * perturbationArray_[iSequence]); |
---|
723 | else if (getStatus(iSequence) == atUpperBound) |
---|
724 | change = CoinMin(-djValue, -10.0 * perturbationArray_[iSequence]); |
---|
725 | cost_[iSequence] = change; |
---|
726 | dj_[iSequence] += change; |
---|
727 | } |
---|
728 | } |
---|
729 | } |
---|
730 | |
---|
731 | // now check solutions |
---|
732 | //checkPrimalSolution( rowActivityWork_, columnActivityWork_); |
---|
733 | //checkDualSolution(); |
---|
734 | checkBothSolutions(); |
---|
735 | objectiveValue_ += objectiveModification / (objectiveScale_ * rhsScale_); |
---|
736 | #if CAN_HAVE_ZERO_OBJ > 1 |
---|
737 | } else { |
---|
738 | checkPrimalSolution(rowActivityWork_, columnActivityWork_); |
---|
739 | #ifndef COIN_REUSE_RANDOM |
---|
740 | memset(dj_, 0, (numberRows_ + numberColumns_) * sizeof(double)); |
---|
741 | #else |
---|
742 | for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { |
---|
743 | double value; |
---|
744 | switch (getStatus(iSequence)) { |
---|
745 | case atLowerBound: |
---|
746 | value = 1.0e-9 * (1.0 + randomNumberGenerator_.randomDouble()); |
---|
747 | break; |
---|
748 | case atUpperBound: |
---|
749 | value = -1.0e-9 * (1.0 + randomNumberGenerator_.randomDouble()); |
---|
750 | break; |
---|
751 | default: |
---|
752 | value = 0.0; |
---|
753 | break; |
---|
754 | } |
---|
755 | } |
---|
756 | #endif |
---|
757 | objectiveValue_ = 0.0; |
---|
758 | } |
---|
759 | #endif |
---|
760 | if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 || largestDualError_ > 1.0e-2)) |
---|
761 | handler_->message(CLP_SIMPLEX_ACCURACY, messages_) |
---|
762 | << largestPrimalError_ |
---|
763 | << largestDualError_ |
---|
764 | << CoinMessageEol; |
---|
765 | if (largestPrimalError_ > 1.0e-1 && numberRows_ > 100 && numberIterations_) { |
---|
766 | // Change factorization tolerance |
---|
767 | if (factorization_->zeroTolerance() > 1.0e-18) |
---|
768 | factorization_->zeroTolerance(1.0e-18); |
---|
769 | } |
---|
770 | int returnCode = 0; |
---|
771 | bool notChanged = true; |
---|
772 | if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_ < 0 || factorization_->pivotTolerance() < 0.9899999999) && (oldLargestDualError || oldLargestPrimalError)) { |
---|
773 | double useOldDualError = oldLargestDualError; |
---|
774 | double useDualError = largestDualError_; |
---|
775 | if (algorithm_ > 0 && nonLinearCost_ && nonLinearCost_->sumInfeasibilities()) { |
---|
776 | double factor = CoinMax(1.0, CoinMin(1.0e3, infeasibilityCost_ * 1.0e-6)); |
---|
777 | useOldDualError /= factor; |
---|
778 | useDualError /= factor; |
---|
779 | } |
---|
780 | if ((largestPrimalError_ > 1.0e3 && oldLargestPrimalError * 1.0e2 < largestPrimalError_) || (useDualError > 1.0e3 && useOldDualError * 1.0e2 < useDualError)) { |
---|
781 | double pivotTolerance = factorization_->pivotTolerance(); |
---|
782 | double factor = (largestPrimalError_ > 1.0e10 || largestDualError_ > 1.0e10) |
---|
783 | ? 2.0 |
---|
784 | : 1.2; |
---|
785 | if (pivotTolerance < 0.1) |
---|
786 | factorization_->pivotTolerance(0.1); |
---|
787 | else if (pivotTolerance < 0.98999999) |
---|
788 | factorization_->pivotTolerance(CoinMin(0.99, pivotTolerance * factor)); |
---|
789 | notChanged = pivotTolerance == factorization_->pivotTolerance(); |
---|
790 | #ifdef CLP_USEFUL_PRINTOUT |
---|
791 | if (pivotTolerance < 0.9899999) { |
---|
792 | printf("Changing pivot tolerance from %g to %g and backtracking\n", |
---|
793 | pivotTolerance, factorization_->pivotTolerance()); |
---|
794 | } |
---|
795 | printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n", |
---|
796 | oldLargestPrimalError, largestPrimalError_, |
---|
797 | oldLargestDualError, largestDualError_, |
---|
798 | pivotTolerance); |
---|
799 | #endif |
---|
800 | if (pivotTolerance < 0.9899999) { |
---|
801 | largestPrimalError_ = 0.0; |
---|
802 | largestDualError_ = 0.0; |
---|
803 | returnCode = -123456789; |
---|
804 | } |
---|
805 | } |
---|
806 | } |
---|
807 | if (progress_.iterationNumber_[0] > 0 && progress_.iterationNumber_[CLP_PROGRESS - 1] - progress_.iterationNumber_[0] < CLP_PROGRESS * 3 && factorization_->pivotTolerance() < 0.25 && notChanged) { |
---|
808 | double pivotTolerance = factorization_->pivotTolerance(); |
---|
809 | factorization_->pivotTolerance(pivotTolerance * 1.5); |
---|
810 | #ifdef CLP_USEFUL_PRINTOUT |
---|
811 | printf("Changing pivot tolerance from %g to %g - inverting too often\n", |
---|
812 | pivotTolerance, factorization_->pivotTolerance()); |
---|
813 | #endif |
---|
814 | } |
---|
815 | // Switch off false values pass indicator |
---|
816 | if (!valuesPass && algorithm_ > 0) |
---|
817 | firstFree_ = -1; |
---|
818 | if (handler_->logLevel() == 63) |
---|
819 | printf("end getsolution algorithm %d status %d npinf %d sum,relaxed %g,%g ndinf %d sum,relaxed %g,%g\n", |
---|
820 | algorithm_, problemStatus_, |
---|
821 | numberPrimalInfeasibilities_, sumPrimalInfeasibilities_, sumOfRelaxedPrimalInfeasibilities_, |
---|
822 | numberDualInfeasibilities_, sumDualInfeasibilities_, sumOfRelaxedDualInfeasibilities_); |
---|
823 | if ((moreSpecialOptions_ & 8388608) != 0) { |
---|
824 | if (algorithm_ < 0) { |
---|
825 | bool doneFiddling = false; |
---|
826 | // Optimization may make exact test iffy |
---|
827 | double testTolerance = minimumPrimalTolerance_ + 1.0e-15; |
---|
828 | while (!doneFiddling) { |
---|
829 | doneFiddling = true; |
---|
830 | while (!sumOfRelaxedPrimalInfeasibilities_ && primalTolerance_ > testTolerance) { |
---|
831 | // feasible - adjust tolerance |
---|
832 | double saveTolerance = primalTolerance_; |
---|
833 | primalTolerance_ = CoinMax(0.25 * primalTolerance_, |
---|
834 | minimumPrimalTolerance_); |
---|
835 | printf("Resetting primal tolerance from %g to %g\n", |
---|
836 | saveTolerance, primalTolerance_); |
---|
837 | dblParam_[ClpPrimalTolerance] = primalTolerance_; |
---|
838 | moreSpecialOptions_ &= ~8388608; |
---|
839 | // redo with switch off |
---|
840 | returnCode = gutsOfSolution(givenDuals, givenPrimals, valuesPass); |
---|
841 | } |
---|
842 | if (primalTolerance_ > testTolerance) |
---|
843 | moreSpecialOptions_ |= 8388608; // back on |
---|
844 | if ((moreSpecialOptions_ & 8388608) != 0) { |
---|
845 | assert(numberPrimalInfeasibilities_); |
---|
846 | // average infeasibility |
---|
847 | double average = sumPrimalInfeasibilities_ / numberPrimalInfeasibilities_; |
---|
848 | double minimum = COIN_DBL_MAX; |
---|
849 | double averageTotal = average; |
---|
850 | bool firstTime = averageInfeasibility_[0] == COIN_DBL_MAX; |
---|
851 | for (int i = 0; i < CLP_INFEAS_SAVE - 1; i++) { |
---|
852 | double value = averageInfeasibility_[i + 1]; |
---|
853 | averageTotal += value; |
---|
854 | averageInfeasibility_[i] = value; |
---|
855 | minimum = CoinMin(minimum, value); |
---|
856 | } |
---|
857 | averageInfeasibility_[CLP_INFEAS_SAVE - 1] = average; |
---|
858 | averageTotal /= CLP_INFEAS_SAVE; |
---|
859 | double oldTolerance = primalTolerance_; |
---|
860 | if (averageInfeasibility_[0] != COIN_DBL_MAX) { |
---|
861 | if (firstTime) { |
---|
862 | primalTolerance_ = CoinMin(0.1, 0.1 * averageTotal); |
---|
863 | primalTolerance_ = CoinMin(primalTolerance_, average); |
---|
864 | } else if (primalTolerance_ > 0.1 * minimum) { |
---|
865 | primalTolerance_ = 0.1 * minimum; |
---|
866 | } |
---|
867 | primalTolerance_ = CoinMax(primalTolerance_, minimumPrimalTolerance_); |
---|
868 | } |
---|
869 | if (primalTolerance_ != oldTolerance) { |
---|
870 | printf("Changing primal tolerance from %g to %g\n", |
---|
871 | oldTolerance, primalTolerance_); |
---|
872 | moreSpecialOptions_ &= ~8388608; |
---|
873 | // redo with switch off |
---|
874 | returnCode = gutsOfSolution(givenDuals, givenPrimals, valuesPass); |
---|
875 | if (primalTolerance_ > testTolerance) |
---|
876 | moreSpecialOptions_ |= 8388608 | 4194304; |
---|
877 | if (!sumOfRelaxedPrimalInfeasibilities_) |
---|
878 | doneFiddling = false; // over done it |
---|
879 | } |
---|
880 | } |
---|
881 | } |
---|
882 | } |
---|
883 | } |
---|
884 | return returnCode; |
---|
885 | } |
---|
886 | void ClpSimplex::computePrimals(const double *rowActivities, |
---|
887 | const double *columnActivities) |
---|
888 | { |
---|
889 | |
---|
890 | //work space |
---|
891 | CoinIndexedVector *workSpace = rowArray_[0]; |
---|
892 | |
---|
893 | CoinIndexedVector *arrayVector = rowArray_[1]; |
---|
894 | arrayVector->clear(); |
---|
895 | CoinIndexedVector *previousVector = rowArray_[2]; |
---|
896 | previousVector->clear(); |
---|
897 | // accumulate non basic stuff |
---|
898 | |
---|
899 | int iRow; |
---|
900 | // order is this way for scaling |
---|
901 | if (columnActivities != columnActivityWork_) |
---|
902 | ClpDisjointCopyN(columnActivities, numberColumns_, columnActivityWork_); |
---|
903 | if (rowActivities != rowActivityWork_) |
---|
904 | ClpDisjointCopyN(rowActivities, numberRows_, rowActivityWork_); |
---|
905 | double *array = arrayVector->denseVector(); |
---|
906 | int *index = arrayVector->getIndices(); |
---|
907 | int number = 0; |
---|
908 | const double *rhsOffset = matrix_->rhsOffset(this, false, true); |
---|
909 | if (!rhsOffset) { |
---|
910 | // Use whole matrix every time to make it easier for ClpMatrixBase |
---|
911 | // So zero out basic |
---|
912 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
913 | int iPivot = pivotVariable_[iRow]; |
---|
914 | assert(iPivot >= 0); |
---|
915 | solution_[iPivot] = 0.0; |
---|
916 | #ifdef CLP_INVESTIGATE |
---|
917 | assert(getStatus(iPivot) == basic); |
---|
918 | #endif |
---|
919 | } |
---|
920 | // Extended solution before "update" |
---|
921 | matrix_->primalExpanded(this, 0); |
---|
922 | times(-1.0, columnActivityWork_, array); |
---|
923 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
924 | double value = array[iRow] + rowActivityWork_[iRow]; |
---|
925 | if (value) { |
---|
926 | array[iRow] = value; |
---|
927 | index[number++] = iRow; |
---|
928 | } else { |
---|
929 | array[iRow] = 0.0; |
---|
930 | } |
---|
931 | } |
---|
932 | } else { |
---|
933 | // we have an effective rhs lying around |
---|
934 | // zero out basic (really just for slacks) |
---|
935 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
936 | int iPivot = pivotVariable_[iRow]; |
---|
937 | solution_[iPivot] = 0.0; |
---|
938 | } |
---|
939 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
940 | double value = rhsOffset[iRow] + rowActivityWork_[iRow]; |
---|
941 | if (value) { |
---|
942 | array[iRow] = value; |
---|
943 | index[number++] = iRow; |
---|
944 | } else { |
---|
945 | array[iRow] = 0.0; |
---|
946 | } |
---|
947 | } |
---|
948 | } |
---|
949 | arrayVector->setNumElements(number); |
---|
950 | #ifdef CLP_DEBUG |
---|
951 | if (numberIterations_ == -3840) { |
---|
952 | int i; |
---|
953 | for (i = 0; i < numberRows_ + numberColumns_; i++) |
---|
954 | printf("%d status %d\n", i, status_[i]); |
---|
955 | printf("xxxxx1\n"); |
---|
956 | for (i = 0; i < numberRows_; i++) |
---|
957 | if (array[i]) |
---|
958 | printf("%d rhs %g\n", i, array[i]); |
---|
959 | printf("xxxxx2\n"); |
---|
960 | for (i = 0; i < numberRows_ + numberColumns_; i++) |
---|
961 | if (getStatus(i) != basic) |
---|
962 | printf("%d non basic %g %g %g\n", i, lower_[i], solution_[i], upper_[i]); |
---|
963 | printf("xxxxx3\n"); |
---|
964 | } |
---|
965 | #endif |
---|
966 | // Ftran adjusted RHS and iterate to improve accuracy |
---|
967 | double lastError = COIN_DBL_MAX; |
---|
968 | int iRefine; |
---|
969 | CoinIndexedVector *thisVector = arrayVector; |
---|
970 | CoinIndexedVector *lastVector = previousVector; |
---|
971 | #if 0 |
---|
972 | static double * xsave=NULL; |
---|
973 | { |
---|
974 | double * xx = thisVector->denseVector(); |
---|
975 | double largest=0.0; |
---|
976 | int iLargest=-1; |
---|
977 | for (int i=0;i<numberRows_;i++) { |
---|
978 | if (fabs(xx[i])>largest) { |
---|
979 | largest=fabs(xx[i]); |
---|
980 | iLargest=i; |
---|
981 | } |
---|
982 | } |
---|
983 | printf("largest incoming rhs %g on row %d\n",largest,iLargest); |
---|
984 | } |
---|
985 | if (numberIterations_<-40722) { |
---|
986 | double * xx = thisVector->denseVector(); |
---|
987 | if (xsave) { |
---|
988 | double * sol = xsave+numberRows_; |
---|
989 | double largest=0.0; |
---|
990 | int iLargest=-1; |
---|
991 | for (int i=0;i<numberRows_;i++) { |
---|
992 | if (fabs(xx[i]-xsave[i])>largest) { |
---|
993 | largest=fabs(xx[i]-xsave[i]); |
---|
994 | iLargest=i; |
---|
995 | } |
---|
996 | } |
---|
997 | printf("error %g on row %d\n",largest,iLargest); |
---|
998 | largest=0.0; |
---|
999 | iLargest=-1; |
---|
1000 | for (int i=0;i<numberColumns_;i++) { |
---|
1001 | if (fabs(solution_[i]-sol[i])>largest) { |
---|
1002 | largest=fabs(solution_[i]-sol[i]); |
---|
1003 | iLargest=i; |
---|
1004 | } |
---|
1005 | } |
---|
1006 | printf("error %g on col %d\n",largest,iLargest); |
---|
1007 | } else { |
---|
1008 | xsave=new double[2*numberRows_+numberColumns_]; |
---|
1009 | } |
---|
1010 | memcpy(xsave,xx,numberRows_*sizeof(double)); |
---|
1011 | memcpy(xsave+numberRows_,solution_,(numberRows_+numberColumns_)*sizeof(double)); |
---|
1012 | } |
---|
1013 | #endif |
---|
1014 | //printf ("ZZ0 n before %d",number); |
---|
1015 | if (number) |
---|
1016 | factorization_->updateColumn(workSpace, thisVector); |
---|
1017 | //printf(" - after %d\n",thisVector->getNumElements()); |
---|
1018 | double *work = workSpace->denseVector(); |
---|
1019 | #ifdef CLP_DEBUG |
---|
1020 | if (numberIterations_ == -3840) { |
---|
1021 | int i; |
---|
1022 | for (i = 0; i < numberRows_; i++) |
---|
1023 | if (array[i]) |
---|
1024 | printf("%d after rhs %g\n", i, array[i]); |
---|
1025 | printf("xxxxx4\n"); |
---|
1026 | } |
---|
1027 | #endif |
---|
1028 | bool goodSolution = true; |
---|
1029 | for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) { |
---|
1030 | |
---|
1031 | int numberIn = thisVector->getNumElements(); |
---|
1032 | int *indexIn = thisVector->getIndices(); |
---|
1033 | double *arrayIn = thisVector->denseVector(); |
---|
1034 | // put solution in correct place |
---|
1035 | if (!rhsOffset) { |
---|
1036 | int j; |
---|
1037 | for (j = 0; j < numberIn; j++) { |
---|
1038 | iRow = indexIn[j]; |
---|
1039 | int iPivot = pivotVariable_[iRow]; |
---|
1040 | solution_[iPivot] = arrayIn[iRow]; |
---|
1041 | //assert (fabs(solution_[iPivot])<1.0e100); |
---|
1042 | } |
---|
1043 | } else { |
---|
1044 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1045 | int iPivot = pivotVariable_[iRow]; |
---|
1046 | solution_[iPivot] = arrayIn[iRow]; |
---|
1047 | //assert (fabs(solution_[iPivot])<1.0e100); |
---|
1048 | } |
---|
1049 | } |
---|
1050 | // Extended solution after "update" |
---|
1051 | matrix_->primalExpanded(this, 1); |
---|
1052 | // check Ax == b (for all) |
---|
1053 | // signal column generated matrix to just do basic (and gub) |
---|
1054 | unsigned int saveOptions = specialOptions(); |
---|
1055 | setSpecialOptions(16); |
---|
1056 | times(-1.0, columnActivityWork_, work); |
---|
1057 | setSpecialOptions(saveOptions); |
---|
1058 | largestPrimalError_ = 0.0; |
---|
1059 | double multiplier = 131072.0; |
---|
1060 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1061 | double value = work[iRow] + rowActivityWork_[iRow]; |
---|
1062 | work[iRow] = value * multiplier; |
---|
1063 | if (fabs(value) > largestPrimalError_) { |
---|
1064 | largestPrimalError_ = fabs(value); |
---|
1065 | } |
---|
1066 | } |
---|
1067 | if (largestPrimalError_ >= lastError) { |
---|
1068 | // restore |
---|
1069 | CoinIndexedVector *temp = thisVector; |
---|
1070 | thisVector = lastVector; |
---|
1071 | lastVector = temp; |
---|
1072 | goodSolution = false; |
---|
1073 | break; |
---|
1074 | } |
---|
1075 | if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) { |
---|
1076 | // try and make better |
---|
1077 | // save this |
---|
1078 | CoinIndexedVector *temp = thisVector; |
---|
1079 | thisVector = lastVector; |
---|
1080 | lastVector = temp; |
---|
1081 | int *indexOut = thisVector->getIndices(); |
---|
1082 | int number = 0; |
---|
1083 | array = thisVector->denseVector(); |
---|
1084 | thisVector->clear(); |
---|
1085 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1086 | double value = work[iRow]; |
---|
1087 | if (value) { |
---|
1088 | array[iRow] = value; |
---|
1089 | indexOut[number++] = iRow; |
---|
1090 | work[iRow] = 0.0; |
---|
1091 | } |
---|
1092 | } |
---|
1093 | thisVector->setNumElements(number); |
---|
1094 | lastError = largestPrimalError_; |
---|
1095 | //printf ("ZZ%d n before %d",iRefine+1,number); |
---|
1096 | factorization_->updateColumn(workSpace, thisVector); |
---|
1097 | //printf(" - after %d\n",thisVector->getNumElements()); |
---|
1098 | multiplier = 1.0 / multiplier; |
---|
1099 | double *previous = lastVector->denseVector(); |
---|
1100 | number = 0; |
---|
1101 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1102 | double value = previous[iRow] + multiplier * array[iRow]; |
---|
1103 | if (value) { |
---|
1104 | array[iRow] = value; |
---|
1105 | indexOut[number++] = iRow; |
---|
1106 | } else { |
---|
1107 | array[iRow] = 0.0; |
---|
1108 | } |
---|
1109 | } |
---|
1110 | thisVector->setNumElements(number); |
---|
1111 | } else { |
---|
1112 | break; |
---|
1113 | } |
---|
1114 | } |
---|
1115 | |
---|
1116 | // solution as accurate as we are going to get |
---|
1117 | ClpFillN(work, numberRows_, 0.0); |
---|
1118 | if (!goodSolution) { |
---|
1119 | array = thisVector->denseVector(); |
---|
1120 | // put solution in correct place |
---|
1121 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1122 | int iPivot = pivotVariable_[iRow]; |
---|
1123 | solution_[iPivot] = array[iRow]; |
---|
1124 | //assert (fabs(solution_[iPivot])<1.0e100); |
---|
1125 | } |
---|
1126 | } |
---|
1127 | arrayVector->clear(); |
---|
1128 | previousVector->clear(); |
---|
1129 | #ifdef CLP_DEBUG |
---|
1130 | if (numberIterations_ == -3840) { |
---|
1131 | exit(77); |
---|
1132 | } |
---|
1133 | #endif |
---|
1134 | } |
---|
1135 | // now dual side |
---|
1136 | void ClpSimplex::computeDuals(double *givenDjs) |
---|
1137 | { |
---|
1138 | #ifndef SLIM_CLP |
---|
1139 | if (objective_->type() == 1 || !objective_->activated()) { |
---|
1140 | #endif |
---|
1141 | // Linear |
---|
1142 | //work space |
---|
1143 | CoinIndexedVector *workSpace = rowArray_[0]; |
---|
1144 | |
---|
1145 | CoinIndexedVector *arrayVector = rowArray_[1]; |
---|
1146 | arrayVector->clear(); |
---|
1147 | CoinIndexedVector *previousVector = rowArray_[2]; |
---|
1148 | previousVector->clear(); |
---|
1149 | int iRow; |
---|
1150 | #ifdef CLP_DEBUG |
---|
1151 | workSpace->checkClear(); |
---|
1152 | #endif |
---|
1153 | double *array = arrayVector->denseVector(); |
---|
1154 | int *index = arrayVector->getIndices(); |
---|
1155 | int number = 0; |
---|
1156 | if (!givenDjs) { |
---|
1157 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1158 | int iPivot = pivotVariable_[iRow]; |
---|
1159 | double value = cost_[iPivot]; |
---|
1160 | if (value) { |
---|
1161 | array[iRow] = value; |
---|
1162 | index[number++] = iRow; |
---|
1163 | } |
---|
1164 | } |
---|
1165 | } else { |
---|
1166 | // dual values pass - djs may not be zero |
---|
1167 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1168 | int iPivot = pivotVariable_[iRow]; |
---|
1169 | // make sure zero if done |
---|
1170 | if (!pivoted(iPivot)) |
---|
1171 | givenDjs[iPivot] = 0.0; |
---|
1172 | double value = cost_[iPivot] - givenDjs[iPivot]; |
---|
1173 | if (value) { |
---|
1174 | array[iRow] = value; |
---|
1175 | index[number++] = iRow; |
---|
1176 | } |
---|
1177 | } |
---|
1178 | } |
---|
1179 | arrayVector->setNumElements(number); |
---|
1180 | // Extended duals before "updateTranspose" |
---|
1181 | matrix_->dualExpanded(this, arrayVector, givenDjs, 0); |
---|
1182 | |
---|
1183 | // Btran basic costs and get as accurate as possible |
---|
1184 | double lastError = COIN_DBL_MAX; |
---|
1185 | int iRefine; |
---|
1186 | double *work = workSpace->denseVector(); |
---|
1187 | CoinIndexedVector *thisVector = arrayVector; |
---|
1188 | CoinIndexedVector *lastVector = previousVector; |
---|
1189 | factorization_->updateColumnTranspose(workSpace, thisVector); |
---|
1190 | |
---|
1191 | for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) { |
---|
1192 | // check basic reduced costs zero |
---|
1193 | largestDualError_ = 0.0; |
---|
1194 | if (!numberExtraRows_) { |
---|
1195 | // Just basic |
---|
1196 | int *index2 = workSpace->getIndices(); |
---|
1197 | // use reduced costs for slacks as work array |
---|
1198 | double *work2 = reducedCostWork_ + numberColumns_; |
---|
1199 | int numberStructurals = 0; |
---|
1200 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1201 | int iPivot = pivotVariable_[iRow]; |
---|
1202 | if (iPivot < numberColumns_) |
---|
1203 | index2[numberStructurals++] = iPivot; |
---|
1204 | } |
---|
1205 | matrix_->listTransposeTimes(this, array, index2, numberStructurals, work2); |
---|
1206 | numberStructurals = 0; |
---|
1207 | if (!givenDjs) { |
---|
1208 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1209 | int iPivot = pivotVariable_[iRow]; |
---|
1210 | double value; |
---|
1211 | if (iPivot >= numberColumns_) { |
---|
1212 | // slack |
---|
1213 | value = rowObjectiveWork_[iPivot - numberColumns_] |
---|
1214 | + array[iPivot - numberColumns_]; |
---|
1215 | } else { |
---|
1216 | // column |
---|
1217 | value = objectiveWork_[iPivot] - work2[numberStructurals++]; |
---|
1218 | } |
---|
1219 | work[iRow] = value; |
---|
1220 | if (fabs(value) > largestDualError_) { |
---|
1221 | largestDualError_ = fabs(value); |
---|
1222 | } |
---|
1223 | } |
---|
1224 | } else { |
---|
1225 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1226 | int iPivot = pivotVariable_[iRow]; |
---|
1227 | if (iPivot >= numberColumns_) { |
---|
1228 | // slack |
---|
1229 | work[iRow] = rowObjectiveWork_[iPivot - numberColumns_] |
---|
1230 | + array[iPivot - numberColumns_] - givenDjs[iPivot]; |
---|
1231 | } else { |
---|
1232 | // column |
---|
1233 | work[iRow] = objectiveWork_[iPivot] - work2[numberStructurals++] |
---|
1234 | - givenDjs[iPivot]; |
---|
1235 | } |
---|
1236 | if (fabs(work[iRow]) > largestDualError_) { |
---|
1237 | largestDualError_ = fabs(work[iRow]); |
---|
1238 | //assert (largestDualError_<1.0e-7); |
---|
1239 | //if (largestDualError_>1.0e-7) |
---|
1240 | //printf("large dual error %g\n",largestDualError_); |
---|
1241 | } |
---|
1242 | } |
---|
1243 | } |
---|
1244 | } else { |
---|
1245 | // extra rows - be more careful |
---|
1246 | #if 1 |
---|
1247 | // would be faster to do just for basic but this reduces code |
---|
1248 | ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_); |
---|
1249 | transposeTimes(-1.0, array, reducedCostWork_); |
---|
1250 | #else |
---|
1251 | // Just basic |
---|
1252 | int *index2 = workSpace->getIndices(); |
---|
1253 | int numberStructurals = 0; |
---|
1254 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1255 | int iPivot = pivotVariable_[iRow]; |
---|
1256 | if (iPivot < numberColumns_) |
---|
1257 | index2[numberStructurals++] = iPivot; |
---|
1258 | } |
---|
1259 | matrix_->listTransposeTimes(this, array, index2, numberStructurals, work); |
---|
1260 | for (iRow = 0; iRow < numberStructurals; iRow++) { |
---|
1261 | int iPivot = index2[iRow]; |
---|
1262 | reducedCostWork_[iPivot] = objectiveWork_[iPivot] - work[iRow]; |
---|
1263 | } |
---|
1264 | #endif |
---|
1265 | // update by duals on sets |
---|
1266 | matrix_->dualExpanded(this, NULL, NULL, 1); |
---|
1267 | if (!givenDjs) { |
---|
1268 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1269 | int iPivot = pivotVariable_[iRow]; |
---|
1270 | double value; |
---|
1271 | if (iPivot >= numberColumns_) { |
---|
1272 | // slack |
---|
1273 | value = rowObjectiveWork_[iPivot - numberColumns_] |
---|
1274 | + array[iPivot - numberColumns_]; |
---|
1275 | } else { |
---|
1276 | // column |
---|
1277 | value = reducedCostWork_[iPivot]; |
---|
1278 | } |
---|
1279 | work[iRow] = value; |
---|
1280 | if (fabs(value) > largestDualError_) { |
---|
1281 | largestDualError_ = fabs(value); |
---|
1282 | } |
---|
1283 | } |
---|
1284 | } else { |
---|
1285 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1286 | int iPivot = pivotVariable_[iRow]; |
---|
1287 | if (iPivot >= numberColumns_) { |
---|
1288 | // slack |
---|
1289 | work[iRow] = rowObjectiveWork_[iPivot - numberColumns_] |
---|
1290 | + array[iPivot - numberColumns_] - givenDjs[iPivot]; |
---|
1291 | } else { |
---|
1292 | // column |
---|
1293 | work[iRow] = reducedCostWork_[iPivot] - givenDjs[iPivot]; |
---|
1294 | } |
---|
1295 | if (fabs(work[iRow]) > largestDualError_) { |
---|
1296 | largestDualError_ = fabs(work[iRow]); |
---|
1297 | //assert (largestDualError_<1.0e-7); |
---|
1298 | //if (largestDualError_>1.0e-7) |
---|
1299 | //printf("large dual error %g\n",largestDualError_); |
---|
1300 | } |
---|
1301 | } |
---|
1302 | } |
---|
1303 | } |
---|
1304 | if (largestDualError_ >= lastError) { |
---|
1305 | // restore |
---|
1306 | CoinIndexedVector *temp = thisVector; |
---|
1307 | thisVector = lastVector; |
---|
1308 | lastVector = temp; |
---|
1309 | break; |
---|
1310 | } |
---|
1311 | if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10 |
---|
1312 | && !givenDjs) { |
---|
1313 | // try and make better |
---|
1314 | // save this |
---|
1315 | CoinIndexedVector *temp = thisVector; |
---|
1316 | thisVector = lastVector; |
---|
1317 | lastVector = temp; |
---|
1318 | int *indexOut = thisVector->getIndices(); |
---|
1319 | int number = 0; |
---|
1320 | array = thisVector->denseVector(); |
---|
1321 | thisVector->clear(); |
---|
1322 | double multiplier = 131072.0; |
---|
1323 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1324 | double value = multiplier * work[iRow]; |
---|
1325 | if (value) { |
---|
1326 | array[iRow] = value; |
---|
1327 | indexOut[number++] = iRow; |
---|
1328 | work[iRow] = 0.0; |
---|
1329 | } |
---|
1330 | work[iRow] = 0.0; |
---|
1331 | } |
---|
1332 | thisVector->setNumElements(number); |
---|
1333 | lastError = largestDualError_; |
---|
1334 | factorization_->updateColumnTranspose(workSpace, thisVector); |
---|
1335 | multiplier = 1.0 / multiplier; |
---|
1336 | double *previous = lastVector->denseVector(); |
---|
1337 | number = 0; |
---|
1338 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1339 | double value = previous[iRow] + multiplier * array[iRow]; |
---|
1340 | if (value) { |
---|
1341 | array[iRow] = value; |
---|
1342 | indexOut[number++] = iRow; |
---|
1343 | } else { |
---|
1344 | array[iRow] = 0.0; |
---|
1345 | } |
---|
1346 | } |
---|
1347 | thisVector->setNumElements(number); |
---|
1348 | } else { |
---|
1349 | break; |
---|
1350 | } |
---|
1351 | } |
---|
1352 | // now look at dual solution |
---|
1353 | array = thisVector->denseVector(); |
---|
1354 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1355 | // slack |
---|
1356 | double value = array[iRow]; |
---|
1357 | dual_[iRow] = value; |
---|
1358 | value += rowObjectiveWork_[iRow]; |
---|
1359 | rowReducedCost_[iRow] = value; |
---|
1360 | } |
---|
1361 | // can use work if problem scaled (for better cache) |
---|
1362 | ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix_); |
---|
1363 | double *saveRowScale = rowScale_; |
---|
1364 | //double * saveColumnScale = columnScale_; |
---|
1365 | if (scaledMatrix_) { |
---|
1366 | rowScale_ = NULL; |
---|
1367 | clpMatrix = scaledMatrix_; |
---|
1368 | } |
---|
1369 | if (clpMatrix && (clpMatrix->flags() & 2) == 0) { |
---|
1370 | CoinIndexedVector *cVector = columnArray_[0]; |
---|
1371 | int *whichColumn = cVector->getIndices(); |
---|
1372 | assert(!cVector->getNumElements()); |
---|
1373 | int n = 0; |
---|
1374 | for (int i = 0; i < numberColumns_; i++) { |
---|
1375 | if (getColumnStatus(i) != basic) { |
---|
1376 | whichColumn[n++] = i; |
---|
1377 | reducedCostWork_[i] = objectiveWork_[i]; |
---|
1378 | } else { |
---|
1379 | reducedCostWork_[i] = 0.0; |
---|
1380 | } |
---|
1381 | } |
---|
1382 | if (numberRows_ > 4000) |
---|
1383 | clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_, |
---|
1384 | rowScale_, columnScale_, work); |
---|
1385 | else |
---|
1386 | clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_, |
---|
1387 | rowScale_, columnScale_, NULL); |
---|
1388 | } else { |
---|
1389 | ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_); |
---|
1390 | if (numberRows_ > 4000) |
---|
1391 | matrix_->transposeTimes(-1.0, dual_, reducedCostWork_, |
---|
1392 | rowScale_, columnScale_, work); |
---|
1393 | else |
---|
1394 | matrix_->transposeTimes(-1.0, dual_, reducedCostWork_, |
---|
1395 | rowScale_, columnScale_, NULL); |
---|
1396 | } |
---|
1397 | rowScale_ = saveRowScale; |
---|
1398 | //columnScale_ = saveColumnScale; |
---|
1399 | ClpFillN(work, numberRows_, 0.0); |
---|
1400 | // Extended duals and check dual infeasibility |
---|
1401 | if (!matrix_->skipDualCheck() || algorithm_ < 0 || problemStatus_ != -2) |
---|
1402 | matrix_->dualExpanded(this, NULL, NULL, 2); |
---|
1403 | // If necessary - override results |
---|
1404 | if (givenDjs) { |
---|
1405 | // restore accurate duals |
---|
1406 | CoinMemcpyN(dj_, (numberRows_ + numberColumns_), givenDjs); |
---|
1407 | } |
---|
1408 | arrayVector->clear(); |
---|
1409 | previousVector->clear(); |
---|
1410 | #ifndef SLIM_CLP |
---|
1411 | } else { |
---|
1412 | // Nonlinear |
---|
1413 | objective_->reducedGradient(this, dj_, false); |
---|
1414 | // get dual_ by moving from reduced costs for slacks |
---|
1415 | CoinMemcpyN(dj_ + numberColumns_, numberRows_, dual_); |
---|
1416 | } |
---|
1417 | #endif |
---|
1418 | } |
---|
1419 | /* Given an existing factorization computes and checks |
---|
1420 | primal and dual solutions. Uses input arrays for variables at |
---|
1421 | bounds. Returns feasibility states */ |
---|
1422 | int ClpSimplex::getSolution(const double * /*rowActivities*/, |
---|
1423 | const double * /*columnActivities*/) |
---|
1424 | { |
---|
1425 | if (!factorization_->status()) { |
---|
1426 | // put in standard form |
---|
1427 | createRim(7 + 8 + 16 + 32, false, -1); |
---|
1428 | if (pivotVariable_[0] < 0) |
---|
1429 | internalFactorize(0); |
---|
1430 | // do work |
---|
1431 | gutsOfSolution(NULL, NULL); |
---|
1432 | // release extra memory |
---|
1433 | deleteRim(0); |
---|
1434 | } |
---|
1435 | return factorization_->status(); |
---|
1436 | } |
---|
1437 | /* Given an existing factorization computes and checks |
---|
1438 | primal and dual solutions. Uses current problem arrays for |
---|
1439 | bounds. Returns feasibility states */ |
---|
1440 | int ClpSimplex::getSolution() |
---|
1441 | { |
---|
1442 | double *rowActivities = new double[numberRows_]; |
---|
1443 | double *columnActivities = new double[numberColumns_]; |
---|
1444 | ClpDisjointCopyN(rowActivityWork_, numberRows_, rowActivities); |
---|
1445 | ClpDisjointCopyN(columnActivityWork_, numberColumns_, columnActivities); |
---|
1446 | int status = getSolution(rowActivities, columnActivities); |
---|
1447 | delete[] rowActivities; |
---|
1448 | delete[] columnActivities; |
---|
1449 | return status; |
---|
1450 | } |
---|
1451 | // Factorizes using current basis. This is for external use |
---|
1452 | // Return codes are as from ClpFactorization |
---|
1453 | int ClpSimplex::factorize() |
---|
1454 | { |
---|
1455 | // put in standard form |
---|
1456 | createRim(7 + 8 + 16 + 32, false); |
---|
1457 | // do work |
---|
1458 | int status = internalFactorize(-1); |
---|
1459 | // release extra memory |
---|
1460 | deleteRim(0); |
---|
1461 | |
---|
1462 | return status; |
---|
1463 | } |
---|
1464 | // Clean up status |
---|
1465 | void ClpSimplex::cleanStatus() |
---|
1466 | { |
---|
1467 | int iRow, iColumn; |
---|
1468 | int numberBasic = 0; |
---|
1469 | // make row activities correct |
---|
1470 | memset(rowActivityWork_, 0, numberRows_ * sizeof(double)); |
---|
1471 | times(1.0, columnActivityWork_, rowActivityWork_); |
---|
1472 | if (!status_) |
---|
1473 | createStatus(); |
---|
1474 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1475 | if (getRowStatus(iRow) == basic) |
---|
1476 | numberBasic++; |
---|
1477 | else { |
---|
1478 | setRowStatus(iRow, superBasic); |
---|
1479 | // but put to bound if close |
---|
1480 | if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow]) |
---|
1481 | <= primalTolerance_) { |
---|
1482 | rowActivityWork_[iRow] = rowLowerWork_[iRow]; |
---|
1483 | setRowStatus(iRow, atLowerBound); |
---|
1484 | } else if (fabs(rowActivityWork_[iRow] - rowUpperWork_[iRow]) |
---|
1485 | <= primalTolerance_) { |
---|
1486 | rowActivityWork_[iRow] = rowUpperWork_[iRow]; |
---|
1487 | setRowStatus(iRow, atUpperBound); |
---|
1488 | } |
---|
1489 | } |
---|
1490 | } |
---|
1491 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1492 | if (getColumnStatus(iColumn) == basic) { |
---|
1493 | if (numberBasic == numberRows_) { |
---|
1494 | // take out of basis |
---|
1495 | setColumnStatus(iColumn, superBasic); |
---|
1496 | // but put to bound if close |
---|
1497 | if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn]) |
---|
1498 | <= primalTolerance_) { |
---|
1499 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1500 | setColumnStatus(iColumn, atLowerBound); |
---|
1501 | } else if (fabs(columnActivityWork_[iColumn] |
---|
1502 | - columnUpperWork_[iColumn]) |
---|
1503 | <= primalTolerance_) { |
---|
1504 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1505 | setColumnStatus(iColumn, atUpperBound); |
---|
1506 | } |
---|
1507 | } else |
---|
1508 | numberBasic++; |
---|
1509 | } else { |
---|
1510 | setColumnStatus(iColumn, superBasic); |
---|
1511 | // but put to bound if close |
---|
1512 | if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn]) |
---|
1513 | <= primalTolerance_) { |
---|
1514 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1515 | setColumnStatus(iColumn, atLowerBound); |
---|
1516 | } else if (fabs(columnActivityWork_[iColumn] |
---|
1517 | - columnUpperWork_[iColumn]) |
---|
1518 | <= primalTolerance_) { |
---|
1519 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1520 | setColumnStatus(iColumn, atUpperBound); |
---|
1521 | } |
---|
1522 | } |
---|
1523 | } |
---|
1524 | } |
---|
1525 | |
---|
1526 | /* Factorizes using current basis. |
---|
1527 | solveType - 1 iterating, 0 initial, -1 external |
---|
1528 | - 2 then iterating but can throw out of basis |
---|
1529 | If 10 added then in primal values pass |
---|
1530 | Return codes are as from ClpFactorization unless initial factorization |
---|
1531 | when total number of singularities is returned. |
---|
1532 | Special case is numberRows_+1 -> all slack basis. |
---|
1533 | */ |
---|
1534 | int ClpSimplex::internalFactorize(int solveType) |
---|
1535 | { |
---|
1536 | int iRow, iColumn; |
---|
1537 | int totalSlacks = numberRows_; |
---|
1538 | if (!status_) |
---|
1539 | createStatus(); |
---|
1540 | |
---|
1541 | bool valuesPass = false; |
---|
1542 | if (solveType >= 10) { |
---|
1543 | valuesPass = true; |
---|
1544 | solveType -= 10; |
---|
1545 | } |
---|
1546 | #ifdef CLP_DEBUG |
---|
1547 | if (solveType > 0) { |
---|
1548 | int numberFreeIn = 0, numberFreeOut = 0; |
---|
1549 | double biggestDj = 0.0; |
---|
1550 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1551 | switch (getColumnStatus(iColumn)) { |
---|
1552 | |
---|
1553 | case basic: |
---|
1554 | if (columnLower_[iColumn] < -largeValue_ |
---|
1555 | && columnUpper_[iColumn] > largeValue_) |
---|
1556 | numberFreeIn++; |
---|
1557 | break; |
---|
1558 | default: |
---|
1559 | if (columnLower_[iColumn] < -largeValue_ |
---|
1560 | && columnUpper_[iColumn] > largeValue_) { |
---|
1561 | numberFreeOut++; |
---|
1562 | biggestDj = CoinMax(fabs(dj_[iColumn]), biggestDj); |
---|
1563 | } |
---|
1564 | break; |
---|
1565 | } |
---|
1566 | } |
---|
1567 | if (numberFreeIn + numberFreeOut) |
---|
1568 | printf("%d in basis, %d out - largest dj %g\n", |
---|
1569 | numberFreeIn, numberFreeOut, biggestDj); |
---|
1570 | } |
---|
1571 | #endif |
---|
1572 | if (solveType <= 0) { |
---|
1573 | // Make sure everything is clean |
---|
1574 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1575 | if (getRowStatus(iRow) == isFixed) { |
---|
1576 | // double check fixed |
---|
1577 | if (rowUpperWork_[iRow] > rowLowerWork_[iRow]) |
---|
1578 | setRowStatus(iRow, atLowerBound); |
---|
1579 | } else if (getRowStatus(iRow) == isFree) { |
---|
1580 | // may not be free after all |
---|
1581 | if (rowLowerWork_[iRow] > -largeValue_ || rowUpperWork_[iRow] < largeValue_) |
---|
1582 | setRowStatus(iRow, superBasic); |
---|
1583 | } |
---|
1584 | } |
---|
1585 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1586 | if (getColumnStatus(iColumn) == isFixed) { |
---|
1587 | // double check fixed |
---|
1588 | if (columnUpperWork_[iColumn] > columnLowerWork_[iColumn]) |
---|
1589 | setColumnStatus(iColumn, atLowerBound); |
---|
1590 | } else if (getColumnStatus(iColumn) == isFree) { |
---|
1591 | // may not be free after all |
---|
1592 | if (columnLowerWork_[iColumn] > -largeValue_ || columnUpperWork_[iColumn] < largeValue_) |
---|
1593 | setColumnStatus(iColumn, superBasic); |
---|
1594 | } |
---|
1595 | } |
---|
1596 | if (!valuesPass) { |
---|
1597 | // not values pass so set to bounds |
---|
1598 | bool allSlack = true; |
---|
1599 | if (status_) { |
---|
1600 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1601 | if (getRowStatus(iRow) != basic) { |
---|
1602 | allSlack = false; |
---|
1603 | break; |
---|
1604 | } |
---|
1605 | } |
---|
1606 | } |
---|
1607 | if (!allSlack) { |
---|
1608 | //#define CLP_INVESTIGATE2 |
---|
1609 | #ifdef CLP_INVESTIGATE3 |
---|
1610 | int numberTotal = numberRows_ + numberColumns_; |
---|
1611 | double *saveSol = valuesPass ? CoinCopyOfArray(solution_, numberTotal) : NULL; |
---|
1612 | #endif |
---|
1613 | // set values from warm start (if sensible) |
---|
1614 | int numberBasic = 0; |
---|
1615 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1616 | switch (getRowStatus(iRow)) { |
---|
1617 | |
---|
1618 | case basic: |
---|
1619 | numberBasic++; |
---|
1620 | break; |
---|
1621 | case atUpperBound: |
---|
1622 | rowActivityWork_[iRow] = rowUpperWork_[iRow]; |
---|
1623 | if (rowActivityWork_[iRow] > largeValue_) { |
---|
1624 | if (rowLowerWork_[iRow] > -largeValue_) { |
---|
1625 | rowActivityWork_[iRow] = rowLowerWork_[iRow]; |
---|
1626 | setRowStatus(iRow, atLowerBound); |
---|
1627 | } else { |
---|
1628 | // say free |
---|
1629 | setRowStatus(iRow, isFree); |
---|
1630 | rowActivityWork_[iRow] = 0.0; |
---|
1631 | } |
---|
1632 | } |
---|
1633 | break; |
---|
1634 | case ClpSimplex::isFixed: |
---|
1635 | case atLowerBound: |
---|
1636 | rowActivityWork_[iRow] = rowLowerWork_[iRow]; |
---|
1637 | if (rowActivityWork_[iRow] < -largeValue_) { |
---|
1638 | if (rowUpperWork_[iRow] < largeValue_) { |
---|
1639 | rowActivityWork_[iRow] = rowUpperWork_[iRow]; |
---|
1640 | setRowStatus(iRow, atUpperBound); |
---|
1641 | } else { |
---|
1642 | // say free |
---|
1643 | setRowStatus(iRow, isFree); |
---|
1644 | rowActivityWork_[iRow] = 0.0; |
---|
1645 | } |
---|
1646 | } |
---|
1647 | break; |
---|
1648 | case isFree: |
---|
1649 | break; |
---|
1650 | // not really free - fall through to superbasic |
---|
1651 | case superBasic: |
---|
1652 | if (rowUpperWork_[iRow] > largeValue_) { |
---|
1653 | if (rowLowerWork_[iRow] > -largeValue_) { |
---|
1654 | rowActivityWork_[iRow] = rowLowerWork_[iRow]; |
---|
1655 | setRowStatus(iRow, atLowerBound); |
---|
1656 | } else { |
---|
1657 | // say free |
---|
1658 | setRowStatus(iRow, isFree); |
---|
1659 | rowActivityWork_[iRow] = 0.0; |
---|
1660 | } |
---|
1661 | } else { |
---|
1662 | if (rowLowerWork_[iRow] > -largeValue_) { |
---|
1663 | // set to nearest |
---|
1664 | if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow]) |
---|
1665 | < fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])) { |
---|
1666 | rowActivityWork_[iRow] = rowLowerWork_[iRow]; |
---|
1667 | setRowStatus(iRow, atLowerBound); |
---|
1668 | } else { |
---|
1669 | rowActivityWork_[iRow] = rowUpperWork_[iRow]; |
---|
1670 | setRowStatus(iRow, atUpperBound); |
---|
1671 | } |
---|
1672 | } else { |
---|
1673 | rowActivityWork_[iRow] = rowUpperWork_[iRow]; |
---|
1674 | setRowStatus(iRow, atUpperBound); |
---|
1675 | } |
---|
1676 | } |
---|
1677 | break; |
---|
1678 | } |
---|
1679 | } |
---|
1680 | totalSlacks = numberBasic; |
---|
1681 | |
---|
1682 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1683 | switch (getColumnStatus(iColumn)) { |
---|
1684 | |
---|
1685 | case basic: |
---|
1686 | if (numberBasic == maximumBasic_) { |
---|
1687 | // take out of basis |
---|
1688 | if (columnLowerWork_[iColumn] > -largeValue_) { |
---|
1689 | if (columnActivityWork_[iColumn] - columnLowerWork_[iColumn] < columnUpperWork_[iColumn] - columnActivityWork_[iColumn]) { |
---|
1690 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1691 | setColumnStatus(iColumn, atLowerBound); |
---|
1692 | } else { |
---|
1693 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1694 | setColumnStatus(iColumn, atUpperBound); |
---|
1695 | } |
---|
1696 | } else if (columnUpperWork_[iColumn] < largeValue_) { |
---|
1697 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1698 | setColumnStatus(iColumn, atUpperBound); |
---|
1699 | } else { |
---|
1700 | columnActivityWork_[iColumn] = 0.0; |
---|
1701 | setColumnStatus(iColumn, isFree); |
---|
1702 | } |
---|
1703 | } else { |
---|
1704 | numberBasic++; |
---|
1705 | } |
---|
1706 | break; |
---|
1707 | case atUpperBound: |
---|
1708 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1709 | if (columnActivityWork_[iColumn] > largeValue_) { |
---|
1710 | if (columnLowerWork_[iColumn] < -largeValue_) { |
---|
1711 | columnActivityWork_[iColumn] = 0.0; |
---|
1712 | setColumnStatus(iColumn, isFree); |
---|
1713 | } else { |
---|
1714 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1715 | setColumnStatus(iColumn, atLowerBound); |
---|
1716 | } |
---|
1717 | } |
---|
1718 | break; |
---|
1719 | case isFixed: |
---|
1720 | case atLowerBound: |
---|
1721 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1722 | if (columnActivityWork_[iColumn] < -largeValue_) { |
---|
1723 | if (columnUpperWork_[iColumn] > largeValue_) { |
---|
1724 | columnActivityWork_[iColumn] = 0.0; |
---|
1725 | setColumnStatus(iColumn, isFree); |
---|
1726 | } else { |
---|
1727 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1728 | setColumnStatus(iColumn, atUpperBound); |
---|
1729 | } |
---|
1730 | } |
---|
1731 | break; |
---|
1732 | case isFree: |
---|
1733 | break; |
---|
1734 | // not really free - fall through to superbasic |
---|
1735 | case superBasic: |
---|
1736 | if (columnUpperWork_[iColumn] > largeValue_) { |
---|
1737 | if (columnLowerWork_[iColumn] > -largeValue_) { |
---|
1738 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1739 | setColumnStatus(iColumn, atLowerBound); |
---|
1740 | } else { |
---|
1741 | // say free |
---|
1742 | setColumnStatus(iColumn, isFree); |
---|
1743 | columnActivityWork_[iColumn] = 0.0; |
---|
1744 | } |
---|
1745 | } else { |
---|
1746 | if (columnLowerWork_[iColumn] > -largeValue_) { |
---|
1747 | // set to nearest |
---|
1748 | if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn]) |
---|
1749 | < fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])) { |
---|
1750 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1751 | setColumnStatus(iColumn, atLowerBound); |
---|
1752 | } else { |
---|
1753 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1754 | setColumnStatus(iColumn, atUpperBound); |
---|
1755 | } |
---|
1756 | } else { |
---|
1757 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1758 | setColumnStatus(iColumn, atUpperBound); |
---|
1759 | } |
---|
1760 | } |
---|
1761 | break; |
---|
1762 | } |
---|
1763 | } |
---|
1764 | #ifdef CLP_INVESTIGATE3 |
---|
1765 | if (saveSol) { |
---|
1766 | int numberChanged = 0; |
---|
1767 | double largestChanged = 0.0; |
---|
1768 | for (int i = 0; i < numberTotal; i++) { |
---|
1769 | double difference = fabs(solution_[i] - saveSol[i]); |
---|
1770 | if (difference > 1.0e-7) { |
---|
1771 | numberChanged++; |
---|
1772 | if (difference > largestChanged) |
---|
1773 | largestChanged = difference; |
---|
1774 | } |
---|
1775 | } |
---|
1776 | if (numberChanged) |
---|
1777 | printf("%d changed, largest %g\n", numberChanged, largestChanged); |
---|
1778 | delete[] saveSol; |
---|
1779 | } |
---|
1780 | #endif |
---|
1781 | #if 0 |
---|
1782 | if (numberBasic < numberRows_) { |
---|
1783 | // add some slacks in case odd warmstart |
---|
1784 | #ifdef CLP_INVESTIGATE |
---|
1785 | printf("BAD %d basic, %d rows %d slacks\n", |
---|
1786 | numberBasic, numberRows_, totalSlacks); |
---|
1787 | #endif |
---|
1788 | int iRow = numberRows_ - 1; |
---|
1789 | while (numberBasic < numberRows_) { |
---|
1790 | if (getRowStatus(iRow) != basic) { |
---|
1791 | setRowStatus(iRow, basic); |
---|
1792 | numberBasic++; |
---|
1793 | totalSlacks++; |
---|
1794 | iRow--; |
---|
1795 | } else { |
---|
1796 | break; |
---|
1797 | } |
---|
1798 | } |
---|
1799 | } |
---|
1800 | #endif |
---|
1801 | } else { |
---|
1802 | // all slack basis |
---|
1803 | int numberBasic = 0; |
---|
1804 | if (!status_) { |
---|
1805 | createStatus(); |
---|
1806 | } |
---|
1807 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1808 | double lower = rowLowerWork_[iRow]; |
---|
1809 | double upper = rowUpperWork_[iRow]; |
---|
1810 | if (lower > -largeValue_ || upper < largeValue_) { |
---|
1811 | if (fabs(lower) <= fabs(upper)) { |
---|
1812 | rowActivityWork_[iRow] = lower; |
---|
1813 | } else { |
---|
1814 | rowActivityWork_[iRow] = upper; |
---|
1815 | } |
---|
1816 | } else { |
---|
1817 | rowActivityWork_[iRow] = 0.0; |
---|
1818 | } |
---|
1819 | setRowStatus(iRow, basic); |
---|
1820 | numberBasic++; |
---|
1821 | } |
---|
1822 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1823 | double lower = columnLowerWork_[iColumn]; |
---|
1824 | double upper = columnUpperWork_[iColumn]; |
---|
1825 | double big_bound = largeValue_; |
---|
1826 | double value = columnActivityWork_[iColumn]; |
---|
1827 | if (lower > -big_bound || upper < big_bound) { |
---|
1828 | if ((getColumnStatus(iColumn) == atLowerBound && columnActivityWork_[iColumn] == lower) || (getColumnStatus(iColumn) == atUpperBound && columnActivityWork_[iColumn] == upper)) { |
---|
1829 | // status looks plausible |
---|
1830 | } else { |
---|
1831 | // set to sensible |
---|
1832 | if (getColumnStatus(iColumn) == atUpperBound |
---|
1833 | && upper < 1.0e20) { |
---|
1834 | columnActivityWork_[iColumn] = upper; |
---|
1835 | } else if (getColumnStatus(iColumn) == atLowerBound |
---|
1836 | && lower > -1.0e20) { |
---|
1837 | columnActivityWork_[iColumn] = lower; |
---|
1838 | } else { |
---|
1839 | if (fabs(lower) <= fabs(upper)) { |
---|
1840 | setColumnStatus(iColumn, atLowerBound); |
---|
1841 | columnActivityWork_[iColumn] = lower; |
---|
1842 | } else { |
---|
1843 | setColumnStatus(iColumn, atUpperBound); |
---|
1844 | columnActivityWork_[iColumn] = upper; |
---|
1845 | } |
---|
1846 | } |
---|
1847 | } |
---|
1848 | } else { |
---|
1849 | setColumnStatus(iColumn, isFree); |
---|
1850 | columnActivityWork_[iColumn] = 0.0; |
---|
1851 | } |
---|
1852 | } |
---|
1853 | } |
---|
1854 | } else { |
---|
1855 | // values pass has less coding |
---|
1856 | // make row activities correct and clean basis a bit |
---|
1857 | cleanStatus(); |
---|
1858 | if (status_) { |
---|
1859 | int numberBasic = 0; |
---|
1860 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1861 | if (getRowStatus(iRow) == basic) |
---|
1862 | numberBasic++; |
---|
1863 | } |
---|
1864 | totalSlacks = numberBasic; |
---|
1865 | #if 0 |
---|
1866 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1867 | if (getColumnStatus(iColumn) == basic) |
---|
1868 | numberBasic++; |
---|
1869 | } |
---|
1870 | #endif |
---|
1871 | } else { |
---|
1872 | // all slack basis |
---|
1873 | int numberBasic = 0; |
---|
1874 | if (!status_) { |
---|
1875 | createStatus(); |
---|
1876 | } |
---|
1877 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1878 | setRowStatus(iRow, basic); |
---|
1879 | numberBasic++; |
---|
1880 | } |
---|
1881 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1882 | setColumnStatus(iColumn, superBasic); |
---|
1883 | // but put to bound if close |
---|
1884 | if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn]) |
---|
1885 | <= primalTolerance_) { |
---|
1886 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1887 | setColumnStatus(iColumn, atLowerBound); |
---|
1888 | } else if (fabs(columnActivityWork_[iColumn] |
---|
1889 | - columnUpperWork_[iColumn]) |
---|
1890 | <= primalTolerance_) { |
---|
1891 | columnActivityWork_[iColumn] = columnUpperWork_[iColumn]; |
---|
1892 | setColumnStatus(iColumn, atUpperBound); |
---|
1893 | } |
---|
1894 | } |
---|
1895 | } |
---|
1896 | } |
---|
1897 | numberRefinements_ = 1; |
---|
1898 | // set fixed if they are |
---|
1899 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1900 | if (getRowStatus(iRow) != basic) { |
---|
1901 | if (rowLowerWork_[iRow] == rowUpperWork_[iRow]) { |
---|
1902 | rowActivityWork_[iRow] = rowLowerWork_[iRow]; |
---|
1903 | setRowStatus(iRow, isFixed); |
---|
1904 | } |
---|
1905 | } |
---|
1906 | } |
---|
1907 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1908 | if (getColumnStatus(iColumn) != basic) { |
---|
1909 | if (columnLowerWork_[iColumn] == columnUpperWork_[iColumn]) { |
---|
1910 | columnActivityWork_[iColumn] = columnLowerWork_[iColumn]; |
---|
1911 | setColumnStatus(iColumn, isFixed); |
---|
1912 | } |
---|
1913 | } |
---|
1914 | } |
---|
1915 | } |
---|
1916 | //for (iRow=0;iRow<numberRows_+numberColumns_;iRow++) { |
---|
1917 | //if (fabs(solution_[iRow])>1.0e10) { |
---|
1918 | // printf("large %g at %d - status %d\n", |
---|
1919 | // solution_[iRow],iRow,status_[iRow]); |
---|
1920 | //} |
---|
1921 | //} |
---|
1922 | #if 0 //ndef _MSC_VER \ |
---|
1923 | // The local static var k is a problem when trying to build a DLL. Since this is \ |
---|
1924 | // just for debugging (likely done on *nix), just hide it from Windows \ |
---|
1925 | // -- lh, 101016 -- |
---|
1926 | if (0) { |
---|
1927 | static int k = 0; |
---|
1928 | printf("start basis\n"); |
---|
1929 | int i; |
---|
1930 | for (i = 0; i < numberRows_; i++) |
---|
1931 | printf ("xx %d %d\n", i, pivotVariable_[i]); |
---|
1932 | for (i = 0; i < numberRows_ + numberColumns_; i++) |
---|
1933 | if (getColumnStatus(i) == basic) |
---|
1934 | printf ("yy %d basic\n", i); |
---|
1935 | if (k > 20) |
---|
1936 | exit(0); |
---|
1937 | k++; |
---|
1938 | } |
---|
1939 | #endif |
---|
1940 | #if 0 //ndef NDEBUG |
---|
1941 | // Make sure everything is clean |
---|
1942 | double sumOutside=0.0; |
---|
1943 | int numberOutside=0; |
---|
1944 | //double sumOutsideLarge=0.0; |
---|
1945 | int numberOutsideLarge=0; |
---|
1946 | double sumInside=0.0; |
---|
1947 | int numberInside=0; |
---|
1948 | //double sumInsideLarge=0.0; |
---|
1949 | int numberInsideLarge=0; |
---|
1950 | int numberTotal=numberRows_+numberColumns_; |
---|
1951 | for (int iSequence = 0; iSequence < numberTotal; iSequence++) { |
---|
1952 | if(getStatus(iSequence) == isFixed) { |
---|
1953 | // double check fixed |
---|
1954 | assert (upper_[iSequence] == lower_[iSequence]); |
---|
1955 | assert (fabs(solution_[iSequence]-lower_[iSequence])<primalTolerance_); |
---|
1956 | } else if (getStatus(iSequence) == isFree) { |
---|
1957 | assert (upper_[iSequence] == COIN_DBL_MAX && lower_[iSequence]==-COIN_DBL_MAX); |
---|
1958 | } else if (getStatus(iSequence) == atLowerBound) { |
---|
1959 | assert (fabs(solution_[iSequence]-lower_[iSequence])<1000.0*primalTolerance_); |
---|
1960 | if (solution_[iSequence]<lower_[iSequence]) { |
---|
1961 | numberOutside++; |
---|
1962 | sumOutside-=solution_[iSequence]-lower_[iSequence]; |
---|
1963 | if (solution_[iSequence]<lower_[iSequence]-primalTolerance_) |
---|
1964 | numberOutsideLarge++; |
---|
1965 | } else if (solution_[iSequence]>lower_[iSequence]) { |
---|
1966 | numberInside++; |
---|
1967 | sumInside+=solution_[iSequence]-lower_[iSequence]; |
---|
1968 | if (solution_[iSequence]>lower_[iSequence]+primalTolerance_) |
---|
1969 | numberInsideLarge++; |
---|
1970 | } |
---|
1971 | } else if (getStatus(iSequence) == atUpperBound) { |
---|
1972 | assert (fabs(solution_[iSequence]-upper_[iSequence])<1000.0*primalTolerance_); |
---|
1973 | if (solution_[iSequence]>upper_[iSequence]) { |
---|
1974 | numberOutside++; |
---|
1975 | sumOutside+=solution_[iSequence]-upper_[iSequence]; |
---|
1976 | if (solution_[iSequence]>upper_[iSequence]+primalTolerance_) |
---|
1977 | numberOutsideLarge++; |
---|
1978 | } else if (solution_[iSequence]<upper_[iSequence]) { |
---|
1979 | numberInside++; |
---|
1980 | sumInside-=solution_[iSequence]-upper_[iSequence]; |
---|
1981 | if (solution_[iSequence]<upper_[iSequence]-primalTolerance_) |
---|
1982 | numberInsideLarge++; |
---|
1983 | } |
---|
1984 | } else if (getStatus(iSequence) == superBasic) { |
---|
1985 | //assert (!valuesPass); |
---|
1986 | } |
---|
1987 | } |
---|
1988 | if (numberInside+numberOutside) |
---|
1989 | printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n", |
---|
1990 | numberOutside,sumOutside,numberOutsideLarge, |
---|
1991 | numberInside,sumInside,numberInsideLarge); |
---|
1992 | #endif |
---|
1993 | int status = factorization_->factorize(this, solveType, valuesPass); |
---|
1994 | if (status) { |
---|
1995 | handler_->message(CLP_SIMPLEX_BADFACTOR, messages_) |
---|
1996 | << status |
---|
1997 | << CoinMessageEol; |
---|
1998 | #ifdef CLP_USEFUL_PRINTOUT |
---|
1999 | printf("Basis singular - pivot tolerance %g\n", |
---|
2000 | factorization_->pivotTolerance()); |
---|
2001 | #endif |
---|
2002 | return -1; |
---|
2003 | } else if (!solveType) { |
---|
2004 | // Initial basis - return number of singularities |
---|
2005 | int numberSlacks = 0; |
---|
2006 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2007 | if (getRowStatus(iRow) == basic) |
---|
2008 | numberSlacks++; |
---|
2009 | } |
---|
2010 | status = CoinMax(numberSlacks - totalSlacks, 0); |
---|
2011 | // special case if all slack |
---|
2012 | if (numberSlacks == numberRows_) { |
---|
2013 | status = numberRows_ + 1; |
---|
2014 | } |
---|
2015 | } |
---|
2016 | |
---|
2017 | // sparse methods |
---|
2018 | //if (factorization_->sparseThreshold()) { |
---|
2019 | // get default value |
---|
2020 | factorization_->sparseThreshold(0); |
---|
2021 | if (!(moreSpecialOptions_ & 1024)) |
---|
2022 | factorization_->goSparse(); |
---|
2023 | //} |
---|
2024 | |
---|
2025 | return status; |
---|
2026 | } |
---|
2027 | /* |
---|
2028 | This does basis housekeeping and does values for in/out variables. |
---|
2029 | Can also decide to re-factorize |
---|
2030 | */ |
---|
2031 | int ClpSimplex::housekeeping(double objectiveChange) |
---|
2032 | { |
---|
2033 | // save value of incoming and outgoing |
---|
2034 | double oldIn = solution_[sequenceIn_]; |
---|
2035 | double oldOut = solution_[sequenceOut_]; |
---|
2036 | numberIterations_++; |
---|
2037 | changeMade_++; // something has happened |
---|
2038 | // incoming variable |
---|
2039 | if (handler_->logLevel() > 7) { |
---|
2040 | //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) { |
---|
2041 | handler_->message(CLP_SIMPLEX_HOUSE1, messages_) |
---|
2042 | << directionOut_ |
---|
2043 | << directionIn_ << theta_ |
---|
2044 | << dualOut_ << dualIn_ << alpha_ |
---|
2045 | << CoinMessageEol; |
---|
2046 | if (getStatus(sequenceIn_) == isFree) { |
---|
2047 | handler_->message(CLP_SIMPLEX_FREEIN, messages_) |
---|
2048 | << sequenceIn_ |
---|
2049 | << CoinMessageEol; |
---|
2050 | } |
---|
2051 | } |
---|
2052 | #if 0 |
---|
2053 | printf("h1 %d %d %g %g %g %g", |
---|
2054 | directionOut_ |
---|
2055 | , directionIn_, theta_ |
---|
2056 | , dualOut_, dualIn_, alpha_); |
---|
2057 | #endif |
---|
2058 | // change of incoming |
---|
2059 | char rowcol[] = { 'R', 'C' }; |
---|
2060 | if (pivotRow_ >= 0) |
---|
2061 | pivotVariable_[pivotRow_] = sequenceIn(); |
---|
2062 | if (upper_[sequenceIn_] > 1.0e20 && lower_[sequenceIn_] < -1.0e20) |
---|
2063 | progressFlag_ |= 2; // making real progress |
---|
2064 | solution_[sequenceIn_] = valueIn_; |
---|
2065 | if (upper_[sequenceOut_] - lower_[sequenceOut_] < 1.0e-12) |
---|
2066 | progressFlag_ |= 1; // making real progress |
---|
2067 | if (sequenceIn_ != sequenceOut_) { |
---|
2068 | if (alphaAccuracy_ > 0.0) { |
---|
2069 | double value = fabs(alpha_); |
---|
2070 | if (value > 1.0) |
---|
2071 | alphaAccuracy_ *= value; |
---|
2072 | else |
---|
2073 | alphaAccuracy_ /= value; |
---|
2074 | } |
---|
2075 | //assert( getStatus(sequenceOut_)== basic); |
---|
2076 | setStatus(sequenceIn_, basic); |
---|
2077 | if (upper_[sequenceOut_] - lower_[sequenceOut_] > 0) { |
---|
2078 | // As Nonlinear costs may have moved bounds (to more feasible) |
---|
2079 | // Redo using value |
---|
2080 | if (fabs(valueOut_ - lower_[sequenceOut_]) < fabs(valueOut_ - upper_[sequenceOut_])) { |
---|
2081 | // going to lower |
---|
2082 | setStatus(sequenceOut_, atLowerBound); |
---|
2083 | oldOut = lower_[sequenceOut_]; |
---|
2084 | } else { |
---|
2085 | // going to upper |
---|
2086 | setStatus(sequenceOut_, atUpperBound); |
---|
2087 | oldOut = upper_[sequenceOut_]; |
---|
2088 | } |
---|
2089 | } else { |
---|
2090 | // fixed |
---|
2091 | setStatus(sequenceOut_, isFixed); |
---|
2092 | } |
---|
2093 | solution_[sequenceOut_] = valueOut_; |
---|
2094 | } else { |
---|
2095 | //if (objective_->type()<2) |
---|
2096 | //assert (fabs(theta_)>1.0e-13); |
---|
2097 | // flip from bound to bound |
---|
2098 | // As Nonlinear costs may have moved bounds (to more feasible) |
---|
2099 | // Redo using value |
---|
2100 | if (fabs(valueIn_ - lower_[sequenceIn_]) < fabs(valueIn_ - upper_[sequenceIn_])) { |
---|
2101 | // as if from upper bound |
---|
2102 | setStatus(sequenceIn_, atLowerBound); |
---|
2103 | } else { |
---|
2104 | // as if from lower bound |
---|
2105 | setStatus(sequenceIn_, atUpperBound); |
---|
2106 | } |
---|
2107 | } |
---|
2108 | |
---|
2109 | // Update hidden stuff e.g. effective RHS and gub |
---|
2110 | int invertNow = matrix_->updatePivot(this, oldIn, oldOut); |
---|
2111 | objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_); |
---|
2112 | if (handler_->logLevel() > 7) { |
---|
2113 | //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) { |
---|
2114 | handler_->message(CLP_SIMPLEX_HOUSE2, messages_) |
---|
2115 | << numberIterations_ << objectiveValue() |
---|
2116 | << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_) |
---|
2117 | << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_); |
---|
2118 | handler_->printing(algorithm_ < 0) << dualOut_ << theta_; |
---|
2119 | handler_->printing(algorithm_ > 0) << dualIn_ << theta_; |
---|
2120 | handler_->message() << CoinMessageEol; |
---|
2121 | } |
---|
2122 | #if 0 |
---|
2123 | if (numberIterations_ > 10000) |
---|
2124 | printf(" it %d %g %c%d %c%d\n" |
---|
2125 | , numberIterations_, objectiveValue() |
---|
2126 | , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_) |
---|
2127 | , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_)); |
---|
2128 | #endif |
---|
2129 | if (trustedUserPointer_ && trustedUserPointer_->typeStruct == 1) { |
---|
2130 | if (algorithm_ > 0 && integerType_ && !nonLinearCost_->numberInfeasibilities()) { |
---|
2131 | if (fabs(theta_) > 1.0e-6 || !numberIterations_) { |
---|
2132 | // For saving solutions |
---|
2133 | typedef struct { |
---|
2134 | int numberSolutions; |
---|
2135 | int maximumSolutions; |
---|
2136 | int numberColumns; |
---|
2137 | double **solution; |
---|
2138 | int *numberUnsatisfied; |
---|
2139 | } clpSolution; |
---|
2140 | clpSolution *solution = reinterpret_cast< clpSolution * >(trustedUserPointer_->data); |
---|
2141 | if (solution->numberSolutions == solution->maximumSolutions) { |
---|
2142 | int n = solution->maximumSolutions; |
---|
2143 | int n2 = (n * 3) / 2 + 10; |
---|
2144 | solution->maximumSolutions = n2; |
---|
2145 | double **temp = new double *[n2]; |
---|
2146 | for (int i = 0; i < n; i++) |
---|
2147 | temp[i] = solution->solution[i]; |
---|
2148 | delete[] solution->solution; |
---|
2149 | solution->solution = temp; |
---|
2150 | int *tempN = new int[n2]; |
---|
2151 | for (int i = 0; i < n; i++) |
---|
2152 | tempN[i] = solution->numberUnsatisfied[i]; |
---|
2153 | delete[] solution->numberUnsatisfied; |
---|
2154 | solution->numberUnsatisfied = tempN; |
---|
2155 | } |
---|
2156 | assert(numberColumns_ == solution->numberColumns); |
---|
2157 | double *sol = new double[numberColumns_]; |
---|
2158 | solution->solution[solution->numberSolutions] = sol; |
---|
2159 | int numberFixed = 0; |
---|
2160 | int numberUnsat = 0; |
---|
2161 | int numberSat = 0; |
---|
2162 | double sumUnsat = 0.0; |
---|
2163 | double tolerance = 10.0 * primalTolerance_; |
---|
2164 | double mostAway = 0.0; |
---|
2165 | for (int i = 0; i < numberColumns_; i++) { |
---|
2166 | // Save anyway |
---|
2167 | sol[i] = columnScale_ ? solution_[i] * columnScale_[i] : solution_[i]; |
---|
2168 | // rest is optional |
---|
2169 | if (upper_[i] > lower_[i]) { |
---|
2170 | double value = solution_[i]; |
---|
2171 | if (value > lower_[i] + tolerance && value < upper_[i] - tolerance && integerType_[i]) { |
---|
2172 | // may have to modify value if scaled |
---|
2173 | if (columnScale_) |
---|
2174 | value *= columnScale_[i]; |
---|
2175 | double closest = floor(value + 0.5); |
---|
2176 | // problem may be perturbed so relax test |
---|
2177 | if (fabs(value - closest) > 1.0e-4) { |
---|
2178 | numberUnsat++; |
---|
2179 | sumUnsat += fabs(value - closest); |
---|
2180 | if (mostAway < fabs(value - closest)) { |
---|
2181 | mostAway = fabs(value - closest); |
---|
2182 | } |
---|
2183 | } else { |
---|
2184 | numberSat++; |
---|
2185 | } |
---|
2186 | } else { |
---|
2187 | numberSat++; |
---|
2188 | } |
---|
2189 | } else { |
---|
2190 | numberFixed++; |
---|
2191 | } |
---|
2192 | } |
---|
2193 | solution->numberUnsatisfied[solution->numberSolutions++] = numberUnsat; |
---|
2194 | COIN_DETAIL_PRINT(printf("iteration %d, %d unsatisfied (%g,%g), %d fixed, %d satisfied\n", |
---|
2195 | numberIterations_, numberUnsat, sumUnsat, mostAway, numberFixed, numberSat)); |
---|
2196 | } |
---|
2197 | } |
---|
2198 | } |
---|
2199 | if (hitMaximumIterations()) |
---|
2200 | return 2; |
---|
2201 | #if 1 |
---|
2202 | //if (numberIterations_>14000) |
---|
2203 | //handler_->setLogLevel(63); |
---|
2204 | //if (numberIterations_>24000) |
---|
2205 | //exit(77); |
---|
2206 | // check for small cycles |
---|
2207 | int in = sequenceIn_; |
---|
2208 | int out = sequenceOut_; |
---|
2209 | matrix_->correctSequence(this, in, out); |
---|
2210 | int cycle = progress_.cycle(in, out, |
---|
2211 | directionIn_, directionOut_); |
---|
2212 | if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) { |
---|
2213 | //if (cycle>0) { |
---|
2214 | if (handler_->logLevel() >= 63) |
---|
2215 | printf("Cycle of %d\n", cycle); |
---|
2216 | // reset |
---|
2217 | progress_.startCheck(); |
---|
2218 | double random = randomNumberGenerator_.randomDouble(); |
---|
2219 | int extra = static_cast< int >(9.999 * random); |
---|
2220 | int off[] = { 1, 1, 1, 1, 2, 2, 2, 3, 3, 4 }; |
---|
2221 | if (factorization_->pivots() > cycle) { |
---|
2222 | forceFactorization_ = CoinMax(1, cycle - off[extra]); |
---|
2223 | } else { |
---|
2224 | /* need to reject something |
---|
2225 | should be better if don't reject incoming |
---|
2226 | as it is in basis */ |
---|
2227 | int iSequence; |
---|
2228 | //if (algorithm_ > 0) |
---|
2229 | // iSequence = sequenceIn_; |
---|
2230 | //else |
---|
2231 | iSequence = sequenceOut_; |
---|
2232 | char x = isColumn(iSequence) ? 'C' : 'R'; |
---|
2233 | if (handler_->logLevel() >= 63) |
---|
2234 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
2235 | << x << sequenceWithin(iSequence) |
---|
2236 | << CoinMessageEol; |
---|
2237 | setFlagged(iSequence); |
---|
2238 | //printf("flagging %d\n",iSequence); |
---|
2239 | } |
---|
2240 | return 1; |
---|
2241 | } |
---|
2242 | #endif |
---|
2243 | // only time to re-factorize if one before real time |
---|
2244 | // this is so user won't be surprised that maximumPivots has exact meaning |
---|
2245 | int numberPivots = factorization_->pivots(); |
---|
2246 | int maximumPivots = factorization_->maximumPivots(); |
---|
2247 | int numberDense = factorization_->numberDense(); |
---|
2248 | bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 > 2 * maximumIterations()); |
---|
2249 | if (numberPivots == maximumPivots || maximumPivots < 2) { |
---|
2250 | // If dense then increase |
---|
2251 | if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots |
---|
2252 | && false) { |
---|
2253 | factorization_->maximumPivots(numberDense); |
---|
2254 | dualRowPivot_->maximumPivotsChanged(); |
---|
2255 | primalColumnPivot_->maximumPivotsChanged(); |
---|
2256 | // and redo arrays |
---|
2257 | for (int iRow = 0; iRow < 4; iRow++) { |
---|
2258 | int length = rowArray_[iRow]->capacity() + numberDense - maximumPivots; |
---|
2259 | rowArray_[iRow]->reserve(length); |
---|
2260 | } |
---|
2261 | } |
---|
2262 | #if CLP_FACTORIZATION_NEW_TIMING > 1 |
---|
2263 | factorization_->statsRefactor('M'); |
---|
2264 | #endif |
---|
2265 | return 1; |
---|
2266 | } else if ((factorization_->timeToRefactorize() && !dontInvert) |
---|
2267 | || invertNow) { |
---|
2268 | //printf("ret after %d pivots\n",factorization_->pivots()); |
---|
2269 | #if CLP_FACTORIZATION_NEW_TIMING > 1 |
---|
2270 | factorization_->statsRefactor('T'); |
---|
2271 | #endif |
---|
2272 | return 1; |
---|
2273 | } else if (forceFactorization_ > 0 && factorization_->pivots() == forceFactorization_) { |
---|
2274 | // relax |
---|
2275 | forceFactorization_ = (3 + 5 * forceFactorization_) / 4; |
---|
2276 | if (forceFactorization_ > factorization_->maximumPivots()) |
---|
2277 | forceFactorization_ = -1; //off |
---|
2278 | #if CLP_FACTORIZATION_NEW_TIMING > 1 |
---|
2279 | factorization_->statsRefactor('F'); |
---|
2280 | #endif |
---|
2281 | return 1; |
---|
2282 | } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type() < 15) { |
---|
2283 | double random = randomNumberGenerator_.randomDouble(); |
---|
2284 | while (random < 0.45) |
---|
2285 | random *= 2.0; |
---|
2286 | int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots); |
---|
2287 | if (factorization_->pivots() >= random * maxNumber) { |
---|
2288 | return 1; |
---|
2289 | } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && numberIterations_ < 1001000 + 10 * (numberRows_ + (numberColumns_ >> 2))) { |
---|
2290 | return 1; |
---|
2291 | } else { |
---|
2292 | // carry on iterating |
---|
2293 | return 0; |
---|
2294 | } |
---|
2295 | } else { |
---|
2296 | // carry on iterating |
---|
2297 | return 0; |
---|
2298 | } |
---|
2299 | } |
---|
2300 | // Copy constructor. |
---|
2301 | ClpSimplex::ClpSimplex(const ClpSimplex &rhs, int scalingMode) |
---|
2302 | : ClpModel(rhs, scalingMode) |
---|
2303 | , bestPossibleImprovement_(0.0) |
---|
2304 | , zeroTolerance_(1.0e-13) |
---|
2305 | , columnPrimalSequence_(-2) |
---|
2306 | , rowPrimalSequence_(-2) |
---|
2307 | , bestObjectiveValue_(rhs.bestObjectiveValue_) |
---|
2308 | , moreSpecialOptions_(2) |
---|
2309 | , baseIteration_(0) |
---|
2310 | , vectorMode_(0) |
---|
2311 | , primalToleranceToGetOptimal_(-1.0) |
---|
2312 | , largeValue_(1.0e15) |
---|
2313 | , largestPrimalError_(0.0) |
---|
2314 | , largestDualError_(0.0) |
---|
2315 | , alphaAccuracy_(-1.0) |
---|
2316 | , dualBound_(1.0e10) |
---|
2317 | , alpha_(0.0) |
---|
2318 | , theta_(0.0) |
---|
2319 | , lowerIn_(0.0) |
---|
2320 | , valueIn_(0.0) |
---|
2321 | , upperIn_(-COIN_DBL_MAX) |
---|
2322 | , dualIn_(0.0) |
---|
2323 | , lowerOut_(-1) |
---|
2324 | , valueOut_(-1) |
---|
2325 | , upperOut_(-1) |
---|
2326 | , dualOut_(-1) |
---|
2327 | , dualTolerance_(1.0e-7) |
---|
2328 | , primalTolerance_(1.0e-7) |
---|
2329 | , sumDualInfeasibilities_(0.0) |
---|
2330 | , sumPrimalInfeasibilities_(0.0) |
---|
2331 | , infeasibilityCost_(1.0e10) |
---|
2332 | , sumOfRelaxedDualInfeasibilities_(0.0) |
---|
2333 | , sumOfRelaxedPrimalInfeasibilities_(0.0) |
---|
2334 | , acceptablePivot_(1.0e-8) |
---|
2335 | , lower_(NULL) |
---|
2336 | , rowLowerWork_(NULL) |
---|
2337 | , columnLowerWork_(NULL) |
---|
2338 | , upper_(NULL) |
---|
2339 | , rowUpperWork_(NULL) |
---|
2340 | , columnUpperWork_(NULL) |
---|
2341 | , cost_(NULL) |
---|
2342 | , rowObjectiveWork_(NULL) |
---|
2343 | , objectiveWork_(NULL) |
---|
2344 | , sequenceIn_(-1) |
---|
2345 | , directionIn_(-1) |
---|
2346 | , sequenceOut_(-1) |
---|
2347 | , directionOut_(-1) |
---|
2348 | , pivotRow_(-1) |
---|
2349 | , lastGoodIteration_(-100) |
---|
2350 | , dj_(NULL) |
---|
2351 | , rowReducedCost_(NULL) |
---|
2352 | , reducedCostWork_(NULL) |
---|
2353 | , solution_(NULL) |
---|
2354 | , rowActivityWork_(NULL) |
---|
2355 | , columnActivityWork_(NULL) |
---|
2356 | , numberDualInfeasibilities_(0) |
---|
2357 | , numberDualInfeasibilitiesWithoutFree_(0) |
---|
2358 | , numberPrimalInfeasibilities_(100) |
---|
2359 | , numberRefinements_(0) |
---|
2360 | , pivotVariable_(NULL) |
---|
2361 | , factorization_(NULL) |
---|
2362 | , savedSolution_(NULL) |
---|
2363 | , numberTimesOptimal_(0) |
---|
2364 | , disasterArea_(NULL) |
---|
2365 | , changeMade_(1) |
---|
2366 | , algorithm_(0) |
---|
2367 | , forceFactorization_(-1) |
---|
2368 | , perturbation_(100) |
---|
2369 | , nonLinearCost_(NULL) |
---|
2370 | , lastBadIteration_(-999999) |
---|
2371 | , lastFlaggedIteration_(-999999) |
---|
2372 | , numberFake_(0) |
---|
2373 | , numberChanged_(0) |
---|
2374 | , progressFlag_(0) |
---|
2375 | , firstFree_(-1) |
---|
2376 | , numberExtraRows_(0) |
---|
2377 | , maximumBasic_(0) |
---|
2378 | , dontFactorizePivots_(0) |
---|
2379 | , incomingInfeasibility_(1.0) |
---|
2380 | , allowedInfeasibility_(10.0) |
---|
2381 | , automaticScale_(0) |
---|
2382 | , maximumPerturbationSize_(0) |
---|
2383 | , perturbationArray_(NULL) |
---|
2384 | , baseModel_(NULL) |
---|
2385 | #ifdef ABC_INHERIT |
---|
2386 | , abcSimplex_(NULL) |
---|
2387 | , abcState_(0) |
---|
2388 | #endif |
---|
2389 | { |
---|
2390 | int i; |
---|
2391 | for (i = 0; i < 6; i++) { |
---|
2392 | rowArray_[i] = NULL; |
---|
2393 | columnArray_[i] = NULL; |
---|
2394 | } |
---|
2395 | for (i = 0; i < 4; i++) { |
---|
2396 | spareIntArray_[i] = 0; |
---|
2397 | spareDoubleArray_[i] = 0.0; |
---|
2398 | } |
---|
2399 | saveStatus_ = NULL; |
---|
2400 | factorization_ = NULL; |
---|
2401 | dualRowPivot_ = NULL; |
---|
2402 | primalColumnPivot_ = NULL; |
---|
2403 | gutsOfDelete(0); |
---|
2404 | delete nonLinearCost_; |
---|
2405 | nonLinearCost_ = NULL; |
---|
2406 | gutsOfCopy(rhs); |
---|
2407 | solveType_ = 1; // say simplex based life form |
---|
2408 | } |
---|
2409 | // Copy constructor from model |
---|
2410 | ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) |
---|
2411 | : ClpModel(rhs, scalingMode) |
---|
2412 | , bestPossibleImprovement_(0.0) |
---|
2413 | , zeroTolerance_(1.0e-13) |
---|
2414 | , columnPrimalSequence_(-2) |
---|
2415 | , rowPrimalSequence_(-2) |
---|
2416 | , bestObjectiveValue_(-COIN_DBL_MAX) |
---|
2417 | , moreSpecialOptions_(2) |
---|
2418 | , baseIteration_(0) |
---|
2419 | , vectorMode_(0) |
---|
2420 | , primalToleranceToGetOptimal_(-1.0) |
---|
2421 | , largeValue_(1.0e15) |
---|
2422 | , largestPrimalError_(0.0) |
---|
2423 | , largestDualError_(0.0) |
---|
2424 | , alphaAccuracy_(-1.0) |
---|
2425 | , dualBound_(1.0e10) |
---|
2426 | , alpha_(0.0) |
---|
2427 | , theta_(0.0) |
---|
2428 | , lowerIn_(0.0) |
---|
2429 | , valueIn_(0.0) |
---|
2430 | , upperIn_(-COIN_DBL_MAX) |
---|
2431 | , dualIn_(0.0) |
---|
2432 | , lowerOut_(-1) |
---|
2433 | , valueOut_(-1) |
---|
2434 | , upperOut_(-1) |
---|
2435 | , dualOut_(-1) |
---|
2436 | , dualTolerance_(1.0e-7) |
---|
2437 | , primalTolerance_(1.0e-7) |
---|
2438 | , sumDualInfeasibilities_(0.0) |
---|
2439 | , sumPrimalInfeasibilities_(0.0) |
---|
2440 | , infeasibilityCost_(1.0e10) |
---|
2441 | , sumOfRelaxedDualInfeasibilities_(0.0) |
---|
2442 | , sumOfRelaxedPrimalInfeasibilities_(0.0) |
---|
2443 | , acceptablePivot_(1.0e-8) |
---|
2444 | , lower_(NULL) |
---|
2445 | , rowLowerWork_(NULL) |
---|
2446 | , columnLowerWork_(NULL) |
---|
2447 | , upper_(NULL) |
---|
2448 | , rowUpperWork_(NULL) |
---|
2449 | , columnUpperWork_(NULL) |
---|
2450 | , cost_(NULL) |
---|
2451 | , rowObjectiveWork_(NULL) |
---|
2452 | , objectiveWork_(NULL) |
---|
2453 | , sequenceIn_(-1) |
---|
2454 | , directionIn_(-1) |
---|
2455 | , sequenceOut_(-1) |
---|
2456 | , directionOut_(-1) |
---|
2457 | , pivotRow_(-1) |
---|
2458 | , lastGoodIteration_(-100) |
---|
2459 | , dj_(NULL) |
---|
2460 | , rowReducedCost_(NULL) |
---|
2461 | , reducedCostWork_(NULL) |
---|
2462 | , solution_(NULL) |
---|
2463 | , rowActivityWork_(NULL) |
---|
2464 | , columnActivityWork_(NULL) |
---|
2465 | , numberDualInfeasibilities_(0) |
---|
2466 | , numberDualInfeasibilitiesWithoutFree_(0) |
---|
2467 | , numberPrimalInfeasibilities_(100) |
---|
2468 | , numberRefinements_(0) |
---|
2469 | , pivotVariable_(NULL) |
---|
2470 | , factorization_(NULL) |
---|
2471 | , savedSolution_(NULL) |
---|
2472 | , numberTimesOptimal_(0) |
---|
2473 | , disasterArea_(NULL) |
---|
2474 | , changeMade_(1) |
---|
2475 | , algorithm_(0) |
---|
2476 | , forceFactorization_(-1) |
---|
2477 | , perturbation_(100) |
---|
2478 | , nonLinearCost_(NULL) |
---|
2479 | , lastBadIteration_(-999999) |
---|
2480 | , lastFlaggedIteration_(-999999) |
---|
2481 | , numberFake_(0) |
---|
2482 | , numberChanged_(0) |
---|
2483 | , progressFlag_(0) |
---|
2484 | , firstFree_(-1) |
---|
2485 | , numberExtraRows_(0) |
---|
2486 | , maximumBasic_(0) |
---|
2487 | , dontFactorizePivots_(0) |
---|
2488 | , incomingInfeasibility_(1.0) |
---|
2489 | , allowedInfeasibility_(10.0) |
---|
2490 | , automaticScale_(0) |
---|
2491 | , maximumPerturbationSize_(0) |
---|
2492 | , perturbationArray_(NULL) |
---|
2493 | , baseModel_(NULL) |
---|
2494 | #ifdef ABC_INHERIT |
---|
2495 | , abcSimplex_(NULL) |
---|
2496 | , abcState_(0) |
---|
2497 | #endif |
---|
2498 | { |
---|
2499 | int i; |
---|
2500 | for (i = 0; i < 6; i++) { |
---|
2501 | rowArray_[i] = NULL; |
---|
2502 | columnArray_[i] = NULL; |
---|
2503 | } |
---|
2504 | for (i = 0; i < 4; i++) { |
---|
2505 | spareIntArray_[i] = 0; |
---|
2506 | spareDoubleArray_[i] = 0.0; |
---|
2507 | } |
---|
2508 | saveStatus_ = NULL; |
---|
2509 | // get an empty factorization so we can set tolerances etc |
---|
2510 | getEmptyFactorization(); |
---|
2511 | // say Steepest pricing |
---|
2512 | dualRowPivot_ = new ClpDualRowSteepest(); |
---|
2513 | // say Steepest pricing |
---|
2514 | primalColumnPivot_ = new ClpPrimalColumnSteepest(); |
---|
2515 | solveType_ = 1; // say simplex based life form |
---|
2516 | } |
---|
2517 | // Assignment operator. This copies the data |
---|
2518 | ClpSimplex & |
---|
2519 | ClpSimplex::operator=(const ClpSimplex &rhs) |
---|
2520 | { |
---|
2521 | if (this != &rhs) { |
---|
2522 | gutsOfDelete(0); |
---|
2523 | delete nonLinearCost_; |
---|
2524 | nonLinearCost_ = NULL; |
---|
2525 | ClpModel::operator=(rhs); |
---|
2526 | gutsOfCopy(rhs); |
---|
2527 | } |
---|
2528 | return *this; |
---|
2529 | } |
---|
2530 | void ClpSimplex::gutsOfCopy(const ClpSimplex &rhs) |
---|
2531 | { |
---|
2532 | assert(numberRows_ == rhs.numberRows_); |
---|
2533 | assert(numberColumns_ == rhs.numberColumns_); |
---|
2534 | numberExtraRows_ = rhs.numberExtraRows_; |
---|
2535 | maximumBasic_ = rhs.maximumBasic_; |
---|
2536 | dontFactorizePivots_ = rhs.dontFactorizePivots_; |
---|
2537 | int numberRows2 = numberRows_ + numberExtraRows_; |
---|
2538 | moreSpecialOptions_ = rhs.moreSpecialOptions_; |
---|
2539 | if ((whatsChanged_ & 1) != 0) { |
---|
2540 | int numberTotal = numberColumns_ + numberRows2; |
---|
2541 | if ((specialOptions_ & 65536) != 0 && maximumRows_ >= 0) { |
---|
2542 | assert(maximumInternalRows_ >= numberRows2); |
---|
2543 | assert(maximumInternalColumns_ >= numberColumns_); |
---|
2544 | numberTotal = 2 * (maximumInternalColumns_ + maximumInternalRows_); |
---|
2545 | } |
---|
2546 | lower_ = ClpCopyOfArray(rhs.lower_, numberTotal); |
---|
2547 | rowLowerWork_ = lower_ + numberColumns_; |
---|
2548 | columnLowerWork_ = lower_; |
---|
2549 | upper_ = ClpCopyOfArray(rhs.upper_, numberTotal); |
---|
2550 | rowUpperWork_ = upper_ + numberColumns_; |
---|
2551 | columnUpperWork_ = upper_; |
---|
2552 | cost_ = ClpCopyOfArray(rhs.cost_, numberTotal); |
---|
2553 | objectiveWork_ = cost_; |
---|
2554 | rowObjectiveWork_ = cost_ + numberColumns_; |
---|
2555 | dj_ = ClpCopyOfArray(rhs.dj_, numberTotal); |
---|
2556 | if (dj_) { |
---|
2557 | reducedCostWork_ = dj_; |
---|
2558 | rowReducedCost_ = dj_ + numberColumns_; |
---|
2559 | } |
---|
2560 | solution_ = ClpCopyOfArray(rhs.solution_, numberTotal); |
---|
2561 | if (solution_) { |
---|
2562 | columnActivityWork_ = solution_; |
---|
2563 | rowActivityWork_ = solution_ + numberColumns_; |
---|
2564 | } |
---|
2565 | if (rhs.pivotVariable_) { |
---|
2566 | pivotVariable_ = new int[numberRows2]; |
---|
2567 | CoinMemcpyN(rhs.pivotVariable_, numberRows2, pivotVariable_); |
---|
2568 | } else { |
---|
2569 | pivotVariable_ = NULL; |
---|
2570 | } |
---|
2571 | savedSolution_ = ClpCopyOfArray(rhs.savedSolution_, numberTotal); |
---|
2572 | int i; |
---|
2573 | for (i = 0; i < 6; i++) { |
---|
2574 | rowArray_[i] = NULL; |
---|
2575 | if (rhs.rowArray_[i]) |
---|
2576 | rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]); |
---|
2577 | columnArray_[i] = NULL; |
---|
2578 | if (rhs.columnArray_[i]) |
---|
2579 | columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]); |
---|
2580 | } |
---|
2581 | if (rhs.saveStatus_) { |
---|
2582 | saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberTotal); |
---|
2583 | } |
---|
2584 | } else { |
---|
2585 | lower_ = NULL; |
---|
2586 | rowLowerWork_ = NULL; |
---|
2587 | columnLowerWork_ = NULL; |
---|
2588 | upper_ = NULL; |
---|
2589 | rowUpperWork_ = NULL; |
---|
2590 | columnUpperWork_ = NULL; |
---|
2591 | cost_ = NULL; |
---|
2592 | objectiveWork_ = NULL; |
---|
2593 | rowObjectiveWork_ = NULL; |
---|
2594 | dj_ = NULL; |
---|
2595 | reducedCostWork_ = NULL; |
---|
2596 | rowReducedCost_ = NULL; |
---|
2597 | solution_ = NULL; |
---|
2598 | columnActivityWork_ = NULL; |
---|
2599 | rowActivityWork_ = NULL; |
---|
2600 | pivotVariable_ = NULL; |
---|
2601 | savedSolution_ = NULL; |
---|
2602 | int i; |
---|
2603 | for (i = 0; i < 6; i++) { |
---|
2604 | rowArray_[i] = NULL; |
---|
2605 | columnArray_[i] = NULL; |
---|
2606 | } |
---|
2607 | saveStatus_ = NULL; |
---|
2608 | } |
---|
2609 | if (rhs.factorization_) { |
---|
2610 | setFactorization(*rhs.factorization_); |
---|
2611 | } else { |
---|
2612 | delete factorization_; |
---|
2613 | factorization_ = NULL; |
---|
2614 | } |
---|
2615 | bestPossibleImprovement_ = rhs.bestPossibleImprovement_; |
---|
2616 | columnPrimalSequence_ = rhs.columnPrimalSequence_; |
---|
2617 | zeroTolerance_ = rhs.zeroTolerance_; |
---|
2618 | rowPrimalSequence_ = rhs.rowPrimalSequence_; |
---|
2619 | bestObjectiveValue_ = rhs.bestObjectiveValue_; |
---|
2620 | baseIteration_ = rhs.baseIteration_; |
---|
2621 | vectorMode_ = rhs.vectorMode_; |
---|
2622 | primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_; |
---|
2623 | largeValue_ = rhs.largeValue_; |
---|
2624 | largestPrimalError_ = rhs.largestPrimalError_; |
---|
2625 | largestDualError_ = rhs.largestDualError_; |
---|
2626 | alphaAccuracy_ = rhs.alphaAccuracy_; |
---|
2627 | dualBound_ = rhs.dualBound_; |
---|
2628 | alpha_ = rhs.alpha_; |
---|
2629 | theta_ = rhs.theta_; |
---|
2630 | lowerIn_ = rhs.lowerIn_; |
---|
2631 | valueIn_ = rhs.valueIn_; |
---|
2632 | upperIn_ = rhs.upperIn_; |
---|
2633 | dualIn_ = rhs.dualIn_; |
---|
2634 | sequenceIn_ = rhs.sequenceIn_; |
---|
2635 | directionIn_ = rhs.directionIn_; |
---|
2636 | lowerOut_ = rhs.lowerOut_; |
---|
2637 | valueOut_ = rhs.valueOut_; |
---|
2638 | upperOut_ = rhs.upperOut_; |
---|
2639 | dualOut_ = rhs.dualOut_; |
---|
2640 | sequenceOut_ = rhs.sequenceOut_; |
---|
2641 | directionOut_ = rhs.directionOut_; |
---|
2642 | pivotRow_ = rhs.pivotRow_; |
---|
2643 | lastGoodIteration_ = rhs.lastGoodIteration_; |
---|
2644 | numberRefinements_ = rhs.numberRefinements_; |
---|
2645 | dualTolerance_ = rhs.dualTolerance_; |
---|
2646 | primalTolerance_ = rhs.primalTolerance_; |
---|
2647 | sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; |
---|
2648 | numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; |
---|
2649 | numberDualInfeasibilitiesWithoutFree_ = rhs.numberDualInfeasibilitiesWithoutFree_; |
---|
2650 | sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; |
---|
2651 | numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; |
---|
2652 | dualRowPivot_ = rhs.dualRowPivot_->clone(true); |
---|
2653 | dualRowPivot_->setModel(this); |
---|
2654 | primalColumnPivot_ = rhs.primalColumnPivot_->clone(true); |
---|
2655 | primalColumnPivot_->setModel(this); |
---|
2656 | numberTimesOptimal_ = rhs.numberTimesOptimal_; |
---|
2657 | disasterArea_ = NULL; |
---|
2658 | changeMade_ = rhs.changeMade_; |
---|
2659 | algorithm_ = rhs.algorithm_; |
---|
2660 | forceFactorization_ = rhs.forceFactorization_; |
---|
2661 | perturbation_ = rhs.perturbation_; |
---|
2662 | infeasibilityCost_ = rhs.infeasibilityCost_; |
---|
2663 | lastBadIteration_ = rhs.lastBadIteration_; |
---|
2664 | lastFlaggedIteration_ = rhs.lastFlaggedIteration_; |
---|
2665 | numberFake_ = rhs.numberFake_; |
---|
2666 | numberChanged_ = rhs.numberChanged_; |
---|
2667 | progressFlag_ = rhs.progressFlag_; |
---|
2668 | firstFree_ = rhs.firstFree_; |
---|
2669 | incomingInfeasibility_ = rhs.incomingInfeasibility_; |
---|
2670 | allowedInfeasibility_ = rhs.allowedInfeasibility_; |
---|
2671 | automaticScale_ = rhs.automaticScale_; |
---|
2672 | #ifdef ABC_INHERIT |
---|
2673 | abcSimplex_ = NULL; |
---|
2674 | abcState_ = rhs.abcState_; |
---|
2675 | #endif |
---|
2676 | maximumPerturbationSize_ = rhs.maximumPerturbationSize_; |
---|
2677 | if (maximumPerturbationSize_ && maximumPerturbationSize_ >= 2 * numberColumns_) { |
---|
2678 | perturbationArray_ = CoinCopyOfArray(rhs.perturbationArray_, |
---|
2679 | maximumPerturbationSize_); |
---|
2680 | } else { |
---|
2681 | maximumPerturbationSize_ = 0; |
---|
2682 | perturbationArray_ = NULL; |
---|
2683 | } |
---|
2684 | if (rhs.baseModel_) { |
---|
2685 | baseModel_ = new ClpSimplex(*rhs.baseModel_); |
---|
2686 | } else { |
---|
2687 | baseModel_ = NULL; |
---|
2688 | } |
---|
2689 | progress_ = rhs.progress_; |
---|
2690 | for (int i = 0; i < 4; i++) { |
---|
2691 | spareIntArray_[i] = rhs.spareIntArray_[i]; |
---|
2692 | spareDoubleArray_[i] = rhs.spareDoubleArray_[i]; |
---|
2693 | } |
---|
2694 | sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; |
---|
2695 | sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; |
---|
2696 | acceptablePivot_ = rhs.acceptablePivot_; |
---|
2697 | if (rhs.nonLinearCost_ != NULL) |
---|
2698 | nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_); |
---|
2699 | else |
---|
2700 | nonLinearCost_ = NULL; |
---|
2701 | solveType_ = rhs.solveType_; |
---|
2702 | eventHandler_->setSimplex(this); |
---|
2703 | } |
---|
2704 | // type == 0 do everything, most + pivot data, 2 factorization data as well |
---|
2705 | void ClpSimplex::gutsOfDelete(int type) |
---|
2706 | { |
---|
2707 | if (!type || (specialOptions_ & 65536) == 0) { |
---|
2708 | maximumInternalColumns_ = -1; |
---|
2709 | maximumInternalRows_ = -1; |
---|
2710 | delete[] lower_; |
---|
2711 | lower_ = NULL; |
---|
2712 | rowLowerWork_ = NULL; |
---|
2713 | columnLowerWork_ = NULL; |
---|
2714 | delete[] upper_; |
---|
2715 | upper_ = NULL; |
---|
2716 | rowUpperWork_ = NULL; |
---|
2717 | columnUpperWork_ = NULL; |
---|
2718 | delete[] cost_; |
---|
2719 | cost_ = NULL; |
---|
2720 | objectiveWork_ = NULL; |
---|
2721 | rowObjectiveWork_ = NULL; |
---|
2722 | delete[] dj_; |
---|
2723 | dj_ = NULL; |
---|
2724 | reducedCostWork_ = NULL; |
---|
2725 | rowReducedCost_ = NULL; |
---|
2726 | delete[] solution_; |
---|
2727 | solution_ = NULL; |
---|
2728 | rowActivityWork_ = NULL; |
---|
2729 | columnActivityWork_ = NULL; |
---|
2730 | delete[] savedSolution_; |
---|
2731 | savedSolution_ = NULL; |
---|
2732 | } |
---|
2733 | if ((specialOptions_ & 2) == 0) { |
---|
2734 | delete nonLinearCost_; |
---|
2735 | nonLinearCost_ = NULL; |
---|
2736 | } |
---|
2737 | int i; |
---|
2738 | if ((specialOptions_ & 65536) == 0) { |
---|
2739 | for (i = 0; i < 6; i++) { |
---|
2740 | delete rowArray_[i]; |
---|
2741 | rowArray_[i] = NULL; |
---|
2742 | delete columnArray_[i]; |
---|
2743 | columnArray_[i] = NULL; |
---|
2744 | } |
---|
2745 | } |
---|
2746 | delete[] saveStatus_; |
---|
2747 | saveStatus_ = NULL; |
---|
2748 | if (type != 1) { |
---|
2749 | delete rowCopy_; |
---|
2750 | rowCopy_ = NULL; |
---|
2751 | } |
---|
2752 | if (!type) { |
---|
2753 | // delete everything |
---|
2754 | setEmptyFactorization(); |
---|
2755 | delete[] pivotVariable_; |
---|
2756 | pivotVariable_ = NULL; |
---|
2757 | delete dualRowPivot_; |
---|
2758 | dualRowPivot_ = NULL; |
---|
2759 | delete primalColumnPivot_; |
---|
2760 | primalColumnPivot_ = NULL; |
---|
2761 | delete baseModel_; |
---|
2762 | baseModel_ = NULL; |
---|
2763 | delete[] perturbationArray_; |
---|
2764 | perturbationArray_ = NULL; |
---|
2765 | maximumPerturbationSize_ = 0; |
---|
2766 | } else { |
---|
2767 | // delete any size information in methods |
---|
2768 | if (type > 1) { |
---|
2769 | //assert (factorization_); |
---|
2770 | if (factorization_) |
---|
2771 | factorization_->clearArrays(); |
---|
2772 | delete[] pivotVariable_; |
---|
2773 | pivotVariable_ = NULL; |
---|
2774 | } |
---|
2775 | dualRowPivot_->clearArrays(); |
---|
2776 | primalColumnPivot_->clearArrays(); |
---|
2777 | } |
---|
2778 | } |
---|
2779 | // This sets largest infeasibility and most infeasible |
---|
2780 | void ClpSimplex::checkPrimalSolution(const double *rowActivities, |
---|
2781 | const double *columnActivities) |
---|
2782 | { |
---|
2783 | double *solution; |
---|
2784 | int iRow, iColumn; |
---|
2785 | |
---|
2786 | objectiveValue_ = 0.0; |
---|
2787 | // now look at primal solution |
---|
2788 | solution = rowActivityWork_; |
---|
2789 | sumPrimalInfeasibilities_ = 0.0; |
---|
2790 | numberPrimalInfeasibilities_ = 0; |
---|
2791 | double primalTolerance = primalTolerance_; |
---|
2792 | double relaxedTolerance = primalTolerance_; |
---|
2793 | // we can't really trust infeasibilities if there is primal error |
---|
2794 | double error = CoinMin(1.0e-2, largestPrimalError_); |
---|
2795 | // allow tolerance at least slightly bigger than standard |
---|
2796 | relaxedTolerance = relaxedTolerance + error; |
---|
2797 | sumOfRelaxedPrimalInfeasibilities_ = 0.0; |
---|
2798 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2799 | //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic); |
---|
2800 | double infeasibility = 0.0; |
---|
2801 | objectiveValue_ += solution[iRow] * rowObjectiveWork_[iRow]; |
---|
2802 | if (solution[iRow] > rowUpperWork_[iRow]) { |
---|
2803 | infeasibility = solution[iRow] - rowUpperWork_[iRow]; |
---|
2804 | } else if (solution[iRow] < rowLowerWork_[iRow]) { |
---|
2805 | infeasibility = rowLowerWork_[iRow] - solution[iRow]; |
---|
2806 | } |
---|
2807 | if (infeasibility > primalTolerance) { |
---|
2808 | sumPrimalInfeasibilities_ += infeasibility - primalTolerance_; |
---|
2809 | if (infeasibility > relaxedTolerance) |
---|
2810 | sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance; |
---|
2811 | numberPrimalInfeasibilities_++; |
---|
2812 | } |
---|
2813 | infeasibility = fabs(rowActivities[iRow] - solution[iRow]); |
---|
2814 | } |
---|
2815 | // Check any infeasibilities from dynamic rows |
---|
2816 | matrix_->primalExpanded(this, 2); |
---|
2817 | solution = columnActivityWork_; |
---|
2818 | if (!matrix_->rhsOffset(this)) { |
---|
2819 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
2820 | //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic); |
---|
2821 | double infeasibility = 0.0; |
---|
2822 | objectiveValue_ += objectiveWork_[iColumn] * solution[iColumn]; |
---|
2823 | if (solution[iColumn] > columnUpperWork_[iColumn]) { |
---|
2824 | infeasibility = solution[iColumn] - columnUpperWork_[iColumn]; |
---|
2825 | } else if (solution[iColumn] < columnLowerWork_[iColumn]) { |
---|
2826 | infeasibility = columnLowerWork_[iColumn] - solution[iColumn]; |
---|
2827 | } |
---|
2828 | if (infeasibility > primalTolerance) { |
---|
2829 | sumPrimalInfeasibilities_ += infeasibility - primalTolerance_; |
---|
2830 | if (infeasibility > relaxedTolerance) |
---|
2831 | sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance; |
---|
2832 | numberPrimalInfeasibilities_++; |
---|
2833 | } |
---|
2834 | infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]); |
---|
2835 | } |
---|
2836 | } else { |
---|
2837 | // as we are using effective rhs we only check basics |
---|
2838 | // But we do need to get objective |
---|
2839 | objectiveValue_ += innerProduct(objectiveWork_, numberColumns_, solution); |
---|
2840 | for (int j = 0; j < numberRows_; j++) { |
---|
2841 | int iColumn = pivotVariable_[j]; |
---|
2842 | //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic); |
---|
2843 | double infeasibility = 0.0; |
---|
2844 | if (solution[iColumn] > columnUpperWork_[iColumn]) { |
---|
2845 | infeasibility = solution[iColumn] - columnUpperWork_[iColumn]; |
---|
2846 | } else if (solution[iColumn] < columnLowerWork_[iColumn]) { |
---|
2847 | infeasibility = columnLowerWork_[iColumn] - solution[iColumn]; |
---|
2848 | } |
---|
2849 | if (infeasibility > primalTolerance) { |
---|
2850 | sumPrimalInfeasibilities_ += infeasibility - primalTolerance_; |
---|
2851 | if (infeasibility > relaxedTolerance) |
---|
2852 | sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance; |
---|
2853 | numberPrimalInfeasibilities_++; |
---|
2854 | } |
---|
2855 | infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]); |
---|
2856 | } |
---|
2857 | } |
---|
2858 | objectiveValue_ += objective_->nonlinearOffset(); |
---|
2859 | objectiveValue_ /= (objectiveScale_ * rhsScale_); |
---|
2860 | } |
---|
2861 | void ClpSimplex::checkDualSolution() |
---|
2862 | { |
---|
2863 | |
---|
2864 | int iRow, iColumn; |
---|
2865 | sumDualInfeasibilities_ = 0.0; |
---|
2866 | numberDualInfeasibilities_ = 0; |
---|
2867 | numberDualInfeasibilitiesWithoutFree_ = 0; |
---|
2868 | if (matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) { |
---|
2869 | // pretend we found dual infeasibilities |
---|
2870 | sumOfRelaxedDualInfeasibilities_ = 1.0; |
---|
2871 | sumDualInfeasibilities_ = 1.0; |
---|
2872 | numberDualInfeasibilities_ = 1; |
---|
2873 | return; |
---|
2874 | } |
---|
2875 | int firstFreePrimal = -1; |
---|
2876 | int firstFreeDual = -1; |
---|
2877 | int numberSuperBasicWithDj = 0; |
---|
2878 | bestPossibleImprovement_ = 0.0; |
---|
2879 | // we can't really trust infeasibilities if there is dual error |
---|
2880 | double error = CoinMin(1.0e-2, largestDualError_); |
---|
2881 | // allow tolerance at least slightly bigger than standard |
---|
2882 | double relaxedTolerance = dualTolerance_ + error; |
---|
2883 | // allow bigger tolerance for possible improvement |
---|
2884 | double possTolerance = 5.0 * relaxedTolerance; |
---|
2885 | sumOfRelaxedDualInfeasibilities_ = 0.0; |
---|
2886 | |
---|
2887 | // Check any djs from dynamic rows |
---|
2888 | matrix_->dualExpanded(this, NULL, NULL, 3); |
---|
2889 | numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_; |
---|
2890 | objectiveValue_ = 0.0; |
---|
2891 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
2892 | objectiveValue_ += objectiveWork_[iColumn] * columnActivityWork_[iColumn]; |
---|
2893 | if (getColumnStatus(iColumn) != basic && !flagged(iColumn)) { |
---|
2894 | // not basic |
---|
2895 | double distanceUp = columnUpperWork_[iColumn] - columnActivityWork_[iColumn]; |
---|
2896 | double distanceDown = columnActivityWork_[iColumn] - columnLowerWork_[iColumn]; |
---|
2897 | if (distanceUp > primalTolerance_) { |
---|
2898 | double value = reducedCostWork_[iColumn]; |
---|
2899 | // Check if "free" |
---|
2900 | if (distanceDown > primalTolerance_) { |
---|
2901 | if (fabs(value) > 1.0e2 * relaxedTolerance) { |
---|
2902 | numberSuperBasicWithDj++; |
---|
2903 | if (firstFreeDual < 0) |
---|
2904 | firstFreeDual = iColumn; |
---|
2905 | } |
---|
2906 | if (firstFreePrimal < 0) |
---|
2907 | firstFreePrimal = iColumn; |
---|
2908 | } |
---|
2909 | // should not be negative |
---|
2910 | if (value < 0.0) { |
---|
2911 | value = -value; |
---|
2912 | if (value > dualTolerance_) { |
---|
2913 | if (getColumnStatus(iColumn) != isFree) { |
---|
2914 | numberDualInfeasibilitiesWithoutFree_++; |
---|
2915 | sumDualInfeasibilities_ += value - dualTolerance_; |
---|
2916 | if (value > possTolerance) |
---|
2917 | bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value; |
---|
2918 | if (value > relaxedTolerance) |
---|
2919 | sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; |
---|
2920 | numberDualInfeasibilities_++; |
---|
2921 | } else { |
---|
2922 | // free so relax a lot |
---|
2923 | value *= 0.01; |
---|
2924 | if (value > dualTolerance_) { |
---|
2925 | sumDualInfeasibilities_ += value - dualTolerance_; |
---|
2926 | if (value > possTolerance) |
---|
2927 | bestPossibleImprovement_ = 1.0e100; |
---|
2928 | if (value > relaxedTolerance) |
---|
2929 | sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; |
---|
2930 | numberDualInfeasibilities_++; |
---|
2931 | } |
---|
2932 | } |
---|
2933 | } |
---|
2934 | } |
---|
2935 | } |
---|
2936 | if (distanceDown > primalTolerance_) { |
---|
2937 | double value = reducedCostWork_[iColumn]; |
---|
2938 | // should not be positive |
---|
2939 | if (value > 0.0) { |
---|
2940 | if (value > dualTolerance_) { |
---|
2941 | sumDualInfeasibilities_ += value - dualTolerance_; |
---|
2942 | if (value > possTolerance) |
---|
2943 | bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10); |
---|
2944 | if (value > relaxedTolerance) |
---|
2945 | sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; |
---|
2946 | numberDualInfeasibilities_++; |
---|
2947 | if (getColumnStatus(iColumn) != isFree) |
---|
2948 | numberDualInfeasibilitiesWithoutFree_++; |
---|
2949 | // maybe we can make feasible by increasing tolerance |
---|
2950 | } |
---|
2951 | } |
---|
2952 | } |
---|
2953 | } |
---|
2954 | } |
---|
2955 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2956 | objectiveValue_ += rowActivityWork_[iRow] * rowObjectiveWork_[iRow]; |
---|
2957 | if (getRowStatus(iRow) != basic && !flagged(iRow + numberColumns_)) { |
---|
2958 | // not basic |
---|
2959 | double distanceUp = rowUpperWork_[iRow] - rowActivityWork_[iRow]; |
---|
2960 | double distanceDown = rowActivityWork_[iRow] - rowLowerWork_[iRow]; |
---|
2961 | if (distanceUp > primalTolerance_) { |
---|
2962 | double value = rowReducedCost_[iRow]; |
---|
2963 | // Check if "free" |
---|
2964 | if (distanceDown > primalTolerance_) { |
---|
2965 | if (fabs(value) > 1.0e2 * relaxedTolerance) { |
---|
2966 | numberSuperBasicWithDj++; |
---|
2967 | if (firstFreeDual < 0) |
---|
2968 | firstFreeDual = iRow + numberColumns_; |
---|
2969 | } |
---|
2970 | if (firstFreePrimal < 0) |
---|
2971 | firstFreePrimal = iRow + numberColumns_; |
---|
2972 | } |
---|
2973 | // should not be negative |
---|
2974 | if (value < 0.0) { |
---|
2975 | value = -value; |
---|
2976 | if (value > dualTolerance_) { |
---|
2977 | sumDualInfeasibilities_ += value - dualTolerance_; |
---|
2978 | if (value > possTolerance) |
---|
2979 | bestPossibleImprovement_ += value * CoinMin(distanceUp, 1.0e10); |
---|
2980 | if (value > relaxedTolerance) |
---|
2981 | sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; |
---|
2982 | numberDualInfeasibilities_++; |
---|
2983 | if (getRowStatus(iRow) != isFree) |
---|
2984 | numberDualInfeasibilitiesWithoutFree_++; |
---|
2985 | } |
---|
2986 | } |
---|
2987 | } |
---|
2988 | if (distanceDown > primalTolerance_) { |
---|
2989 | double value = rowReducedCost_[iRow]; |
---|
2990 | // should not be positive |
---|
2991 | if (value > 0.0) { |
---|
2992 | if (value > dualTolerance_) { |
---|
2993 | sumDualInfeasibilities_ += value - dualTolerance_; |
---|
2994 | if (value > possTolerance) |
---|
2995 | bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10); |
---|
2996 | if (value > relaxedTolerance) |
---|
2997 | sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; |
---|
2998 | numberDualInfeasibilities_++; |
---|
2999 | if (getRowStatus(iRow) != isFree) |
---|
3000 | numberDualInfeasibilitiesWithoutFree_++; |
---|
3001 | // maybe we can make feasible by increasing tolerance |
---|
3002 | } |
---|
3003 | } |
---|
3004 | } |
---|
3005 | } |
---|
3006 | } |
---|
3007 | if (algorithm_ < 0 && firstFreeDual >= 0) { |
---|
3008 | // dual |
---|
3009 | firstFree_ = firstFreeDual; |
---|
3010 | } else if (numberSuperBasicWithDj || (progress_.lastIterationNumber(0) <= 0)) { |
---|
3011 | firstFree_ = firstFreePrimal; |
---|
3012 | } |
---|
3013 | objectiveValue_ += objective_->nonlinearOffset(); |
---|
3014 | objectiveValue_ /= (objectiveScale_ * rhsScale_); |
---|
3015 | } |
---|
3016 | /* This sets sum and number of infeasibilities (Dual and Primal) */ |
---|
3017 | void ClpSimplex::checkBothSolutions() |
---|
3018 | { |
---|
3019 | if ((matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) || matrix_->rhsOffset(this)) { |
---|
3020 | // Say may be free or superbasic |
---|
3021 | moreSpecialOptions_ &= ~8; |
---|
3022 | // old way |
---|
3023 | checkPrimalSolution(rowActivityWork_, columnActivityWork_); |
---|
3024 | checkDualSolution(); |
---|
3025 | return; |
---|
3026 | } |
---|
3027 | int iSequence; |
---|
3028 | assert(dualTolerance_ > 0.0 && dualTolerance_ < 1.0e10); |
---|
3029 | assert(primalTolerance_ > 0.0 && primalTolerance_ < 1.0e10); |
---|
3030 | objectiveValue_ = 0.0; |
---|
3031 | sumPrimalInfeasibilities_ = 0.0; |
---|
3032 | numberPrimalInfeasibilities_ = 0; |
---|
3033 | double primalTolerance = primalTolerance_; |
---|
3034 | double relaxedToleranceP = primalTolerance_; |
---|
3035 | // we can't really trust infeasibilities if there is primal error |
---|
3036 | double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_, 0.0 * primalTolerance_)); |
---|
3037 | // allow tolerance at least slightly bigger than standard |
---|
3038 | relaxedToleranceP = relaxedToleranceP + error; |
---|
3039 | sumOfRelaxedPrimalInfeasibilities_ = 0.0; |
---|
3040 | sumDualInfeasibilities_ = 0.0; |
---|
3041 | numberDualInfeasibilities_ = 0; |
---|
3042 | double dualTolerance = dualTolerance_; |
---|
3043 | double relaxedToleranceD = dualTolerance; |
---|
3044 | // we can't really trust infeasibilities if there is dual error |
---|
3045 | error = CoinMin(1.0e-2, CoinMax(largestDualError_, 5.0 * dualTolerance_)); |
---|
3046 | // allow tolerance at least slightly bigger than standard |
---|
3047 | relaxedToleranceD = relaxedToleranceD + error; |
---|
3048 | // allow bigger tolerance for possible improvement |
---|
3049 | double possTolerance = 5.0 * relaxedToleranceD; |
---|
3050 | sumOfRelaxedDualInfeasibilities_ = 0.0; |
---|
3051 | bestPossibleImprovement_ = 0.0; |
---|
3052 | |
---|
3053 | // Check any infeasibilities from dynamic rows |
---|
3054 | matrix_->primalExpanded(this, 2); |
---|
3055 | // Check any djs from dynamic rows |
---|
3056 | matrix_->dualExpanded(this, NULL, NULL, 3); |
---|
3057 | int numberDualInfeasibilitiesFree = 0; |
---|
3058 | int firstFreePrimal = -1; |
---|
3059 | int firstFreeDual = -1; |
---|
3060 | int numberSuperBasicWithDj = 0; |
---|
3061 | |
---|
3062 | int numberTotal = numberRows_ + numberColumns_; |
---|
3063 | // Say no free or superbasic |
---|
3064 | moreSpecialOptions_ |= 8; |
---|
3065 | //#define PRINT_INFEAS |
---|
3066 | #ifdef PRINT_INFEAS |
---|
3067 | int seqInf[10]; |
---|
3068 | #endif |
---|
3069 | for (iSequence = 0; iSequence < numberTotal; iSequence++) { |
---|
3070 | //#define LIKELY_SUPERBASIC |
---|
3071 | #ifdef LIKELY_SUPERBASIC |
---|
3072 | if (getStatus(iSequence) == isFree || getStatus(iSequence) == superBasic) |
---|
3073 | moreSpecialOptions_ &= ~8; // Say superbasic variables exist |
---|
3074 | #endif |
---|
3075 | double value = solution_[iSequence]; |
---|
3076 | #ifdef COIN_DEBUG |
---|
3077 | if (fabs(value) > 1.0e20) |
---|
3078 | printf("%d values %g %g %g - status %d\n", iSequence, lower_[iSequence], |
---|
3079 | solution_[iSequence], upper_[iSequence], status_[iSequence]); |
---|
3080 | #endif |
---|
3081 | objectiveValue_ += value * cost_[iSequence]; |
---|
3082 | double distanceUp = upper_[iSequence] - value; |
---|
3083 | double distanceDown = value - lower_[iSequence]; |
---|
3084 | if (distanceUp < -primalTolerance) { |
---|
3085 | double infeasibility = -distanceUp; |
---|
3086 | #ifndef LIKELY_SUPERBASIC |
---|
3087 | if (getStatus(iSequence) != basic) |
---|
3088 | moreSpecialOptions_ &= ~8; // Say superbasic variables exist |
---|
3089 | #endif |
---|
3090 | sumPrimalInfeasibilities_ += infeasibility - primalTolerance_; |
---|
3091 | if (infeasibility > relaxedToleranceP) |
---|
3092 | sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP; |
---|
3093 | #ifdef PRINT_INFEAS |
---|
3094 | if (numberPrimalInfeasibilities_ < 10) { |
---|
3095 | seqInf[numberPrimalInfeasibilities_] = iSequence; |
---|
3096 | } |
---|
3097 | #endif |
---|
3098 | numberPrimalInfeasibilities_++; |
---|
3099 | } else if (distanceDown < -primalTolerance) { |
---|
3100 | double infeasibility = -distanceDown; |
---|
3101 | #ifndef LIKELY_SUPERBASIC |
---|
3102 | if (getStatus(iSequence) != basic) |
---|
3103 | moreSpecialOptions_ &= ~8; // Say superbasic variables exist |
---|
3104 | #endif |
---|
3105 | sumPrimalInfeasibilities_ += infeasibility - primalTolerance_; |
---|
3106 | if (infeasibility > relaxedToleranceP) |
---|
3107 | sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP; |
---|
3108 | #ifdef PRINT_INFEAS |
---|
3109 | if (numberPrimalInfeasibilities_ < 10) { |
---|
3110 | seqInf[numberPrimalInfeasibilities_] = iSequence; |
---|
3111 | } |
---|
3112 | #endif |
---|
3113 | numberPrimalInfeasibilities_++; |
---|
3114 | } else { |
---|
3115 | // feasible (so could be free) |
---|
3116 | if (getStatus(iSequence) != basic && !flagged(iSequence)) { |
---|
3117 | // not basic |
---|
3118 | double djValue = dj_[iSequence]; |
---|
3119 | if (distanceDown < primalTolerance) { |
---|
3120 | if (distanceUp > primalTolerance && djValue < -dualTolerance) { |
---|
3121 | sumDualInfeasibilities_ -= djValue + dualTolerance; |
---|
3122 | if (djValue < -possTolerance) |
---|
3123 | bestPossibleImprovement_ -= distanceUp * djValue; |
---|
3124 | if (djValue < -relaxedToleranceD) |
---|
3125 | sumOfRelaxedDualInfeasibilities_ -= djValue + relaxedToleranceD; |
---|
3126 | numberDualInfeasibilities_++; |
---|
3127 | } |
---|
3128 | } else if (distanceUp < primalTolerance) { |
---|
3129 | if (djValue > dualTolerance) { |
---|
3130 | sumDualInfeasibilities_ += djValue - dualTolerance; |
---|
3131 | if (djValue > possTolerance) |
---|
3132 | bestPossibleImprovement_ += distanceDown * djValue; |
---|
3133 | if (djValue > relaxedToleranceD) |
---|
3134 | sumOfRelaxedDualInfeasibilities_ += djValue - relaxedToleranceD; |
---|
3135 | numberDualInfeasibilities_++; |
---|
3136 | } |
---|
3137 | } else { |
---|
3138 | // may be free |
---|
3139 | // Say free or superbasic |
---|
3140 | moreSpecialOptions_ &= ~8; |
---|
3141 | djValue *= 0.01; |
---|
3142 | if (fabs(djValue) > dualTolerance) { |
---|
3143 | if (getStatus(iSequence) == isFree) |
---|
3144 | numberDualInfeasibilitiesFree++; |
---|
3145 | sumDualInfeasibilities_ += fabs(djValue) - dualTolerance; |
---|
3146 | bestPossibleImprovement_ = 1.0e100; |
---|
3147 | numberDualInfeasibilities_++; |
---|
3148 | if (fabs(djValue) > relaxedToleranceD) { |
---|
3149 | sumOfRelaxedDualInfeasibilities_ += value - relaxedToleranceD; |
---|
3150 | numberSuperBasicWithDj++; |
---|
3151 | if (firstFreeDual < 0) |
---|
3152 | firstFreeDual = iSequence; |
---|
3153 | } |
---|
3154 | } |
---|
3155 | if (firstFreePrimal < 0) |
---|
3156 | firstFreePrimal = iSequence; |
---|
3157 | } |
---|
3158 | } |
---|
3159 | } |
---|
3160 | } |
---|
3161 | objectiveValue_ += objective_->nonlinearOffset(); |
---|
3162 | objectiveValue_ /= (objectiveScale_ * rhsScale_); |
---|
3163 | numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberDualInfeasibilitiesFree; |
---|
3164 | #ifdef PRINT_INFEAS |
---|
3165 | if (numberPrimalInfeasibilities_ <= 10) { |
---|
3166 | printf("---------------start-----------\n"); |
---|
3167 | if (!rowScale_) { |
---|
3168 | for (int i = 0; i < numberPrimalInfeasibilities_; i++) { |
---|
3169 | int iSeq = seqInf[i]; |
---|
3170 | double infeas; |
---|
3171 | if (solution_[iSeq] < lower_[iSeq]) |
---|
3172 | infeas = lower_[iSeq] - solution_[iSeq]; |
---|
3173 | else |
---|
3174 | infeas = solution_[iSeq] - upper_[iSeq]; |
---|
3175 | if (iSeq < numberColumns_) { |
---|
3176 | printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g\n", |
---|
3177 | iSeq, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas); |
---|
3178 | } else { |
---|
3179 | printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g\n", |
---|
3180 | iSeq - numberColumns_, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas); |
---|
3181 | } |
---|
3182 | } |
---|
3183 | } else { |
---|
3184 | for (int i = 0; i < numberPrimalInfeasibilities_; i++) { |
---|
3185 | int iSeq = seqInf[i]; |
---|
3186 | double infeas; |
---|
3187 | if (solution_[iSeq] < lower_[iSeq]) |
---|
3188 | infeas = lower_[iSeq] - solution_[iSeq]; |
---|
3189 | else |
---|
3190 | infeas = solution_[iSeq] - upper_[iSeq]; |
---|
3191 | double unscaled = infeas; |
---|
3192 | if (iSeq < numberColumns_) { |
---|
3193 | unscaled *= columnScale_[iSeq]; |
---|
3194 | printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n", |
---|
3195 | iSeq, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas, unscaled); |
---|
3196 | } else { |
---|
3197 | unscaled /= rowScale_[iSeq - numberColumns_]; |
---|
3198 | printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n", |
---|
3199 | iSeq - numberColumns_, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas, unscaled); |
---|
3200 | } |
---|
3201 | } |
---|
3202 | } |
---|
3203 | } |
---|
3204 | #endif |
---|
3205 | if (algorithm_ < 0 && firstFreeDual >= 0) { |
---|
3206 | // dual |
---|
3207 | firstFree_ = firstFreeDual; |
---|
3208 | } else if (numberSuperBasicWithDj || (progress_.lastIterationNumber(0) <= 0)) { |
---|
3209 | firstFree_ = firstFreePrimal; |
---|
3210 | } |
---|
3211 | } |
---|
3212 | /* Adds multiple of a column into an array */ |
---|
3213 | void ClpSimplex::add(double *array, |
---|
3214 | int sequence, double multiplier) const |
---|
3215 | { |
---|
3216 | if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) { |
---|
3217 | //slack |
---|
3218 | array[sequence - numberColumns_] -= multiplier; |
---|
3219 | } else { |
---|
3220 | // column |
---|
3221 | matrix_->add(this, array, sequence, multiplier); |
---|
3222 | } |
---|
3223 | } |
---|
3224 | /* |
---|
3225 | Unpacks one column of the matrix into indexed array |
---|
3226 | */ |
---|
3227 | void ClpSimplex::unpack(CoinIndexedVector *rowArray) const |
---|
3228 | { |
---|
3229 | rowArray->clear(); |
---|
3230 | if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) { |
---|
3231 | //slack |
---|
3232 | rowArray->insert(sequenceIn_ - numberColumns_, -1.0); |
---|
3233 | } else { |
---|
3234 | // column |
---|
3235 | matrix_->unpack(this, rowArray, sequenceIn_); |
---|
3236 | } |
---|
3237 | } |
---|
3238 | void ClpSimplex::unpack(CoinIndexedVector *rowArray, int sequence) const |
---|
3239 | { |
---|
3240 | rowArray->clear(); |
---|
3241 | if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) { |
---|
3242 | //slack |
---|
3243 | rowArray->insert(sequence - numberColumns_, -1.0); |
---|
3244 | } else { |
---|
3245 | // column |
---|
3246 | matrix_->unpack(this, rowArray, sequence); |
---|
3247 | } |
---|
3248 | } |
---|
3249 | /* |
---|
3250 | Unpacks one column of the matrix into indexed array |
---|
3251 | */ |
---|
3252 | void ClpSimplex::unpackPacked(CoinIndexedVector *rowArray) |
---|
3253 | { |
---|
3254 | rowArray->clear(); |
---|
3255 | if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) { |
---|
3256 | //slack |
---|
3257 | int *index = rowArray->getIndices(); |
---|
3258 | double *array = rowArray->denseVector(); |
---|
3259 | array[0] = -1.0; |
---|
3260 | index[0] = sequenceIn_ - numberColumns_; |
---|
3261 | rowArray->setNumElements(1); |
---|
3262 | rowArray->setPackedMode(true); |
---|
3263 | } else { |
---|
3264 | // column |
---|
3265 | matrix_->unpackPacked(this, rowArray, sequenceIn_); |
---|
3266 | } |
---|
3267 | } |
---|
3268 | void ClpSimplex::unpackPacked(CoinIndexedVector *rowArray, int sequence) |
---|
3269 | { |
---|
3270 | rowArray->clear(); |
---|
3271 | if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) { |
---|
3272 | //slack |
---|
3273 | int *index = rowArray->getIndices(); |
---|
3274 | double *array = rowArray->denseVector(); |
---|
3275 | array[0] = -1.0; |
---|
3276 | index[0] = sequence - numberColumns_; |
---|
3277 | rowArray->setNumElements(1); |
---|
3278 | rowArray->setPackedMode(true); |
---|
3279 | } else { |
---|
3280 | // column |
---|
3281 | matrix_->unpackPacked(this, rowArray, sequence); |
---|
3282 | } |
---|
3283 | } |
---|
3284 | //static int x_gaps[4]={0,0,0,0}; |
---|
3285 | //static int scale_times[]={0,0,0,0}; |
---|
3286 | bool ClpSimplex::createRim(int what, bool makeRowCopy, int startFinishOptions) |
---|
3287 | { |
---|
3288 | bool goodMatrix = true; |
---|
3289 | int saveLevel = handler_->logLevel(); |
---|
3290 | spareIntArray_[0] = 0; |
---|
3291 | if (!matrix_->canGetRowCopy()) |
---|
3292 | makeRowCopy = false; // switch off row copy if can't produce |
---|
3293 | // Arrays will be there and correct size unless what is 63 |
---|
3294 | bool newArrays = (what == 63); |
---|
3295 | // We may be restarting with same size |
---|
3296 | bool keepPivots = false; |
---|
3297 | if (startFinishOptions == -1) { |
---|
3298 | startFinishOptions = 0; |
---|
3299 | keepPivots = true; |
---|
3300 | } |
---|
3301 | bool oldMatrix = ((startFinishOptions & 4) != 0 && (whatsChanged_ & 1) != 0); |
---|
3302 | if (what == 63) { |
---|
3303 | pivotRow_ = -1; |
---|
3304 | if (!status_) |
---|
3305 | createStatus(); |
---|
3306 | if (oldMatrix) |
---|
3307 | newArrays = false; |
---|
3308 | if (problemStatus_ == 10) { |
---|
3309 | handler_->setLogLevel(0); // switch off messages |
---|
3310 | if (rowArray_[0]) { |
---|
3311 | // stuff is still there |
---|
3312 | oldMatrix = true; |
---|
3313 | newArrays = false; |
---|
3314 | keepPivots = true; |
---|
3315 | for (int iRow = 0; iRow < 4; iRow++) { |
---|
3316 | rowArray_[iRow]->clear(); |
---|
3317 | } |
---|
3318 | for (int iColumn = 0; iColumn < SHORT_REGION; iColumn++) { |
---|
3319 | columnArray_[iColumn]->clear(); |
---|
3320 | } |
---|
3321 | } |
---|
3322 | } else if (factorization_) { |
---|
3323 | // match up factorization messages |
---|
3324 | if (handler_->logLevel() < 3) |
---|
3325 | factorization_->messageLevel(0); |
---|
3326 | else |
---|
3327 | factorization_->messageLevel(CoinMax(3, factorization_->messageLevel())); |
---|
3328 | /* Faster to keep pivots rather than re-scan matrix. Matrix may have changed |
---|
3329 | i.e. oldMatrix false but okay as long as same number rows and status array exists |
---|
3330 | */ |
---|
3331 | if ((startFinishOptions & 2) != 0 && factorization_->numberRows() == numberRows_ && status_) |
---|
3332 | keepPivots = true; |
---|
3333 | } |
---|
3334 | numberExtraRows_ = matrix_->generalExpanded(this, 2, maximumBasic_); |
---|
3335 | if (numberExtraRows_ && newArrays) { |
---|
3336 | // make sure status array large enough |
---|
3337 | assert(status_); |
---|
3338 | int numberOld = numberRows_ + numberColumns_; |
---|
3339 | int numberNew = numberRows_ + numberColumns_ + numberExtraRows_; |
---|
3340 | unsigned char *newStatus = new unsigned char[numberNew]; |
---|
3341 | memset(newStatus + numberOld, 0, numberExtraRows_); |
---|
3342 | CoinMemcpyN(status_, numberOld, newStatus); |
---|
3343 | delete[] status_; |
---|
3344 | status_ = newStatus; |
---|
3345 | } |
---|
3346 | } |
---|
3347 | int numberRows2 = numberRows_ + numberExtraRows_; |
---|
3348 | int numberTotal = numberRows2 + numberColumns_; |
---|
3349 | if ((specialOptions_ & 65536) != 0) { |
---|
3350 | assert(!numberExtraRows_); |
---|
3351 | if (!cost_ || numberRows2 > maximumInternalRows_ || numberColumns_ > maximumInternalColumns_) { |
---|
3352 | newArrays = true; |
---|
3353 | keepPivots = false; |
---|
3354 | COIN_DETAIL_PRINT(printf("createrim a %d rows, %d maximum rows %d maxinternal\n", |
---|
3355 | numberRows_, maximumRows_, maximumInternalRows_)); |
---|
3356 | int oldMaximumRows = maximumInternalRows_; |
---|
3357 | int oldMaximumColumns = maximumInternalColumns_; |
---|
3358 | if (cost_) { |
---|
3359 | if (numberRows2 > maximumInternalRows_) |
---|
3360 | maximumInternalRows_ = numberRows2; |
---|
3361 | if (numberColumns_ > maximumInternalColumns_) |
---|
3362 | maximumInternalColumns_ = numberColumns_; |
---|
3363 | } else { |
---|
3364 | maximumInternalRows_ = numberRows2; |
---|
3365 | maximumInternalColumns_ = numberColumns_; |
---|
3366 | } |
---|
3367 | //maximumRows_=CoinMax(maximumInternalRows_,maximumRows_); |
---|
3368 | //maximumColumns_=CoinMax(maximumInternalColumns_,maximumColumns_); |
---|
3369 | assert(maximumInternalRows_ == maximumRows_); |
---|
3370 | assert(maximumInternalColumns_ == maximumColumns_); |
---|
3371 | COIN_DETAIL_PRINT(printf("createrim b %d rows, %d maximum rows, %d maxinternal\n", |
---|
3372 | numberRows_, maximumRows_, maximumInternalRows_)); |
---|
3373 | int numberTotal2 = (maximumInternalRows_ + maximumInternalColumns_) * 2; |
---|
3374 | delete[] cost_; |
---|
3375 | cost_ = new double[numberTotal2]; |
---|
3376 | delete[] lower_; |
---|
3377 | delete[] upper_; |
---|
3378 | lower_ = new double[numberTotal2]; |
---|
3379 | upper_ = new double[numberTotal2]; |
---|
3380 | delete[] dj_; |
---|
3381 | dj_ = new double[numberTotal2]; |
---|
3382 | delete[] solution_; |
---|
3383 | solution_ = new double[numberTotal2]; |
---|
3384 | // ***** should be non NULL but seems to be too much |
---|
3385 | //printf("resize %d savedRowScale %x\n",maximumRows_,savedRowScale_); |
---|
3386 | if (savedRowScale_) { |
---|
3387 | assert(oldMaximumRows > 0); |
---|
3388 | double *temp; |
---|
3389 | temp = new double[4 * maximumRows_]; |
---|
3390 | CoinFillN(temp, 4 * maximumRows_, 1.0); |
---|
3391 | CoinMemcpyN(savedRowScale_, numberRows_, temp); |
---|
3392 | CoinMemcpyN(savedRowScale_ + oldMaximumRows, numberRows_, temp + maximumRows_); |
---|
3393 | CoinMemcpyN(savedRowScale_ + 2 * oldMaximumRows, numberRows_, temp + 2 * maximumRows_); |
---|
3394 | CoinMemcpyN(savedRowScale_ + 3 * oldMaximumRows, numberRows_, temp + 3 * maximumRows_); |
---|
3395 | delete[] savedRowScale_; |
---|
3396 | savedRowScale_ = temp; |
---|
3397 | temp = new double[4 * maximumColumns_]; |
---|
3398 | CoinFillN(temp, 4 * maximumColumns_, 1.0); |
---|
3399 | CoinMemcpyN(savedColumnScale_, numberColumns_, temp); |
---|
3400 | CoinMemcpyN(savedColumnScale_ + oldMaximumColumns, numberColumns_, temp + maximumColumns_); |
---|
3401 | CoinMemcpyN(savedColumnScale_ + 2 * oldMaximumColumns, numberColumns_, temp + 2 * maximumColumns_); |
---|
3402 | CoinMemcpyN(savedColumnScale_ + 3 * oldMaximumColumns, numberColumns_, temp + 3 * maximumColumns_); |
---|
3403 | delete[] savedColumnScale_; |
---|
3404 | savedColumnScale_ = temp; |
---|
3405 | } |
---|
3406 | } |
---|
3407 | } |
---|
3408 | int i; |
---|
3409 | bool doSanityCheck = true; |
---|
3410 | if (what == 63) { |
---|
3411 | // We may want to switch stuff off for speed |
---|
3412 | if ((specialOptions_ & 256) != 0) |
---|
3413 | makeRowCopy = false; // no row copy |
---|
3414 | if ((specialOptions_ & 128) != 0) |
---|
3415 | doSanityCheck = false; // no sanity check |
---|
3416 | //check matrix |
---|
3417 | if (!matrix_) |
---|
3418 | matrix_ = new ClpPackedMatrix(); |
---|
3419 | int checkType = (doSanityCheck) ? 15 : 14; |
---|
3420 | if (oldMatrix) |
---|
3421 | checkType = 14; |
---|
3422 | bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0; |
---|
3423 | if (inCbcOrOther) |
---|
3424 | checkType -= 4; // don't check for duplicates |
---|
3425 | if (!matrix_->allElementsInRange(this, smallElement_, 1.0e20, checkType)) { |
---|
3426 | problemStatus_ = 4; |
---|
3427 | secondaryStatus_ = 8; |
---|
3428 | //goodMatrix= false; |
---|
3429 | return false; |
---|
3430 | } |
---|
3431 | bool rowCopyIsScaled; |
---|
3432 | if (makeRowCopy) { |
---|
3433 | if (!oldMatrix || !rowCopy_) { |
---|
3434 | delete rowCopy_; |
---|
3435 | // may return NULL if can't give row copy |
---|
3436 | rowCopy_ = matrix_->reverseOrderedCopy(); |
---|
3437 | rowCopyIsScaled = false; |
---|
3438 | } else { |
---|
3439 | rowCopyIsScaled = true; |
---|
3440 | } |
---|
3441 | } |
---|
3442 | #if 0 |
---|
3443 | if (what == 63) { |
---|
3444 | int k = rowScale_ ? 1 : 0; |
---|
3445 | if (oldMatrix) |
---|
3446 | k += 2; |
---|
3447 | scale_times[k]++; |
---|
3448 | if ((scale_times[0] + scale_times[1] + scale_times[2] + scale_times[3]) % 1000 == 0) |
---|
3449 | printf("scale counts %d %d %d %d\n", |
---|
3450 | scale_times[0], scale_times[1], scale_times[2], scale_times[3]); |
---|
3451 | } |
---|
3452 | #endif |
---|
3453 | // do scaling if needed |
---|
3454 | if (!oldMatrix && scalingFlag_ < 0) { |
---|
3455 | if (scalingFlag_ < 0 && rowScale_) { |
---|
3456 | //if (handler_->logLevel()>0) |
---|
3457 | printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n", |
---|
3458 | scalingFlag_); |
---|
3459 | scalingFlag_ = 0; |
---|
3460 | } |
---|
3461 | delete[] rowScale_; |
---|
3462 | delete[] columnScale_; |
---|
3463 | rowScale_ = NULL; |
---|
3464 | columnScale_ = NULL; |
---|
3465 | } |
---|
3466 | inverseRowScale_ = NULL; |
---|
3467 | inverseColumnScale_ = NULL; |
---|
3468 | if (scalingFlag_ > 0 && (specialOptions_ & 65536) != 0 && rowScale_ && rowScale_ == savedRowScale_) |
---|
3469 | rowScale_ = NULL; |
---|
3470 | if (scalingFlag_ > 0 && !rowScale_) { |
---|
3471 | if ((specialOptions_ & 65536) != 0) { |
---|
3472 | assert(!rowScale_); |
---|
3473 | rowScale_ = savedRowScale_; |
---|
3474 | columnScale_ = savedColumnScale_; |
---|
3475 | // put back original |
---|
3476 | if (savedRowScale_) { |
---|
3477 | inverseRowScale_ = savedRowScale_ + maximumInternalRows_; |
---|
3478 | inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_; |
---|
3479 | CoinMemcpyN(savedRowScale_ + 2 * maximumInternalRows_, |
---|
3480 | numberRows2, savedRowScale_); |
---|
3481 | CoinMemcpyN(savedRowScale_ + 3 * maximumInternalRows_, |
---|
3482 | numberRows2, inverseRowScale_); |
---|
3483 | CoinMemcpyN(savedColumnScale_ + 2 * maximumColumns_, |
---|
3484 | numberColumns_, savedColumnScale_); |
---|
3485 | CoinMemcpyN(savedColumnScale_ + 3 * maximumColumns_, |
---|
3486 | numberColumns_, inverseColumnScale_); |
---|
3487 | } |
---|
3488 | } |
---|
3489 | if (matrix_->scale(this, this)) |
---|
3490 | scalingFlag_ = -scalingFlag_; // not scaled after all |
---|
3491 | if (rowScale_ && automaticScale_) { |
---|
3492 | if (!savedRowScale_) { |
---|
3493 | inverseRowScale_ = rowScale_ + numberRows2; |
---|
3494 | inverseColumnScale_ = columnScale_ + numberColumns_; |
---|
3495 | } |
---|
3496 | // try automatic scaling |
---|
3497 | double smallestObj = 1.0e100; |
---|
3498 | double largestObj = 0.0; |
---|
3499 | double largestRhs = 0.0; |
---|
3500 | const double *obj = objective(); |
---|
3501 | for (i = 0; i < numberColumns_; i++) { |
---|
3502 | double value = fabs(obj[i]); |
---|
3503 | value *= columnScale_[i]; |
---|
3504 | if (value && columnLower_[i] != columnUpper_[i]) { |
---|
3505 | smallestObj = CoinMin(smallestObj, value); |
---|
3506 | largestObj = CoinMax(largestObj, value); |
---|
3507 | } |
---|
3508 | if (columnLower_[i] > 0.0 || columnUpper_[i] < 0.0) { |
---|
3509 | double scale = 1.0 * inverseColumnScale_[i]; |
---|
3510 | //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs); |
---|
3511 | if (columnLower_[i] > 0) |
---|
3512 | largestRhs = CoinMax(largestRhs, columnLower_[i] * scale); |
---|
3513 | if (columnUpper_[i] < 0.0) |
---|
3514 | largestRhs = CoinMax(largestRhs, -columnUpper_[i] * scale); |
---|
3515 | } |
---|
3516 | } |
---|
3517 | for (i = 0; i < numberRows_; i++) { |
---|
3518 | if (rowLower_[i] > 0.0 || rowUpper_[i] < 0.0) { |
---|
3519 | double scale = rowScale_[i]; |
---|
3520 | //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs); |
---|
3521 | if (rowLower_[i] > 0) |
---|
3522 | largestRhs = CoinMax(largestRhs, rowLower_[i] * scale); |
---|
3523 | if (rowUpper_[i] < 0.0) |
---|
3524 | largestRhs = CoinMax(largestRhs, -rowUpper_[i] * scale); |
---|
3525 | } |
---|
3526 | } |
---|
3527 | COIN_DETAIL_PRINT(printf("small obj %g, large %g - rhs %g\n", smallestObj, largestObj, largestRhs)); |
---|
3528 | bool scalingDone = false; |
---|
3529 | // look at element range |
---|
3530 | double smallestNegative; |
---|
3531 | double largestNegative; |
---|
3532 | double smallestPositive; |
---|
3533 | double largestPositive; |
---|
3534 | matrix_->rangeOfElements(smallestNegative, largestNegative, |
---|
3535 | smallestPositive, largestPositive); |
---|
3536 | smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive); |
---|
3537 | largestPositive = CoinMax(fabs(largestNegative), largestPositive); |
---|
3538 | if (largestObj) { |
---|
3539 | double ratio = largestObj / smallestObj; |
---|
3540 | double scale = 1.0; |
---|
3541 | if (ratio < 1.0e8) { |
---|
3542 | // reasonable |
---|
3543 | if (smallestObj < 1.0e-4) { |
---|
3544 | // may as well scale up |
---|
3545 | scalingDone = true; |
---|
3546 | scale = 1.0e-3 / smallestObj; |
---|
3547 | } else if (largestObj < 1.0e6 || (algorithm_ > 0 && largestObj < 1.0e-4 * infeasibilityCost_)) { |
---|
3548 | //done=true; |
---|
3549 | } else { |
---|
3550 | scalingDone = true; |
---|
3551 | if (algorithm_ < 0) { |
---|
3552 | scale = 1.0e6 / largestObj; |
---|
3553 | } else { |
---|
3554 | scale = CoinMax(1.0e6, 1.0e-4 * infeasibilityCost_) / largestObj; |
---|
3555 | } |
---|
3556 | } |
---|
3557 | } else if (ratio < 1.0e12) { |
---|
3558 | // not so good |
---|
3559 | if (smallestObj < 1.0e-7) { |
---|
3560 | // may as well scale up |
---|
3561 | scalingDone = true; |
---|
3562 | scale = 1.0e-6 / smallestObj; |
---|
3563 | } else if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) { |
---|
3564 | //done=true; |
---|
3565 | } else { |
---|
3566 | scalingDone = true; |
---|
3567 | if (algorithm_ < 0) { |
---|
3568 | scale = 1.0e7 / largestObj; |
---|
3569 | } else { |
---|
3570 | scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj; |
---|
3571 | } |
---|
3572 | } |
---|
3573 | } else { |
---|
3574 | // Really nasty problem |
---|
3575 | if (smallestObj < 1.0e-8) { |
---|
3576 | // may as well scale up |
---|
3577 | scalingDone = true; |
---|
3578 | scale = 1.0e-7 / smallestObj; |
---|
3579 | largestObj *= scale; |
---|
3580 | } |
---|
3581 | if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) { |
---|
3582 | //done=true; |
---|
3583 | } else { |
---|
3584 | scalingDone = true; |
---|
3585 | if (algorithm_ < 0) { |
---|
3586 | scale = 1.0e7 / largestObj; |
---|
3587 | } else { |
---|
3588 | scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj; |
---|
3589 | } |
---|
3590 | } |
---|
3591 | } |
---|
3592 | objectiveScale_ = scale; |
---|
3593 | } |
---|
3594 | if (largestRhs > 1.0e12) { |
---|
3595 | scalingDone = true; |
---|
3596 | rhsScale_ = 1.0e9 / largestRhs; |
---|
3597 | } else if (largestPositive > 1.0e-14 * smallestPositive && largestRhs > 1.0e6) { |
---|
3598 | scalingDone = true; |
---|
3599 | rhsScale_ = 1.0e6 / largestRhs; |
---|
3600 | } else { |
---|
3601 | rhsScale_ = 1.0; |
---|
3602 | } |
---|
3603 | if (scalingDone) { |
---|
3604 | handler_->message(CLP_RIM_SCALE, messages_) |
---|
3605 | << objectiveScale_ << rhsScale_ |
---|
3606 | << CoinMessageEol; |
---|
3607 | } |
---|
3608 | } |
---|
3609 | } else if (makeRowCopy && scalingFlag_ > 0 && !rowCopyIsScaled) { |
---|
3610 | matrix_->scaleRowCopy(this); |
---|
3611 | } |
---|
3612 | if (rowScale_ && !savedRowScale_) { |
---|
3613 | inverseRowScale_ = rowScale_ + numberRows2; |
---|
3614 | inverseColumnScale_ = columnScale_ + numberColumns_; |
---|
3615 | } |
---|
3616 | // See if we can try for faster row copy |
---|
3617 | if (makeRowCopy && !oldMatrix) { |
---|
3618 | ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix_); |
---|
3619 | if (clpMatrix && numberThreads_) |
---|
3620 | clpMatrix->specialRowCopy(this, rowCopy_); |
---|
3621 | if (clpMatrix) |
---|
3622 | clpMatrix->specialColumnCopy(this); |
---|
3623 | } |
---|
3624 | } |
---|
3625 | if (what == 63) { |
---|
3626 | #if 0 |
---|
3627 | { |
---|
3628 | x_gaps[0]++; |
---|
3629 | ClpPackedMatrix* clpMatrix = |
---|
3630 | dynamic_cast< ClpPackedMatrix*>(matrix_); |
---|
3631 | if (clpMatrix) { |
---|
3632 | if (!clpMatrix->getPackedMatrix()->hasGaps()) |
---|
3633 | x_gaps[1]++; |
---|
3634 | if ((clpMatrix->flags() & 2) == 0) |
---|
3635 | x_gaps[3]++; |
---|
3636 | } else { |
---|
3637 | x_gaps[2]++; |
---|
3638 | } |
---|
3639 | if ((x_gaps[0] % 1000) == 0) |
---|
3640 | printf("create %d times, no gaps %d times - not clp %d times - flagged %d\n", |
---|
3641 | x_gaps[0], x_gaps[1], x_gaps[2], x_gaps[3]); |
---|
3642 | } |
---|
3643 | #endif |
---|
3644 | if (newArrays && (specialOptions_ & 65536) == 0) { |
---|
3645 | delete[] cost_; |
---|
3646 | cost_ = new double[2 * numberTotal]; |
---|
3647 | delete[] lower_; |
---|
3648 | delete[] upper_; |
---|
3649 | lower_ = new double[numberTotal]; |
---|
3650 | upper_ = new double[numberTotal]; |
---|
3651 | delete[] dj_; |
---|
3652 | dj_ = new double[numberTotal]; |
---|
3653 | delete[] solution_; |
---|
3654 | solution_ = new double[numberTotal]; |
---|
3655 | } |
---|
3656 | reducedCostWork_ = dj_; |
---|
3657 | rowReducedCost_ = dj_ + numberColumns_; |
---|
3658 | columnActivityWork_ = solution_; |
---|
3659 | rowActivityWork_ = solution_ + numberColumns_; |
---|
3660 | objectiveWork_ = cost_; |
---|
3661 | rowObjectiveWork_ = cost_ + numberColumns_; |
---|
3662 | rowLowerWork_ = lower_ + numberColumns_; |
---|
3663 | columnLowerWork_ = lower_; |
---|
3664 | rowUpperWork_ = upper_ + numberColumns_; |
---|
3665 | columnUpperWork_ = upper_; |
---|
3666 | } |
---|
3667 | if ((what & 4) != 0) { |
---|
3668 | double direction = optimizationDirection_ * objectiveScale_; |
---|
3669 | const double *obj = objective(); |
---|
3670 | const double *rowScale = rowScale_; |
---|
3671 | const double *columnScale = columnScale_; |
---|
3672 | // and also scale by scale factors |
---|
3673 | if (rowScale) { |
---|
3674 | if (rowObjective_) { |
---|
3675 | for (i = 0; i < numberRows_; i++) |
---|
3676 | rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i]; |
---|
3677 | } else { |
---|
3678 | memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double)); |
---|
3679 | } |
---|
3680 | // If scaled then do all columns later in one loop |
---|
3681 | if (what != 63) { |
---|
3682 | for (i = 0; i < numberColumns_; i++) { |
---|
3683 | CoinAssert(fabs(obj[i]) < 1.0e25); |
---|
3684 | objectiveWork_[i] = obj[i] * direction * columnScale[i]; |
---|
3685 | } |
---|
3686 | } |
---|
3687 | } else { |
---|
3688 | if (rowObjective_) { |
---|
3689 | for (i = 0; i < numberRows_; i++) |
---|
3690 | rowObjectiveWork_[i] = rowObjective_[i] * direction; |
---|
3691 | } else { |
---|
3692 | memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double)); |
---|
3693 | } |
---|
3694 | for (i = 0; i < numberColumns_; i++) { |
---|
3695 | CoinAssert(fabs(obj[i]) < 1.0e25); |
---|
3696 | objectiveWork_[i] = obj[i] * direction; |
---|
3697 | } |
---|
3698 | } |
---|
3699 | } |
---|
3700 | if ((what & 1) != 0) { |
---|
3701 | const double *rowScale = rowScale_; |
---|
3702 | // clean up any mismatches on infinity |
---|
3703 | // and fix any variables with tiny gaps |
---|
3704 | double primalTolerance = dblParam_[ClpPrimalTolerance]; |
---|
3705 | if (rowScale) { |
---|
3706 | // If scaled then do all columns later in one loop |
---|
3707 | if (what != 63) { |
---|
3708 | const double *inverseScale = inverseColumnScale_; |
---|
3709 | for (i = 0; i < numberColumns_; i++) { |
---|
3710 | double multiplier = rhsScale_ * inverseScale[i]; |
---|
3711 | double lowerValue = columnLower_[i]; |
---|
3712 | double upperValue = columnUpper_[i]; |
---|
3713 | if (lowerValue > -1.0e20) { |
---|
3714 | columnLowerWork_[i] = lowerValue * multiplier; |
---|
3715 | if (upperValue >= 1.0e20) { |
---|
3716 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
3717 | } else { |
---|
3718 | columnUpperWork_[i] = upperValue * multiplier; |
---|
3719 | if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) { |
---|
3720 | if (columnLowerWork_[i] >= 0.0) { |
---|
3721 | columnUpperWork_[i] = columnLowerWork_[i]; |
---|
3722 | } else if (columnUpperWork_[i] <= 0.0) { |
---|
3723 | columnLowerWork_[i] = columnUpperWork_[i]; |
---|
3724 | } else { |
---|
3725 | columnUpperWork_[i] = 0.0; |
---|
3726 | columnLowerWork_[i] = 0.0; |
---|
3727 | } |
---|
3728 | } |
---|
3729 | } |
---|
3730 | } else if (upperValue < 1.0e20) { |
---|
3731 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
3732 | columnUpperWork_[i] = upperValue * multiplier; |
---|
3733 | } else { |
---|
3734 | // free |
---|
3735 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
3736 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
3737 | } |
---|
3738 | } |
---|
3739 | } |
---|
3740 | for (i = 0; i < numberRows_; i++) { |
---|
3741 | double multiplier = rhsScale_ * rowScale[i]; |
---|
3742 | double lowerValue = rowLower_[i]; |
---|
3743 | double upperValue = rowUpper_[i]; |
---|
3744 | if (lowerValue > -1.0e20) { |
---|
3745 | rowLowerWork_[i] = lowerValue * multiplier; |
---|
3746 | if (upperValue >= 1.0e20) { |
---|
3747 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
3748 | } else { |
---|
3749 | rowUpperWork_[i] = upperValue * multiplier; |
---|
3750 | if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) { |
---|
3751 | if (rowLowerWork_[i] >= 0.0) { |
---|
3752 | rowUpperWork_[i] = rowLowerWork_[i]; |
---|
3753 | } else if (rowUpperWork_[i] <= 0.0) { |
---|
3754 | rowLowerWork_[i] = rowUpperWork_[i]; |
---|
3755 | } else { |
---|
3756 | rowUpperWork_[i] = 0.0; |
---|
3757 | rowLowerWork_[i] = 0.0; |
---|
3758 | } |
---|
3759 | } |
---|
3760 | } |
---|
3761 | } else if (upperValue < 1.0e20) { |
---|
3762 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
3763 | rowUpperWork_[i] = upperValue * multiplier; |
---|
3764 | } else { |
---|
3765 | // free |
---|
3766 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
3767 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
3768 | } |
---|
3769 | } |
---|
3770 | } else if (rhsScale_ != 1.0) { |
---|
3771 | for (i = 0; i < numberColumns_; i++) { |
---|
3772 | double lowerValue = columnLower_[i]; |
---|
3773 | double upperValue = columnUpper_[i]; |
---|
3774 | if (lowerValue > -1.0e20) { |
---|
3775 | columnLowerWork_[i] = lowerValue * rhsScale_; |
---|
3776 | if (upperValue >= 1.0e20) { |
---|
3777 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
3778 | } else { |
---|
3779 | columnUpperWork_[i] = upperValue * rhsScale_; |
---|
3780 | if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) { |
---|
3781 | if (columnLowerWork_[i] >= 0.0) { |
---|
3782 | columnUpperWork_[i] = columnLowerWork_[i]; |
---|
3783 | } else if (columnUpperWork_[i] <= 0.0) { |
---|
3784 | columnLowerWork_[i] = columnUpperWork_[i]; |
---|
3785 | } else { |
---|
3786 | columnUpperWork_[i] = 0.0; |
---|
3787 | columnLowerWork_[i] = 0.0; |
---|
3788 | } |
---|
3789 | } |
---|
3790 | } |
---|
3791 | } else if (upperValue < 1.0e20) { |
---|
3792 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
3793 | columnUpperWork_[i] = upperValue * rhsScale_; |
---|
3794 | } else { |
---|
3795 | // free |
---|
3796 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
3797 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
3798 | } |
---|
3799 | } |
---|
3800 | for (i = 0; i < numberRows_; i++) { |
---|
3801 | double lowerValue = rowLower_[i]; |
---|
3802 | double upperValue = rowUpper_[i]; |
---|
3803 | if (lowerValue > -1.0e20) { |
---|
3804 | rowLowerWork_[i] = lowerValue * rhsScale_; |
---|
3805 | if (upperValue >= 1.0e20) { |
---|
3806 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
3807 | } else { |
---|
3808 | rowUpperWork_[i] = upperValue * rhsScale_; |
---|
3809 | if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) { |
---|
3810 | if (rowLowerWork_[i] >= 0.0) { |
---|
3811 | rowUpperWork_[i] = rowLowerWork_[i]; |
---|
3812 | } else if (rowUpperWork_[i] <= 0.0) { |
---|
3813 | rowLowerWork_[i] = rowUpperWork_[i]; |
---|
3814 | } else { |
---|
3815 | rowUpperWork_[i] = 0.0; |
---|
3816 | rowLowerWork_[i] = 0.0; |
---|
3817 | } |
---|
3818 | } |
---|
3819 | } |
---|
3820 | } else if (upperValue < 1.0e20) { |
---|
3821 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
3822 | rowUpperWork_[i] = upperValue * rhsScale_; |
---|
3823 | } else { |
---|
3824 | // free |
---|
3825 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
3826 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
3827 | } |
---|
3828 | } |
---|
3829 | } else { |
---|
3830 | for (i = 0; i < numberColumns_; i++) { |
---|
3831 | double lowerValue = columnLower_[i]; |
---|
3832 | double upperValue = columnUpper_[i]; |
---|
3833 | if (lowerValue > -1.0e20) { |
---|
3834 | columnLowerWork_[i] = lowerValue; |
---|
3835 | if (upperValue >= 1.0e20) { |
---|
3836 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
3837 | } else { |
---|
3838 | columnUpperWork_[i] = upperValue; |
---|
3839 | if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) { |
---|
3840 | if (columnLowerWork_[i] >= 0.0) { |
---|
3841 | columnUpperWork_[i] = columnLowerWork_[i]; |
---|
3842 | } else if (columnUpperWork_[i] <= 0.0) { |
---|
3843 | columnLowerWork_[i] = columnUpperWork_[i]; |
---|
3844 | } else { |
---|
3845 | columnUpperWork_[i] = 0.0; |
---|
3846 | columnLowerWork_[i] = 0.0; |
---|
3847 | } |
---|
3848 | } |
---|
3849 | } |
---|
3850 | } else if (upperValue < 1.0e20) { |
---|
3851 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
3852 | columnUpperWork_[i] = upperValue; |
---|
3853 | } else { |
---|
3854 | // free |
---|
3855 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
3856 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
3857 | } |
---|
3858 | } |
---|
3859 | for (i = 0; i < numberRows_; i++) { |
---|
3860 | double lowerValue = rowLower_[i]; |
---|
3861 | double upperValue = rowUpper_[i]; |
---|
3862 | if (lowerValue > -1.0e20) { |
---|
3863 | rowLowerWork_[i] = lowerValue; |
---|
3864 | if (upperValue >= 1.0e20) { |
---|
3865 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
3866 | } else { |
---|
3867 | rowUpperWork_[i] = upperValue; |
---|
3868 | if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) { |
---|
3869 | if (rowLowerWork_[i] >= 0.0) { |
---|
3870 | rowUpperWork_[i] = rowLowerWork_[i]; |
---|
3871 | } else if (rowUpperWork_[i] <= 0.0) { |
---|
3872 | rowLowerWork_[i] = rowUpperWork_[i]; |
---|
3873 | } else { |
---|
3874 | rowUpperWork_[i] = 0.0; |
---|
3875 | rowLowerWork_[i] = 0.0; |
---|
3876 | } |
---|
3877 | } |
---|
3878 | } |
---|
3879 | } else if (upperValue < 1.0e20) { |
---|
3880 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
3881 | rowUpperWork_[i] = upperValue; |
---|
3882 | } else { |
---|
3883 | // free |
---|
3884 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
3885 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
3886 | } |
---|
3887 | } |
---|
3888 | } |
---|
3889 | } |
---|
3890 | if (what == 63) { |
---|
3891 | // move information to work arrays |
---|
3892 | double direction = optimizationDirection_; |
---|
3893 | // direction is actually scale out not scale in |
---|
3894 | if (direction) |
---|
3895 | direction = 1.0 / direction; |
---|
3896 | if (direction != 1.0) { |
---|
3897 | // reverse all dual signs |
---|
3898 | for (i = 0; i < numberColumns_; i++) |
---|
3899 | reducedCost_[i] *= direction; |
---|
3900 | for (i = 0; i < numberRows_; i++) |
---|
3901 | dual_[i] *= direction; |
---|
3902 | } |
---|
3903 | for (i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
3904 | setFakeBound(i, noFake); |
---|
3905 | } |
---|
3906 | if (rowScale_) { |
---|
3907 | const double *obj = objective(); |
---|
3908 | double direction = optimizationDirection_ * objectiveScale_; |
---|
3909 | // clean up any mismatches on infinity |
---|
3910 | // and fix any variables with tiny gaps |
---|
3911 | double primalTolerance = dblParam_[ClpPrimalTolerance]; |
---|
3912 | // on entry |
---|
3913 | const double *inverseScale = inverseColumnScale_; |
---|
3914 | for (i = 0; i < numberColumns_; i++) { |
---|
3915 | CoinAssert(fabs(obj[i]) < 1.0e25); |
---|
3916 | double scaleFactor = columnScale_[i]; |
---|
3917 | double multiplier = rhsScale_ * inverseScale[i]; |
---|
3918 | scaleFactor *= direction; |
---|
3919 | objectiveWork_[i] = obj[i] * scaleFactor; |
---|
3920 | reducedCostWork_[i] = reducedCost_[i] * scaleFactor; |
---|
3921 | double lowerValue = columnLower_[i]; |
---|
3922 | double upperValue = columnUpper_[i]; |
---|
3923 | if (lowerValue > -1.0e20) { |
---|
3924 | columnLowerWork_[i] = lowerValue * multiplier; |
---|
3925 | if (upperValue >= 1.0e20) { |
---|
3926 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
3927 | } else { |
---|
3928 | columnUpperWork_[i] = upperValue * multiplier; |
---|
3929 | if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) { |
---|
3930 | if (columnLowerWork_[i] >= 0.0) { |
---|
3931 | columnUpperWork_[i] = columnLowerWork_[i]; |
---|
3932 | } else if (columnUpperWork_[i] <= 0.0) { |
---|
3933 | columnLowerWork_[i] = columnUpperWork_[i]; |
---|
3934 | } else { |
---|
3935 | columnUpperWork_[i] = 0.0; |
---|
3936 | columnLowerWork_[i] = 0.0; |
---|
3937 | } |
---|
3938 | } |
---|
3939 | } |
---|
3940 | } else if (upperValue < 1.0e20) { |
---|
3941 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
3942 | columnUpperWork_[i] = upperValue * multiplier; |
---|
3943 | } else { |
---|
3944 | // free |
---|
3945 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
3946 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
3947 | } |
---|
3948 | double value = columnActivity_[i] * multiplier; |
---|
3949 | if (fabs(value) > 1.0e20) { |
---|
3950 | //printf("bad value of %g for column %d\n",value,i); |
---|
3951 | setColumnStatus(i, superBasic); |
---|
3952 | if (columnUpperWork_[i] < 0.0) { |
---|
3953 | value = columnUpperWork_[i]; |
---|
3954 | } else if (columnLowerWork_[i] > 0.0) { |
---|
3955 | value = columnLowerWork_[i]; |
---|
3956 | } else { |
---|
3957 | value = 0.0; |
---|
3958 | } |
---|
3959 | } |
---|
3960 | columnActivityWork_[i] = value; |
---|
3961 | } |
---|
3962 | inverseScale = inverseRowScale_; |
---|
3963 | for (i = 0; i < numberRows_; i++) { |
---|
3964 | dual_[i] *= inverseScale[i]; |
---|
3965 | dual_[i] *= objectiveScale_; |
---|
3966 | rowReducedCost_[i] = dual_[i]; |
---|
3967 | double multiplier = rhsScale_ * rowScale_[i]; |
---|
3968 | double value = rowActivity_[i] * multiplier; |
---|
3969 | if (fabs(value) > 1.0e20) { |
---|
3970 | //printf("bad value of %g for row %d\n",value,i); |
---|
3971 | setRowStatus(i, superBasic); |
---|
3972 | if (rowUpperWork_[i] < 0.0) { |
---|
3973 | value = rowUpperWork_[i]; |
---|
3974 | } else if (rowLowerWork_[i] > 0.0) { |
---|
3975 | value = rowLowerWork_[i]; |
---|
3976 | } else { |
---|
3977 | value = 0.0; |
---|
3978 | } |
---|
3979 | } |
---|
3980 | rowActivityWork_[i] = value; |
---|
3981 | } |
---|
3982 | } else if (objectiveScale_ != 1.0 || rhsScale_ != 1.0) { |
---|
3983 | // on entry |
---|
3984 | for (i = 0; i < numberColumns_; i++) { |
---|
3985 | double value = columnActivity_[i]; |
---|
3986 | value *= rhsScale_; |
---|
3987 | if (fabs(value) > 1.0e20) { |
---|
3988 | //printf("bad value of %g for column %d\n",value,i); |
---|
3989 | setColumnStatus(i, superBasic); |
---|
3990 | if (columnUpperWork_[i] < 0.0) { |
---|
3991 | value = columnUpperWork_[i]; |
---|
3992 | } else if (columnLowerWork_[i] > 0.0) { |
---|
3993 | value = columnLowerWork_[i]; |
---|
3994 | } else { |
---|
3995 | value = 0.0; |
---|
3996 | } |
---|
3997 | } |
---|
3998 | columnActivityWork_[i] = value; |
---|
3999 | reducedCostWork_[i] = reducedCost_[i] * objectiveScale_; |
---|
4000 | } |
---|
4001 | for (i = 0; i < numberRows_; i++) { |
---|
4002 | double value = rowActivity_[i]; |
---|
4003 | value *= rhsScale_; |
---|
4004 | if (fabs(value) > 1.0e20) { |
---|
4005 | //printf("bad value of %g for row %d\n",value,i); |
---|
4006 | setRowStatus(i, superBasic); |
---|
4007 | if (rowUpperWork_[i] < 0.0) { |
---|
4008 | value = rowUpperWork_[i]; |
---|
4009 | } else if (rowLowerWork_[i] > 0.0) { |
---|
4010 | value = rowLowerWork_[i]; |
---|
4011 | } else { |
---|
4012 | value = 0.0; |
---|
4013 | } |
---|
4014 | } |
---|
4015 | rowActivityWork_[i] = value; |
---|
4016 | dual_[i] *= objectiveScale_; |
---|
4017 | rowReducedCost_[i] = dual_[i]; |
---|
4018 | } |
---|
4019 | } else { |
---|
4020 | // on entry |
---|
4021 | for (i = 0; i < numberColumns_; i++) { |
---|
4022 | double value = columnActivity_[i]; |
---|
4023 | if (fabs(value) > 1.0e20) { |
---|
4024 | //printf("bad value of %g for column %d\n",value,i); |
---|
4025 | setColumnStatus(i, superBasic); |
---|
4026 | if (columnUpperWork_[i] < 0.0) { |
---|
4027 | value = columnUpperWork_[i]; |
---|
4028 | } else if (columnLowerWork_[i] > 0.0) { |
---|
4029 | value = columnLowerWork_[i]; |
---|
4030 | } else { |
---|
4031 | value = 0.0; |
---|
4032 | } |
---|
4033 | } |
---|
4034 | columnActivityWork_[i] = value; |
---|
4035 | reducedCostWork_[i] = reducedCost_[i]; |
---|
4036 | } |
---|
4037 | for (i = 0; i < numberRows_; i++) { |
---|
4038 | double value = rowActivity_[i]; |
---|
4039 | if (fabs(value) > 1.0e20) { |
---|
4040 | //printf("bad value of %g for row %d\n",value,i); |
---|
4041 | setRowStatus(i, superBasic); |
---|
4042 | if (rowUpperWork_[i] < 0.0) { |
---|
4043 | value = rowUpperWork_[i]; |
---|
4044 | } else if (rowLowerWork_[i] > 0.0) { |
---|
4045 | value = rowLowerWork_[i]; |
---|
4046 | } else { |
---|
4047 | value = 0.0; |
---|
4048 | } |
---|
4049 | } |
---|
4050 | rowActivityWork_[i] = value; |
---|
4051 | rowReducedCost_[i] = dual_[i]; |
---|
4052 | } |
---|
4053 | } |
---|
4054 | } |
---|
4055 | |
---|
4056 | if (what == 63 && doSanityCheck) { |
---|
4057 | // check rim of problem okay |
---|
4058 | if (!sanityCheck()) |
---|
4059 | goodMatrix = false; |
---|
4060 | } |
---|
4061 | // we need to treat matrix as if each element by rowScaleIn and columnScaleout?? |
---|
4062 | // maybe we need to move scales to SimplexModel for factorization? |
---|
4063 | if ((what == 63 && !pivotVariable_) || (newArrays && !keepPivots)) { |
---|
4064 | delete[] pivotVariable_; |
---|
4065 | pivotVariable_ = new int[numberRows2 + 1]; |
---|
4066 | for (int i = 0; i < numberRows2 + 1; i++) |
---|
4067 | pivotVariable_[i] = -1; |
---|
4068 | } else if (what == 63 && !keepPivots) { |
---|
4069 | // just reset |
---|
4070 | for (int i = 0; i < numberRows2; i++) |
---|
4071 | pivotVariable_[i] = -1; |
---|
4072 | } else if (what == 63) { |
---|
4073 | // check pivots |
---|
4074 | for (int i = 0; i < numberRows2; i++) { |
---|
4075 | int iSequence = pivotVariable_[i]; |
---|
4076 | if (iSequence < 0 || (iSequence < numberRows_ + numberColumns_ && getStatus(iSequence) != basic)) { |
---|
4077 | keepPivots = false; |
---|
4078 | break; |
---|
4079 | } |
---|
4080 | } |
---|
4081 | if (!keepPivots) { |
---|
4082 | // reset |
---|
4083 | for (int i = 0; i < numberRows2; i++) |
---|
4084 | pivotVariable_[i] = -1; |
---|
4085 | } else { |
---|
4086 | // clean |
---|
4087 | for (int i = 0; i < numberColumns_ + numberRows_; i++) { |
---|
4088 | Status status = getStatus(i); |
---|
4089 | if (status != basic) { |
---|
4090 | if (upper_[i] == lower_[i]) { |
---|
4091 | setStatus(i, isFixed); |
---|
4092 | solution_[i] = lower_[i]; |
---|
4093 | } else if (status == atLowerBound) { |
---|
4094 | if (lower_[i] > -1.0e20) { |
---|
4095 | solution_[i] = lower_[i]; |
---|
4096 | } else { |
---|
4097 | //printf("seq %d at lower of %g\n",i,lower_[i]); |
---|
4098 | if (upper_[i] < 1.0e20) { |
---|
4099 | solution_[i] = upper_[i]; |
---|
4100 | setStatus(i, atUpperBound); |
---|
4101 | } else { |
---|
4102 | setStatus(i, isFree); |
---|
4103 | } |
---|
4104 | } |
---|
4105 | } else if (status == atUpperBound) { |
---|
4106 | if (upper_[i] < 1.0e20) { |
---|
4107 | solution_[i] = upper_[i]; |
---|
4108 | } else { |
---|
4109 | //printf("seq %d at upper of %g\n",i,upper_[i]); |
---|
4110 | if (lower_[i] > -1.0e20) { |
---|
4111 | solution_[i] = lower_[i]; |
---|
4112 | setStatus(i, atLowerBound); |
---|
4113 | } else { |
---|
4114 | setStatus(i, isFree); |
---|
4115 | } |
---|
4116 | } |
---|
4117 | } else if (status == isFixed && upper_[i] > lower_[i]) { |
---|
4118 | // was fixed - not now |
---|
4119 | if (solution_[i] <= lower_[i] + primalTolerance_) { |
---|
4120 | setStatus(i, atLowerBound); |
---|
4121 | } else if (solution_[i] >= upper_[i] - primalTolerance_) { |
---|
4122 | setStatus(i, atUpperBound); |
---|
4123 | } else { |
---|
4124 | setStatus(i, superBasic); |
---|
4125 | } |
---|
4126 | } |
---|
4127 | } |
---|
4128 | } |
---|
4129 | } |
---|
4130 | } |
---|
4131 | |
---|
4132 | if (what == 63) { |
---|
4133 | // Safer to get new arrays anyway so test on newArrays removed |
---|
4134 | // get some arrays |
---|
4135 | int iRow, iColumn; |
---|
4136 | // these are "indexed" arrays so we always know where nonzeros are |
---|
4137 | /********************************************************** |
---|
4138 | rowArray_[3] is long enough for rows+columns (2 also maybe) |
---|
4139 | rowArray_[1] is long enough for max(rows,columns) |
---|
4140 | *********************************************************/ |
---|
4141 | for (iRow = 0; iRow < 4; iRow++) { |
---|
4142 | int length = numberRows2 + factorization_->maximumPivots(); |
---|
4143 | if (iRow > SHORT_REGION || objective_->type() > 1) |
---|
4144 | length += numberColumns_; |
---|
4145 | else if (iRow == 1) |
---|
4146 | length = CoinMax(length, numberColumns_); |
---|
4147 | if ((specialOptions_ & 65536) == 0 || !rowArray_[iRow]) { |
---|
4148 | delete rowArray_[iRow]; |
---|
4149 | rowArray_[iRow] = new CoinIndexedVector(); |
---|
4150 | } |
---|
4151 | rowArray_[iRow]->reserve(length); |
---|
4152 | } |
---|
4153 | |
---|
4154 | for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) { |
---|
4155 | if ((specialOptions_ & 65536) == 0 || !columnArray_[iColumn]) { |
---|
4156 | delete columnArray_[iColumn]; |
---|
4157 | columnArray_[iColumn] = new CoinIndexedVector(); |
---|
4158 | } |
---|
4159 | columnArray_[iColumn]->reserve(numberColumns_ + numberRows2); |
---|
4160 | } |
---|
4161 | } |
---|
4162 | if (problemStatus_ == 10) { |
---|
4163 | problemStatus_ = -1; |
---|
4164 | handler_->setLogLevel(saveLevel); // switch back messages |
---|
4165 | } |
---|
4166 | if ((what & 5) != 0) |
---|
4167 | matrix_->generalExpanded(this, 9, what); // update costs and bounds if necessary |
---|
4168 | if (goodMatrix && (specialOptions_ & 65536) != 0) { |
---|
4169 | int save = maximumColumns_ + maximumRows_; |
---|
4170 | CoinMemcpyN(cost_, numberTotal, cost_ + save); |
---|
4171 | CoinMemcpyN(lower_, numberTotal, lower_ + save); |
---|
4172 | CoinMemcpyN(upper_, numberTotal, upper_ + save); |
---|
4173 | CoinMemcpyN(dj_, numberTotal, dj_ + save); |
---|
4174 | CoinMemcpyN(solution_, numberTotal, solution_ + save); |
---|
4175 | if (rowScale_ && !savedRowScale_) { |
---|
4176 | double *temp; |
---|
4177 | temp = new double[4 * maximumRows_]; |
---|
4178 | CoinFillN(temp, 4 * maximumRows_, 1.0); |
---|
4179 | CoinMemcpyN(rowScale_, numberRows2, temp); |
---|
4180 | CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + maximumRows_); |
---|
4181 | CoinMemcpyN(rowScale_, numberRows2, temp + 2 * maximumRows_); |
---|
4182 | CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + 3 * maximumRows_); |
---|
4183 | delete[] rowScale_; |
---|
4184 | savedRowScale_ = temp; |
---|
4185 | rowScale_ = savedRowScale_; |
---|
4186 | inverseRowScale_ = savedRowScale_ + maximumInternalRows_; |
---|
4187 | temp = new double[4 * maximumColumns_]; |
---|
4188 | CoinFillN(temp, 4 * maximumColumns_, 1.0); |
---|
4189 | CoinMemcpyN(columnScale_, numberColumns_, temp); |
---|
4190 | CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + maximumColumns_); |
---|
4191 | CoinMemcpyN(columnScale_, numberColumns_, temp + 2 * maximumColumns_); |
---|
4192 | CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + 3 * maximumColumns_); |
---|
4193 | delete[] columnScale_; |
---|
4194 | savedColumnScale_ = temp; |
---|
4195 | columnScale_ = savedColumnScale_; |
---|
4196 | inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_; |
---|
4197 | } |
---|
4198 | } |
---|
4199 | #ifdef CLP_USER_DRIVEN |
---|
4200 | eventHandler_->event(ClpEventHandler::endOfCreateRim); |
---|
4201 | #endif |
---|
4202 | return goodMatrix; |
---|
4203 | } |
---|
4204 | // Does rows and columns |
---|
4205 | void ClpSimplex::createRim1(bool initial) |
---|
4206 | { |
---|
4207 | int i; |
---|
4208 | int numberRows2 = numberRows_ + numberExtraRows_; |
---|
4209 | int numberTotal = numberRows2 + numberColumns_; |
---|
4210 | if ((specialOptions_ & 65536) != 0 && true) { |
---|
4211 | assert(!initial); |
---|
4212 | int save = maximumColumns_ + maximumRows_; |
---|
4213 | CoinMemcpyN(lower_ + save, numberTotal, lower_); |
---|
4214 | CoinMemcpyN(upper_ + save, numberTotal, upper_); |
---|
4215 | return; |
---|
4216 | } |
---|
4217 | const double *rowScale = rowScale_; |
---|
4218 | // clean up any mismatches on infinity |
---|
4219 | // and fix any variables with tiny gaps |
---|
4220 | double primalTolerance = dblParam_[ClpPrimalTolerance]; |
---|
4221 | if (rowScale) { |
---|
4222 | // If scaled then do all columns later in one loop |
---|
4223 | if (!initial) { |
---|
4224 | const double *inverseScale = inverseColumnScale_; |
---|
4225 | for (i = 0; i < numberColumns_; i++) { |
---|
4226 | double multiplier = rhsScale_ * inverseScale[i]; |
---|
4227 | double lowerValue = columnLower_[i]; |
---|
4228 | double upperValue = columnUpper_[i]; |
---|
4229 | if (lowerValue > -1.0e20) { |
---|
4230 | columnLowerWork_[i] = lowerValue * multiplier; |
---|
4231 | if (upperValue >= 1.0e20) { |
---|
4232 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
4233 | } else { |
---|
4234 | columnUpperWork_[i] = upperValue * multiplier; |
---|
4235 | if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) { |
---|
4236 | if (columnLowerWork_[i] >= 0.0) { |
---|
4237 | columnUpperWork_[i] = columnLowerWork_[i]; |
---|
4238 | } else if (columnUpperWork_[i] <= 0.0) { |
---|
4239 | columnLowerWork_[i] = columnUpperWork_[i]; |
---|
4240 | } else { |
---|
4241 | columnUpperWork_[i] = 0.0; |
---|
4242 | columnLowerWork_[i] = 0.0; |
---|
4243 | } |
---|
4244 | } |
---|
4245 | } |
---|
4246 | } else if (upperValue < 1.0e20) { |
---|
4247 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
4248 | columnUpperWork_[i] = upperValue * multiplier; |
---|
4249 | } else { |
---|
4250 | // free |
---|
4251 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
4252 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
4253 | } |
---|
4254 | } |
---|
4255 | } |
---|
4256 | for (i = 0; i < numberRows_; i++) { |
---|
4257 | double multiplier = rhsScale_ * rowScale[i]; |
---|
4258 | double lowerValue = rowLower_[i]; |
---|
4259 | double upperValue = rowUpper_[i]; |
---|
4260 | if (lowerValue > -1.0e20) { |
---|
4261 | rowLowerWork_[i] = lowerValue * multiplier; |
---|
4262 | if (upperValue >= 1.0e20) { |
---|
4263 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
4264 | } else { |
---|
4265 | rowUpperWork_[i] = upperValue * multiplier; |
---|
4266 | if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) { |
---|
4267 | if (rowLowerWork_[i] >= 0.0) { |
---|
4268 | rowUpperWork_[i] = rowLowerWork_[i]; |
---|
4269 | } else if (rowUpperWork_[i] <= 0.0) { |
---|
4270 | rowLowerWork_[i] = rowUpperWork_[i]; |
---|
4271 | } else { |
---|
4272 | rowUpperWork_[i] = 0.0; |
---|
4273 | rowLowerWork_[i] = 0.0; |
---|
4274 | } |
---|
4275 | } |
---|
4276 | } |
---|
4277 | } else if (upperValue < 1.0e20) { |
---|
4278 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
4279 | rowUpperWork_[i] = upperValue * multiplier; |
---|
4280 | } else { |
---|
4281 | // free |
---|
4282 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
4283 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
4284 | } |
---|
4285 | } |
---|
4286 | } else if (rhsScale_ != 1.0) { |
---|
4287 | for (i = 0; i < numberColumns_; i++) { |
---|
4288 | double lowerValue = columnLower_[i]; |
---|
4289 | double upperValue = columnUpper_[i]; |
---|
4290 | if (lowerValue > -1.0e20) { |
---|
4291 | columnLowerWork_[i] = lowerValue * rhsScale_; |
---|
4292 | if (upperValue >= 1.0e20) { |
---|
4293 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
4294 | } else { |
---|
4295 | columnUpperWork_[i] = upperValue * rhsScale_; |
---|
4296 | if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) { |
---|
4297 | if (columnLowerWork_[i] >= 0.0) { |
---|
4298 | columnUpperWork_[i] = columnLowerWork_[i]; |
---|
4299 | } else if (columnUpperWork_[i] <= 0.0) { |
---|
4300 | columnLowerWork_[i] = columnUpperWork_[i]; |
---|
4301 | } else { |
---|
4302 | columnUpperWork_[i] = 0.0; |
---|
4303 | columnLowerWork_[i] = 0.0; |
---|
4304 | } |
---|
4305 | } |
---|
4306 | } |
---|
4307 | } else if (upperValue < 1.0e20) { |
---|
4308 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
4309 | columnUpperWork_[i] = upperValue * rhsScale_; |
---|
4310 | } else { |
---|
4311 | // free |
---|
4312 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
4313 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
4314 | } |
---|
4315 | } |
---|
4316 | for (i = 0; i < numberRows_; i++) { |
---|
4317 | double lowerValue = rowLower_[i]; |
---|
4318 | double upperValue = rowUpper_[i]; |
---|
4319 | if (lowerValue > -1.0e20) { |
---|
4320 | rowLowerWork_[i] = lowerValue * rhsScale_; |
---|
4321 | if (upperValue >= 1.0e20) { |
---|
4322 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
4323 | } else { |
---|
4324 | rowUpperWork_[i] = upperValue * rhsScale_; |
---|
4325 | if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) { |
---|
4326 | if (rowLowerWork_[i] >= 0.0) { |
---|
4327 | rowUpperWork_[i] = rowLowerWork_[i]; |
---|
4328 | } else if (rowUpperWork_[i] <= 0.0) { |
---|
4329 | rowLowerWork_[i] = rowUpperWork_[i]; |
---|
4330 | } else { |
---|
4331 | rowUpperWork_[i] = 0.0; |
---|
4332 | rowLowerWork_[i] = 0.0; |
---|
4333 | } |
---|
4334 | } |
---|
4335 | } |
---|
4336 | } else if (upperValue < 1.0e20) { |
---|
4337 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
4338 | rowUpperWork_[i] = upperValue * rhsScale_; |
---|
4339 | } else { |
---|
4340 | // free |
---|
4341 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
4342 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
4343 | } |
---|
4344 | } |
---|
4345 | } else { |
---|
4346 | for (i = 0; i < numberColumns_; i++) { |
---|
4347 | double lowerValue = columnLower_[i]; |
---|
4348 | double upperValue = columnUpper_[i]; |
---|
4349 | if (lowerValue > -1.0e20) { |
---|
4350 | columnLowerWork_[i] = lowerValue; |
---|
4351 | if (upperValue >= 1.0e20) { |
---|
4352 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
4353 | } else { |
---|
4354 | columnUpperWork_[i] = upperValue; |
---|
4355 | if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) { |
---|
4356 | if (columnLowerWork_[i] >= 0.0) { |
---|
4357 | columnUpperWork_[i] = columnLowerWork_[i]; |
---|
4358 | } else if (columnUpperWork_[i] <= 0.0) { |
---|
4359 | columnLowerWork_[i] = columnUpperWork_[i]; |
---|
4360 | } else { |
---|
4361 | columnUpperWork_[i] = 0.0; |
---|
4362 | columnLowerWork_[i] = 0.0; |
---|
4363 | } |
---|
4364 | } |
---|
4365 | } |
---|
4366 | } else if (upperValue < 1.0e20) { |
---|
4367 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
4368 | columnUpperWork_[i] = upperValue; |
---|
4369 | } else { |
---|
4370 | // free |
---|
4371 | columnLowerWork_[i] = -COIN_DBL_MAX; |
---|
4372 | columnUpperWork_[i] = COIN_DBL_MAX; |
---|
4373 | } |
---|
4374 | } |
---|
4375 | for (i = 0; i < numberRows_; i++) { |
---|
4376 | double lowerValue = rowLower_[i]; |
---|
4377 | double upperValue = rowUpper_[i]; |
---|
4378 | if (lowerValue > -1.0e20) { |
---|
4379 | rowLowerWork_[i] = lowerValue; |
---|
4380 | if (upperValue >= 1.0e20) { |
---|
4381 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
4382 | } else { |
---|
4383 | rowUpperWork_[i] = upperValue; |
---|
4384 | if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) { |
---|
4385 | if (rowLowerWork_[i] >= 0.0) { |
---|
4386 | rowUpperWork_[i] = rowLowerWork_[i]; |
---|
4387 | } else if (rowUpperWork_[i] <= 0.0) { |
---|
4388 | rowLowerWork_[i] = rowUpperWork_[i]; |
---|
4389 | } else { |
---|
4390 | rowUpperWork_[i] = 0.0; |
---|
4391 | rowLowerWork_[i] = 0.0; |
---|
4392 | } |
---|
4393 | } |
---|
4394 | } |
---|
4395 | } else if (upperValue < 1.0e20) { |
---|
4396 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
4397 | rowUpperWork_[i] = upperValue; |
---|
4398 | } else { |
---|
4399 | // free |
---|
4400 | rowLowerWork_[i] = -COIN_DBL_MAX; |
---|
4401 | rowUpperWork_[i] = COIN_DBL_MAX; |
---|
4402 | } |
---|
4403 | } |
---|
4404 | } |
---|
4405 | #ifndef NDEBUG |
---|
4406 | if ((specialOptions_ & 65536) != 0 && false) { |
---|
4407 | assert(!initial); |
---|
4408 | int save = maximumColumns_ + maximumRows_; |
---|
4409 | for (int i = 0; i < numberTotal; i++) { |
---|
4410 | assert(fabs(lower_[i] - lower_[i + save]) < 1.0e-5); |
---|
4411 | assert(fabs(upper_[i] - upper_[i + save]) < 1.0e-5); |
---|
4412 | } |
---|
4413 | } |
---|
4414 | #endif |
---|
4415 | } |
---|
4416 | // Does objective |
---|
4417 | void ClpSimplex::createRim4(bool initial) |
---|
4418 | { |
---|
4419 | int i; |
---|
4420 | int numberRows2 = numberRows_ + numberExtraRows_; |
---|
4421 | int numberTotal = numberRows2 + numberColumns_; |
---|
4422 | if ((specialOptions_ & 65536) != 0 && true) { |
---|
4423 | assert(!initial); |
---|
4424 | int save = maximumColumns_ + maximumRows_; |
---|
4425 | CoinMemcpyN(cost_ + save, numberTotal, cost_); |
---|
4426 | return; |
---|
4427 | } |
---|
4428 | double direction = optimizationDirection_ * objectiveScale_; |
---|
4429 | const double *obj = objective(); |
---|
4430 | const double *rowScale = rowScale_; |
---|
4431 | const double *columnScale = columnScale_; |
---|
4432 | // and also scale by scale factors |
---|
4433 | if (rowScale) { |
---|
4434 | if (rowObjective_) { |
---|
4435 | for (i = 0; i < numberRows_; i++) |
---|
4436 | rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i]; |
---|
4437 | } else { |
---|
4438 | memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double)); |
---|
4439 | } |
---|
4440 | // If scaled then do all columns later in one loop |
---|
4441 | if (!initial) { |
---|
4442 | for (i = 0; i < numberColumns_; i++) { |
---|
4443 | CoinAssert(fabs(obj[i]) < 1.0e25); |
---|
4444 | objectiveWork_[i] = obj[i] * direction * columnScale[i]; |
---|
4445 | } |
---|
4446 | } |
---|
4447 | } else { |
---|
4448 | if (rowObjective_) { |
---|
4449 | for (i = 0; i < numberRows_; i++) |
---|
4450 | rowObjectiveWork_[i] = rowObjective_[i] * direction; |
---|
4451 | } else { |
---|
4452 | memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double)); |
---|
4453 | } |
---|
4454 | for (i = 0; i < numberColumns_; i++) { |
---|
4455 | CoinAssert(fabs(obj[i]) < 1.0e25); |
---|
4456 | objectiveWork_[i] = obj[i] * direction; |
---|
4457 | } |
---|
4458 | } |
---|
4459 | } |
---|
4460 | // Does rows and columns and objective |
---|
4461 | void ClpSimplex::createRim5(bool initial) |
---|
4462 | { |
---|
4463 | createRim4(initial); |
---|
4464 | createRim1(initial); |
---|
4465 | } |
---|
4466 | void ClpSimplex::deleteRim(int getRidOfFactorizationData) |
---|
4467 | { |
---|
4468 | #ifdef CLP_USER_DRIVEN |
---|
4469 | eventHandler_->event(ClpEventHandler::beforeDeleteRim); |
---|
4470 | #endif |
---|
4471 | // Just possible empty problem |
---|
4472 | int numberRows = numberRows_; |
---|
4473 | int numberColumns = numberColumns_; |
---|
4474 | if (!numberRows || !numberColumns) { |
---|
4475 | numberRows = 0; |
---|
4476 | if (objective_->type() < 2) |
---|
4477 | numberColumns = 0; |
---|
4478 | } |
---|
4479 | int i; |
---|
4480 | if (problemStatus_ != 1 && problemStatus_ != 2) { |
---|
4481 | delete[] ray_; |
---|
4482 | ray_ = NULL; |
---|
4483 | } |
---|
4484 | // set upperOut_ to furthest away from bound so can use in dual for dualBound_ |
---|
4485 | upperOut_ = 1.0; |
---|
4486 | #if 0 |
---|
4487 | { |
---|
4488 | int nBad = 0; |
---|
4489 | for (i = 0; i < numberColumns; i++) { |
---|
4490 | if (lower_[i] == upper_[i] && getColumnStatus(i) == basic) |
---|
4491 | nBad++; |
---|
4492 | } |
---|
4493 | if (nBad) |
---|
4494 | printf("yy %d basic fixed\n", nBad); |
---|
4495 | } |
---|
4496 | #endif |
---|
4497 | if ((moreSpecialOptions_ & 4194304) != 0) { |
---|
4498 | // preset tolerances were changed |
---|
4499 | moreSpecialOptions_ &= ~4194304; |
---|
4500 | primalTolerance_ = 1.0e-7; |
---|
4501 | dblParam_[ClpPrimalTolerance] = primalTolerance_; |
---|
4502 | dualTolerance_ = 1.0e-7; |
---|
4503 | dblParam_[ClpDualTolerance] = dualTolerance_; |
---|
4504 | } |
---|
4505 | // ray may be null if in branch and bound |
---|
4506 | if (rowScale_ && solution_) { |
---|
4507 | // Collect infeasibilities |
---|
4508 | int numberPrimalScaled = 0; |
---|
4509 | int numberPrimalUnscaled = 0; |
---|
4510 | int numberDualScaled = 0; |
---|
4511 | int numberDualUnscaled = 0; |
---|
4512 | double scaleC = 1.0 / objectiveScale_; |
---|
4513 | double scaleR = 1.0 / rhsScale_; |
---|
4514 | const double *inverseScale = inverseColumnScale_; |
---|
4515 | for (i = 0; i < numberColumns; i++) { |
---|
4516 | double scaleFactor = columnScale_[i]; |
---|
4517 | double valueScaled = columnActivityWork_[i]; |
---|
4518 | double lowerScaled = columnLowerWork_[i]; |
---|
4519 | double upperScaled = columnUpperWork_[i]; |
---|
4520 | if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) { |
---|
4521 | if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_) |
---|
4522 | numberPrimalScaled++; |
---|
4523 | else |
---|
4524 | upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled)); |
---|
4525 | } |
---|
4526 | columnActivity_[i] = valueScaled * scaleFactor * scaleR; |
---|
4527 | double value = columnActivity_[i]; |
---|
4528 | if (value < columnLower_[i] - primalTolerance_) |
---|
4529 | numberPrimalUnscaled++; |
---|
4530 | else if (value > columnUpper_[i] + primalTolerance_) |
---|
4531 | numberPrimalUnscaled++; |
---|
4532 | double valueScaledDual = reducedCostWork_[i]; |
---|
4533 | if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_) |
---|
4534 | numberDualScaled++; |
---|
4535 | if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_) |
---|
4536 | numberDualScaled++; |
---|
4537 | reducedCost_[i] = (valueScaledDual * scaleC) * inverseScale[i]; |
---|
4538 | double valueDual = reducedCost_[i]; |
---|
4539 | if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_) |
---|
4540 | numberDualUnscaled++; |
---|
4541 | if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_) |
---|
4542 | numberDualUnscaled++; |
---|
4543 | } |
---|
4544 | inverseScale = inverseRowScale_; |
---|
4545 | for (i = 0; i < numberRows; i++) { |
---|
4546 | double scaleFactor = rowScale_[i]; |
---|
4547 | double valueScaled = rowActivityWork_[i]; |
---|
4548 | double lowerScaled = rowLowerWork_[i]; |
---|
4549 | double upperScaled = rowUpperWork_[i]; |
---|
4550 | if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) { |
---|
4551 | if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_) |
---|
4552 | numberPrimalScaled++; |
---|
4553 | else |
---|
4554 | upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled)); |
---|
4555 | } |
---|
4556 | rowActivity_[i] = (valueScaled * scaleR) * inverseScale[i]; |
---|
4557 | double value = rowActivity_[i]; |
---|
4558 | if (value < rowLower_[i] - primalTolerance_) |
---|
4559 | numberPrimalUnscaled++; |
---|
4560 | else if (value > rowUpper_[i] + primalTolerance_) |
---|
4561 | numberPrimalUnscaled++; |
---|
4562 | double valueScaledDual = dual_[i] + rowObjectiveWork_[i]; |
---|
4563 | ; |
---|
4564 | if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_) |
---|
4565 | numberDualScaled++; |
---|
4566 | if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_) |
---|
4567 | numberDualScaled++; |
---|
4568 | dual_[i] *= scaleFactor * scaleC; |
---|
4569 | double valueDual = dual_[i]; |
---|
4570 | if (rowObjective_) |
---|
4571 | valueDual += rowObjective_[i]; |
---|
4572 | if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_) |
---|
4573 | numberDualUnscaled++; |
---|
4574 | if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_) |
---|
4575 | numberDualUnscaled++; |
---|
4576 | } |
---|
4577 | if (!problemStatus_ && !secondaryStatus_) { |
---|
4578 | // See if we need to set secondary status |
---|
4579 | if (numberPrimalUnscaled) { |
---|
4580 | if (numberDualUnscaled) |
---|
4581 | secondaryStatus_ = 4; |
---|
4582 | else |
---|
4583 | secondaryStatus_ = 2; |
---|
4584 | } else { |
---|
4585 | if (numberDualUnscaled) |
---|
4586 | secondaryStatus_ = 3; |
---|
4587 | } |
---|
4588 | } |
---|
4589 | if (problemStatus_ == 2 && ray_) { |
---|
4590 | for (i = 0; i < numberColumns; i++) { |
---|
4591 | ray_[i] *= columnScale_[i]; |
---|
4592 | } |
---|
4593 | } else if (problemStatus_ == 1 && ray_) { |
---|
4594 | for (i = 0; i < numberRows; i++) { |
---|
4595 | ray_[i] *= rowScale_[i]; |
---|
4596 | } |
---|
4597 | } |
---|
4598 | } else if (rhsScale_ != 1.0 || objectiveScale_ != 1.0) { |
---|
4599 | // Collect infeasibilities |
---|
4600 | int numberPrimalScaled = 0; |
---|
4601 | int numberPrimalUnscaled = 0; |
---|
4602 | int numberDualScaled = 0; |
---|
4603 | int numberDualUnscaled = 0; |
---|
4604 | double scaleC = 1.0 / objectiveScale_; |
---|
4605 | double scaleR = 1.0 / rhsScale_; |
---|
4606 | for (i = 0; i < numberColumns; i++) { |
---|
4607 | double valueScaled = columnActivityWork_[i]; |
---|
4608 | double lowerScaled = columnLowerWork_[i]; |
---|
4609 | double upperScaled = columnUpperWork_[i]; |
---|
4610 | if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) { |
---|
4611 | if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_) |
---|
4612 | numberPrimalScaled++; |
---|
4613 | else |
---|
4614 | upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled)); |
---|
4615 | } |
---|
4616 | columnActivity_[i] = valueScaled * scaleR; |
---|
4617 | double value = columnActivity_[i]; |
---|
4618 | if (value < columnLower_[i] - primalTolerance_) |
---|
4619 | numberPrimalUnscaled++; |
---|
4620 | else if (value > columnUpper_[i] + primalTolerance_) |
---|
4621 | numberPrimalUnscaled++; |
---|
4622 | double valueScaledDual = reducedCostWork_[i]; |
---|
4623 | if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_) |
---|
4624 | numberDualScaled++; |
---|
4625 | if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_) |
---|
4626 | numberDualScaled++; |
---|
4627 | reducedCost_[i] = valueScaledDual * scaleC; |
---|
4628 | double valueDual = reducedCost_[i]; |
---|
4629 | if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_) |
---|
4630 | numberDualUnscaled++; |
---|
4631 | if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_) |
---|
4632 | numberDualUnscaled++; |
---|
4633 | } |
---|
4634 | for (i = 0; i < numberRows; i++) { |
---|
4635 | double valueScaled = rowActivityWork_[i]; |
---|
4636 | double lowerScaled = rowLowerWork_[i]; |
---|
4637 | double upperScaled = rowUpperWork_[i]; |
---|
4638 | if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) { |
---|
4639 | if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_) |
---|
4640 | numberPrimalScaled++; |
---|
4641 | else |
---|
4642 | upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled)); |
---|
4643 | } |
---|
4644 | rowActivity_[i] = valueScaled * scaleR; |
---|
4645 | double value = rowActivity_[i]; |
---|
4646 | if (value < rowLower_[i] - primalTolerance_) |
---|
4647 | numberPrimalUnscaled++; |
---|
4648 | else if (value > rowUpper_[i] + primalTolerance_) |
---|
4649 | numberPrimalUnscaled++; |
---|
4650 | double valueScaledDual = dual_[i] + rowObjectiveWork_[i]; |
---|
4651 | ; |
---|
4652 | if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_) |
---|
4653 | numberDualScaled++; |
---|
4654 | if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_) |
---|
4655 | numberDualScaled++; |
---|
4656 | dual_[i] *= scaleC; |
---|
4657 | double valueDual = dual_[i]; |
---|
4658 | if (rowObjective_) |
---|
4659 | valueDual += rowObjective_[i]; |
---|
4660 | if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_) |
---|
4661 | numberDualUnscaled++; |
---|
4662 | if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_) |
---|
4663 | numberDualUnscaled++; |
---|
4664 | } |
---|
4665 | if (!problemStatus_ && !secondaryStatus_) { |
---|
4666 | // See if we need to set secondary status |
---|
4667 | if (numberPrimalUnscaled) { |
---|
4668 | if (numberDualUnscaled) |
---|
4669 | secondaryStatus_ = 4; |
---|
4670 | else |
---|
4671 | secondaryStatus_ = 2; |
---|
4672 | } else { |
---|
4673 | if (numberDualUnscaled) |
---|
4674 | secondaryStatus_ = 3; |
---|
4675 | } |
---|
4676 | } |
---|
4677 | } else { |
---|
4678 | if (columnActivityWork_) { |
---|
4679 | for (i = 0; i < numberColumns; i++) { |
---|
4680 | double value = columnActivityWork_[i]; |
---|
4681 | double lower = columnLowerWork_[i]; |
---|
4682 | double upper = columnUpperWork_[i]; |
---|
4683 | if (lower > -1.0e20 || upper < 1.0e20) { |
---|
4684 | if (value > lower && value < upper) |
---|
4685 | upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value)); |
---|
4686 | } |
---|
4687 | columnActivity_[i] = columnActivityWork_[i]; |
---|
4688 | reducedCost_[i] = reducedCostWork_[i]; |
---|
4689 | } |
---|
4690 | for (i = 0; i < numberRows; i++) { |
---|
4691 | double value = rowActivityWork_[i]; |
---|
4692 | double lower = rowLowerWork_[i]; |
---|
4693 | double upper = rowUpperWork_[i]; |
---|
4694 | if (lower > -1.0e20 || upper < 1.0e20) { |
---|
4695 | if (value > lower && value < upper) |
---|
4696 | upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value)); |
---|
4697 | } |
---|
4698 | rowActivity_[i] = rowActivityWork_[i]; |
---|
4699 | } |
---|
4700 | } |
---|
4701 | } |
---|
4702 | // switch off scalefactor if auto |
---|
4703 | if (automaticScale_) { |
---|
4704 | rhsScale_ = 1.0; |
---|
4705 | objectiveScale_ = 1.0; |
---|
4706 | } |
---|
4707 | if (optimizationDirection_ != 1.0) { |
---|
4708 | // and modify all dual signs |
---|
4709 | for (i = 0; i < numberColumns; i++) |
---|
4710 | reducedCost_[i] *= optimizationDirection_; |
---|
4711 | for (i = 0; i < numberRows; i++) |
---|
4712 | dual_[i] *= optimizationDirection_; |
---|
4713 | } |
---|
4714 | // scaling may have been turned off |
---|
4715 | scalingFlag_ = abs(scalingFlag_); |
---|
4716 | if (getRidOfFactorizationData > 0) { |
---|
4717 | gutsOfDelete(getRidOfFactorizationData + 1); |
---|
4718 | } else { |
---|
4719 | // at least get rid of nonLinearCost_ |
---|
4720 | delete nonLinearCost_; |
---|
4721 | nonLinearCost_ = NULL; |
---|
4722 | } |
---|
4723 | if (!rowObjective_ && problemStatus_ == 0 && objective_->type() == 1 && numberRows && numberColumns) { |
---|
4724 | // Redo objective value |
---|
4725 | double objectiveValue = 0.0; |
---|
4726 | const double *cost = objective(); |
---|
4727 | for (int i = 0; i < numberColumns; i++) { |
---|
4728 | double value = columnActivity_[i]; |
---|
4729 | objectiveValue += value * cost[i]; |
---|
4730 | } |
---|
4731 | //if (fabs(objectiveValue_ -objectiveValue*optimizationDirection())>1.0e-5) |
---|
4732 | //printf("old obj %g new %g\n",objectiveValue_, objectiveValue*optimizationDirection()); |
---|
4733 | objectiveValue_ = objectiveValue * optimizationDirection(); |
---|
4734 | } |
---|
4735 | // get rid of data |
---|
4736 | matrix_->generalExpanded(this, 13, scalingFlag_); |
---|
4737 | } |
---|
4738 | void ClpSimplex::setDualBound(double value) |
---|
4739 | { |
---|
4740 | if (value > 0.0) |
---|
4741 | dualBound_ = value; |
---|
4742 | } |
---|
4743 | void ClpSimplex::setInfeasibilityCost(double value) |
---|
4744 | { |
---|
4745 | if (value > 0.0) |
---|
4746 | infeasibilityCost_ = value; |
---|
4747 | } |
---|
4748 | void ClpSimplex::setNumberRefinements(int value) |
---|
4749 | { |
---|
4750 | if (value >= 0 && value < 10) |
---|
4751 | numberRefinements_ = value; |
---|
4752 | } |
---|
4753 | // Sets row pivot choice algorithm in dual |
---|
4754 | void ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot &choice) |
---|
4755 | { |
---|
4756 | delete dualRowPivot_; |
---|
4757 | dualRowPivot_ = choice.clone(true); |
---|
4758 | dualRowPivot_->setModel(this); |
---|
4759 | } |
---|
4760 | // Sets column pivot choice algorithm in primal |
---|
4761 | void ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot &choice) |
---|
4762 | { |
---|
4763 | delete primalColumnPivot_; |
---|
4764 | primalColumnPivot_ = choice.clone(true); |
---|
4765 | primalColumnPivot_->setModel(this); |
---|
4766 | } |
---|
4767 | void ClpSimplex::setFactorization(ClpFactorization &factorization) |
---|
4768 | { |
---|
4769 | if (factorization_) |
---|
4770 | factorization_->setFactorization(factorization); |
---|
4771 | else |
---|
4772 | factorization_ = new ClpFactorization(factorization, |
---|
4773 | numberRows_); |
---|
4774 | } |
---|
4775 | |
---|
4776 | // Swaps factorization |
---|
4777 | ClpFactorization * |
---|
4778 | ClpSimplex::swapFactorization(ClpFactorization *factorization) |
---|
4779 | { |
---|
4780 | ClpFactorization *swap = factorization_; |
---|
4781 | factorization_ = factorization; |
---|
4782 | return swap; |
---|
4783 | } |
---|
4784 | // Copies in factorization to existing one |
---|
4785 | void ClpSimplex::copyFactorization(ClpFactorization &factorization) |
---|
4786 | { |
---|
4787 | *factorization_ = factorization; |
---|
4788 | } |
---|
4789 | /* Perturbation: |
---|
4790 | -50 to +50 - perturb by this power of ten (-6 sounds good) |
---|
4791 | 100 - auto perturb if takes too long (1.0e-6 largest nonzero) |
---|
4792 | 101 - we are perturbed |
---|
4793 | 102 - don't try perturbing again |
---|
4794 | default is 100 |
---|
4795 | */ |
---|
4796 | void ClpSimplex::setPerturbation(int value) |
---|
4797 | { |
---|
4798 | if (value <= 102 && value >= -1000) { |
---|
4799 | perturbation_ = value; |
---|
4800 | } |
---|
4801 | } |
---|
4802 | // Sparsity on or off |
---|
4803 | bool ClpSimplex::sparseFactorization() const |
---|
4804 | { |
---|
4805 | return factorization_->sparseThreshold() != 0; |
---|
4806 | } |
---|
4807 | void ClpSimplex::setSparseFactorization(bool value) |
---|
4808 | { |
---|
4809 | if (value) { |
---|
4810 | if (!factorization_->sparseThreshold()) |
---|
4811 | factorization_->goSparse(); |
---|
4812 | } else { |
---|
4813 | factorization_->sparseThreshold(0); |
---|
4814 | } |
---|
4815 | } |
---|
4816 | void checkCorrect(ClpSimplex * /*model*/, int iRow, |
---|
4817 | const double *element, const int *rowStart, const int *rowLength, |
---|
4818 | const int *column, |
---|
4819 | const double *columnLower_, const double *columnUpper_, |
---|
4820 | int /*infiniteUpperC*/, |
---|
4821 | int /*infiniteLowerC*/, |
---|
4822 | double &maximumUpC, |
---|
4823 | double &maximumDownC) |
---|
4824 | { |
---|
4825 | int infiniteUpper = 0; |
---|
4826 | int infiniteLower = 0; |
---|
4827 | double maximumUp = 0.0; |
---|
4828 | double maximumDown = 0.0; |
---|
4829 | CoinBigIndex rStart = rowStart[iRow]; |
---|
4830 | CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow]; |
---|
4831 | CoinBigIndex j; |
---|
4832 | double large = 1.0e15; |
---|
4833 | int iColumn; |
---|
4834 | // Compute possible lower and upper ranges |
---|
4835 | |
---|
4836 | for (j = rStart; j < rEnd; ++j) { |
---|
4837 | double value = element[j]; |
---|
4838 | iColumn = column[j]; |
---|
4839 | if (value > 0.0) { |
---|
4840 | if (columnUpper_[iColumn] >= large) { |
---|
4841 | ++infiniteUpper; |
---|
4842 | } else { |
---|
4843 | maximumUp += columnUpper_[iColumn] * value; |
---|
4844 | } |
---|
4845 | if (columnLower_[iColumn] <= -large) { |
---|
4846 | ++infiniteLower; |
---|
4847 | } else { |
---|
4848 | maximumDown += columnLower_[iColumn] * value; |
---|
4849 | } |
---|
4850 | } else if (value < 0.0) { |
---|
4851 | if (columnUpper_[iColumn] >= large) { |
---|
4852 | ++infiniteLower; |
---|
4853 | } else { |
---|
4854 | maximumDown += columnUpper_[iColumn] * value; |
---|
4855 | } |
---|
4856 | if (columnLower_[iColumn] <= -large) { |
---|
4857 | ++infiniteUpper; |
---|
4858 | } else { |
---|
4859 | maximumUp += columnLower_[iColumn] * value; |
---|
4860 | } |
---|
4861 | } |
---|
4862 | } |
---|
4863 | //assert (infiniteLowerC==infiniteLower); |
---|
4864 | //assert (infiniteUpperC==infiniteUpper); |
---|
4865 | if (fabs(maximumUp - maximumUpC) > 1.0e-12 * CoinMax(fabs(maximumUp), fabs(maximumUpC))) |
---|
4866 | COIN_DETAIL_PRINT(printf("row %d comp up %g, true up %g\n", iRow, |
---|
4867 | maximumUpC, maximumUp)); |
---|
4868 | if (fabs(maximumDown - maximumDownC) > 1.0e-12 * CoinMax(fabs(maximumDown), fabs(maximumDownC))) |
---|
4869 | COIN_DETAIL_PRINT(printf("row %d comp down %g, true down %g\n", iRow, |
---|
4870 | maximumDownC, maximumDown)); |
---|
4871 | maximumUpC = maximumUp; |
---|
4872 | maximumDownC = maximumDown; |
---|
4873 | } |
---|
4874 | |
---|
4875 | /* Tightens primal bounds to make dual faster. Unless |
---|
4876 | fixed, bounds are slightly looser than they could be. |
---|
4877 | This is to make dual go faster and is probably not needed |
---|
4878 | with a presolve. Returns non-zero if problem infeasible |
---|
4879 | |
---|
4880 | Fudge for branch and bound - put bounds on columns of factor * |
---|
4881 | largest value (at continuous) - should improve stability |
---|
4882 | in branch and bound on infeasible branches (0.0 is off) |
---|
4883 | */ |
---|
4884 | int ClpSimplex::tightenPrimalBounds(double factor, int doTight, bool tightIntegers) |
---|
4885 | { |
---|
4886 | |
---|
4887 | // Get a row copy in standard format |
---|
4888 | CoinPackedMatrix copy; |
---|
4889 | copy.setExtraGap(0.0); |
---|
4890 | copy.setExtraMajor(0.0); |
---|
4891 | copy.reverseOrderedCopyOf(*matrix()); |
---|
4892 | // Matrix may have been created so get rid of it |
---|
4893 | matrix_->releasePackedMatrix(); |
---|
4894 | // get matrix data pointers |
---|
4895 | const int *column = copy.getIndices(); |
---|
4896 | const CoinBigIndex *rowStart = copy.getVectorStarts(); |
---|
4897 | const int *rowLength = copy.getVectorLengths(); |
---|
4898 | const double *element = copy.getElements(); |
---|
4899 | int numberChanged = 1, iPass = 0; |
---|
4900 | double large = largeValue(); // treat bounds > this as infinite |
---|
4901 | #ifndef NDEBUG |
---|
4902 | double large2 = 1.0e10 * large; |
---|
4903 | #endif |
---|
4904 | int numberInfeasible = 0; |
---|
4905 | int totalTightened = 0; |
---|
4906 | |
---|
4907 | double tolerance = primalTolerance(); |
---|
4908 | |
---|
4909 | // Save column bounds |
---|
4910 | double *saveLower = new double[numberColumns_]; |
---|
4911 | CoinMemcpyN(columnLower_, numberColumns_, saveLower); |
---|
4912 | double *saveUpper = new double[numberColumns_]; |
---|
4913 | CoinMemcpyN(columnUpper_, numberColumns_, saveUpper); |
---|
4914 | |
---|
4915 | int iRow, iColumn; |
---|
4916 | // If wanted compute a reasonable dualBound_ |
---|
4917 | if (factor == COIN_DBL_MAX) { |
---|
4918 | factor = 0.0; |
---|
4919 | if (dualBound_ == 1.0e10) { |
---|
4920 | // get largest scaled away from bound |
---|
4921 | double largest = 1.0e-12; |
---|
4922 | double largestScaled = 1.0e-12; |
---|
4923 | int iRow; |
---|
4924 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
4925 | double value = rowActivity_[iRow]; |
---|
4926 | double above = value - rowLower_[iRow]; |
---|
4927 | double below = rowUpper_[iRow] - value; |
---|
4928 | if (above < 1.0e12) { |
---|
4929 | largest = CoinMax(largest, above); |
---|
4930 | } |
---|
4931 | if (below < 1.0e12) { |
---|
4932 | largest = CoinMax(largest, below); |
---|
4933 | } |
---|
4934 | if (rowScale_) { |
---|
4935 | double multiplier = rowScale_[iRow]; |
---|
4936 | above *= multiplier; |
---|
4937 | below *= multiplier; |
---|
4938 | } |
---|
4939 | if (above < 1.0e12) { |
---|
4940 | largestScaled = CoinMax(largestScaled, above); |
---|
4941 | } |
---|
4942 | if (below < 1.0e12) { |
---|
4943 | largestScaled = CoinMax(largestScaled, below); |
---|
4944 | } |
---|
4945 | } |
---|
4946 | |
---|
4947 | int iColumn; |
---|
4948 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
4949 | double value = columnActivity_[iColumn]; |
---|
4950 | double above = value - columnLower_[iColumn]; |
---|
4951 | double below = columnUpper_[iColumn] - value; |
---|
4952 | if (above < 1.0e12) { |
---|
4953 | largest = CoinMax(largest, above); |
---|
4954 | } |
---|
4955 | if (below < 1.0e12) { |
---|
4956 | largest = CoinMax(largest, below); |
---|
4957 | } |
---|
4958 | if (columnScale_) { |
---|
4959 | double multiplier = 1.0 / columnScale_[iColumn]; |
---|
4960 | above *= multiplier; |
---|
4961 | below *= multiplier; |
---|
4962 | } |
---|
4963 | if (above < 1.0e12) { |
---|
4964 | largestScaled = CoinMax(largestScaled, above); |
---|
4965 | } |
---|
4966 | if (below < 1.0e12) { |
---|
4967 | largestScaled = CoinMax(largestScaled, below); |
---|
4968 | } |
---|
4969 | } |
---|
4970 | std::cout << "Largest (scaled) away from bound " << largestScaled |
---|
4971 | << " unscaled " << largest << std::endl; |
---|
4972 | dualBound_ = CoinMax(1.0001e7, CoinMin(100.0 * largest, 1.00001e10)); |
---|
4973 | } |
---|
4974 | } |
---|
4975 | |
---|
4976 | // If wanted - tighten column bounds using solution |
---|
4977 | if (factor) { |
---|
4978 | double largest = 0.0; |
---|
4979 | if (factor > 0.0) { |
---|
4980 | assert(factor > 1.0); |
---|
4981 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
4982 | if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) { |
---|
4983 | largest = CoinMax(largest, fabs(columnActivity_[iColumn])); |
---|
4984 | } |
---|
4985 | } |
---|
4986 | largest *= factor; |
---|
4987 | } else { |
---|
4988 | // absolute |
---|
4989 | largest = -factor; |
---|
4990 | } |
---|
4991 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
4992 | if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) { |
---|
4993 | columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn], largest); |
---|
4994 | columnLower_[iColumn] = CoinMax(columnLower_[iColumn], -largest); |
---|
4995 | } |
---|
4996 | } |
---|
4997 | } |
---|
4998 | #define MAXPASS 10 |
---|
4999 | |
---|
5000 | // Loop round seeing if we can tighten bounds |
---|
5001 | // Would be faster to have a stack of possible rows |
---|
5002 | // and we put altered rows back on stack |
---|
5003 | int numberCheck = -1; |
---|
5004 | while (numberChanged > numberCheck) { |
---|
5005 | |
---|
5006 | numberChanged = 0; // Bounds tightened this pass |
---|
5007 | |
---|
5008 | if (iPass == MAXPASS) |
---|
5009 | break; |
---|
5010 | iPass++; |
---|
5011 | |
---|
5012 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
5013 | |
---|
5014 | if (rowLower_[iRow] > -large || rowUpper_[iRow] < large) { |
---|
5015 | |
---|
5016 | // possible row |
---|
5017 | int infiniteUpper = 0; |
---|
5018 | int infiniteLower = 0; |
---|
5019 | double maximumUp = 0.0; |
---|
5020 | double maximumDown = 0.0; |
---|
5021 | double newBound; |
---|
5022 | CoinBigIndex rStart = rowStart[iRow]; |
---|
5023 | CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow]; |
---|
5024 | CoinBigIndex j; |
---|
5025 | // Compute possible lower and upper ranges |
---|
5026 | |
---|
5027 | for (j = rStart; j < rEnd; ++j) { |
---|
5028 | double value = element[j]; |
---|
5029 | iColumn = column[j]; |
---|
5030 | if (value > 0.0) { |
---|
5031 | if (columnUpper_[iColumn] >= large) { |
---|
5032 | ++infiniteUpper; |
---|
5033 | } else { |
---|
5034 | maximumUp += columnUpper_[iColumn] * value; |
---|
5035 | } |
---|
5036 | if (columnLower_[iColumn] <= -large) { |
---|
5037 | ++infiniteLower; |
---|
5038 | } else { |
---|
5039 | maximumDown += columnLower_[iColumn] * value; |
---|
5040 | } |
---|
5041 | } else if (value < 0.0) { |
---|
5042 | if (columnUpper_[iColumn] >= large) { |
---|
5043 | ++infiniteLower; |
---|
5044 | } else { |
---|
5045 | maximumDown += columnUpper_[iColumn] * value; |
---|
5046 | } |
---|
5047 | if (columnLower_[iColumn] <= -large) { |
---|
5048 | ++infiniteUpper; |
---|
5049 | } else { |
---|
5050 | maximumUp += columnLower_[iColumn] * value; |
---|
5051 | } |
---|
5052 | } |
---|
5053 | } |
---|
5054 | // Build in a margin of error |
---|
5055 | maximumUp += 1.0e-8 * fabs(maximumUp); |
---|
5056 | maximumDown -= 1.0e-8 * fabs(maximumDown); |
---|
5057 | double maxUp = maximumUp + infiniteUpper * 1.0e31; |
---|
5058 | double maxDown = maximumDown - infiniteLower * 1.0e31; |
---|
5059 | if (maxUp <= rowUpper_[iRow] + tolerance && maxDown >= rowLower_[iRow] - tolerance) { |
---|
5060 | |
---|
5061 | // Row is redundant - make totally free |
---|
5062 | // NO - can't do this for postsolve |
---|
5063 | // rowLower_[iRow]=-COIN_DBL_MAX; |
---|
5064 | // rowUpper_[iRow]=COIN_DBL_MAX; |
---|
5065 | //printf("Redundant row in presolveX %d\n",iRow); |
---|
5066 | |
---|
5067 | } else { |
---|
5068 | if (maxUp < rowLower_[iRow] - 100.0 * tolerance || maxDown > rowUpper_[iRow] + 100.0 * tolerance) { |
---|
5069 | // problem is infeasible - exit at once |
---|
5070 | numberInfeasible++; |
---|
5071 | break; |
---|
5072 | } |
---|
5073 | double lower = rowLower_[iRow]; |
---|
5074 | double upper = rowUpper_[iRow]; |
---|
5075 | for (j = rStart; j < rEnd; ++j) { |
---|
5076 | double value = element[j]; |
---|
5077 | iColumn = column[j]; |
---|
5078 | double nowLower = columnLower_[iColumn]; |
---|
5079 | double nowUpper = columnUpper_[iColumn]; |
---|
5080 | if (value > 0.0) { |
---|
5081 | // positive value |
---|
5082 | if (lower > -large) { |
---|
5083 | if (!infiniteUpper) { |
---|
5084 | assert(nowUpper < large2); |
---|
5085 | newBound = nowUpper + (lower - maximumUp) / value; |
---|
5086 | // relax if original was large |
---|
5087 | if (fabs(maximumUp) > 1.0e8) |
---|
5088 | newBound -= 1.0e-12 * fabs(maximumUp); |
---|
5089 | } else if (infiniteUpper == 1 && nowUpper > large) { |
---|
5090 | newBound = (lower - maximumUp) / value; |
---|
5091 | // relax if original was large |
---|
5092 | if (fabs(maximumUp) > 1.0e8) |
---|
5093 | newBound -= 1.0e-12 * fabs(maximumUp); |
---|
5094 | } else { |
---|
5095 | newBound = -COIN_DBL_MAX; |
---|
5096 | } |
---|
5097 | if (newBound > nowLower + 1.0e-12 && newBound > -large) { |
---|
5098 | // Tighten the lower bound |
---|
5099 | numberChanged++; |
---|
5100 | // check infeasible (relaxed) |
---|
5101 | if (nowUpper < newBound) { |
---|
5102 | if (nowUpper - newBound < -100.0 * tolerance) |
---|
5103 | numberInfeasible++; |
---|
5104 | else |
---|
5105 | newBound = nowUpper; |
---|
5106 | } |
---|
5107 | columnLower_[iColumn] = newBound; |
---|
5108 | // adjust |
---|
5109 | double now; |
---|
5110 | if (nowLower < -large) { |
---|
5111 | now = 0.0; |
---|
5112 | infiniteLower--; |
---|
5113 | } else { |
---|
5114 | now = nowLower; |
---|
5115 | } |
---|
5116 | maximumDown += (newBound - now) * value; |
---|
5117 | nowLower = newBound; |
---|
5118 | #ifdef DEBUG |
---|
5119 | checkCorrect(this, iRow, |
---|
5120 | element, rowStart, rowLength, |
---|
5121 | column, |
---|
5122 | columnLower_, columnUpper_, |
---|
5123 | infiniteUpper, |
---|
5124 | infiniteLower, |
---|
5125 | maximumUp, |
---|
5126 | maximumDown); |
---|
5127 | #endif |
---|
5128 | } |
---|
5129 | } |
---|
5130 | if (upper < large) { |
---|
5131 | if (!infiniteLower) { |
---|
5132 | assert(nowLower > -large2); |
---|
5133 | newBound = nowLower + (upper - maximumDown) / value; |
---|
5134 | // relax if original was large |
---|
5135 | if (fabs(maximumDown) > 1.0e8) |
---|
5136 | newBound += 1.0e-12 * fabs(maximumDown); |
---|
5137 | } else if (infiniteLower == 1 && nowLower < -large) { |
---|
5138 | newBound = (upper - maximumDown) / value; |
---|
5139 | // relax if original was large |
---|
5140 | if (fabs(maximumDown) > 1.0e8) |
---|
5141 | newBound += 1.0e-12 * fabs(maximumDown); |
---|
5142 | } else { |
---|
5143 | newBound = COIN_DBL_MAX; |
---|
5144 | } |
---|
5145 | if (newBound < nowUpper - 1.0e-12 && newBound < large) { |
---|
5146 | // Tighten the upper bound |
---|
5147 | numberChanged++; |
---|
5148 | // check infeasible (relaxed) |
---|
5149 | if (nowLower > newBound) { |
---|
5150 | if (newBound - nowLower < -100.0 * tolerance) |
---|
5151 | numberInfeasible++; |
---|
5152 | else |
---|
5153 | newBound = nowLower; |
---|
5154 | } |
---|
5155 | columnUpper_[iColumn] = newBound; |
---|
5156 | // adjust |
---|
5157 | double now; |
---|
5158 | if (nowUpper > large) { |
---|
5159 | now = 0.0; |
---|
5160 | infiniteUpper--; |
---|
5161 | } else { |
---|
5162 | now = nowUpper; |
---|
5163 | } |
---|
5164 | maximumUp += (newBound - now) * value; |
---|
5165 | nowUpper = newBound; |
---|
5166 | #ifdef DEBUG |
---|
5167 | checkCorrect(this, iRow, |
---|
5168 | element, rowStart, rowLength, |
---|
5169 | column, |
---|
5170 | columnLower_, columnUpper_, |
---|
5171 | infiniteUpper, |
---|
5172 | infiniteLower, |
---|
5173 | maximumUp, |
---|
5174 | maximumDown); |
---|
5175 | #endif |
---|
5176 | } |
---|
5177 | } |
---|
5178 | } else { |
---|
5179 | // negative value |
---|
5180 | if (lower > -large) { |
---|
5181 | if (!infiniteUpper) { |
---|
5182 | assert(nowLower < large2); |
---|
5183 | newBound = nowLower + (lower - maximumUp) / value; |
---|
5184 | // relax if original was large |
---|
5185 | if (fabs(maximumUp) > 1.0e8) |
---|
5186 | newBound += 1.0e-12 * fabs(maximumUp); |
---|
5187 | } else if (infiniteUpper == 1 && nowLower < -large) { |
---|
5188 | newBound = (lower - maximumUp) / value; |
---|
5189 | // relax if original was large |
---|
5190 | if (fabs(maximumUp) > 1.0e8) |
---|
5191 | newBound += 1.0e-12 * fabs(maximumUp); |
---|
5192 | } else { |
---|
5193 | newBound = COIN_DBL_MAX; |
---|
5194 | } |
---|
5195 | if (newBound < nowUpper - 1.0e-12 && newBound < large) { |
---|
5196 | // Tighten the upper bound |
---|
5197 | numberChanged++; |
---|
5198 | // check infeasible (relaxed) |
---|
5199 | if (nowLower > newBound) { |
---|
5200 | if (newBound - nowLower < -100.0 * tolerance) |
---|
5201 | numberInfeasible++; |
---|
5202 | else |
---|
5203 | newBound = nowLower; |
---|
5204 | } |
---|
5205 | columnUpper_[iColumn] = newBound; |
---|
5206 | // adjust |
---|
5207 | double now; |
---|
5208 | if (nowUpper > large) { |
---|
5209 | now = 0.0; |
---|
5210 | infiniteLower--; |
---|
5211 | } else { |
---|
5212 | now = nowUpper; |
---|
5213 | } |
---|
5214 | maximumDown += (newBound - now) * value; |
---|
5215 | nowUpper = newBound; |
---|
5216 | #ifdef DEBUG |
---|
5217 | checkCorrect(this, iRow, |
---|
5218 | element, rowStart, rowLength, |
---|
5219 | column, |
---|
5220 | columnLower_, columnUpper_, |
---|
5221 | infiniteUpper, |
---|
5222 | infiniteLower, |
---|
5223 | maximumUp, |
---|
5224 | maximumDown); |
---|
5225 | #endif |
---|
5226 | } |
---|
5227 | } |
---|
5228 | if (upper < large) { |
---|
5229 | if (!infiniteLower) { |
---|
5230 | assert(nowUpper < large2); |
---|
5231 | newBound = nowUpper + (upper - maximumDown) / value; |
---|
5232 | // relax if original was large |
---|
5233 | if (fabs(maximumDown) > 1.0e8) |
---|
5234 | newBound -= 1.0e-12 * fabs(maximumDown); |
---|
5235 | } else if (infiniteLower == 1 && nowUpper > large) { |
---|
5236 | newBound = (upper - maximumDown) / value; |
---|
5237 | // relax if original was large |
---|
5238 | if (fabs(maximumDown) > 1.0e8) |
---|
5239 | newBound -= 1.0e-12 * fabs(maximumDown); |
---|
5240 | } else { |
---|
5241 | newBound = -COIN_DBL_MAX; |
---|
5242 | } |
---|
5243 | if (newBound > nowLower + 1.0e-12 && newBound > -large) { |
---|
5244 | // Tighten the lower bound |
---|
5245 | numberChanged++; |
---|
5246 | // check infeasible (relaxed) |
---|
5247 | if (nowUpper < newBound) { |
---|
5248 | if (nowUpper - newBound < -100.0 * tolerance) |
---|
5249 | numberInfeasible++; |
---|
5250 | else |
---|
5251 | newBound = nowUpper; |
---|
5252 | } |
---|
5253 | columnLower_[iColumn] = newBound; |
---|
5254 | // adjust |
---|
5255 | double now; |
---|
5256 | if (nowLower < -large) { |
---|
5257 | now = 0.0; |
---|
5258 | infiniteUpper--; |
---|
5259 | } else { |
---|
5260 | now = nowLower; |
---|
5261 | } |
---|
5262 | maximumUp += (newBound - now) * value; |
---|
5263 | nowLower = newBound; |
---|
5264 | #ifdef DEBUG |
---|
5265 | checkCorrect(this, iRow, |
---|
5266 | element, rowStart, rowLength, |
---|
5267 | column, |
---|
5268 | columnLower_, columnUpper_, |
---|
5269 | infiniteUpper, |
---|
5270 | infiniteLower, |
---|
5271 | maximumUp, |
---|
5272 | maximumDown); |
---|
5273 | #endif |
---|
5274 | } |
---|
5275 | } |
---|
5276 | } |
---|
5277 | } |
---|
5278 | } |
---|
5279 | } |
---|
5280 | } |
---|
5281 | totalTightened += numberChanged; |
---|
5282 | if (iPass == 1) |
---|
5283 | numberCheck = numberChanged >> 4; |
---|
5284 | if (numberInfeasible) |
---|
5285 | break; |
---|
5286 | } |
---|
5287 | if (!numberInfeasible) { |
---|
5288 | handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN, messages_) |
---|
5289 | << totalTightened |
---|
5290 | << CoinMessageEol; |
---|
5291 | // Set bounds slightly loose |
---|
5292 | double useTolerance = 1.0e-3; |
---|
5293 | if (doTight > 0) { |
---|
5294 | if (doTight > 10) { |
---|
5295 | useTolerance = 0.0; |
---|
5296 | } else { |
---|
5297 | while (doTight) { |
---|
5298 | useTolerance *= 0.1; |
---|
5299 | doTight--; |
---|
5300 | } |
---|
5301 | } |
---|
5302 | } |
---|
5303 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
5304 | if (saveUpper[iColumn] > saveLower[iColumn] + tolerance) { |
---|
5305 | // Make large bounds stay infinite |
---|
5306 | if (saveUpper[iColumn] > 1.0e30 && columnUpper_[iColumn] > 1.0e10) { |
---|
5307 | columnUpper_[iColumn] = COIN_DBL_MAX; |
---|
5308 | } |
---|
5309 | if (saveLower[iColumn] < -1.0e30 && columnLower_[iColumn] < -1.0e10) { |
---|
5310 | columnLower_[iColumn] = -COIN_DBL_MAX; |
---|
5311 | } |
---|
5312 | #ifdef KEEP_GOING_IF_FIXED |
---|
5313 | double multiplier = 5.0e-3 * floor(100.0 * randomNumberGenerator_.randomDouble()) + 1.0; |
---|
5314 | multiplier *= 100.0; |
---|
5315 | #else |
---|
5316 | double multiplier = 100.0; |
---|
5317 | #endif |
---|
5318 | if (columnUpper_[iColumn] - columnLower_[iColumn] < useTolerance + 1.0e-8) { |
---|
5319 | // relax enough so will have correct dj |
---|
5320 | #if 1 |
---|
5321 | columnLower_[iColumn] = CoinMax(saveLower[iColumn], |
---|
5322 | columnLower_[iColumn] - multiplier * useTolerance); |
---|
5323 | columnUpper_[iColumn] = CoinMin(saveUpper[iColumn], |
---|
5324 | columnUpper_[iColumn] + multiplier * useTolerance); |
---|
5325 | #else |
---|
5326 | if (fabs(columnUpper_[iColumn]) < fabs(columnLower_[iColumn])) { |
---|
5327 | if (columnUpper_[iColumn] - multiplier * useTolerance > saveLower[iColumn]) { |
---|
5328 | columnLower_[iColumn] = columnUpper_[iColumn] - multiplier * useTolerance; |
---|
5329 | } else { |
---|
5330 | columnLower_[iColumn] = saveLower[iColumn]; |
---|
5331 | columnUpper_[iColumn] = CoinMin(saveUpper[iColumn], |
---|
5332 | saveLower[iColumn] + multiplier * useTolerance); |
---|
5333 | } |
---|
5334 | } else { |
---|
5335 | if (columnLower_[iColumn] + multiplier * useTolerance < saveUpper[iColumn]) { |
---|
5336 | columnUpper_[iColumn] = columnLower_[iColumn] + multiplier * useTolerance; |
---|
5337 | } else { |
---|
5338 | columnUpper_[iColumn] = saveUpper[iColumn]; |
---|
5339 | columnLower_[iColumn] = CoinMax(saveLower[iColumn], |
---|
5340 | saveUpper[iColumn] - multiplier * useTolerance); |
---|
5341 | } |
---|
5342 | } |
---|
5343 | #endif |
---|
5344 | } else { |
---|
5345 | if (columnUpper_[iColumn] < saveUpper[iColumn]) { |
---|
5346 | // relax a bit |
---|
5347 | columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn] + multiplier * useTolerance, |
---|
5348 | saveUpper[iColumn]); |
---|
5349 | } |
---|
5350 | if (columnLower_[iColumn] > saveLower[iColumn]) { |
---|
5351 | // relax a bit |
---|
5352 | columnLower_[iColumn] = CoinMax(columnLower_[iColumn] - multiplier * useTolerance, |
---|
5353 | saveLower[iColumn]); |
---|
5354 | } |
---|
5355 | } |
---|
5356 | } |
---|
5357 | } |
---|
5358 | if (tightIntegers && integerType_) { |
---|
5359 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
5360 | if (integerType_[iColumn]) { |
---|
5361 | double value; |
---|
5362 | value = floor(columnLower_[iColumn] + 0.5); |
---|
5363 | if (fabs(value - columnLower_[iColumn]) > primalTolerance_) |
---|
5364 | value = ceil(columnLower_[iColumn]); |
---|
5365 | columnLower_[iColumn] = value; |
---|
5366 | value = floor(columnUpper_[iColumn] + 0.5); |
---|
5367 | if (fabs(value - columnUpper_[iColumn]) > primalTolerance_) |
---|
5368 | value = floor(columnUpper_[iColumn]); |
---|
5369 | columnUpper_[iColumn] = value; |
---|
5370 | if (columnLower_[iColumn] > columnUpper_[iColumn]) |
---|
5371 | numberInfeasible++; |
---|
5372 | } |
---|
5373 | } |
---|
5374 | if (numberInfeasible) { |
---|
5375 | handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_) |
---|
5376 | << numberInfeasible |
---|
5377 | << CoinMessageEol; |
---|
5378 | // restore column bounds |
---|
5379 | CoinMemcpyN(saveLower, numberColumns_, columnLower_); |
---|
5380 | CoinMemcpyN(saveUpper, numberColumns_, columnUpper_); |
---|
5381 | } |
---|
5382 | } |
---|
5383 | } else { |
---|
5384 | handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_) |
---|
5385 | << numberInfeasible |
---|
5386 | << CoinMessageEol; |
---|
5387 | // restore column bounds |
---|
5388 | CoinMemcpyN(saveLower, numberColumns_, columnLower_); |
---|
5389 | CoinMemcpyN(saveUpper, numberColumns_, columnUpper_); |
---|
5390 | } |
---|
5391 | delete[] saveLower; |
---|
5392 | delete[] saveUpper; |
---|
5393 | return (numberInfeasible); |
---|
5394 | } |
---|
5395 | //#define SAVE_AND_RESTORE |
---|
5396 | // dual |
---|
5397 | #include "ClpSimplexDual.hpp" |
---|
5398 | #include "ClpSimplexPrimal.hpp" |
---|
5399 | #ifndef SAVE_AND_RESTORE |
---|
5400 | int ClpSimplex::dual(int ifValuesPass, int startFinishOptions) |
---|
5401 | #else |
---|
5402 | int ClpSimplex::dual(int ifValuesPass, int startFinishOptions) |
---|
5403 | { |
---|
5404 | // May be empty problem |
---|
5405 | if (numberRows_ && numberColumns_) { |
---|
5406 | // Save on file for debug |
---|
5407 | int returnCode; |
---|
5408 | returnCode = saveModel("debug.sav"); |
---|
5409 | if (returnCode) |
---|