1 | /* $Id: ClpSimplexOther.cpp 2128 2015-03-16 08:57:46Z forrest $ */ |
---|
2 | // Copyright (C) 2004, 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 | #include "CoinPragma.hpp" |
---|
7 | |
---|
8 | #include <math.h> |
---|
9 | |
---|
10 | #include "CoinHelperFunctions.hpp" |
---|
11 | #include "ClpSimplexOther.hpp" |
---|
12 | #include "ClpSimplexDual.hpp" |
---|
13 | #include "ClpSimplexPrimal.hpp" |
---|
14 | #include "ClpEventHandler.hpp" |
---|
15 | #include "ClpHelperFunctions.hpp" |
---|
16 | #include "ClpFactorization.hpp" |
---|
17 | #include "ClpDualRowDantzig.hpp" |
---|
18 | #include "ClpNonLinearCost.hpp" |
---|
19 | #include "ClpDynamicMatrix.hpp" |
---|
20 | #include "CoinPackedMatrix.hpp" |
---|
21 | #include "CoinIndexedVector.hpp" |
---|
22 | #include "CoinBuild.hpp" |
---|
23 | #include "CoinMpsIO.hpp" |
---|
24 | #include "CoinFloatEqual.hpp" |
---|
25 | #include "ClpMessage.hpp" |
---|
26 | #include <cfloat> |
---|
27 | #include <cassert> |
---|
28 | #include <string> |
---|
29 | #include <stdio.h> |
---|
30 | #include <iostream> |
---|
31 | #ifdef INT_IS_8 |
---|
32 | #define COIN_ANY_BITS_PER_INT 64 |
---|
33 | #define COIN_ANY_SHIFT_PER_INT 6 |
---|
34 | #define COIN_ANY_MASK_PER_INT 0x3f |
---|
35 | #else |
---|
36 | #define COIN_ANY_BITS_PER_INT 32 |
---|
37 | #define COIN_ANY_SHIFT_PER_INT 5 |
---|
38 | #define COIN_ANY_MASK_PER_INT 0x1f |
---|
39 | #endif |
---|
40 | /* Dual ranging. |
---|
41 | This computes increase/decrease in cost for each given variable and corresponding |
---|
42 | sequence numbers which would change basis. Sequence numbers are 0..numberColumns |
---|
43 | and numberColumns.. for artificials/slacks. |
---|
44 | For non-basic variables the sequence number will be that of the non-basic variables. |
---|
45 | |
---|
46 | Up to user to provide correct length arrays. |
---|
47 | |
---|
48 | */ |
---|
49 | void ClpSimplexOther::dualRanging(int numberCheck, const int * which, |
---|
50 | double * costIncreased, int * sequenceIncreased, |
---|
51 | double * costDecreased, int * sequenceDecreased, |
---|
52 | double * valueIncrease, double * valueDecrease) |
---|
53 | { |
---|
54 | rowArray_[1]->clear(); |
---|
55 | columnArray_[1]->clear(); |
---|
56 | // long enough for rows+columns |
---|
57 | assert(rowArray_[3]->capacity() >= numberRows_ + numberColumns_); |
---|
58 | rowArray_[3]->clear(); |
---|
59 | int * backPivot = rowArray_[3]->getIndices(); |
---|
60 | int i; |
---|
61 | for ( i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
62 | backPivot[i] = -1; |
---|
63 | } |
---|
64 | for (i = 0; i < numberRows_; i++) { |
---|
65 | int iSequence = pivotVariable_[i]; |
---|
66 | backPivot[iSequence] = i; |
---|
67 | } |
---|
68 | // dualTolerance may be zero if from CBC. In fact use that fact |
---|
69 | bool inCBC = !dualTolerance_; |
---|
70 | if (inCBC) |
---|
71 | assert (integerType_); |
---|
72 | dualTolerance_ = dblParam_[ClpDualTolerance]; |
---|
73 | double * arrayX = rowArray_[0]->denseVector(); |
---|
74 | for ( i = 0; i < numberCheck; i++) { |
---|
75 | rowArray_[0]->clear(); |
---|
76 | //rowArray_[0]->checkClear(); |
---|
77 | //rowArray_[1]->checkClear(); |
---|
78 | //columnArray_[1]->checkClear(); |
---|
79 | columnArray_[0]->clear(); |
---|
80 | //columnArray_[0]->checkClear(); |
---|
81 | int iSequence = which[i]; |
---|
82 | if (iSequence < 0) { |
---|
83 | costIncreased[i] = 0.0; |
---|
84 | sequenceIncreased[i] = -1; |
---|
85 | costDecreased[i] = 0.0; |
---|
86 | sequenceDecreased[i] = -1; |
---|
87 | continue; |
---|
88 | } |
---|
89 | double costIncrease = COIN_DBL_MAX; |
---|
90 | double costDecrease = COIN_DBL_MAX; |
---|
91 | int sequenceIncrease = -1; |
---|
92 | int sequenceDecrease = -1; |
---|
93 | if (valueIncrease) { |
---|
94 | assert (valueDecrease); |
---|
95 | valueIncrease[i] = iSequence < numberColumns_ ? columnActivity_[iSequence] : rowActivity_[iSequence-numberColumns_]; |
---|
96 | valueDecrease[i] = valueIncrease[i]; |
---|
97 | } |
---|
98 | |
---|
99 | switch(getStatus(iSequence)) { |
---|
100 | |
---|
101 | case basic: { |
---|
102 | // non-trvial |
---|
103 | // Get pivot row |
---|
104 | int iRow = backPivot[iSequence]; |
---|
105 | assert (iRow >= 0); |
---|
106 | #ifndef COIN_FAC_NEW |
---|
107 | double plusOne = 1.0; |
---|
108 | rowArray_[0]->createPacked(1, &iRow, &plusOne); |
---|
109 | #else |
---|
110 | rowArray_[0]->createOneUnpackedElement( iRow, 1.0); |
---|
111 | #endif |
---|
112 | factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]); |
---|
113 | // put row of tableau in rowArray[0] and columnArray[0] |
---|
114 | matrix_->transposeTimes(this, -1.0, |
---|
115 | rowArray_[0], columnArray_[1], columnArray_[0]); |
---|
116 | #ifdef COIN_FAC_NEW |
---|
117 | assert (!rowArray_[0]->packedMode()); |
---|
118 | #endif |
---|
119 | double alphaIncrease; |
---|
120 | double alphaDecrease; |
---|
121 | // do ratio test up and down |
---|
122 | checkDualRatios(rowArray_[0], columnArray_[0], costIncrease, sequenceIncrease, alphaIncrease, |
---|
123 | costDecrease, sequenceDecrease, alphaDecrease); |
---|
124 | if (!inCBC) { |
---|
125 | if (valueIncrease) { |
---|
126 | if (sequenceIncrease >= 0) |
---|
127 | valueIncrease[i] = primalRanging1(sequenceIncrease, iSequence); |
---|
128 | if (sequenceDecrease >= 0) |
---|
129 | valueDecrease[i] = primalRanging1(sequenceDecrease, iSequence); |
---|
130 | } |
---|
131 | } else { |
---|
132 | int number = rowArray_[0]->getNumElements(); |
---|
133 | #ifdef COIN_FAC_NEW |
---|
134 | const int * index = rowArray_[0]->getIndices(); |
---|
135 | #endif |
---|
136 | double scale2 = 0.0; |
---|
137 | int j; |
---|
138 | for (j = 0; j < number; j++) { |
---|
139 | #ifndef COIN_FAC_NEW |
---|
140 | scale2 += arrayX[j] * arrayX[j]; |
---|
141 | #else |
---|
142 | int iRow=index[j]; |
---|
143 | scale2 += arrayX[iRow] * arrayX[iRow]; |
---|
144 | #endif |
---|
145 | } |
---|
146 | scale2 = 1.0 / sqrt(scale2); |
---|
147 | //valueIncrease[i] = scale2; |
---|
148 | if (sequenceIncrease >= 0) { |
---|
149 | double djValue = dj_[sequenceIncrease]; |
---|
150 | if (fabs(djValue) > 10.0 * dualTolerance_) { |
---|
151 | // we are going to use for cutoff so be exact |
---|
152 | costIncrease = fabs(djValue / alphaIncrease); |
---|
153 | /* Not sure this is good idea as I don't think correct e.g. |
---|
154 | suppose a continuous variable has dj slightly greater. */ |
---|
155 | if(false && sequenceIncrease < numberColumns_ && integerType_[sequenceIncrease]) { |
---|
156 | // can improve |
---|
157 | double movement = (columnScale_ == NULL) ? 1.0 : |
---|
158 | rhsScale_ * inverseColumnScale_[sequenceIncrease]; |
---|
159 | costIncrease = CoinMax(fabs(djValue * movement), costIncrease); |
---|
160 | } |
---|
161 | } else { |
---|
162 | costIncrease = 0.0; |
---|
163 | } |
---|
164 | } |
---|
165 | if (sequenceDecrease >= 0) { |
---|
166 | double djValue = dj_[sequenceDecrease]; |
---|
167 | if (fabs(djValue) > 10.0 * dualTolerance_) { |
---|
168 | // we are going to use for cutoff so be exact |
---|
169 | costDecrease = fabs(djValue / alphaDecrease); |
---|
170 | if(sequenceDecrease < numberColumns_ && integerType_[sequenceDecrease]) { |
---|
171 | // can improve |
---|
172 | double movement = (columnScale_ == NULL) ? 1.0 : |
---|
173 | rhsScale_ * inverseColumnScale_[sequenceDecrease]; |
---|
174 | costDecrease = CoinMax(fabs(djValue * movement), costDecrease); |
---|
175 | } |
---|
176 | } else { |
---|
177 | costDecrease = 0.0; |
---|
178 | } |
---|
179 | } |
---|
180 | costIncrease *= scale2; |
---|
181 | costDecrease *= scale2; |
---|
182 | } |
---|
183 | } |
---|
184 | break; |
---|
185 | case isFixed: |
---|
186 | break; |
---|
187 | case isFree: |
---|
188 | case superBasic: |
---|
189 | costIncrease = 0.0; |
---|
190 | costDecrease = 0.0; |
---|
191 | sequenceIncrease = iSequence; |
---|
192 | sequenceDecrease = iSequence; |
---|
193 | break; |
---|
194 | case atUpperBound: |
---|
195 | costIncrease = CoinMax(0.0, -dj_[iSequence]); |
---|
196 | sequenceIncrease = iSequence; |
---|
197 | if (valueIncrease) |
---|
198 | valueIncrease[i] = primalRanging1(iSequence, iSequence); |
---|
199 | break; |
---|
200 | case atLowerBound: |
---|
201 | costDecrease = CoinMax(0.0, dj_[iSequence]); |
---|
202 | sequenceDecrease = iSequence; |
---|
203 | if (valueIncrease) |
---|
204 | valueDecrease[i] = primalRanging1(iSequence, iSequence); |
---|
205 | break; |
---|
206 | } |
---|
207 | double scaleFactor; |
---|
208 | if (rowScale_) { |
---|
209 | if (iSequence < numberColumns_) |
---|
210 | scaleFactor = 1.0 / (objectiveScale_ * columnScale_[iSequence]); |
---|
211 | else |
---|
212 | scaleFactor = rowScale_[iSequence-numberColumns_] / objectiveScale_; |
---|
213 | } else { |
---|
214 | scaleFactor = 1.0 / objectiveScale_; |
---|
215 | } |
---|
216 | if (costIncrease < 1.0e30) |
---|
217 | costIncrease *= scaleFactor; |
---|
218 | if (costDecrease < 1.0e30) |
---|
219 | costDecrease *= scaleFactor; |
---|
220 | if (optimizationDirection_ == 1.0) { |
---|
221 | costIncreased[i] = costIncrease; |
---|
222 | sequenceIncreased[i] = sequenceIncrease; |
---|
223 | costDecreased[i] = costDecrease; |
---|
224 | sequenceDecreased[i] = sequenceDecrease; |
---|
225 | } else if (optimizationDirection_ == -1.0) { |
---|
226 | costIncreased[i] = costDecrease; |
---|
227 | sequenceIncreased[i] = sequenceDecrease; |
---|
228 | costDecreased[i] = costIncrease; |
---|
229 | sequenceDecreased[i] = sequenceIncrease; |
---|
230 | if (valueIncrease) { |
---|
231 | double temp = valueIncrease[i]; |
---|
232 | valueIncrease[i] = valueDecrease[i]; |
---|
233 | valueDecrease[i] = temp; |
---|
234 | } |
---|
235 | } else if (optimizationDirection_ == 0.0) { |
---|
236 | // !!!!!! ??? |
---|
237 | costIncreased[i] = COIN_DBL_MAX; |
---|
238 | sequenceIncreased[i] = -1; |
---|
239 | costDecreased[i] = COIN_DBL_MAX; |
---|
240 | sequenceDecreased[i] = -1; |
---|
241 | } else { |
---|
242 | abort(); |
---|
243 | } |
---|
244 | } |
---|
245 | rowArray_[0]->clear(); |
---|
246 | //rowArray_[1]->clear(); |
---|
247 | //columnArray_[1]->clear(); |
---|
248 | columnArray_[0]->clear(); |
---|
249 | //rowArray_[3]->clear(); |
---|
250 | if (!optimizationDirection_) |
---|
251 | printf("*** ????? Ranging with zero optimization costs\n"); |
---|
252 | } |
---|
253 | /* |
---|
254 | Row array has row part of pivot row |
---|
255 | Column array has column part. |
---|
256 | This is used in dual ranging |
---|
257 | */ |
---|
258 | void |
---|
259 | ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray, |
---|
260 | CoinIndexedVector * columnArray, |
---|
261 | double & costIncrease, int & sequenceIncrease, double & alphaIncrease, |
---|
262 | double & costDecrease, int & sequenceDecrease, double & alphaDecrease) |
---|
263 | { |
---|
264 | double acceptablePivot = 1.0e-9; |
---|
265 | double * work; |
---|
266 | int number; |
---|
267 | int * which; |
---|
268 | int iSection; |
---|
269 | |
---|
270 | double thetaDown = 1.0e31; |
---|
271 | double thetaUp = 1.0e31; |
---|
272 | int sequenceDown = -1; |
---|
273 | int sequenceUp = -1; |
---|
274 | double alphaDown = 0.0; |
---|
275 | double alphaUp = 0.0; |
---|
276 | |
---|
277 | int addSequence; |
---|
278 | |
---|
279 | for (iSection = 0; iSection < 2; iSection++) { |
---|
280 | |
---|
281 | int i; |
---|
282 | if (!iSection) { |
---|
283 | work = rowArray->denseVector(); |
---|
284 | number = rowArray->getNumElements(); |
---|
285 | which = rowArray->getIndices(); |
---|
286 | addSequence = numberColumns_; |
---|
287 | } else { |
---|
288 | work = columnArray->denseVector(); |
---|
289 | number = columnArray->getNumElements(); |
---|
290 | which = columnArray->getIndices(); |
---|
291 | addSequence = 0; |
---|
292 | } |
---|
293 | |
---|
294 | for (i = 0; i < number; i++) { |
---|
295 | int iSequence = which[i]; |
---|
296 | int iSequence2 = iSequence + addSequence; |
---|
297 | #ifndef COIN_FAC_NEW |
---|
298 | double alpha = work[i]; |
---|
299 | #else |
---|
300 | double alpha = !addSequence ? work[i] : work[iSequence]; |
---|
301 | #endif |
---|
302 | if (fabs(alpha) < acceptablePivot) |
---|
303 | continue; |
---|
304 | double oldValue = dj_[iSequence2]; |
---|
305 | |
---|
306 | switch(getStatus(iSequence2)) { |
---|
307 | |
---|
308 | case basic: |
---|
309 | break; |
---|
310 | case ClpSimplex::isFixed: |
---|
311 | break; |
---|
312 | case isFree: |
---|
313 | case superBasic: |
---|
314 | // treat dj as if zero |
---|
315 | thetaDown = 0.0; |
---|
316 | thetaUp = 0.0; |
---|
317 | sequenceDown = iSequence2; |
---|
318 | sequenceUp = iSequence2; |
---|
319 | break; |
---|
320 | case atUpperBound: |
---|
321 | if (alpha > 0.0) { |
---|
322 | // test up |
---|
323 | if (oldValue + thetaUp * alpha > dualTolerance_) { |
---|
324 | thetaUp = (dualTolerance_ - oldValue) / alpha; |
---|
325 | sequenceUp = iSequence2; |
---|
326 | alphaUp = alpha; |
---|
327 | } |
---|
328 | } else { |
---|
329 | // test down |
---|
330 | if (oldValue - thetaDown * alpha > dualTolerance_) { |
---|
331 | thetaDown = -(dualTolerance_ - oldValue) / alpha; |
---|
332 | sequenceDown = iSequence2; |
---|
333 | alphaDown = alpha; |
---|
334 | } |
---|
335 | } |
---|
336 | break; |
---|
337 | case atLowerBound: |
---|
338 | if (alpha < 0.0) { |
---|
339 | // test up |
---|
340 | if (oldValue + thetaUp * alpha < - dualTolerance_) { |
---|
341 | thetaUp = -(dualTolerance_ + oldValue) / alpha; |
---|
342 | sequenceUp = iSequence2; |
---|
343 | alphaUp = alpha; |
---|
344 | } |
---|
345 | } else { |
---|
346 | // test down |
---|
347 | if (oldValue - thetaDown * alpha < -dualTolerance_) { |
---|
348 | thetaDown = (dualTolerance_ + oldValue) / alpha; |
---|
349 | sequenceDown = iSequence2; |
---|
350 | alphaDown = alpha; |
---|
351 | } |
---|
352 | } |
---|
353 | break; |
---|
354 | } |
---|
355 | } |
---|
356 | } |
---|
357 | if (sequenceUp >= 0) { |
---|
358 | costIncrease = thetaUp; |
---|
359 | sequenceIncrease = sequenceUp; |
---|
360 | alphaIncrease = alphaUp; |
---|
361 | } |
---|
362 | if (sequenceDown >= 0) { |
---|
363 | costDecrease = thetaDown; |
---|
364 | sequenceDecrease = sequenceDown; |
---|
365 | alphaDecrease = alphaDown; |
---|
366 | } |
---|
367 | } |
---|
368 | /** Primal ranging. |
---|
369 | This computes increase/decrease in value for each given variable and corresponding |
---|
370 | sequence numbers which would change basis. Sequence numbers are 0..numberColumns |
---|
371 | and numberColumns.. for artificials/slacks. |
---|
372 | For basic variables the sequence number will be that of the basic variables. |
---|
373 | |
---|
374 | Up to user to provide correct length arrays. |
---|
375 | |
---|
376 | When here - guaranteed optimal |
---|
377 | */ |
---|
378 | void |
---|
379 | ClpSimplexOther::primalRanging(int numberCheck, const int * which, |
---|
380 | double * valueIncreased, int * sequenceIncreased, |
---|
381 | double * valueDecreased, int * sequenceDecreased) |
---|
382 | { |
---|
383 | rowArray_[0]->clear(); |
---|
384 | rowArray_[1]->clear(); |
---|
385 | lowerIn_ = -COIN_DBL_MAX; |
---|
386 | upperIn_ = COIN_DBL_MAX; |
---|
387 | valueIn_ = 0.0; |
---|
388 | for ( int i = 0; i < numberCheck; i++) { |
---|
389 | int iSequence = which[i]; |
---|
390 | double valueIncrease = COIN_DBL_MAX; |
---|
391 | double valueDecrease = COIN_DBL_MAX; |
---|
392 | int sequenceIncrease = -1; |
---|
393 | int sequenceDecrease = -1; |
---|
394 | |
---|
395 | switch(getStatus(iSequence)) { |
---|
396 | |
---|
397 | case basic: |
---|
398 | case isFree: |
---|
399 | case superBasic: |
---|
400 | // Easy |
---|
401 | valueDecrease = CoinMax(0.0, upper_[iSequence] - solution_[iSequence]); |
---|
402 | valueIncrease = CoinMax(0.0, solution_[iSequence] - lower_[iSequence]); |
---|
403 | sequenceDecrease = iSequence; |
---|
404 | sequenceIncrease = iSequence; |
---|
405 | break; |
---|
406 | case isFixed: |
---|
407 | case atUpperBound: |
---|
408 | case atLowerBound: { |
---|
409 | // Non trivial |
---|
410 | // Other bound is ignored |
---|
411 | #ifndef COIN_FAC_NEW |
---|
412 | unpackPacked(rowArray_[1], iSequence); |
---|
413 | #else |
---|
414 | unpack(rowArray_[1], iSequence); |
---|
415 | #endif |
---|
416 | factorization_->updateColumn(rowArray_[2], rowArray_[1]); |
---|
417 | // Get extra rows |
---|
418 | matrix_->extendUpdated(this, rowArray_[1], 0); |
---|
419 | // do ratio test |
---|
420 | checkPrimalRatios(rowArray_[1], 1); |
---|
421 | if (pivotRow_ >= 0) { |
---|
422 | valueIncrease = theta_; |
---|
423 | sequenceIncrease = pivotVariable_[pivotRow_]; |
---|
424 | } |
---|
425 | checkPrimalRatios(rowArray_[1], -1); |
---|
426 | if (pivotRow_ >= 0) { |
---|
427 | valueDecrease = theta_; |
---|
428 | sequenceDecrease = pivotVariable_[pivotRow_]; |
---|
429 | } |
---|
430 | rowArray_[1]->clear(); |
---|
431 | } |
---|
432 | break; |
---|
433 | } |
---|
434 | double scaleFactor; |
---|
435 | if (rowScale_) { |
---|
436 | if (iSequence < numberColumns_) |
---|
437 | scaleFactor = columnScale_[iSequence] / rhsScale_; |
---|
438 | else |
---|
439 | scaleFactor = 1.0 / (rowScale_[iSequence-numberColumns_] * rhsScale_); |
---|
440 | } else { |
---|
441 | scaleFactor = 1.0 / rhsScale_; |
---|
442 | } |
---|
443 | if (valueIncrease < 1.0e30) |
---|
444 | valueIncrease *= scaleFactor; |
---|
445 | else |
---|
446 | valueIncrease = COIN_DBL_MAX; |
---|
447 | if (valueDecrease < 1.0e30) |
---|
448 | valueDecrease *= scaleFactor; |
---|
449 | else |
---|
450 | valueDecrease = COIN_DBL_MAX; |
---|
451 | valueIncreased[i] = valueIncrease; |
---|
452 | sequenceIncreased[i] = sequenceIncrease; |
---|
453 | valueDecreased[i] = valueDecrease; |
---|
454 | sequenceDecreased[i] = sequenceDecrease; |
---|
455 | } |
---|
456 | } |
---|
457 | // Returns new value of whichOther when whichIn enters basis |
---|
458 | double |
---|
459 | ClpSimplexOther::primalRanging1(int whichIn, int whichOther) |
---|
460 | { |
---|
461 | rowArray_[0]->clear(); |
---|
462 | rowArray_[1]->clear(); |
---|
463 | int iSequence = whichIn; |
---|
464 | double newValue = solution_[whichOther]; |
---|
465 | double alphaOther = 0.0; |
---|
466 | Status status = getStatus(iSequence); |
---|
467 | assert (status == atLowerBound || status == atUpperBound); |
---|
468 | int wayIn = (status == atLowerBound) ? 1 : -1; |
---|
469 | |
---|
470 | switch(getStatus(iSequence)) { |
---|
471 | |
---|
472 | case basic: |
---|
473 | case isFree: |
---|
474 | case superBasic: |
---|
475 | assert (whichIn == whichOther); |
---|
476 | // Easy |
---|
477 | newValue = wayIn > 0 ? upper_[iSequence] : lower_[iSequence]; |
---|
478 | break; |
---|
479 | case isFixed: |
---|
480 | case atUpperBound: |
---|
481 | case atLowerBound: |
---|
482 | // Non trivial |
---|
483 | { |
---|
484 | // Other bound is ignored |
---|
485 | #ifndef COIN_FAC_NEW |
---|
486 | unpackPacked(rowArray_[1], iSequence); |
---|
487 | #else |
---|
488 | unpack(rowArray_[1], iSequence); |
---|
489 | #endif |
---|
490 | factorization_->updateColumn(rowArray_[2], rowArray_[1]); |
---|
491 | // Get extra rows |
---|
492 | matrix_->extendUpdated(this, rowArray_[1], 0); |
---|
493 | // do ratio test |
---|
494 | double acceptablePivot = 1.0e-7; |
---|
495 | double * work = rowArray_[1]->denseVector(); |
---|
496 | int number = rowArray_[1]->getNumElements(); |
---|
497 | int * which = rowArray_[1]->getIndices(); |
---|
498 | |
---|
499 | // we may need to swap sign |
---|
500 | double way = wayIn; |
---|
501 | double theta = 1.0e30; |
---|
502 | for (int iIndex = 0; iIndex < number; iIndex++) { |
---|
503 | |
---|
504 | int iRow = which[iIndex]; |
---|
505 | #ifndef COIN_FAC_NEW |
---|
506 | double alpha = work[iIndex] * way; |
---|
507 | #else |
---|
508 | double alpha = work[iRow] * way; |
---|
509 | #endif |
---|
510 | int iPivot = pivotVariable_[iRow]; |
---|
511 | if (iPivot == whichOther) { |
---|
512 | alphaOther = alpha; |
---|
513 | continue; |
---|
514 | } |
---|
515 | double oldValue = solution_[iPivot]; |
---|
516 | if (fabs(alpha) > acceptablePivot) { |
---|
517 | if (alpha > 0.0) { |
---|
518 | // basic variable going towards lower bound |
---|
519 | double bound = lower_[iPivot]; |
---|
520 | oldValue -= bound; |
---|
521 | if (oldValue - theta * alpha < 0.0) { |
---|
522 | theta = CoinMax(0.0, oldValue / alpha); |
---|
523 | } |
---|
524 | } else { |
---|
525 | // basic variable going towards upper bound |
---|
526 | double bound = upper_[iPivot]; |
---|
527 | oldValue = oldValue - bound; |
---|
528 | if (oldValue - theta * alpha > 0.0) { |
---|
529 | theta = CoinMax(0.0, oldValue / alpha); |
---|
530 | } |
---|
531 | } |
---|
532 | } |
---|
533 | } |
---|
534 | if (whichIn != whichOther) { |
---|
535 | if (theta < 1.0e30) |
---|
536 | newValue -= theta * alphaOther; |
---|
537 | else |
---|
538 | newValue = alphaOther > 0.0 ? -1.0e30 : 1.0e30; |
---|
539 | } else { |
---|
540 | newValue += theta * wayIn; |
---|
541 | } |
---|
542 | } |
---|
543 | rowArray_[1]->clear(); |
---|
544 | break; |
---|
545 | } |
---|
546 | double scaleFactor; |
---|
547 | if (rowScale_) { |
---|
548 | if (whichOther < numberColumns_) |
---|
549 | scaleFactor = columnScale_[whichOther] / rhsScale_; |
---|
550 | else |
---|
551 | scaleFactor = 1.0 / (rowScale_[whichOther-numberColumns_] * rhsScale_); |
---|
552 | } else { |
---|
553 | scaleFactor = 1.0 / rhsScale_; |
---|
554 | } |
---|
555 | if (newValue < 1.0e29) |
---|
556 | if (newValue > -1.0e29) |
---|
557 | newValue *= scaleFactor; |
---|
558 | else |
---|
559 | newValue = -COIN_DBL_MAX; |
---|
560 | else |
---|
561 | newValue = COIN_DBL_MAX; |
---|
562 | return newValue; |
---|
563 | } |
---|
564 | /* |
---|
565 | Row array has pivot column |
---|
566 | This is used in primal ranging |
---|
567 | */ |
---|
568 | void |
---|
569 | ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray, |
---|
570 | int direction) |
---|
571 | { |
---|
572 | // sequence stays as row number until end |
---|
573 | pivotRow_ = -1; |
---|
574 | double acceptablePivot = 1.0e-7; |
---|
575 | double * work = rowArray->denseVector(); |
---|
576 | int number = rowArray->getNumElements(); |
---|
577 | int * which = rowArray->getIndices(); |
---|
578 | |
---|
579 | // we need to swap sign if going down |
---|
580 | double way = direction; |
---|
581 | theta_ = 1.0e30; |
---|
582 | for (int iIndex = 0; iIndex < number; iIndex++) { |
---|
583 | |
---|
584 | int iRow = which[iIndex]; |
---|
585 | #ifndef COIN_FAC_NEW |
---|
586 | double alpha = work[iIndex] * way; |
---|
587 | #else |
---|
588 | double alpha = work[iRow] * way; |
---|
589 | #endif |
---|
590 | int iPivot = pivotVariable_[iRow]; |
---|
591 | double oldValue = solution_[iPivot]; |
---|
592 | if (fabs(alpha) > acceptablePivot) { |
---|
593 | if (alpha > 0.0) { |
---|
594 | // basic variable going towards lower bound |
---|
595 | double bound = lower_[iPivot]; |
---|
596 | oldValue -= bound; |
---|
597 | if (oldValue - theta_ * alpha < 0.0) { |
---|
598 | pivotRow_ = iRow; |
---|
599 | theta_ = CoinMax(0.0, oldValue / alpha); |
---|
600 | } |
---|
601 | } else { |
---|
602 | // basic variable going towards upper bound |
---|
603 | double bound = upper_[iPivot]; |
---|
604 | oldValue = oldValue - bound; |
---|
605 | if (oldValue - theta_ * alpha > 0.0) { |
---|
606 | pivotRow_ = iRow; |
---|
607 | theta_ = CoinMax(0.0, oldValue / alpha); |
---|
608 | } |
---|
609 | } |
---|
610 | } |
---|
611 | } |
---|
612 | } |
---|
613 | /* Write the basis in MPS format to the specified file. |
---|
614 | If writeValues true writes values of structurals |
---|
615 | (and adds VALUES to end of NAME card) |
---|
616 | |
---|
617 | Row and column names may be null. |
---|
618 | formatType is |
---|
619 | <ul> |
---|
620 | <li> 0 - normal |
---|
621 | <li> 1 - extra accuracy |
---|
622 | <li> 2 - IEEE hex (later) |
---|
623 | </ul> |
---|
624 | |
---|
625 | Returns non-zero on I/O error |
---|
626 | |
---|
627 | This is based on code contributed by Thorsten Koch |
---|
628 | */ |
---|
629 | int |
---|
630 | ClpSimplexOther::writeBasis(const char *filename, |
---|
631 | bool writeValues, |
---|
632 | int formatType) const |
---|
633 | { |
---|
634 | formatType = CoinMax(0, formatType); |
---|
635 | formatType = CoinMin(2, formatType); |
---|
636 | if (!writeValues) |
---|
637 | formatType = 0; |
---|
638 | // See if INTEL if IEEE |
---|
639 | if (formatType == 2) { |
---|
640 | // test intel here and add 1 if not intel |
---|
641 | double value = 1.0; |
---|
642 | char x[8]; |
---|
643 | memcpy(x, &value, 8); |
---|
644 | if (x[0] == 63) { |
---|
645 | formatType ++; // not intel |
---|
646 | } else { |
---|
647 | assert (x[0] == 0); |
---|
648 | } |
---|
649 | } |
---|
650 | |
---|
651 | char number[20]; |
---|
652 | FILE * fp = fopen(filename, "w"); |
---|
653 | if (!fp) |
---|
654 | return -1; |
---|
655 | |
---|
656 | // NAME card |
---|
657 | |
---|
658 | if (strcmp(strParam_[ClpProbName].c_str(), "") == 0) { |
---|
659 | fprintf(fp, "NAME BLANK "); |
---|
660 | } else { |
---|
661 | fprintf(fp, "NAME %s ", strParam_[ClpProbName].c_str()); |
---|
662 | } |
---|
663 | if (formatType >= 2) |
---|
664 | fprintf(fp, "FREEIEEE"); |
---|
665 | else if (writeValues) |
---|
666 | fprintf(fp, "VALUES"); |
---|
667 | // finish off name |
---|
668 | fprintf(fp, "\n"); |
---|
669 | int iRow = 0; |
---|
670 | for(int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
671 | bool printit = false; |
---|
672 | if( getColumnStatus(iColumn) == ClpSimplex::basic) { |
---|
673 | printit = true; |
---|
674 | // Find non basic row |
---|
675 | for(; iRow < numberRows_; iRow++) { |
---|
676 | if (getRowStatus(iRow) != ClpSimplex::basic) |
---|
677 | break; |
---|
678 | } |
---|
679 | if (lengthNames_) { |
---|
680 | if (iRow != numberRows_) { |
---|
681 | fprintf(fp, " %s %-8s %s", |
---|
682 | getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL", |
---|
683 | columnNames_[iColumn].c_str(), |
---|
684 | rowNames_[iRow].c_str()); |
---|
685 | iRow++; |
---|
686 | } else { |
---|
687 | // Allow for too many basics! |
---|
688 | fprintf(fp, " BS %-8s ", |
---|
689 | columnNames_[iColumn].c_str()); |
---|
690 | // Dummy row name if values |
---|
691 | if (writeValues) |
---|
692 | fprintf(fp, " _dummy_"); |
---|
693 | } |
---|
694 | } else { |
---|
695 | // no names |
---|
696 | if (iRow != numberRows_) { |
---|
697 | fprintf(fp, " %s C%7.7d R%7.7d", |
---|
698 | getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL", |
---|
699 | iColumn, iRow); |
---|
700 | iRow++; |
---|
701 | } else { |
---|
702 | // Allow for too many basics! |
---|
703 | fprintf(fp, " BS C%7.7d", iColumn); |
---|
704 | // Dummy row name if values |
---|
705 | if (writeValues) |
---|
706 | fprintf(fp, " _dummy_"); |
---|
707 | } |
---|
708 | } |
---|
709 | } else { |
---|
710 | if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) { |
---|
711 | printit = true; |
---|
712 | if (lengthNames_) |
---|
713 | fprintf(fp, " UL %s", columnNames_[iColumn].c_str()); |
---|
714 | else |
---|
715 | fprintf(fp, " UL C%7.7d", iColumn); |
---|
716 | // Dummy row name if values |
---|
717 | if (writeValues) |
---|
718 | fprintf(fp, " _dummy_"); |
---|
719 | } else if( (getColumnStatus(iColumn) == ClpSimplex::superBasic|| |
---|
720 | getColumnStatus(iColumn) == ClpSimplex::isFree)&& |
---|
721 | writeValues) { |
---|
722 | printit = true; |
---|
723 | if (lengthNames_) |
---|
724 | fprintf(fp, " BS %s", columnNames_[iColumn].c_str()); |
---|
725 | else |
---|
726 | fprintf(fp, " BS C%7.7d", iColumn); |
---|
727 | // Dummy row name if values |
---|
728 | if (writeValues) |
---|
729 | fprintf(fp, " _dummy_"); |
---|
730 | } |
---|
731 | } |
---|
732 | if (printit && writeValues) { |
---|
733 | // add value |
---|
734 | CoinConvertDouble(0, formatType, columnActivity_[iColumn], number); |
---|
735 | fprintf(fp, " %s", number); |
---|
736 | } |
---|
737 | if (printit) |
---|
738 | fprintf(fp, "\n"); |
---|
739 | } |
---|
740 | fprintf(fp, "ENDATA\n"); |
---|
741 | fclose(fp); |
---|
742 | return 0; |
---|
743 | } |
---|
744 | // Read a basis from the given filename |
---|
745 | int |
---|
746 | ClpSimplexOther::readBasis(const char *fileName) |
---|
747 | { |
---|
748 | int status = 0; |
---|
749 | if (strcmp(fileName, "-") != 0 && strcmp(fileName, "stdin") != 0) { |
---|
750 | FILE *fp = fopen(fileName, "r"); |
---|
751 | if (fp) { |
---|
752 | // can open - lets go for it |
---|
753 | fclose(fp); |
---|
754 | } else { |
---|
755 | handler_->message(CLP_UNABLE_OPEN, messages_) |
---|
756 | << fileName << CoinMessageEol; |
---|
757 | return -1; |
---|
758 | } |
---|
759 | } |
---|
760 | CoinMpsIO m; |
---|
761 | m.passInMessageHandler(handler_); |
---|
762 | *m.messagesPointer() = coinMessages(); |
---|
763 | bool savePrefix = m.messageHandler()->prefix(); |
---|
764 | m.messageHandler()->setPrefix(handler_->prefix()); |
---|
765 | status = m.readBasis(fileName, "", columnActivity_, status_ + numberColumns_, |
---|
766 | status_, |
---|
767 | columnNames_, numberColumns_, |
---|
768 | rowNames_, numberRows_); |
---|
769 | m.messageHandler()->setPrefix(savePrefix); |
---|
770 | if (status >= 0) { |
---|
771 | if (!status) { |
---|
772 | // set values |
---|
773 | int iColumn, iRow; |
---|
774 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
775 | if (getRowStatus(iRow) == atLowerBound) |
---|
776 | rowActivity_[iRow] = rowLower_[iRow]; |
---|
777 | else if (getRowStatus(iRow) == atUpperBound) |
---|
778 | rowActivity_[iRow] = rowUpper_[iRow]; |
---|
779 | } |
---|
780 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
781 | if (getColumnStatus(iColumn) == atLowerBound) |
---|
782 | columnActivity_[iColumn] = columnLower_[iColumn]; |
---|
783 | else if (getColumnStatus(iColumn) == atUpperBound) |
---|
784 | columnActivity_[iColumn] = columnUpper_[iColumn]; |
---|
785 | } |
---|
786 | } else { |
---|
787 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
788 | matrix_->times(-1.0, columnActivity_, rowActivity_); |
---|
789 | } |
---|
790 | } else { |
---|
791 | // errors |
---|
792 | handler_->message(CLP_IMPORT_ERRORS, messages_) |
---|
793 | << status << fileName << CoinMessageEol; |
---|
794 | } |
---|
795 | return status; |
---|
796 | } |
---|
797 | /* Creates dual of a problem if looks plausible |
---|
798 | (defaults will always create model) |
---|
799 | fractionRowRanges is fraction of rows allowed to have ranges |
---|
800 | fractionColumnRanges is fraction of columns allowed to have ranges |
---|
801 | */ |
---|
802 | ClpSimplex * |
---|
803 | ClpSimplexOther::dualOfModel(double fractionRowRanges, double fractionColumnRanges) const |
---|
804 | { |
---|
805 | const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this); |
---|
806 | bool changed = false; |
---|
807 | int numberChanged = 0; |
---|
808 | int numberFreeColumnsInPrimal=0; |
---|
809 | int iColumn; |
---|
810 | // check if we need to change bounds to rows |
---|
811 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
812 | if (columnUpper_[iColumn] < 1.0e20) { |
---|
813 | if (columnLower_[iColumn] > -1.0e20) { |
---|
814 | changed = true; |
---|
815 | numberChanged++; |
---|
816 | } |
---|
817 | } else if (columnLower_[iColumn] < -1.0e20) { |
---|
818 | numberFreeColumnsInPrimal++; |
---|
819 | } |
---|
820 | } |
---|
821 | int iRow; |
---|
822 | int numberExtraRows = 0; |
---|
823 | int numberFreeColumnsInDual=0; |
---|
824 | if (numberChanged <= fractionColumnRanges * numberColumns_) { |
---|
825 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
826 | if (rowLower_[iRow] > -1.0e20 && |
---|
827 | rowUpper_[iRow] < 1.0e20) { |
---|
828 | if (rowUpper_[iRow] != rowLower_[iRow]) |
---|
829 | numberExtraRows++; |
---|
830 | else |
---|
831 | numberFreeColumnsInDual++; |
---|
832 | } |
---|
833 | } |
---|
834 | if (numberExtraRows > fractionRowRanges * numberRows_) |
---|
835 | return NULL; |
---|
836 | } else { |
---|
837 | return NULL; |
---|
838 | } |
---|
839 | printf("would have %d free columns in primal, %d in dual\n", |
---|
840 | numberFreeColumnsInPrimal,numberFreeColumnsInDual); |
---|
841 | if (4*(numberFreeColumnsInDual-numberFreeColumnsInPrimal)> |
---|
842 | numberColumns_&&fractionRowRanges<1.0) |
---|
843 | return NULL; //dangerous (well anyway in dual) |
---|
844 | if (changed) { |
---|
845 | ClpSimplex * model3 = new ClpSimplex(*model2); |
---|
846 | CoinBuild build; |
---|
847 | double one = 1.0; |
---|
848 | int numberColumns = model3->numberColumns(); |
---|
849 | const double * columnLower = model3->columnLower(); |
---|
850 | const double * columnUpper = model3->columnUpper(); |
---|
851 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
852 | if (columnUpper[iColumn] < 1.0e20 && |
---|
853 | columnLower[iColumn] > -1.0e20) { |
---|
854 | if (fabs(columnLower[iColumn]) < fabs(columnUpper[iColumn])) { |
---|
855 | double value = columnUpper[iColumn]; |
---|
856 | model3->setColumnUpper(iColumn, COIN_DBL_MAX); |
---|
857 | build.addRow(1, &iColumn, &one, -COIN_DBL_MAX, value); |
---|
858 | } else { |
---|
859 | double value = columnLower[iColumn]; |
---|
860 | model3->setColumnLower(iColumn, -COIN_DBL_MAX); |
---|
861 | build.addRow(1, &iColumn, &one, value, COIN_DBL_MAX); |
---|
862 | } |
---|
863 | } |
---|
864 | } |
---|
865 | model3->addRows(build); |
---|
866 | model2 = model3; |
---|
867 | } |
---|
868 | int numberColumns = model2->numberColumns(); |
---|
869 | const double * columnLower = model2->columnLower(); |
---|
870 | const double * columnUpper = model2->columnUpper(); |
---|
871 | int numberRows = model2->numberRows(); |
---|
872 | double * rowLower = CoinCopyOfArray(model2->rowLower(), numberRows); |
---|
873 | double * rowUpper = CoinCopyOfArray(model2->rowUpper(), numberRows); |
---|
874 | |
---|
875 | const double * objective = model2->objective(); |
---|
876 | CoinPackedMatrix * matrix = model2->matrix(); |
---|
877 | // get transpose |
---|
878 | CoinPackedMatrix rowCopy = *matrix; |
---|
879 | const int * row = matrix->getIndices(); |
---|
880 | const int * columnLength = matrix->getVectorLengths(); |
---|
881 | const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
---|
882 | const double * elementByColumn = matrix->getElements(); |
---|
883 | double objOffset = 0.0; |
---|
884 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
885 | double offset = 0.0; |
---|
886 | double objValue = optimizationDirection_ * objective[iColumn]; |
---|
887 | if (columnUpper[iColumn] > 1.0e20) { |
---|
888 | if (columnLower[iColumn] > -1.0e20) |
---|
889 | offset = columnLower[iColumn]; |
---|
890 | } else if (columnLower[iColumn] < -1.0e20) { |
---|
891 | offset = columnUpper[iColumn]; |
---|
892 | } else { |
---|
893 | // taken care of before |
---|
894 | abort(); |
---|
895 | } |
---|
896 | if (offset) { |
---|
897 | objOffset += offset * objValue; |
---|
898 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
899 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
900 | int iRow = row[j]; |
---|
901 | if (rowLower[iRow] > -1.0e20) |
---|
902 | rowLower[iRow] -= offset * elementByColumn[j]; |
---|
903 | if (rowUpper[iRow] < 1.0e20) |
---|
904 | rowUpper[iRow] -= offset * elementByColumn[j]; |
---|
905 | } |
---|
906 | } |
---|
907 | } |
---|
908 | int * which = new int[numberRows+numberExtraRows]; |
---|
909 | rowCopy.reverseOrdering(); |
---|
910 | rowCopy.transpose(); |
---|
911 | double * fromRowsLower = new double[numberRows+numberExtraRows]; |
---|
912 | double * fromRowsUpper = new double[numberRows+numberExtraRows]; |
---|
913 | double * newObjective = new double[numberRows+numberExtraRows]; |
---|
914 | double * fromColumnsLower = new double[numberColumns]; |
---|
915 | double * fromColumnsUpper = new double[numberColumns]; |
---|
916 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
917 | double objValue = optimizationDirection_ * objective[iColumn]; |
---|
918 | // Offset is already in |
---|
919 | if (columnUpper[iColumn] > 1.0e20) { |
---|
920 | if (columnLower[iColumn] > -1.0e20) { |
---|
921 | fromColumnsLower[iColumn] = -COIN_DBL_MAX; |
---|
922 | fromColumnsUpper[iColumn] = objValue; |
---|
923 | } else { |
---|
924 | // free |
---|
925 | fromColumnsLower[iColumn] = objValue; |
---|
926 | fromColumnsUpper[iColumn] = objValue; |
---|
927 | } |
---|
928 | } else if (columnLower[iColumn] < -1.0e20) { |
---|
929 | fromColumnsLower[iColumn] = objValue; |
---|
930 | fromColumnsUpper[iColumn] = COIN_DBL_MAX; |
---|
931 | } else { |
---|
932 | abort(); |
---|
933 | } |
---|
934 | } |
---|
935 | int kRow = 0; |
---|
936 | int kExtraRow = numberRows; |
---|
937 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
938 | if (rowLower[iRow] < -1.0e20) { |
---|
939 | assert (rowUpper[iRow] < 1.0e20); |
---|
940 | newObjective[kRow] = -rowUpper[iRow]; |
---|
941 | fromRowsLower[kRow] = -COIN_DBL_MAX; |
---|
942 | fromRowsUpper[kRow] = 0.0; |
---|
943 | which[kRow] = iRow; |
---|
944 | kRow++; |
---|
945 | } else if (rowUpper[iRow] > 1.0e20) { |
---|
946 | newObjective[kRow] = -rowLower[iRow]; |
---|
947 | fromRowsLower[kRow] = 0.0; |
---|
948 | fromRowsUpper[kRow] = COIN_DBL_MAX; |
---|
949 | which[kRow] = iRow; |
---|
950 | kRow++; |
---|
951 | } else { |
---|
952 | if (rowUpper[iRow] == rowLower[iRow]) { |
---|
953 | newObjective[kRow] = -rowLower[iRow]; |
---|
954 | fromRowsLower[kRow] = -COIN_DBL_MAX;; |
---|
955 | fromRowsUpper[kRow] = COIN_DBL_MAX; |
---|
956 | which[kRow] = iRow; |
---|
957 | kRow++; |
---|
958 | } else { |
---|
959 | // range |
---|
960 | newObjective[kRow] = -rowUpper[iRow]; |
---|
961 | fromRowsLower[kRow] = -COIN_DBL_MAX; |
---|
962 | fromRowsUpper[kRow] = 0.0; |
---|
963 | which[kRow] = iRow; |
---|
964 | kRow++; |
---|
965 | newObjective[kExtraRow] = -rowLower[iRow]; |
---|
966 | fromRowsLower[kExtraRow] = 0.0; |
---|
967 | fromRowsUpper[kExtraRow] = COIN_DBL_MAX; |
---|
968 | which[kExtraRow] = iRow; |
---|
969 | kExtraRow++; |
---|
970 | } |
---|
971 | } |
---|
972 | } |
---|
973 | if (numberExtraRows) { |
---|
974 | CoinPackedMatrix newCopy; |
---|
975 | newCopy.setExtraGap(0.0); |
---|
976 | newCopy.setExtraMajor(0.0); |
---|
977 | newCopy.submatrixOfWithDuplicates(rowCopy, kExtraRow, which); |
---|
978 | rowCopy = newCopy; |
---|
979 | } |
---|
980 | ClpSimplex * modelDual = new ClpSimplex(); |
---|
981 | modelDual->loadProblem(rowCopy, fromRowsLower, fromRowsUpper, newObjective, |
---|
982 | fromColumnsLower, fromColumnsUpper); |
---|
983 | modelDual->setObjectiveOffset(objOffset); |
---|
984 | modelDual->setDualBound(model2->dualBound()); |
---|
985 | modelDual->setInfeasibilityCost(model2->infeasibilityCost()); |
---|
986 | modelDual->setDualTolerance(model2->dualTolerance()); |
---|
987 | modelDual->setPrimalTolerance(model2->primalTolerance()); |
---|
988 | modelDual->setPerturbation(model2->perturbation()); |
---|
989 | modelDual->setSpecialOptions(model2->specialOptions()); |
---|
990 | modelDual->setMoreSpecialOptions(model2->moreSpecialOptions()); |
---|
991 | modelDual->setMaximumIterations(model2->maximumIterations()); |
---|
992 | modelDual->setFactorizationFrequency(model2->factorizationFrequency()); |
---|
993 | modelDual->setLogLevel(model2->logLevel()); |
---|
994 | delete [] fromRowsLower; |
---|
995 | delete [] fromRowsUpper; |
---|
996 | delete [] fromColumnsLower; |
---|
997 | delete [] fromColumnsUpper; |
---|
998 | delete [] newObjective; |
---|
999 | delete [] which; |
---|
1000 | delete [] rowLower; |
---|
1001 | delete [] rowUpper; |
---|
1002 | if (changed) |
---|
1003 | delete model2; |
---|
1004 | modelDual->createStatus(); |
---|
1005 | return modelDual; |
---|
1006 | } |
---|
1007 | // Restores solution from dualized problem |
---|
1008 | int |
---|
1009 | ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem, |
---|
1010 | bool checkAccuracy) |
---|
1011 | { |
---|
1012 | int returnCode = 0;; |
---|
1013 | createStatus(); |
---|
1014 | // Number of rows in dual problem was original number of columns |
---|
1015 | assert (numberColumns_ == dualProblem->numberRows()); |
---|
1016 | // If slack on d-row basic then column at bound otherwise column basic |
---|
1017 | // If d-column basic then rhs tight |
---|
1018 | int numberBasic = 0; |
---|
1019 | int iRow, iColumn = 0; |
---|
1020 | // Get number of extra rows from ranges |
---|
1021 | int numberExtraRows = 0; |
---|
1022 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1023 | if (rowLower_[iRow] > -1.0e20 && |
---|
1024 | rowUpper_[iRow] < 1.0e20) { |
---|
1025 | if (rowUpper_[iRow] != rowLower_[iRow]) |
---|
1026 | numberExtraRows++; |
---|
1027 | } |
---|
1028 | } |
---|
1029 | const double * objective = this->objective(); |
---|
1030 | const double * dualDual = dualProblem->dualRowSolution(); |
---|
1031 | const double * dualDj = dualProblem->dualColumnSolution(); |
---|
1032 | const double * dualSol = dualProblem->primalColumnSolution(); |
---|
1033 | const double * dualActs = dualProblem->primalRowSolution(); |
---|
1034 | #if 0 |
---|
1035 | ClpSimplex thisCopy = *this; |
---|
1036 | thisCopy.dual(); // for testing |
---|
1037 | const double * primalDual = thisCopy.dualRowSolution(); |
---|
1038 | const double * primalDj = thisCopy.dualColumnSolution(); |
---|
1039 | const double * primalSol = thisCopy.primalColumnSolution(); |
---|
1040 | const double * primalActs = thisCopy.primalRowSolution(); |
---|
1041 | char ss[] = {'F', 'B', 'U', 'L', 'S', 'F'}; |
---|
1042 | printf ("Dual problem row info %d rows\n", dualProblem->numberRows()); |
---|
1043 | for (iRow = 0; iRow < dualProblem->numberRows(); iRow++) |
---|
1044 | printf("%d at %c primal %g dual %g\n", |
---|
1045 | iRow, ss[dualProblem->getRowStatus(iRow)], |
---|
1046 | dualActs[iRow], dualDual[iRow]); |
---|
1047 | printf ("Dual problem column info %d columns\n", dualProblem->numberColumns()); |
---|
1048 | for (iColumn = 0; iColumn < dualProblem->numberColumns(); iColumn++) |
---|
1049 | printf("%d at %c primal %g dual %g\n", |
---|
1050 | iColumn, ss[dualProblem->getColumnStatus(iColumn)], |
---|
1051 | dualSol[iColumn], dualDj[iColumn]); |
---|
1052 | printf ("Primal problem row info %d rows\n", thisCopy.numberRows()); |
---|
1053 | for (iRow = 0; iRow < thisCopy.numberRows(); iRow++) |
---|
1054 | printf("%d at %c primal %g dual %g\n", |
---|
1055 | iRow, ss[thisCopy.getRowStatus(iRow)], |
---|
1056 | primalActs[iRow], primalDual[iRow]); |
---|
1057 | printf ("Primal problem column info %d columns\n", thisCopy.numberColumns()); |
---|
1058 | for (iColumn = 0; iColumn < thisCopy.numberColumns(); iColumn++) |
---|
1059 | printf("%d at %c primal %g dual %g\n", |
---|
1060 | iColumn, ss[thisCopy.getColumnStatus(iColumn)], |
---|
1061 | primalSol[iColumn], primalDj[iColumn]); |
---|
1062 | #endif |
---|
1063 | // position at bound information |
---|
1064 | int jColumn = numberRows_; |
---|
1065 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1066 | double objValue = optimizationDirection_ * objective[iColumn]; |
---|
1067 | Status status = dualProblem->getRowStatus(iColumn); |
---|
1068 | double otherValue = COIN_DBL_MAX; |
---|
1069 | if (columnUpper_[iColumn] < 1.0e20 && |
---|
1070 | columnLower_[iColumn] > -1.0e20) { |
---|
1071 | if (fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])) { |
---|
1072 | otherValue = columnUpper_[iColumn] + dualDj[jColumn]; |
---|
1073 | } else { |
---|
1074 | otherValue = columnLower_[iColumn] + dualDj[jColumn]; |
---|
1075 | } |
---|
1076 | jColumn++; |
---|
1077 | } |
---|
1078 | if (status == basic) { |
---|
1079 | // column is at bound |
---|
1080 | if (otherValue == COIN_DBL_MAX) { |
---|
1081 | reducedCost_[iColumn] = objValue - dualActs[iColumn]; |
---|
1082 | if (columnUpper_[iColumn] > 1.0e20) { |
---|
1083 | if (columnLower_[iColumn] > -1.0e20) { |
---|
1084 | if (columnUpper_[iColumn] > columnLower_[iColumn]) |
---|
1085 | setColumnStatus(iColumn, atLowerBound); |
---|
1086 | else |
---|
1087 | setColumnStatus(iColumn, isFixed); |
---|
1088 | columnActivity_[iColumn] = columnLower_[iColumn]; |
---|
1089 | } else { |
---|
1090 | // free |
---|
1091 | setColumnStatus(iColumn, isFree); |
---|
1092 | columnActivity_[iColumn] = 0.0; |
---|
1093 | } |
---|
1094 | } else { |
---|
1095 | setColumnStatus(iColumn, atUpperBound); |
---|
1096 | columnActivity_[iColumn] = columnUpper_[iColumn]; |
---|
1097 | } |
---|
1098 | } else { |
---|
1099 | reducedCost_[iColumn] = objValue - dualActs[iColumn]; |
---|
1100 | //printf("other dual sol %g\n",otherValue); |
---|
1101 | if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) { |
---|
1102 | if (columnUpper_[iColumn] > columnLower_[iColumn]) |
---|
1103 | setColumnStatus(iColumn, atLowerBound); |
---|
1104 | else |
---|
1105 | setColumnStatus(iColumn, isFixed); |
---|
1106 | columnActivity_[iColumn] = columnLower_[iColumn]; |
---|
1107 | } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) { |
---|
1108 | if (columnUpper_[iColumn] > columnLower_[iColumn]) |
---|
1109 | setColumnStatus(iColumn, atUpperBound); |
---|
1110 | else |
---|
1111 | setColumnStatus(iColumn, isFixed); |
---|
1112 | columnActivity_[iColumn] = columnUpper_[iColumn]; |
---|
1113 | } else { |
---|
1114 | setColumnStatus(iColumn, superBasic); |
---|
1115 | columnActivity_[iColumn] = otherValue; |
---|
1116 | } |
---|
1117 | } |
---|
1118 | } else { |
---|
1119 | if (otherValue == COIN_DBL_MAX) { |
---|
1120 | // column basic |
---|
1121 | setColumnStatus(iColumn, basic); |
---|
1122 | numberBasic++; |
---|
1123 | if (columnLower_[iColumn] > -1.0e20) { |
---|
1124 | columnActivity_[iColumn] = -dualDual[iColumn] + columnLower_[iColumn]; |
---|
1125 | } else if (columnUpper_[iColumn] < 1.0e20) { |
---|
1126 | columnActivity_[iColumn] = -dualDual[iColumn] + columnUpper_[iColumn]; |
---|
1127 | } else { |
---|
1128 | columnActivity_[iColumn] = -dualDual[iColumn]; |
---|
1129 | } |
---|
1130 | reducedCost_[iColumn] = 0.0; |
---|
1131 | } else { |
---|
1132 | // may be at other bound |
---|
1133 | //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1); |
---|
1134 | if (dualProblem->getColumnStatus(jColumn - 1) != basic) { |
---|
1135 | // column basic |
---|
1136 | setColumnStatus(iColumn, basic); |
---|
1137 | numberBasic++; |
---|
1138 | //printf("Col %d otherV %g dualDual %g\n",iColumn, |
---|
1139 | // otherValue,dualDual[iColumn]); |
---|
1140 | columnActivity_[iColumn] = -dualDual[iColumn]; |
---|
1141 | columnActivity_[iColumn] = otherValue; |
---|
1142 | reducedCost_[iColumn] = 0.0; |
---|
1143 | } else { |
---|
1144 | reducedCost_[iColumn] = objValue - dualActs[iColumn]; |
---|
1145 | if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) { |
---|
1146 | if (columnUpper_[iColumn] > columnLower_[iColumn]) |
---|
1147 | setColumnStatus(iColumn, atLowerBound); |
---|
1148 | else |
---|
1149 | setColumnStatus(iColumn, isFixed); |
---|
1150 | columnActivity_[iColumn] = columnLower_[iColumn]; |
---|
1151 | } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) { |
---|
1152 | if (columnUpper_[iColumn] > columnLower_[iColumn]) |
---|
1153 | setColumnStatus(iColumn, atUpperBound); |
---|
1154 | else |
---|
1155 | setColumnStatus(iColumn, isFixed); |
---|
1156 | columnActivity_[iColumn] = columnUpper_[iColumn]; |
---|
1157 | } else { |
---|
1158 | setColumnStatus(iColumn, superBasic); |
---|
1159 | columnActivity_[iColumn] = otherValue; |
---|
1160 | } |
---|
1161 | } |
---|
1162 | } |
---|
1163 | } |
---|
1164 | } |
---|
1165 | // now rows |
---|
1166 | int kExtraRow = jColumn; |
---|
1167 | int numberRanges = 0; |
---|
1168 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1169 | Status status = dualProblem->getColumnStatus(iRow); |
---|
1170 | if (status == basic) { |
---|
1171 | // row is at bound |
---|
1172 | dual_[iRow] = dualSol[iRow];; |
---|
1173 | } else { |
---|
1174 | // row basic |
---|
1175 | setRowStatus(iRow, basic); |
---|
1176 | numberBasic++; |
---|
1177 | dual_[iRow] = 0.0; |
---|
1178 | } |
---|
1179 | if (rowLower_[iRow] < -1.0e20) { |
---|
1180 | if (status == basic) { |
---|
1181 | rowActivity_[iRow] = rowUpper_[iRow]; |
---|
1182 | setRowStatus(iRow, atUpperBound); |
---|
1183 | } else { |
---|
1184 | // might be stopped assert (dualDj[iRow] < 1.0e-5); |
---|
1185 | rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow]; |
---|
1186 | } |
---|
1187 | } else if (rowUpper_[iRow] > 1.0e20) { |
---|
1188 | if (status == basic) { |
---|
1189 | rowActivity_[iRow] = rowLower_[iRow]; |
---|
1190 | setRowStatus(iRow, atLowerBound); |
---|
1191 | } else { |
---|
1192 | rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow]; |
---|
1193 | // might be stopped assert (dualDj[iRow] > -1.0e-5); |
---|
1194 | } |
---|
1195 | } else { |
---|
1196 | if (rowUpper_[iRow] == rowLower_[iRow]) { |
---|
1197 | rowActivity_[iRow] = rowLower_[iRow]; |
---|
1198 | if (status == basic) { |
---|
1199 | setRowStatus(iRow, isFixed); |
---|
1200 | } |
---|
1201 | } else { |
---|
1202 | // range |
---|
1203 | numberRanges++; |
---|
1204 | Status statusL = dualProblem->getColumnStatus(kExtraRow); |
---|
1205 | //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n", |
---|
1206 | // iRow,status,kExtraRow,statusL, dualSol[iRow], |
---|
1207 | // dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]); |
---|
1208 | if (status == basic) { |
---|
1209 | // might be stopped assert (statusL != basic); |
---|
1210 | rowActivity_[iRow] = rowUpper_[iRow]; |
---|
1211 | setRowStatus(iRow, atUpperBound); |
---|
1212 | } else if (statusL == basic) { |
---|
1213 | numberBasic--; // already counted |
---|
1214 | rowActivity_[iRow] = rowLower_[iRow]; |
---|
1215 | setRowStatus(iRow, atLowerBound); |
---|
1216 | dual_[iRow] = dualSol[kExtraRow];; |
---|
1217 | } else { |
---|
1218 | rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow]; |
---|
1219 | // might be stopped assert (dualDj[iRow] < 1.0e-5); |
---|
1220 | // row basic |
---|
1221 | //setRowStatus(iRow,basic); |
---|
1222 | //numberBasic++; |
---|
1223 | dual_[iRow] = 0.0; |
---|
1224 | } |
---|
1225 | kExtraRow++; |
---|
1226 | } |
---|
1227 | } |
---|
1228 | } |
---|
1229 | if (numberBasic != numberRows_) { |
---|
1230 | printf("Bad basis - ranges - coding needed\n"); |
---|
1231 | assert (numberRanges); |
---|
1232 | abort(); |
---|
1233 | } |
---|
1234 | if (optimizationDirection_ < 0.0) { |
---|
1235 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1236 | dual_[iRow] = -dual_[iRow]; |
---|
1237 | } |
---|
1238 | } |
---|
1239 | // redo row activities |
---|
1240 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
1241 | matrix_->times(1.0, columnActivity_, rowActivity_); |
---|
1242 | // redo reduced costs |
---|
1243 | memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double)); |
---|
1244 | matrix_->transposeTimes(-1.0, dual_, reducedCost_); |
---|
1245 | checkSolutionInternal(); |
---|
1246 | if (sumDualInfeasibilities_ > 1.0e-5 || sumPrimalInfeasibilities_ > 1.0e-5) { |
---|
1247 | returnCode = 1; |
---|
1248 | #ifdef CLP_INVESTIGATE |
---|
1249 | printf("There are %d dual infeasibilities summing to %g ", |
---|
1250 | numberDualInfeasibilities_, sumDualInfeasibilities_); |
---|
1251 | printf("and %d primal infeasibilities summing to %g\n", |
---|
1252 | numberPrimalInfeasibilities_, sumPrimalInfeasibilities_); |
---|
1253 | #endif |
---|
1254 | } |
---|
1255 | // Below will go to ..DEBUG later |
---|
1256 | #if 1 //ndef NDEBUG |
---|
1257 | if (checkAccuracy) { |
---|
1258 | // Check if correct |
---|
1259 | double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_); |
---|
1260 | double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_); |
---|
1261 | double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_); |
---|
1262 | double * dual = CoinCopyOfArray(dual_, numberRows_); |
---|
1263 | this->dual(); //primal(); |
---|
1264 | CoinRelFltEq eq(1.0e-5); |
---|
1265 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1266 | assert(eq(dual[iRow], dual_[iRow])); |
---|
1267 | } |
---|
1268 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1269 | assert(eq(columnActivity[iColumn], columnActivity_[iColumn])); |
---|
1270 | } |
---|
1271 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1272 | assert(eq(rowActivity[iRow], rowActivity_[iRow])); |
---|
1273 | } |
---|
1274 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1275 | assert(eq(reducedCost[iColumn], reducedCost_[iColumn])); |
---|
1276 | } |
---|
1277 | delete [] columnActivity; |
---|
1278 | delete [] rowActivity; |
---|
1279 | delete [] reducedCost; |
---|
1280 | delete [] dual; |
---|
1281 | } |
---|
1282 | #endif |
---|
1283 | return returnCode; |
---|
1284 | } |
---|
1285 | /* Sets solution in dualized problem |
---|
1286 | non-zero return code indicates minor problems |
---|
1287 | */ |
---|
1288 | int |
---|
1289 | ClpSimplexOther::setInDual(ClpSimplex * dualProblem) |
---|
1290 | { |
---|
1291 | // Number of rows in dual problem was original number of columns |
---|
1292 | assert (numberColumns_ == dualProblem->numberRows()); |
---|
1293 | // out If slack on d-row basic then column at bound otherwise column basic |
---|
1294 | // out If d-column basic then rhs tight |
---|
1295 | // if column at bound then slack on d-row basic |
---|
1296 | // if column basic then slack on d-row at bound |
---|
1297 | // if rhs non-basic then d-column basic |
---|
1298 | // if rhs basic then d-column ? |
---|
1299 | int numberBasic = 0; |
---|
1300 | int iRow, iColumn = 0; |
---|
1301 | //int numberExtraRows = dualProblem->numberColumns()-numberRows_; |
---|
1302 | //const double * objective = this->objective(); |
---|
1303 | //double * dualDual = dualProblem->dualRowSolution(); |
---|
1304 | //double * dualDj = dualProblem->dualColumnSolution(); |
---|
1305 | double * dualSol = dualProblem->primalColumnSolution(); |
---|
1306 | //double * dualActs = dualProblem->primalRowSolution(); |
---|
1307 | const double * lower = dualProblem->columnLower(); |
---|
1308 | const double * upper = dualProblem->columnUpper(); |
---|
1309 | // position at bound information |
---|
1310 | int jColumn = numberRows_; |
---|
1311 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1312 | Status status = getColumnStatus(iColumn); |
---|
1313 | Status statusD = dualProblem->getRowStatus(iColumn); |
---|
1314 | Status statusDJ = dualProblem->getColumnStatus(jColumn); |
---|
1315 | if (status==atLowerBound|| |
---|
1316 | status==isFixed|| |
---|
1317 | status==atUpperBound) { |
---|
1318 | dualProblem->setRowStatus(iColumn,basic); |
---|
1319 | numberBasic++; |
---|
1320 | if (columnUpper_[iColumn] < 1.0e20 && |
---|
1321 | columnLower_[iColumn] > -1.0e20) { |
---|
1322 | bool mainLower =(fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])); |
---|
1323 | // fix this |
---|
1324 | if (mainLower) { |
---|
1325 | if (status==atUpperBound) { |
---|
1326 | dualProblem->setColumnStatus(jColumn,atUpperBound); |
---|
1327 | } else { |
---|
1328 | dualProblem->setColumnStatus(jColumn,atUpperBound); |
---|
1329 | } |
---|
1330 | } else { |
---|
1331 | if (status==atUpperBound) { |
---|
1332 | dualProblem->setColumnStatus(jColumn,atLowerBound); |
---|
1333 | } else { |
---|
1334 | dualProblem->setColumnStatus(jColumn,atLowerBound); |
---|
1335 | } |
---|
1336 | } |
---|
1337 | assert(statusDJ == dualProblem->getColumnStatus(jColumn)); |
---|
1338 | jColumn++; |
---|
1339 | } |
---|
1340 | } else if (status==isFree) { |
---|
1341 | dualProblem->setRowStatus(iColumn,basic); |
---|
1342 | numberBasic++; |
---|
1343 | } else { |
---|
1344 | assert (status==basic); |
---|
1345 | //numberBasic++; |
---|
1346 | } |
---|
1347 | assert(statusD == dualProblem->getRowStatus(iColumn)); |
---|
1348 | } |
---|
1349 | // now rows (no ranges at first) |
---|
1350 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1351 | Status status = getRowStatus(iRow); |
---|
1352 | Status statusD = dualProblem->getColumnStatus(iRow); |
---|
1353 | if (status == basic) { |
---|
1354 | // dual variable is at bound |
---|
1355 | if (!lower[iRow]) { |
---|
1356 | dualProblem->setColumnStatus(iRow,atLowerBound); |
---|
1357 | } else if (!upper[iRow]) { |
---|
1358 | dualProblem->setColumnStatus(iRow,atUpperBound); |
---|
1359 | } else { |
---|
1360 | dualProblem->setColumnStatus(iRow,isFree); |
---|
1361 | dualSol[iRow]=0.0; |
---|
1362 | } |
---|
1363 | } else { |
---|
1364 | // dual variable is basic |
---|
1365 | dualProblem->setColumnStatus(iRow,basic); |
---|
1366 | numberBasic++; |
---|
1367 | } |
---|
1368 | if (rowLower_[iRow] < -1.0e20 && |
---|
1369 | rowUpper_[iRow] > 1.0e20) { |
---|
1370 | if (rowUpper_[iRow] != rowLower_[iRow]) { |
---|
1371 | printf("can't handle ranges yet\n"); |
---|
1372 | abort(); |
---|
1373 | } |
---|
1374 | } |
---|
1375 | assert(statusD == dualProblem->getColumnStatus(iRow)); |
---|
1376 | } |
---|
1377 | if (numberBasic != numberColumns_) { |
---|
1378 | printf("Bad basis - ranges - coding needed ??\n"); |
---|
1379 | abort(); |
---|
1380 | } |
---|
1381 | return 0; |
---|
1382 | } |
---|
1383 | /* Does very cursory presolve. |
---|
1384 | rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns |
---|
1385 | */ |
---|
1386 | ClpSimplex * |
---|
1387 | ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn, |
---|
1388 | int & nBound, bool moreBounds, bool tightenBounds) |
---|
1389 | { |
---|
1390 | //#define CHECK_STATUS |
---|
1391 | #ifdef CHECK_STATUS |
---|
1392 | { |
---|
1393 | int n = 0; |
---|
1394 | int i; |
---|
1395 | for (i = 0; i < numberColumns_; i++) |
---|
1396 | if (getColumnStatus(i) == ClpSimplex::basic) |
---|
1397 | n++; |
---|
1398 | for (i = 0; i < numberRows_; i++) |
---|
1399 | if (getRowStatus(i) == ClpSimplex::basic) |
---|
1400 | n++; |
---|
1401 | assert (n == numberRows_); |
---|
1402 | } |
---|
1403 | #endif |
---|
1404 | |
---|
1405 | const double * element = matrix_->getElements(); |
---|
1406 | const int * row = matrix_->getIndices(); |
---|
1407 | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
---|
1408 | const int * columnLength = matrix_->getVectorLengths(); |
---|
1409 | |
---|
1410 | CoinZeroN(rhs, numberRows_); |
---|
1411 | int iColumn; |
---|
1412 | int iRow; |
---|
1413 | CoinZeroN(whichRow, numberRows_); |
---|
1414 | int * backColumn = whichColumn + numberColumns_; |
---|
1415 | int numberRows2 = 0; |
---|
1416 | int numberColumns2 = 0; |
---|
1417 | double offset = 0.0; |
---|
1418 | const double * objective = this->objective(); |
---|
1419 | double * solution = columnActivity_; |
---|
1420 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1421 | double lower = columnLower_[iColumn]; |
---|
1422 | double upper = columnUpper_[iColumn]; |
---|
1423 | if (upper > lower || getColumnStatus(iColumn) == ClpSimplex::basic) { |
---|
1424 | backColumn[iColumn] = numberColumns2; |
---|
1425 | whichColumn[numberColumns2++] = iColumn; |
---|
1426 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1427 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1428 | int iRow = row[j]; |
---|
1429 | int n = whichRow[iRow]; |
---|
1430 | if (n == 0 && element[j]) |
---|
1431 | whichRow[iRow] = -iColumn - 1; |
---|
1432 | else if (n < 0) |
---|
1433 | whichRow[iRow] = 2; |
---|
1434 | } |
---|
1435 | } else { |
---|
1436 | // fixed |
---|
1437 | backColumn[iColumn] = -1; |
---|
1438 | solution[iColumn] = upper; |
---|
1439 | if (upper) { |
---|
1440 | offset += objective[iColumn] * upper; |
---|
1441 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1442 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1443 | int iRow = row[j]; |
---|
1444 | double value = element[j]; |
---|
1445 | rhs[iRow] += upper * value; |
---|
1446 | } |
---|
1447 | } |
---|
1448 | } |
---|
1449 | } |
---|
1450 | int returnCode = 0; |
---|
1451 | double tolerance = primalTolerance(); |
---|
1452 | nBound = 2 * numberRows_; |
---|
1453 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1454 | int n = whichRow[iRow]; |
---|
1455 | if (n > 0) { |
---|
1456 | whichRow[numberRows2++] = iRow; |
---|
1457 | } else if (n < 0) { |
---|
1458 | //whichRow[numberRows2++]=iRow; |
---|
1459 | //continue; |
---|
1460 | // Can only do in certain circumstances as we don't know current value |
---|
1461 | if (rowLower_[iRow] == rowUpper_[iRow] || getRowStatus(iRow) == ClpSimplex::basic) { |
---|
1462 | // save row and column for bound |
---|
1463 | whichRow[--nBound] = iRow; |
---|
1464 | whichRow[nBound+numberRows_] = -n - 1; |
---|
1465 | } else if (moreBounds) { |
---|
1466 | // save row and column for bound |
---|
1467 | whichRow[--nBound] = iRow; |
---|
1468 | whichRow[nBound+numberRows_] = -n - 1; |
---|
1469 | } else { |
---|
1470 | whichRow[numberRows2++] = iRow; |
---|
1471 | } |
---|
1472 | } else { |
---|
1473 | // empty |
---|
1474 | double rhsValue = rhs[iRow]; |
---|
1475 | if (rhsValue < rowLower_[iRow] - tolerance || rhsValue > rowUpper_[iRow] + tolerance) { |
---|
1476 | returnCode = 1; // infeasible |
---|
1477 | } |
---|
1478 | } |
---|
1479 | } |
---|
1480 | ClpSimplex * small = NULL; |
---|
1481 | if (!returnCode) { |
---|
1482 | //printf("CRUNCH from (%d,%d) to (%d,%d)\n", |
---|
1483 | // numberRows_,numberColumns_,numberRows2,numberColumns2); |
---|
1484 | small = new ClpSimplex(this, numberRows2, whichRow, |
---|
1485 | numberColumns2, whichColumn, true, false); |
---|
1486 | #if 0 |
---|
1487 | ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_); |
---|
1488 | if (rowCopy) { |
---|
1489 | assert(!small->rowCopy()); |
---|
1490 | small->setNewRowCopy(new ClpPackedMatrix(*rowCopy, numberRows2, whichRow, |
---|
1491 | numberColumns2, whichColumn)); |
---|
1492 | } |
---|
1493 | #endif |
---|
1494 | // Set some stuff |
---|
1495 | small->setDualBound(dualBound_); |
---|
1496 | small->setInfeasibilityCost(infeasibilityCost_); |
---|
1497 | small->setSpecialOptions(specialOptions_); |
---|
1498 | small->setPerturbation(perturbation_); |
---|
1499 | small->defaultFactorizationFrequency(); |
---|
1500 | small->setAlphaAccuracy(alphaAccuracy_); |
---|
1501 | // If no rows left then no tightening! |
---|
1502 | if (!numberRows2 || !numberColumns2) |
---|
1503 | tightenBounds = false; |
---|
1504 | |
---|
1505 | int numberElements = getNumElements(); |
---|
1506 | int numberElements2 = small->getNumElements(); |
---|
1507 | small->setObjectiveOffset(objectiveOffset() - offset); |
---|
1508 | handler_->message(CLP_CRUNCH_STATS, messages_) |
---|
1509 | << numberRows2 << -(numberRows_ - numberRows2) |
---|
1510 | << numberColumns2 << -(numberColumns_ - numberColumns2) |
---|
1511 | << numberElements2 << -(numberElements - numberElements2) |
---|
1512 | << CoinMessageEol; |
---|
1513 | // And set objective value to match |
---|
1514 | small->setObjectiveValue(this->objectiveValue()); |
---|
1515 | double * rowLower2 = small->rowLower(); |
---|
1516 | double * rowUpper2 = small->rowUpper(); |
---|
1517 | int jRow; |
---|
1518 | for (jRow = 0; jRow < numberRows2; jRow++) { |
---|
1519 | iRow = whichRow[jRow]; |
---|
1520 | if (rowLower2[jRow] > -1.0e20) |
---|
1521 | rowLower2[jRow] -= rhs[iRow]; |
---|
1522 | if (rowUpper2[jRow] < 1.0e20) |
---|
1523 | rowUpper2[jRow] -= rhs[iRow]; |
---|
1524 | } |
---|
1525 | // and bounds |
---|
1526 | double * columnLower2 = small->columnLower(); |
---|
1527 | double * columnUpper2 = small->columnUpper(); |
---|
1528 | const char * integerInformation = integerType_; |
---|
1529 | for (jRow = nBound; jRow < 2 * numberRows_; jRow++) { |
---|
1530 | iRow = whichRow[jRow]; |
---|
1531 | iColumn = whichRow[jRow+numberRows_]; |
---|
1532 | double lowerRow = rowLower_[iRow]; |
---|
1533 | if (lowerRow > -1.0e20) |
---|
1534 | lowerRow -= rhs[iRow]; |
---|
1535 | double upperRow = rowUpper_[iRow]; |
---|
1536 | if (upperRow < 1.0e20) |
---|
1537 | upperRow -= rhs[iRow]; |
---|
1538 | int jColumn = backColumn[iColumn]; |
---|
1539 | double lower = columnLower2[jColumn]; |
---|
1540 | double upper = columnUpper2[jColumn]; |
---|
1541 | double value = 0.0; |
---|
1542 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1543 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1544 | if (iRow == row[j]) { |
---|
1545 | value = element[j]; |
---|
1546 | break; |
---|
1547 | } |
---|
1548 | } |
---|
1549 | assert (value); |
---|
1550 | // convert rowLower and Upper to implied bounds on column |
---|
1551 | double newLower = -COIN_DBL_MAX; |
---|
1552 | double newUpper = COIN_DBL_MAX; |
---|
1553 | if (value > 0.0) { |
---|
1554 | if (lowerRow > -1.0e20) |
---|
1555 | newLower = lowerRow / value; |
---|
1556 | if (upperRow < 1.0e20) |
---|
1557 | newUpper = upperRow / value; |
---|
1558 | } else { |
---|
1559 | if (upperRow < 1.0e20) |
---|
1560 | newLower = upperRow / value; |
---|
1561 | if (lowerRow > -1.0e20) |
---|
1562 | newUpper = lowerRow / value; |
---|
1563 | } |
---|
1564 | if (integerInformation && integerInformation[iColumn]) { |
---|
1565 | if (newLower - floor(newLower) < 10.0 * tolerance) |
---|
1566 | newLower = floor(newLower); |
---|
1567 | else |
---|
1568 | newLower = ceil(newLower); |
---|
1569 | if (ceil(newUpper) - newUpper < 10.0 * tolerance) |
---|
1570 | newUpper = ceil(newUpper); |
---|
1571 | else |
---|
1572 | newUpper = floor(newUpper); |
---|
1573 | } |
---|
1574 | newLower = CoinMax(lower, newLower); |
---|
1575 | newUpper = CoinMin(upper, newUpper); |
---|
1576 | if (newLower > newUpper + tolerance) { |
---|
1577 | //printf("XXYY inf on bound\n"); |
---|
1578 | returnCode = 1; |
---|
1579 | } |
---|
1580 | columnLower2[jColumn] = newLower; |
---|
1581 | columnUpper2[jColumn] = CoinMax(newLower, newUpper); |
---|
1582 | if (getRowStatus(iRow) != ClpSimplex::basic) { |
---|
1583 | if (getColumnStatus(iColumn) == ClpSimplex::basic) { |
---|
1584 | if (columnLower2[jColumn] == columnUpper2[jColumn]) { |
---|
1585 | // can only get here if will be fixed |
---|
1586 | small->setColumnStatus(jColumn, ClpSimplex::isFixed); |
---|
1587 | } else { |
---|
1588 | // solution is valid |
---|
1589 | if (fabs(columnActivity_[iColumn] - columnLower2[jColumn]) < |
---|
1590 | fabs(columnActivity_[iColumn] - columnUpper2[jColumn])) |
---|
1591 | small->setColumnStatus(jColumn, ClpSimplex::atLowerBound); |
---|
1592 | else |
---|
1593 | small->setColumnStatus(jColumn, ClpSimplex::atUpperBound); |
---|
1594 | } |
---|
1595 | } else { |
---|
1596 | //printf("what now neither basic\n"); |
---|
1597 | } |
---|
1598 | } |
---|
1599 | } |
---|
1600 | if (returnCode) { |
---|
1601 | delete small; |
---|
1602 | small = NULL; |
---|
1603 | } else if (tightenBounds && integerInformation) { |
---|
1604 | // See if we can tighten any bounds |
---|
1605 | // use rhs for upper and small duals for lower |
---|
1606 | double * up = rhs; |
---|
1607 | double * lo = small->dualRowSolution(); |
---|
1608 | const double * element = small->clpMatrix()->getElements(); |
---|
1609 | const int * row = small->clpMatrix()->getIndices(); |
---|
1610 | const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts(); |
---|
1611 | //const int * columnLength = small->clpMatrix()->getVectorLengths(); |
---|
1612 | CoinZeroN(lo, numberRows2); |
---|
1613 | CoinZeroN(up, numberRows2); |
---|
1614 | for (int iColumn = 0; iColumn < numberColumns2; iColumn++) { |
---|
1615 | double upper = columnUpper2[iColumn]; |
---|
1616 | double lower = columnLower2[iColumn]; |
---|
1617 | //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]); |
---|
1618 | for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) { |
---|
1619 | int iRow = row[j]; |
---|
1620 | double value = element[j]; |
---|
1621 | if (value > 0.0) { |
---|
1622 | if (upper < 1.0e20) |
---|
1623 | up[iRow] += upper * value; |
---|
1624 | else |
---|
1625 | up[iRow] = COIN_DBL_MAX; |
---|
1626 | if (lower > -1.0e20) |
---|
1627 | lo[iRow] += lower * value; |
---|
1628 | else |
---|
1629 | lo[iRow] = -COIN_DBL_MAX; |
---|
1630 | } else { |
---|
1631 | if (upper < 1.0e20) |
---|
1632 | lo[iRow] += upper * value; |
---|
1633 | else |
---|
1634 | lo[iRow] = -COIN_DBL_MAX; |
---|
1635 | if (lower > -1.0e20) |
---|
1636 | up[iRow] += lower * value; |
---|
1637 | else |
---|
1638 | up[iRow] = COIN_DBL_MAX; |
---|
1639 | } |
---|
1640 | } |
---|
1641 | } |
---|
1642 | double * rowLower2 = small->rowLower(); |
---|
1643 | double * rowUpper2 = small->rowUpper(); |
---|
1644 | bool feasible = true; |
---|
1645 | // make safer |
---|
1646 | for (int iRow = 0; iRow < numberRows2; iRow++) { |
---|
1647 | double lower = lo[iRow]; |
---|
1648 | if (lower > rowUpper2[iRow] + tolerance) { |
---|
1649 | feasible = false; |
---|
1650 | break; |
---|
1651 | } else { |
---|
1652 | lo[iRow] = CoinMin(lower - rowUpper2[iRow], 0.0) - tolerance; |
---|
1653 | } |
---|
1654 | double upper = up[iRow]; |
---|
1655 | if (upper < rowLower2[iRow] - tolerance) { |
---|
1656 | feasible = false; |
---|
1657 | break; |
---|
1658 | } else { |
---|
1659 | up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance; |
---|
1660 | } |
---|
1661 | // tighten row bounds |
---|
1662 | if (lower>-1.0e10) |
---|
1663 | rowLower2[iRow] = CoinMax(rowLower2[iRow], |
---|
1664 | lower - 1.0e-6*(1.0+fabs(lower))); |
---|
1665 | if (upper<1.0e10) |
---|
1666 | rowUpper2[iRow] = CoinMin(rowUpper2[iRow], |
---|
1667 | upper + 1.0e-6*(1.0+fabs(upper))); |
---|
1668 | } |
---|
1669 | if (!feasible) { |
---|
1670 | delete small; |
---|
1671 | small = NULL; |
---|
1672 | } else { |
---|
1673 | // and tighten |
---|
1674 | for (int iColumn = 0; iColumn < numberColumns2; iColumn++) { |
---|
1675 | if (integerInformation[whichColumn[iColumn]]) { |
---|
1676 | double upper = columnUpper2[iColumn]; |
---|
1677 | double lower = columnLower2[iColumn]; |
---|
1678 | double newUpper = upper; |
---|
1679 | double newLower = lower; |
---|
1680 | double difference = upper - lower; |
---|
1681 | if (lower > -1000.0 && upper < 1000.0) { |
---|
1682 | for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) { |
---|
1683 | int iRow = row[j]; |
---|
1684 | double value = element[j]; |
---|
1685 | if (value > 0.0) { |
---|
1686 | double upWithOut = up[iRow] - value * difference; |
---|
1687 | if (upWithOut < 0.0) { |
---|
1688 | newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value); |
---|
1689 | } |
---|
1690 | double lowWithOut = lo[iRow] + value * difference; |
---|
1691 | if (lowWithOut > 0.0) { |
---|
1692 | newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value); |
---|
1693 | } |
---|
1694 | } else { |
---|
1695 | double upWithOut = up[iRow] + value * difference; |
---|
1696 | if (upWithOut < 0.0) { |
---|
1697 | newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value); |
---|
1698 | } |
---|
1699 | double lowWithOut = lo[iRow] - value * difference; |
---|
1700 | if (lowWithOut > 0.0) { |
---|
1701 | newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value); |
---|
1702 | } |
---|
1703 | } |
---|
1704 | } |
---|
1705 | if (newLower > lower || newUpper < upper) { |
---|
1706 | if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6) |
---|
1707 | newUpper = floor(newUpper); |
---|
1708 | else |
---|
1709 | newUpper = floor(newUpper + 0.5); |
---|
1710 | if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6) |
---|
1711 | newLower = ceil(newLower); |
---|
1712 | else |
---|
1713 | newLower = ceil(newLower - 0.5); |
---|
1714 | // change may be too small - check |
---|
1715 | if (newLower > lower || newUpper < upper) { |
---|
1716 | if (newUpper >= newLower) { |
---|
1717 | // Could also tighten in this |
---|
1718 | //printf("%d bounds %g %g tightened to %g %g\n", |
---|
1719 | // iColumn,columnLower2[iColumn],columnUpper2[iColumn], |
---|
1720 | // newLower,newUpper); |
---|
1721 | #if 1 |
---|
1722 | columnUpper2[iColumn] = newUpper; |
---|
1723 | columnLower2[iColumn] = newLower; |
---|
1724 | columnUpper_[whichColumn[iColumn]] = newUpper; |
---|
1725 | columnLower_[whichColumn[iColumn]] = newLower; |
---|
1726 | #endif |
---|
1727 | // and adjust bounds on rows |
---|
1728 | newUpper -= upper; |
---|
1729 | newLower -= lower; |
---|
1730 | for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) { |
---|
1731 | int iRow = row[j]; |
---|
1732 | double value = element[j]; |
---|
1733 | if (value > 0.0) { |
---|
1734 | up[iRow] += newUpper * value; |
---|
1735 | lo[iRow] += newLower * value; |
---|
1736 | } else { |
---|
1737 | lo[iRow] += newUpper * value; |
---|
1738 | up[iRow] += newLower * value; |
---|
1739 | } |
---|
1740 | } |
---|
1741 | } else { |
---|
1742 | // infeasible |
---|
1743 | //printf("%d bounds infeasible %g %g tightened to %g %g\n", |
---|
1744 | // iColumn,columnLower2[iColumn],columnUpper2[iColumn], |
---|
1745 | // newLower,newUpper); |
---|
1746 | #if 1 |
---|
1747 | delete small; |
---|
1748 | small = NULL; |
---|
1749 | break; |
---|
1750 | #endif |
---|
1751 | } |
---|
1752 | } |
---|
1753 | } |
---|
1754 | } |
---|
1755 | } |
---|
1756 | } |
---|
1757 | } |
---|
1758 | } |
---|
1759 | } |
---|
1760 | #if 0 |
---|
1761 | if (small) { |
---|
1762 | static int which = 0; |
---|
1763 | which++; |
---|
1764 | char xxxx[20]; |
---|
1765 | sprintf(xxxx, "bad%d.mps", which); |
---|
1766 | small->writeMps(xxxx, 0, 1); |
---|
1767 | sprintf(xxxx, "largebad%d.mps", which); |
---|
1768 | writeMps(xxxx, 0, 1); |
---|
1769 | printf("bad%d %x old size %d %d new %d %d\n", which, small, |
---|
1770 | numberRows_, numberColumns_, small->numberRows(), small->numberColumns()); |
---|
1771 | #if 0 |
---|
1772 | for (int i = 0; i < numberColumns_; i++) |
---|
1773 | printf("Bound %d %g %g\n", i, columnLower_[i], columnUpper_[i]); |
---|
1774 | for (int i = 0; i < numberRows_; i++) |
---|
1775 | printf("Row bound %d %g %g\n", i, rowLower_[i], rowUpper_[i]); |
---|
1776 | #endif |
---|
1777 | } |
---|
1778 | #endif |
---|
1779 | #ifdef CHECK_STATUS |
---|
1780 | { |
---|
1781 | int n = 0; |
---|
1782 | int i; |
---|
1783 | for (i = 0; i < small->numberColumns(); i++) |
---|
1784 | if (small->getColumnStatus(i) == ClpSimplex::basic) |
---|
1785 | n++; |
---|
1786 | for (i = 0; i < small->numberRows(); i++) |
---|
1787 | if (small->getRowStatus(i) == ClpSimplex::basic) |
---|
1788 | n++; |
---|
1789 | assert (n == small->numberRows()); |
---|
1790 | } |
---|
1791 | #endif |
---|
1792 | return small; |
---|
1793 | } |
---|
1794 | /* After very cursory presolve. |
---|
1795 | rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns. |
---|
1796 | */ |
---|
1797 | void |
---|
1798 | ClpSimplexOther::afterCrunch(const ClpSimplex & small, |
---|
1799 | const int * whichRow, |
---|
1800 | const int * whichColumn, int nBound) |
---|
1801 | { |
---|
1802 | #ifndef NDEBUG |
---|
1803 | for (int i = 0; i < small.numberRows(); i++) |
---|
1804 | assert (whichRow[i] >= 0 && whichRow[i] < numberRows_); |
---|
1805 | for (int i = 0; i < small.numberColumns(); i++) |
---|
1806 | assert (whichColumn[i] >= 0 && whichColumn[i] < numberColumns_); |
---|
1807 | #endif |
---|
1808 | getbackSolution(small, whichRow, whichColumn); |
---|
1809 | // and deal with status for bounds |
---|
1810 | const double * element = matrix_->getElements(); |
---|
1811 | const int * row = matrix_->getIndices(); |
---|
1812 | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
---|
1813 | const int * columnLength = matrix_->getVectorLengths(); |
---|
1814 | double tolerance = primalTolerance(); |
---|
1815 | double djTolerance = dualTolerance(); |
---|
1816 | for (int jRow = nBound; jRow < 2 * numberRows_; jRow++) { |
---|
1817 | int iRow = whichRow[jRow]; |
---|
1818 | int iColumn = whichRow[jRow+numberRows_]; |
---|
1819 | if (getColumnStatus(iColumn) != ClpSimplex::basic) { |
---|
1820 | double lower = columnLower_[iColumn]; |
---|
1821 | double upper = columnUpper_[iColumn]; |
---|
1822 | double value = columnActivity_[iColumn]; |
---|
1823 | double djValue = reducedCost_[iColumn]; |
---|
1824 | dual_[iRow] = 0.0; |
---|
1825 | if (upper > lower) { |
---|
1826 | if (value < lower + tolerance && djValue > -djTolerance) { |
---|
1827 | setColumnStatus(iColumn, ClpSimplex::atLowerBound); |
---|
1828 | setRowStatus(iRow, ClpSimplex::basic); |
---|
1829 | } else if (value > upper - tolerance && djValue < djTolerance) { |
---|
1830 | setColumnStatus(iColumn, ClpSimplex::atUpperBound); |
---|
1831 | setRowStatus(iRow, ClpSimplex::basic); |
---|
1832 | } else { |
---|
1833 | // has to be basic |
---|
1834 | setColumnStatus(iColumn, ClpSimplex::basic); |
---|
1835 | reducedCost_[iColumn] = 0.0; |
---|
1836 | double value = 0.0; |
---|
1837 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1838 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1839 | if (iRow == row[j]) { |
---|
1840 | value = element[j]; |
---|
1841 | break; |
---|
1842 | } |
---|
1843 | } |
---|
1844 | dual_[iRow] = djValue / value; |
---|
1845 | if (rowUpper_[iRow] > rowLower_[iRow]) { |
---|
1846 | if (fabs(rowActivity_[iRow] - rowLower_[iRow]) < |
---|
1847 | fabs(rowActivity_[iRow] - rowUpper_[iRow])) |
---|
1848 | setRowStatus(iRow, ClpSimplex::atLowerBound); |
---|
1849 | else |
---|
1850 | setRowStatus(iRow, ClpSimplex::atUpperBound); |
---|
1851 | } else { |
---|
1852 | setRowStatus(iRow, ClpSimplex::isFixed); |
---|
1853 | } |
---|
1854 | } |
---|
1855 | } else { |
---|
1856 | // row can always be basic |
---|
1857 | setRowStatus(iRow, ClpSimplex::basic); |
---|
1858 | } |
---|
1859 | } else { |
---|
1860 | // row can always be basic |
---|
1861 | setRowStatus(iRow, ClpSimplex::basic); |
---|
1862 | } |
---|
1863 | } |
---|
1864 | //#ifndef NDEBUG |
---|
1865 | #if 0 |
---|
1866 | if (small.status() == 0) { |
---|
1867 | int n = 0; |
---|
1868 | int i; |
---|
1869 | for (i = 0; i < numberColumns; i++) |
---|
1870 | if (getColumnStatus(i) == ClpSimplex::basic) |
---|
1871 | n++; |
---|
1872 | for (i = 0; i < numberRows; i++) |
---|
1873 | if (getRowStatus(i) == ClpSimplex::basic) |
---|
1874 | n++; |
---|
1875 | assert (n == numberRows); |
---|
1876 | } |
---|
1877 | #endif |
---|
1878 | } |
---|
1879 | /* Tightens integer bounds - returns number tightened or -1 if infeasible |
---|
1880 | */ |
---|
1881 | int |
---|
1882 | ClpSimplexOther::tightenIntegerBounds(double * rhsSpace) |
---|
1883 | { |
---|
1884 | // See if we can tighten any bounds |
---|
1885 | // use rhs for upper and small duals for lower |
---|
1886 | double * up = rhsSpace; |
---|
1887 | double * lo = dual_; |
---|
1888 | const double * element = matrix_->getElements(); |
---|
1889 | const int * row = matrix_->getIndices(); |
---|
1890 | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
---|
1891 | const int * columnLength = matrix_->getVectorLengths(); |
---|
1892 | CoinZeroN(lo, numberRows_); |
---|
1893 | CoinZeroN(up, numberRows_); |
---|
1894 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1895 | double upper = columnUpper_[iColumn]; |
---|
1896 | double lower = columnLower_[iColumn]; |
---|
1897 | //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]); |
---|
1898 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1899 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1900 | int iRow = row[j]; |
---|
1901 | double value = element[j]; |
---|
1902 | if (value > 0.0) { |
---|
1903 | if (upper < 1.0e20) |
---|
1904 | up[iRow] += upper * value; |
---|
1905 | else |
---|
1906 | up[iRow] = COIN_DBL_MAX; |
---|
1907 | if (lower > -1.0e20) |
---|
1908 | lo[iRow] += lower * value; |
---|
1909 | else |
---|
1910 | lo[iRow] = -COIN_DBL_MAX; |
---|
1911 | } else { |
---|
1912 | if (upper < 1.0e20) |
---|
1913 | lo[iRow] += upper * value; |
---|
1914 | else |
---|
1915 | lo[iRow] = -COIN_DBL_MAX; |
---|
1916 | if (lower > -1.0e20) |
---|
1917 | up[iRow] += lower * value; |
---|
1918 | else |
---|
1919 | up[iRow] = COIN_DBL_MAX; |
---|
1920 | } |
---|
1921 | } |
---|
1922 | } |
---|
1923 | bool feasible = true; |
---|
1924 | // make safer |
---|
1925 | double tolerance = primalTolerance(); |
---|
1926 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
1927 | double lower = lo[iRow]; |
---|
1928 | if (lower > rowUpper_[iRow] + tolerance) { |
---|
1929 | feasible = false; |
---|
1930 | break; |
---|
1931 | } else { |
---|
1932 | lo[iRow] = CoinMin(lower - rowUpper_[iRow], 0.0) - tolerance; |
---|
1933 | } |
---|
1934 | double upper = up[iRow]; |
---|
1935 | if (upper < rowLower_[iRow] - tolerance) { |
---|
1936 | feasible = false; |
---|
1937 | break; |
---|
1938 | } else { |
---|
1939 | up[iRow] = CoinMax(upper - rowLower_[iRow], 0.0) + tolerance; |
---|
1940 | } |
---|
1941 | } |
---|
1942 | int numberTightened = 0; |
---|
1943 | if (!feasible) { |
---|
1944 | return -1; |
---|
1945 | } else if (integerType_) { |
---|
1946 | // and tighten |
---|
1947 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
1948 | if (integerType_[iColumn]) { |
---|
1949 | double upper = columnUpper_[iColumn]; |
---|
1950 | double lower = columnLower_[iColumn]; |
---|
1951 | double newUpper = upper; |
---|
1952 | double newLower = lower; |
---|
1953 | double difference = upper - lower; |
---|
1954 | if (lower > -1000.0 && upper < 1000.0) { |
---|
1955 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1956 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1957 | int iRow = row[j]; |
---|
1958 | double value = element[j]; |
---|
1959 | if (value > 0.0) { |
---|
1960 | double upWithOut = up[iRow] - value * difference; |
---|
1961 | if (upWithOut < 0.0) { |
---|
1962 | newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value); |
---|
1963 | } |
---|
1964 | double lowWithOut = lo[iRow] + value * difference; |
---|
1965 | if (lowWithOut > 0.0) { |
---|
1966 | newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value); |
---|
1967 | } |
---|
1968 | } else { |
---|
1969 | double upWithOut = up[iRow] + value * difference; |
---|
1970 | if (upWithOut < 0.0) { |
---|
1971 | newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value); |
---|
1972 | } |
---|
1973 | double lowWithOut = lo[iRow] - value * difference; |
---|
1974 | if (lowWithOut > 0.0) { |
---|
1975 | newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value); |
---|
1976 | } |
---|
1977 | } |
---|
1978 | } |
---|
1979 | if (newLower > lower || newUpper < upper) { |
---|
1980 | if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6) |
---|
1981 | newUpper = floor(newUpper); |
---|
1982 | else |
---|
1983 | newUpper = floor(newUpper + 0.5); |
---|
1984 | if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6) |
---|
1985 | newLower = ceil(newLower); |
---|
1986 | else |
---|
1987 | newLower = ceil(newLower - 0.5); |
---|
1988 | // change may be too small - check |
---|
1989 | if (newLower > lower || newUpper < upper) { |
---|
1990 | if (newUpper >= newLower) { |
---|
1991 | numberTightened++; |
---|
1992 | //printf("%d bounds %g %g tightened to %g %g\n", |
---|
1993 | // iColumn,columnLower_[iColumn],columnUpper_[iColumn], |
---|
1994 | // newLower,newUpper); |
---|
1995 | columnUpper_[iColumn] = newUpper; |
---|
1996 | columnLower_[iColumn] = newLower; |
---|
1997 | // and adjust bounds on rows |
---|
1998 | newUpper -= upper; |
---|
1999 | newLower -= lower; |
---|
2000 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
2001 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2002 | int iRow = row[j]; |
---|
2003 | double value = element[j]; |
---|
2004 | if (value > 0.0) { |
---|
2005 | up[iRow] += newUpper * value; |
---|
2006 | lo[iRow] += newLower * value; |
---|
2007 | } else { |
---|
2008 | lo[iRow] += newUpper * value; |
---|
2009 | up[iRow] += newLower * value; |
---|
2010 | } |
---|
2011 | } |
---|
2012 | } else { |
---|
2013 | // infeasible |
---|
2014 | //printf("%d bounds infeasible %g %g tightened to %g %g\n", |
---|
2015 | // iColumn,columnLower_[iColumn],columnUpper_[iColumn], |
---|
2016 | // newLower,newUpper); |
---|
2017 | return -1; |
---|
2018 | } |
---|
2019 | } |
---|
2020 | } |
---|
2021 | } |
---|
2022 | } |
---|
2023 | } |
---|
2024 | } |
---|
2025 | return numberTightened; |
---|
2026 | } |
---|
2027 | /* Parametrics |
---|
2028 | This is an initial slow version. |
---|
2029 | The code uses current bounds + theta * change (if change array not NULL) |
---|
2030 | and similarly for objective. |
---|
2031 | It starts at startingTheta and returns ending theta in endingTheta. |
---|
2032 | If reportIncrement 0.0 it will report on any movement |
---|
2033 | If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement. |
---|
2034 | If it can not reach input endingTheta return code will be 1 for infeasible, |
---|
2035 | 2 for unbounded, if error on ranges -1, otherwise 0. |
---|
2036 | Normal report is just theta and objective but |
---|
2037 | if event handler exists it may do more |
---|
2038 | On exit endingTheta is maximum reached (can be used for next startingTheta) |
---|
2039 | */ |
---|
2040 | int |
---|
2041 | ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, double reportIncrement, |
---|
2042 | const double * lowerChangeBound, const double * upperChangeBound, |
---|
2043 | const double * lowerChangeRhs, const double * upperChangeRhs, |
---|
2044 | const double * changeObjective) |
---|
2045 | { |
---|
2046 | bool needToDoSomething = true; |
---|
2047 | bool canTryQuick = (reportIncrement) ? true : false; |
---|
2048 | // Save copy of model |
---|
2049 | ClpSimplex copyModel = *this; |
---|
2050 | int savePerturbation = perturbation_; |
---|
2051 | perturbation_ = 102; // switch off |
---|
2052 | while (needToDoSomething) { |
---|
2053 | needToDoSomething = false; |
---|
2054 | algorithm_ = -1; |
---|
2055 | |
---|
2056 | // save data |
---|
2057 | ClpDataSave data = saveData(); |
---|
2058 | // Dantzig |
---|
2059 | ClpDualRowPivot * savePivot = dualRowPivot_; |
---|
2060 | dualRowPivot_ = new ClpDualRowDantzig(); |
---|
2061 | dualRowPivot_->setModel(this); |
---|
2062 | int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0); |
---|
2063 | int iRow, iColumn; |
---|
2064 | double * chgUpper = NULL; |
---|
2065 | double * chgLower = NULL; |
---|
2066 | double * chgObjective = NULL; |
---|
2067 | |
---|
2068 | |
---|
2069 | if (!returnCode) { |
---|
2070 | // Find theta when bounds will cross over and create arrays |
---|
2071 | int numberTotal = numberRows_ + numberColumns_; |
---|
2072 | chgLower = new double[numberTotal]; |
---|
2073 | memset(chgLower, 0, numberTotal * sizeof(double)); |
---|
2074 | chgUpper = new double[numberTotal]; |
---|
2075 | memset(chgUpper, 0, numberTotal * sizeof(double)); |
---|
2076 | chgObjective = new double[numberTotal]; |
---|
2077 | memset(chgObjective, 0, numberTotal * sizeof(double)); |
---|
2078 | assert (!rowScale_); |
---|
2079 | double maxTheta = 1.0e50; |
---|
2080 | if (lowerChangeRhs || upperChangeRhs) { |
---|
2081 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2082 | double lower = rowLower_[iRow]; |
---|
2083 | double upper = rowUpper_[iRow]; |
---|
2084 | if (lower > upper) { |
---|
2085 | maxTheta = -1.0; |
---|
2086 | break; |
---|
2087 | } |
---|
2088 | double lowerChange = (lowerChangeRhs) ? lowerChangeRhs[iRow] : 0.0; |
---|
2089 | double upperChange = (upperChangeRhs) ? upperChangeRhs[iRow] : 0.0; |
---|
2090 | if (lower > -1.0e20 && upper < 1.0e20) { |
---|
2091 | if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) { |
---|
2092 | maxTheta = (upper - lower) / (lowerChange - upperChange); |
---|
2093 | } |
---|
2094 | } |
---|
2095 | if (lower > -1.0e20) { |
---|
2096 | lower_[numberColumns_+iRow] += startingTheta * lowerChange; |
---|
2097 | chgLower[numberColumns_+iRow] = lowerChange; |
---|
2098 | } |
---|
2099 | if (upper < 1.0e20) { |
---|
2100 | upper_[numberColumns_+iRow] += startingTheta * upperChange; |
---|
2101 | chgUpper[numberColumns_+iRow] = upperChange; |
---|
2102 | } |
---|
2103 | } |
---|
2104 | } |
---|
2105 | if (maxTheta > 0.0) { |
---|
2106 | if (lowerChangeBound || upperChangeBound) { |
---|
2107 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
2108 | double lower = columnLower_[iColumn]; |
---|
2109 | double upper = columnUpper_[iColumn]; |
---|
2110 | if (lower > upper) { |
---|
2111 | maxTheta = -1.0; |
---|
2112 | break; |
---|
2113 | } |
---|
2114 | double lowerChange = (lowerChangeBound) ? lowerChangeBound[iColumn] : 0.0; |
---|
2115 | double upperChange = (upperChangeBound) ? upperChangeBound[iColumn] : 0.0; |
---|
2116 | if (lower > -1.0e20 && upper < 1.0e20) { |
---|
2117 | if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) { |
---|
2118 | maxTheta = (upper - lower) / (lowerChange - upperChange); |
---|
2119 | } |
---|
2120 | } |
---|
2121 | if (lower > -1.0e20) { |
---|
2122 | lower_[iColumn] += startingTheta * lowerChange; |
---|
2123 | chgLower[iColumn] = lowerChange; |
---|
2124 | } |
---|
2125 | if (upper < 1.0e20) { |
---|
2126 | upper_[iColumn] += startingTheta * upperChange; |
---|
2127 | chgUpper[iColumn] = upperChange; |
---|
2128 | } |
---|
2129 | } |
---|
2130 | } |
---|
2131 | if (maxTheta == 1.0e50) |
---|
2132 | maxTheta = COIN_DBL_MAX; |
---|
2133 | } |
---|
2134 | if (maxTheta < 0.0) { |
---|
2135 | // bad ranges or initial |
---|
2136 | returnCode = -1; |
---|
2137 | } |
---|
2138 | if (maxTheta < endingTheta) { |
---|
2139 | char line[100]; |
---|
2140 | sprintf(line,"Crossover considerations reduce ending theta from %g to %g\n", |
---|
2141 | endingTheta,maxTheta); |
---|
2142 | handler_->message(CLP_GENERAL,messages_) |
---|
2143 | << line << CoinMessageEol; |
---|
2144 | endingTheta = maxTheta; |
---|
2145 | } |
---|
2146 | if (endingTheta < startingTheta) { |
---|
2147 | // bad initial |
---|
2148 | returnCode = -2; |
---|
2149 | } |
---|
2150 | } |
---|
2151 | double saveEndingTheta = endingTheta; |
---|
2152 | if (!returnCode) { |
---|
2153 | if (changeObjective) { |
---|
2154 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
2155 | chgObjective[iColumn] = changeObjective[iColumn]; |
---|
2156 | cost_[iColumn] += startingTheta * changeObjective[iColumn]; |
---|
2157 | } |
---|
2158 | } |
---|
2159 | double * saveDuals = NULL; |
---|
2160 | reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data); |
---|
2161 | assert (!problemStatus_); |
---|
2162 | for (int i=0;i<numberRows_+numberColumns_;i++) |
---|
2163 | setFakeBound(i, noFake); |
---|
2164 | // Now do parametrics |
---|
2165 | handler_->message(CLP_PARAMETRICS_STATS, messages_) |
---|
2166 | << startingTheta << objectiveValue() << CoinMessageEol; |
---|
2167 | while (!returnCode) { |
---|
2168 | //assert (reportIncrement); |
---|
2169 | parametricsData paramData; |
---|
2170 | paramData.startingTheta=startingTheta; |
---|
2171 | paramData.endingTheta=endingTheta; |
---|
2172 | paramData.maxTheta=COIN_DBL_MAX; |
---|
2173 | paramData.lowerChange = chgLower; |
---|
2174 | paramData.upperChange = chgUpper; |
---|
2175 | returnCode = parametricsLoop(paramData, reportIncrement, |
---|
2176 | chgLower, chgUpper, chgObjective, data, |
---|
2177 | canTryQuick); |
---|
2178 | startingTheta=paramData.startingTheta; |
---|
2179 | endingTheta=paramData.endingTheta; |
---|
2180 | if (!returnCode) { |
---|
2181 | //double change = endingTheta-startingTheta; |
---|
2182 | startingTheta = endingTheta; |
---|
2183 | endingTheta = saveEndingTheta; |
---|
2184 | //for (int i=0;i<numberTotal;i++) { |
---|
2185 | //lower_[i] += change*chgLower[i]; |
---|
2186 | //upper_[i] += change*chgUpper[i]; |
---|
2187 | //cost_[i] += change*chgObjective[i]; |
---|
2188 | //} |
---|
2189 | handler_->message(CLP_PARAMETRICS_STATS, messages_) |
---|
2190 | << startingTheta << objectiveValue() << CoinMessageEol; |
---|
2191 | if (startingTheta >= endingTheta) |
---|
2192 | break; |
---|
2193 | } else if (returnCode == -1) { |
---|
2194 | // trouble - do external solve |
---|
2195 | needToDoSomething = true; |
---|
2196 | } else if (problemStatus_==1) { |
---|
2197 | // can't move any further |
---|
2198 | if (!canTryQuick) { |
---|
2199 | handler_->message(CLP_PARAMETRICS_STATS, messages_) |
---|
2200 | << endingTheta << objectiveValue() << CoinMessageEol; |
---|
2201 | problemStatus_=0; |
---|
2202 | } |
---|
2203 | } else { |
---|
2204 | abort(); |
---|
2205 | } |
---|
2206 | } |
---|
2207 | } |
---|
2208 | reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0); |
---|
2209 | |
---|
2210 | delete dualRowPivot_; |
---|
2211 | dualRowPivot_ = savePivot; |
---|
2212 | // Restore any saved stuff |
---|
2213 | restoreData(data); |
---|
2214 | if (needToDoSomething) { |
---|
2215 | double saveStartingTheta = startingTheta; // known to be feasible |
---|
2216 | int cleanedUp = 1; |
---|
2217 | while (cleanedUp) { |
---|
2218 | // tweak |
---|
2219 | if (cleanedUp == 1) { |
---|
2220 | if (!reportIncrement) |
---|
2221 | startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta); |
---|
2222 | else |
---|
2223 | startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta); |
---|
2224 | } else { |
---|
2225 | // restoring to go slowly |
---|
2226 | startingTheta = saveStartingTheta; |
---|
2227 | } |
---|
2228 | // only works if not scaled |
---|
2229 | int i; |
---|
2230 | const double * obj1 = objective(); |
---|
2231 | double * obj2 = copyModel.objective(); |
---|
2232 | const double * lower1 = columnLower_; |
---|
2233 | double * lower2 = copyModel.columnLower(); |
---|
2234 | const double * upper1 = columnUpper_; |
---|
2235 | double * upper2 = copyModel.columnUpper(); |
---|
2236 | for (i = 0; i < numberColumns_; i++) { |
---|
2237 | obj2[i] = obj1[i] + startingTheta * chgObjective[i]; |
---|
2238 | lower2[i] = lower1[i] + startingTheta * chgLower[i]; |
---|
2239 | upper2[i] = upper1[i] + startingTheta * chgUpper[i]; |
---|
2240 | } |
---|
2241 | lower1 = rowLower_; |
---|
2242 | lower2 = copyModel.rowLower(); |
---|
2243 | upper1 = rowUpper_; |
---|
2244 | upper2 = copyModel.rowUpper(); |
---|
2245 | for (i = 0; i < numberRows_; i++) { |
---|
2246 | lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_]; |
---|
2247 | upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_]; |
---|
2248 | } |
---|
2249 | copyModel.dual(); |
---|
2250 | if (copyModel.problemStatus()) { |
---|
2251 | char line[100]; |
---|
2252 | sprintf(line,"Can not get to theta of %g\n", startingTheta); |
---|
2253 | handler_->message(CLP_GENERAL,messages_) |
---|
2254 | << line << CoinMessageEol; |
---|
2255 | canTryQuick = false; // do slowly to get exact amount |
---|
2256 | // back to last known good |
---|
2257 | if (cleanedUp == 1) |
---|
2258 | cleanedUp = 2; |
---|
2259 | else |
---|
2260 | abort(); |
---|
2261 | } else { |
---|
2262 | // and move stuff back |
---|
2263 | int numberTotal = numberRows_ + numberColumns_; |
---|
2264 | CoinMemcpyN(copyModel.statusArray(), numberTotal, status_); |
---|
2265 | CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_); |
---|
2266 | CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_); |
---|
2267 | cleanedUp = 0; |
---|
2268 | } |
---|
2269 | } |
---|
2270 | } |
---|
2271 | delete [] chgLower; |
---|
2272 | delete [] chgUpper; |
---|
2273 | delete [] chgObjective; |
---|
2274 | } |
---|
2275 | perturbation_ = savePerturbation; |
---|
2276 | char line[100]; |
---|
2277 | sprintf(line,"Ending theta %g\n", endingTheta); |
---|
2278 | handler_->message(CLP_GENERAL,messages_) |
---|
2279 | << line << CoinMessageEol; |
---|
2280 | return problemStatus_; |
---|
2281 | } |
---|
2282 | /* Version of parametrics which reads from file |
---|
2283 | See CbcClpParam.cpp for details of format |
---|
2284 | Returns -2 if unable to open file */ |
---|
2285 | int |
---|
2286 | ClpSimplexOther::parametrics(const char * dataFile) |
---|
2287 | { |
---|
2288 | int returnCode=-2; |
---|
2289 | FILE *fp = fopen(dataFile, "r"); |
---|
2290 | char line[200]; |
---|
2291 | if (!fp) { |
---|
2292 | handler_->message(CLP_UNABLE_OPEN, messages_) |
---|
2293 | << dataFile << CoinMessageEol; |
---|
2294 | return -2; |
---|
2295 | } |
---|
2296 | |
---|
2297 | if (!fgets(line, 200, fp)) { |
---|
2298 | sprintf(line,"Empty parametrics file %s?",dataFile); |
---|
2299 | handler_->message(CLP_GENERAL,messages_) |
---|
2300 | << line << CoinMessageEol; |
---|
2301 | fclose(fp); |
---|
2302 | return -2; |
---|
2303 | } |
---|
2304 | char * pos = line; |
---|
2305 | char * put = line; |
---|
2306 | while (*pos >= ' ' && *pos != '\n') { |
---|
2307 | if (*pos != ' ' && *pos != '\t') { |
---|
2308 | *put = static_cast<char>(tolower(*pos)); |
---|
2309 | put++; |
---|
2310 | } |
---|
2311 | pos++; |
---|
2312 | } |
---|
2313 | *put = '\0'; |
---|
2314 | pos = line; |
---|
2315 | double startTheta=0.0; |
---|
2316 | double endTheta=0.0; |
---|
2317 | double intervalTheta=COIN_DBL_MAX; |
---|
2318 | int detail=0; |
---|
2319 | bool good = true; |
---|
2320 | while (good) { |
---|
2321 | good=false; |
---|
2322 | // check ROWS |
---|
2323 | char * comma = strchr(pos, ','); |
---|
2324 | if (!comma) |
---|
2325 | break; |
---|
2326 | *comma = '\0'; |
---|
2327 | if (strcmp(pos,"rows")) |
---|
2328 | break; |
---|
2329 | *comma = ','; |
---|
2330 | pos = comma+1; |
---|
2331 | // check lower theta |
---|
2332 | comma = strchr(pos, ','); |
---|
2333 | if (!comma) |
---|
2334 | break; |
---|
2335 | *comma = '\0'; |
---|
2336 | startTheta = atof(pos); |
---|
2337 | *comma = ','; |
---|
2338 | pos = comma+1; |
---|
2339 | // check upper theta |
---|
2340 | comma = strchr(pos, ','); |
---|
2341 | good=true; |
---|
2342 | if (comma) |
---|
2343 | *comma = '\0'; |
---|
2344 | endTheta = atof(pos); |
---|
2345 | if (comma) { |
---|
2346 | *comma = ','; |
---|
2347 | pos = comma+1; |
---|
2348 | comma = strchr(pos, ','); |
---|
2349 | if (comma) |
---|
2350 | *comma = '\0'; |
---|
2351 | intervalTheta = atof(pos); |
---|
2352 | if (comma) { |
---|
2353 | *comma = ','; |
---|
2354 | pos = comma+1; |
---|
2355 | comma = strchr(pos, ','); |
---|
2356 | if (comma) |
---|
2357 | *comma = '\0'; |
---|
2358 | detail = atoi(pos); |
---|
2359 | if (comma) |
---|
2360 | *comma = ','; |
---|
2361 | } |
---|
2362 | } |
---|
2363 | break; |
---|
2364 | } |
---|
2365 | if (good) { |
---|
2366 | if (startTheta<0.0|| |
---|
2367 | startTheta>endTheta|| |
---|
2368 | intervalTheta<0.0) |
---|
2369 | good=false; |
---|
2370 | if (detail<0||detail>1) |
---|
2371 | good=false; |
---|
2372 | } |
---|
2373 | if (intervalTheta>=endTheta) |
---|
2374 | intervalTheta=0.0; |
---|
2375 | if (!good) { |
---|
2376 | sprintf(line,"Odd first line %s on file %s?",line,dataFile); |
---|
2377 | handler_->message(CLP_GENERAL,messages_) |
---|
2378 | << line << CoinMessageEol; |
---|
2379 | fclose(fp); |
---|
2380 | return -2; |
---|
2381 | } |
---|
2382 | if (!fgets(line, 200, fp)) { |
---|
2383 | sprintf(line,"Not enough records on parametrics file %s?",dataFile); |
---|
2384 | handler_->message(CLP_GENERAL,messages_) |
---|
2385 | << line << CoinMessageEol; |
---|
2386 | fclose(fp); |
---|
2387 | return -2; |
---|
2388 | } |
---|
2389 | double * lowerRowMove = NULL; |
---|
2390 | double * upperRowMove = NULL; |
---|
2391 | double * lowerColumnMove = NULL; |
---|
2392 | double * upperColumnMove = NULL; |
---|
2393 | double * objectiveMove = NULL; |
---|
2394 | char saveLine[200]; |
---|
2395 | saveLine[0]='\0'; |
---|
2396 | std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"}; |
---|
2397 | int gotRow[] = { -1, -1, -1, -1, -1}; |
---|
2398 | int orderRow[5]; |
---|
2399 | assert(sizeof(gotRow) == sizeof(orderRow)); |
---|
2400 | int nAcross = 0; |
---|
2401 | pos = line; |
---|
2402 | put = line; |
---|
2403 | while (*pos >= ' ' && *pos != '\n') { |
---|
2404 | if (*pos != ' ' && *pos != '\t') { |
---|
2405 | *put = static_cast<char>(tolower(*pos)); |
---|
2406 | put++; |
---|
2407 | } |
---|
2408 | pos++; |
---|
2409 | } |
---|
2410 | *put = '\0'; |
---|
2411 | pos = line; |
---|
2412 | int i; |
---|
2413 | good = true; |
---|
2414 | if (strncmp(line,"column",6)) { |
---|
2415 | while (pos) { |
---|
2416 | char * comma = strchr(pos, ','); |
---|
2417 | if (comma) |
---|
2418 | *comma = '\0'; |
---|
2419 | for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) { |
---|
2420 | if (headingsRow[i] == pos) { |
---|
2421 | if (gotRow[i] < 0) { |
---|
2422 | orderRow[nAcross] = i; |
---|
2423 | gotRow[i] = nAcross++; |
---|
2424 | } else { |
---|
2425 | // duplicate |
---|
2426 | good = false; |
---|
2427 | } |
---|
2428 | break; |
---|
2429 | } |
---|
2430 | } |
---|
2431 | if (i == static_cast<int> (sizeof(gotRow) / sizeof(int))) |
---|
2432 | good = false; |
---|
2433 | if (comma) { |
---|
2434 | *comma = ','; |
---|
2435 | pos = comma + 1; |
---|
2436 | } else { |
---|
2437 | break; |
---|
2438 | } |
---|
2439 | } |
---|
2440 | if (gotRow[0] < 0 && gotRow[1] < 0) |
---|
2441 | good = false; |
---|
2442 | if (gotRow[0] >= 0 && gotRow[1] >= 0) |
---|
2443 | good = false; |
---|
2444 | if (gotRow[0] >= 0 && !lengthNames()) |
---|
2445 | good = false; |
---|
2446 | if (gotRow[4]<0) { |
---|
2447 | if (gotRow[2] < 0 && gotRow[3] >= 0) |
---|
2448 | good = false; |
---|
2449 | else if (gotRow[3] < 0 && gotRow[2] >= 0) |
---|
2450 | good = false; |
---|
2451 | } else if (gotRow[2]>=0||gotRow[3]>=0) { |
---|
2452 | good = false; |
---|
2453 | } |
---|
2454 | if (good) { |
---|
2455 | char ** rowNames = new char * [numberRows_]; |
---|
2456 | int iRow; |
---|
2457 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2458 | rowNames[iRow] = |
---|
2459 | CoinStrdup(rowName(iRow).c_str()); |
---|
2460 | } |
---|
2461 | lowerRowMove = new double [numberRows_]; |
---|
2462 | memset(lowerRowMove,0,numberRows_*sizeof(double)); |
---|
2463 | upperRowMove = new double [numberRows_]; |
---|
2464 | memset(upperRowMove,0,numberRows_*sizeof(double)); |
---|
2465 | int nLine = 0; |
---|
2466 | int nBadLine = 0; |
---|
2467 | int nBadName = 0; |
---|
2468 | while (fgets(line, 200, fp)) { |
---|
2469 | if (!strncmp(line, "ENDATA", 6)|| |
---|
2470 | !strncmp(line, "COLUMN",6)) |
---|
2471 | break; |
---|
2472 | nLine++; |
---|
2473 | iRow = -1; |
---|
2474 | double upper = 0.0; |
---|
2475 | double lower = 0.0; |
---|
2476 | char * pos = line; |
---|
2477 | char * put = line; |
---|
2478 | while (*pos >= ' ' && *pos != '\n') { |
---|
2479 | if (*pos != ' ' && *pos != '\t') { |
---|
2480 | *put = *pos; |
---|
2481 | put++; |
---|
2482 | } |
---|
2483 | pos++; |
---|
2484 | } |
---|
2485 | *put = '\0'; |
---|
2486 | pos = line; |
---|
2487 | for (int i = 0; i < nAcross; i++) { |
---|
2488 | char * comma = strchr(pos, ','); |
---|
2489 | if (comma) { |
---|
2490 | *comma = '\0'; |
---|
2491 | } else if (i < nAcross - 1) { |
---|
2492 | nBadLine++; |
---|
2493 | break; |
---|
2494 | } |
---|
2495 | switch (orderRow[i]) { |
---|
2496 | // name |
---|
2497 | case 0: |
---|
2498 | // For large problems this could be slow |
---|
2499 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2500 | if (!strcmp(rowNames[iRow], pos)) |
---|
2501 | break; |
---|
2502 | } |
---|
2503 | if (iRow == numberRows_) |
---|
2504 | iRow = -1; |
---|
2505 | break; |
---|
2506 | // number |
---|
2507 | case 1: |
---|
2508 | iRow = atoi(pos); |
---|
2509 | if (iRow < 0 || iRow >= numberRows_) |
---|
2510 | iRow = -1; |
---|
2511 | break; |
---|
2512 | // lower |
---|
2513 | case 2: |
---|
2514 | upper = atof(pos); |
---|
2515 | break; |
---|
2516 | // upper |
---|
2517 | case 3: |
---|
2518 | lower = atof(pos); |
---|
2519 | break; |
---|
2520 | // rhs |
---|
2521 | case 4: |
---|
2522 | lower = atof(pos); |
---|
2523 | upper = lower; |
---|
2524 | break; |
---|
2525 | } |
---|
2526 | if (comma) { |
---|
2527 | *comma = ','; |
---|
2528 | pos = comma + 1; |
---|
2529 | } |
---|
2530 | } |
---|
2531 | if (iRow >= 0) { |
---|
2532 | if (rowLower_[iRow]>-1.0e20) |
---|
2533 | lowerRowMove[iRow] = lower; |
---|
2534 | else |
---|
2535 | lowerRowMove[iRow]=0.0; |
---|
2536 | if (rowUpper_[iRow]<1.0e20) |
---|
2537 | upperRowMove[iRow] = upper; |
---|
2538 | else |
---|
2539 | upperRowMove[iRow] = lower; |
---|
2540 | } else { |
---|
2541 | nBadName++; |
---|
2542 | if(saveLine[0]=='\0') |
---|
2543 | strcpy(saveLine,line); |
---|
2544 | } |
---|
2545 | } |
---|
2546 | sprintf(line,"%d Row fields and %d records", nAcross, nLine); |
---|
2547 | handler_->message(CLP_GENERAL,messages_) |
---|
2548 | << line << CoinMessageEol; |
---|
2549 | if (nBadName) { |
---|
2550 | sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine); |
---|
2551 | handler_->message(CLP_GENERAL,messages_) |
---|
2552 | << line << CoinMessageEol; |
---|
2553 | returnCode=-1; |
---|
2554 | good=false; |
---|
2555 | } |
---|
2556 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2557 | free(rowNames[iRow]); |
---|
2558 | } |
---|
2559 | delete [] rowNames; |
---|
2560 | } else { |
---|
2561 | sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong"); |
---|
2562 | handler_->message(CLP_GENERAL,messages_) |
---|
2563 | << line << CoinMessageEol; |
---|
2564 | returnCode=-1; |
---|
2565 | good=false; |
---|
2566 | } |
---|
2567 | } |
---|
2568 | if (good&&(!strncmp(line, "COLUMN",6)||!strncmp(line, "column",6))) { |
---|
2569 | if (!fgets(line, 200, fp)) { |
---|
2570 | sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile); |
---|
2571 | handler_->message(CLP_GENERAL,messages_) |
---|
2572 | << line << CoinMessageEol; |
---|
2573 | fclose(fp); |
---|
2574 | return -2; |
---|
2575 | } |
---|
2576 | std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"}; |
---|
2577 | saveLine[0]='\0'; |
---|
2578 | int gotColumn[] = { -1, -1, -1, -1, -1}; |
---|
2579 | int orderColumn[5]; |
---|
2580 | assert(sizeof(gotColumn) == sizeof(orderColumn)); |
---|
2581 | nAcross = 0; |
---|
2582 | pos = line; |
---|
2583 | put = line; |
---|
2584 | while (*pos >= ' ' && *pos != '\n') { |
---|
2585 | if (*pos != ' ' && *pos != '\t') { |
---|
2586 | *put = static_cast<char>(tolower(*pos)); |
---|
2587 | put++; |
---|
2588 | } |
---|
2589 | pos++; |
---|
2590 | } |
---|
2591 | *put = '\0'; |
---|
2592 | pos = line; |
---|
2593 | int i; |
---|
2594 | if (strncmp(line,"endata",6)&&good) { |
---|
2595 | while (pos) { |
---|
2596 | char * comma = strchr(pos, ','); |
---|
2597 | if (comma) |
---|
2598 | *comma = '\0'; |
---|
2599 | for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) { |
---|
2600 | if (headingsColumn[i] == pos) { |
---|
2601 | if (gotColumn[i] < 0) { |
---|
2602 | orderColumn[nAcross] = i; |
---|
2603 | gotColumn[i] = nAcross++; |
---|
2604 | } else { |
---|
2605 | // duplicate |
---|
2606 | good = false; |
---|
2607 | } |
---|
2608 | break; |
---|
2609 | } |
---|
2610 | } |
---|
2611 | if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int))) |
---|
2612 | good = false; |
---|
2613 | if (comma) { |
---|
2614 | *comma = ','; |
---|
2615 | pos = comma + 1; |
---|
2616 | } else { |
---|
2617 | break; |
---|
2618 | } |
---|
2619 | } |
---|
2620 | if (gotColumn[0] < 0 && gotColumn[1] < 0) |
---|
2621 | good = false; |
---|
2622 | if (gotColumn[0] >= 0 && gotColumn[1] >= 0) |
---|
2623 | good = false; |
---|
2624 | if (gotColumn[0] >= 0 && !lengthNames()) |
---|
2625 | good = false; |
---|
2626 | if (good) { |
---|
2627 | char ** columnNames = new char * [numberColumns_]; |
---|
2628 | int iColumn; |
---|
2629 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
2630 | columnNames[iColumn] = |
---|
2631 | CoinStrdup(columnName(iColumn).c_str()); |
---|
2632 | } |
---|
2633 | lowerColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double))); |
---|
2634 | memset(lowerColumnMove,0,numberColumns_*sizeof(double)); |
---|
2635 | upperColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double))); |
---|
2636 | memset(upperColumnMove,0,numberColumns_*sizeof(double)); |
---|
2637 | objectiveMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double))); |
---|
2638 | memset(objectiveMove,0,numberColumns_*sizeof(double)); |
---|
2639 | int nLine = 0; |
---|
2640 | int nBadLine = 0; |
---|
2641 | int nBadName = 0; |
---|
2642 | while (fgets(line, 200, fp)) { |
---|
2643 | if (!strncmp(line, "ENDATA", 6)) |
---|
2644 | break; |
---|
2645 | nLine++; |
---|
2646 | iColumn = -1; |
---|
2647 | double upper = 0.0; |
---|
2648 | double lower = 0.0; |
---|
2649 | double obj =0.0; |
---|
2650 | char * pos = line; |
---|
2651 | char * put = line; |
---|
2652 | while (*pos >= ' ' && *pos != '\n') { |
---|
2653 | if (*pos != ' ' && *pos != '\t') { |
---|
2654 | *put = *pos; |
---|
2655 | put++; |
---|
2656 | } |
---|
2657 | pos++; |
---|
2658 | } |
---|
2659 | *put = '\0'; |
---|
2660 | pos = line; |
---|
2661 | for (int i = 0; i < nAcross; i++) { |
---|
2662 | char * comma = strchr(pos, ','); |
---|
2663 | if (comma) { |
---|
2664 | *comma = '\0'; |
---|
2665 | } else if (i < nAcross - 1) { |
---|
2666 | nBadLine++; |
---|
2667 | break; |
---|
2668 | } |
---|
2669 | switch (orderColumn[i]) { |
---|
2670 | // name |
---|
2671 | case 0: |
---|
2672 | // For large problems this could be slow |
---|
2673 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
2674 | if (!strcmp(columnNames[iColumn], pos)) |
---|
2675 | break; |
---|
2676 | } |
---|
2677 | if (iColumn == numberColumns_) |
---|
2678 | iColumn = -1; |
---|
2679 | break; |
---|
2680 | // number |
---|
2681 | case 1: |
---|
2682 | iColumn = atoi(pos); |
---|
2683 | if (iColumn < 0 || iColumn >= numberColumns_) |
---|
2684 | iColumn = -1; |
---|
2685 | break; |
---|
2686 | // lower |
---|
2687 | case 2: |
---|
2688 | upper = atof(pos); |
---|
2689 | break; |
---|
2690 | // upper |
---|
2691 | case 3: |
---|
2692 | lower = atof(pos); |
---|
2693 | break; |
---|
2694 | // objective |
---|
2695 | case 4: |
---|
2696 | obj = atof(pos); |
---|
2697 | upper = lower; |
---|
2698 | break; |
---|
2699 | } |
---|
2700 | if (comma) { |
---|
2701 | *comma = ','; |
---|
2702 | pos = comma + 1; |
---|
2703 | } |
---|
2704 | } |
---|
2705 | if (iColumn >= 0) { |
---|
2706 | if (columnLower_[iColumn]>-1.0e20) |
---|
2707 | lowerColumnMove[iColumn] = lower; |
---|
2708 | else |
---|
2709 | lowerColumnMove[iColumn]=0.0; |
---|
2710 | if (columnUpper_[iColumn]<1.0e20) |
---|
2711 | upperColumnMove[iColumn] = upper; |
---|
2712 | else |
---|
2713 | upperColumnMove[iColumn] = lower; |
---|
2714 | objectiveMove[iColumn] = obj; |
---|
2715 | } else { |
---|
2716 | nBadName++; |
---|
2717 | if(saveLine[0]=='\0') |
---|
2718 | strcpy(saveLine,line); |
---|
2719 | } |
---|
2720 | } |
---|
2721 | sprintf(line,"%d Column fields and %d records", nAcross, nLine); |
---|
2722 | handler_->message(CLP_GENERAL,messages_) |
---|
2723 | << line << CoinMessageEol; |
---|
2724 | if (nBadName) { |
---|
2725 | sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine); |
---|
2726 | handler_->message(CLP_GENERAL,messages_) |
---|
2727 | << line << CoinMessageEol; |
---|
2728 | returnCode=-1; |
---|
2729 | good=false; |
---|
2730 | } |
---|
2731 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
2732 | free(columnNames[iColumn]); |
---|
2733 | } |
---|
2734 | delete [] columnNames; |
---|
2735 | } else { |
---|
2736 | sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong"); |
---|
2737 | handler_->message(CLP_GENERAL,messages_) |
---|
2738 | << line << CoinMessageEol; |
---|
2739 | returnCode=-1; |
---|
2740 | good=false; |
---|
2741 | } |
---|
2742 | } |
---|
2743 | } |
---|
2744 | returnCode=-1; |
---|
2745 | if (good) { |
---|
2746 | // clean arrays |
---|
2747 | if (lowerRowMove) { |
---|
2748 | bool empty=true; |
---|
2749 | for (int i=0;i<numberRows_;i++) { |
---|
2750 | if (lowerRowMove[i]) { |
---|
2751 | empty=false; |
---|
2752 | break; |
---|
2753 | } |
---|
2754 | } |
---|
2755 | if (empty) { |
---|
2756 | delete [] lowerRowMove; |
---|
2757 | lowerRowMove=NULL; |
---|
2758 | } |
---|
2759 | } |
---|
2760 | if (upperRowMove) { |
---|
2761 | bool empty=true; |
---|
2762 | for (int i=0;i<numberRows_;i++) { |
---|
2763 | if (upperRowMove[i]) { |
---|
2764 | empty=false; |
---|
2765 | break; |
---|
2766 | } |
---|
2767 | } |
---|
2768 | if (empty) { |
---|
2769 | delete [] upperRowMove; |
---|
2770 | upperRowMove=NULL; |
---|
2771 | } |
---|
2772 | } |
---|
2773 | if (lowerColumnMove) { |
---|
2774 | bool empty=true; |
---|
2775 | for (int i=0;i<numberColumns_;i++) { |
---|
2776 | if (lowerColumnMove[i]) { |
---|
2777 | empty=false; |
---|
2778 | break; |
---|
2779 | } |
---|
2780 | } |
---|
2781 | if (empty) { |
---|
2782 | delete [] lowerColumnMove; |
---|
2783 | lowerColumnMove=NULL; |
---|
2784 | } |
---|
2785 | } |
---|
2786 | if (upperColumnMove) { |
---|
2787 | bool empty=true; |
---|
2788 | for (int i=0;i<numberColumns_;i++) { |
---|
2789 | if (upperColumnMove[i]) { |
---|
2790 | empty=false; |
---|
2791 | break; |
---|
2792 | } |
---|
2793 | } |
---|
2794 | if (empty) { |
---|
2795 | delete [] upperColumnMove; |
---|
2796 | upperColumnMove=NULL; |
---|
2797 | } |
---|
2798 | } |
---|
2799 | if (objectiveMove) { |
---|
2800 | bool empty=true; |
---|
2801 | for (int i=0;i<numberColumns_;i++) { |
---|
2802 | if (objectiveMove[i]) { |
---|
2803 | empty=false; |
---|
2804 | break; |
---|
2805 | } |
---|
2806 | } |
---|
2807 | if (empty) { |
---|
2808 | delete [] objectiveMove; |
---|
2809 | objectiveMove=NULL; |
---|
2810 | } |
---|
2811 | } |
---|
2812 | int saveScaling = scalingFlag_; |
---|
2813 | scalingFlag_ = 0; |
---|
2814 | int saveLogLevel = handler_->logLevel(); |
---|
2815 | if (detail>0&&!intervalTheta) |
---|
2816 | handler_->setLogLevel(3); |
---|
2817 | else |
---|
2818 | handler_->setLogLevel(1); |
---|
2819 | returnCode = parametrics(startTheta,endTheta,intervalTheta, |
---|
2820 | lowerColumnMove,upperColumnMove, |
---|
2821 | lowerRowMove,upperRowMove, |
---|
2822 | objectiveMove); |
---|
2823 | scalingFlag_ = saveScaling; |
---|
2824 | handler_->setLogLevel(saveLogLevel); |
---|
2825 | } |
---|
2826 | delete [] lowerRowMove; |
---|
2827 | delete [] upperRowMove; |
---|
2828 | delete [] lowerColumnMove; |
---|
2829 | delete [] upperColumnMove; |
---|
2830 | delete [] objectiveMove; |
---|
2831 | fclose(fp); |
---|
2832 | return returnCode; |
---|
2833 | } |
---|
2834 | int |
---|
2835 | ClpSimplexOther::parametricsLoop(parametricsData & paramData,double reportIncrement, |
---|
2836 | const double * lowerChange, const double * upperChange, |
---|
2837 | const double * changeObjective, ClpDataSave & data, |
---|
2838 | bool canTryQuick) |
---|
2839 | { |
---|
2840 | double startingTheta = paramData.startingTheta; |
---|
2841 | double & endingTheta = paramData.endingTheta; |
---|
2842 | // stuff is already at starting |
---|
2843 | // For this crude version just try and go to end |
---|
2844 | double change = 0.0; |
---|
2845 | if (reportIncrement && canTryQuick) { |
---|
2846 | endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement); |
---|
2847 | change = endingTheta - startingTheta; |
---|
2848 | } |
---|
2849 | int numberTotal = numberRows_ + numberColumns_; |
---|
2850 | int i; |
---|
2851 | for ( i = 0; i < numberTotal; i++) { |
---|
2852 | lower_[i] += change * lowerChange[i]; |
---|
2853 | upper_[i] += change * upperChange[i]; |
---|
2854 | switch(getStatus(i)) { |
---|
2855 | |
---|
2856 | case basic: |
---|
2857 | case isFree: |
---|
2858 | case superBasic: |
---|
2859 | break; |
---|
2860 | case isFixed: |
---|
2861 | case atUpperBound: |
---|
2862 | solution_[i] = upper_[i]; |
---|
2863 | break; |
---|
2864 | case atLowerBound: |
---|
2865 | solution_[i] = lower_[i]; |
---|
2866 | break; |
---|
2867 | } |
---|
2868 | cost_[i] += change * changeObjective[i]; |
---|
2869 | } |
---|
2870 | problemStatus_ = -1; |
---|
2871 | |
---|
2872 | // This says whether to restore things etc |
---|
2873 | // startup will have factorized so can skip |
---|
2874 | int factorType = 0; |
---|
2875 | // Start check for cycles |
---|
2876 | progress_.startCheck(); |
---|
2877 | // Say change made on first iteration |
---|
2878 | changeMade_ = 1; |
---|
2879 | /* |
---|
2880 | Status of problem: |
---|
2881 | 0 - optimal |
---|
2882 | 1 - infeasible |
---|
2883 | 2 - unbounded |
---|
2884 | -1 - iterating |
---|
2885 | -2 - factorization wanted |
---|
2886 | -3 - redo checking without factorization |
---|
2887 | -4 - looks infeasible |
---|
2888 | */ |
---|
2889 | while (problemStatus_ < 0) { |
---|
2890 | int iRow, iColumn; |
---|
2891 | // clear |
---|
2892 | for (iRow = 0; iRow < 4; iRow++) { |
---|
2893 | rowArray_[iRow]->clear(); |
---|
2894 | } |
---|
2895 | |
---|
2896 | for (iColumn = 0; iColumn < 2; iColumn++) { |
---|
2897 | columnArray_[iColumn]->clear(); |
---|
2898 | } |
---|
2899 | |
---|
2900 | // give matrix (and model costs and bounds a chance to be |
---|
2901 | // refreshed (normally null) |
---|
2902 | matrix_->refresh(this); |
---|
2903 | // may factorize, checks if problem finished |
---|
2904 | statusOfProblemInParametrics(factorType, data); |
---|
2905 | // Say good factorization |
---|
2906 | factorType = 1; |
---|
2907 | if (data.sparseThreshold_) { |
---|
2908 | // use default at present |
---|
2909 | factorization_->sparseThreshold(0); |
---|
2910 | factorization_->goSparse(); |
---|
2911 | } |
---|
2912 | |
---|
2913 | // exit if victory declared |
---|
2914 | if (problemStatus_ >= 0 && |
---|
2915 | (canTryQuick || startingTheta>=endingTheta-1.0e-7) ) |
---|
2916 | break; |
---|
2917 | |
---|
2918 | // test for maximum iterations |
---|
2919 | if (hitMaximumIterations()) { |
---|
2920 | problemStatus_ = 3; |
---|
2921 | break; |
---|
2922 | } |
---|
2923 | // Check event |
---|
2924 | { |
---|
2925 | int status = eventHandler_->event(ClpEventHandler::endOfFactorization); |
---|
2926 | if (status >= 0) { |
---|
2927 | problemStatus_ = 5; |
---|
2928 | secondaryStatus_ = ClpEventHandler::endOfFactorization; |
---|
2929 | break; |
---|
2930 | } |
---|
2931 | } |
---|
2932 | // Do iterations |
---|
2933 | problemStatus_=-1; |
---|
2934 | if (canTryQuick) { |
---|
2935 | double * saveDuals = NULL; |
---|
2936 | reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0); |
---|
2937 | } else { |
---|
2938 | whileIterating(paramData, reportIncrement, |
---|
2939 | changeObjective); |
---|
2940 | startingTheta = endingTheta; |
---|
2941 | } |
---|
2942 | } |
---|
2943 | if (!problemStatus_) { |
---|
2944 | theta_ = change + startingTheta; |
---|
2945 | eventHandler_->event(ClpEventHandler::theta); |
---|
2946 | return 0; |
---|
2947 | } else if (problemStatus_ == 10) { |
---|
2948 | return -1; |
---|
2949 | } else { |
---|
2950 | return problemStatus_; |
---|
2951 | } |
---|
2952 | } |
---|
2953 | /* Parametrics |
---|
2954 | The code uses current bounds + theta * change (if change array not NULL) |
---|
2955 | It starts at startingTheta and returns ending theta in endingTheta. |
---|
2956 | If it can not reach input endingTheta return code will be 1 for infeasible, |
---|
2957 | 2 for unbounded, if error on ranges -1, otherwise 0. |
---|
2958 | Event handler may do more |
---|
2959 | On exit endingTheta is maximum reached (can be used for next startingTheta) |
---|
2960 | */ |
---|
2961 | int |
---|
2962 | ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, |
---|
2963 | const double * lowerChangeBound, const double * upperChangeBound, |
---|
2964 | const double * lowerChangeRhs, const double * upperChangeRhs) |
---|
2965 | { |
---|
2966 | int savePerturbation = perturbation_; |
---|
2967 | perturbation_ = 102; // switch off |
---|
2968 | algorithm_ = -1; |
---|
2969 | // extra region |
---|
2970 | int maximumPivots = factorization_->maximumPivots(); |
---|
2971 | int numberDense = factorization_->numberDense(); |
---|
2972 | int length = numberRows_ + numberDense + maximumPivots; |
---|
2973 | assert (!rowArray_[4]); |
---|
2974 | rowArray_[4]=new CoinIndexedVector(length); |
---|
2975 | assert (!rowArray_[5]); |
---|
2976 | rowArray_[5]=new CoinIndexedVector(length); |
---|
2977 | |
---|
2978 | // save data |
---|
2979 | ClpDataSave data = saveData(); |
---|
2980 | int numberTotal = numberRows_ + numberColumns_; |
---|
2981 | int ratio = static_cast<int>((2*sizeof(int))/sizeof(double)); |
---|
2982 | assert (ratio==1||ratio==2); |
---|
2983 | // allow for unscaled - even if not needed |
---|
2984 | int lengthArrays = 4*numberTotal+(3*numberTotal+2)*ratio+2*numberRows_+1; |
---|
2985 | int unscaledChangesOffset=lengthArrays; // need extra copy of change |
---|
2986 | lengthArrays += numberTotal; |
---|
2987 | |
---|
2988 | /* |
---|
2989 | Save information and modify |
---|
2990 | */ |
---|
2991 | double * saveLower = new double [lengthArrays]; |
---|
2992 | double * saveUpper = new double [lengthArrays]; |
---|
2993 | double * lowerCopy = saveLower+2*numberTotal; |
---|
2994 | double * upperCopy = saveUpper+2*numberTotal; |
---|
2995 | double * lowerChange = saveLower+numberTotal; |
---|
2996 | double * upperChange = saveUpper+numberTotal; |
---|
2997 | double * lowerGap = saveLower+4*numberTotal; |
---|
2998 | double * lowerCoefficient = lowerGap+numberRows_; |
---|
2999 | double * upperGap = saveUpper+4*numberTotal; |
---|
3000 | double * upperCoefficient = upperGap+numberRows_; |
---|
3001 | int * lowerList = (reinterpret_cast<int *>(saveLower+4*numberTotal+2*numberRows_))+2; |
---|
3002 | int * upperList = (reinterpret_cast<int *>(saveUpper+4*numberTotal+2*numberRows_))+2; |
---|
3003 | int * lowerActive = lowerList+numberTotal+1; |
---|
3004 | int * upperActive = upperList+numberTotal+1; |
---|
3005 | // To mark as odd |
---|
3006 | char * markDone = reinterpret_cast<char *>(lowerActive+numberTotal); |
---|
3007 | //memset(markDone,0,numberTotal); |
---|
3008 | int * backwardBasic = upperActive+numberTotal; |
---|
3009 | parametricsData paramData; |
---|
3010 | paramData.lowerChange = lowerChange; |
---|
3011 | paramData.lowerList=lowerList; |
---|
3012 | paramData.upperChange = upperChange; |
---|
3013 | paramData.upperList=upperList; |
---|
3014 | paramData.markDone=markDone; |
---|
3015 | paramData.backwardBasic=backwardBasic; |
---|
3016 | paramData.lowerActive = lowerActive; |
---|
3017 | paramData.lowerGap = lowerGap; |
---|
3018 | paramData.lowerCoefficient = lowerCoefficient; |
---|
3019 | paramData.upperActive = upperActive; |
---|
3020 | paramData.upperGap = upperGap; |
---|
3021 | paramData.upperCoefficient = upperCoefficient; |
---|
3022 | paramData.unscaledChangesOffset = unscaledChangesOffset-numberTotal; |
---|
3023 | paramData.firstIteration=true; |
---|
3024 | // Find theta when bounds will cross over and create arrays |
---|
3025 | memset(lowerChange, 0, numberTotal * sizeof(double)); |
---|
3026 | memset(upperChange, 0, numberTotal * sizeof(double)); |
---|
3027 | if (lowerChangeBound) |
---|
3028 | memcpy(lowerChange,lowerChangeBound,numberColumns_*sizeof(double)); |
---|
3029 | if (upperChangeBound) |
---|
3030 | memcpy(upperChange,upperChangeBound,numberColumns_*sizeof(double)); |
---|
3031 | if (lowerChangeRhs) |
---|
3032 | memcpy(lowerChange+numberColumns_, |
---|
3033 | lowerChangeRhs,numberRows_*sizeof(double)); |
---|
3034 | if (upperChangeRhs) |
---|
3035 | memcpy(upperChange+numberColumns_, |
---|
3036 | upperChangeRhs,numberRows_*sizeof(double)); |
---|
3037 | // clean |
---|
3038 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
3039 | double lower = rowLower_[iRow]; |
---|
3040 | double upper = rowUpper_[iRow]; |
---|
3041 | if (lower<-1.0e30) |
---|
3042 | lowerChange[numberColumns_+iRow]=0.0; |
---|
3043 | if (upper>1.0e30) |
---|
3044 | upperChange[numberColumns_+iRow]=0.0; |
---|
3045 | } |
---|
3046 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
3047 | double lower = columnLower_[iColumn]; |
---|
3048 | double upper = columnUpper_[iColumn]; |
---|
3049 | if (lower<-1.0e30) |
---|
3050 | lowerChange[iColumn]=0.0; |
---|
3051 | if (upper>1.0e30) |
---|
3052 | upperChange[iColumn]=0.0; |
---|
3053 | } |
---|
3054 | // save unscaled version of changes |
---|
3055 | memcpy(saveLower+unscaledChangesOffset,lowerChange,numberTotal*sizeof(double)); |
---|
3056 | memcpy(saveUpper+unscaledChangesOffset,upperChange,numberTotal*sizeof(double)); |
---|
3057 | int nLowerChange=0; |
---|
3058 | int nUpperChange=0; |
---|
3059 | for (int i=0;i<numberColumns_;i++) { |
---|
3060 | if (lowerChange[i]) { |
---|
3061 | lowerList[nLowerChange++]=i; |
---|
3062 | } |
---|
3063 | if (upperChange[i]) { |
---|
3064 | upperList[nUpperChange++]=i; |
---|
3065 | } |
---|
3066 | } |
---|
3067 | lowerList[-2]=nLowerChange; |
---|
3068 | upperList[-2]=nUpperChange; |
---|
3069 | for (int i=numberColumns_;i<numberTotal;i++) { |
---|
3070 | if (lowerChange[i]) { |
---|
3071 | lowerList[nLowerChange++]=i; |
---|
3072 | } |
---|
3073 | if (upperChange[i]) { |
---|
3074 | upperList[nUpperChange++]=i; |
---|
3075 | } |
---|
3076 | } |
---|
3077 | lowerList[-1]=nLowerChange; |
---|
3078 | upperList[-1]=nUpperChange; |
---|
3079 | memcpy(lowerCopy,columnLower_,numberColumns_*sizeof(double)); |
---|
3080 | memcpy(upperCopy,columnUpper_,numberColumns_*sizeof(double)); |
---|
3081 | memcpy(lowerCopy+numberColumns_, |
---|
3082 | rowLower_,numberRows_*sizeof(double)); |
---|
3083 | memcpy(upperCopy+numberColumns_, |
---|
3084 | rowUpper_,numberRows_*sizeof(double)); |
---|
3085 | { |
---|
3086 | // extra for unscaled |
---|
3087 | double * unscaledCopy; |
---|
3088 | unscaledCopy = lowerCopy + numberTotal; |
---|
3089 | memcpy(unscaledCopy,columnLower_,numberColumns_*sizeof(double)); |
---|
3090 | memcpy(unscaledCopy+numberColumns_, |
---|
3091 | rowLower_,numberRows_*sizeof(double)); |
---|
3092 | unscaledCopy = upperCopy + numberTotal; |
---|
3093 | memcpy(unscaledCopy,columnUpper_,numberColumns_*sizeof(double)); |
---|
3094 | memcpy(unscaledCopy+numberColumns_, |
---|
3095 | rowUpper_,numberRows_*sizeof(double)); |
---|
3096 | } |
---|
3097 | int returnCode=0; |
---|
3098 | paramData.startingTheta=startingTheta; |
---|
3099 | paramData.endingTheta=endingTheta; |
---|
3100 | paramData.maxTheta=endingTheta; |
---|
3101 | computeRhsEtc(paramData); |
---|
3102 | bool swapped=false; |
---|
3103 | // Dantzig |
---|
3104 | #define ALL_DANTZIG |
---|
3105 | #ifdef ALL_DANTZIG |
---|
3106 | ClpDualRowPivot * savePivot = dualRowPivot_; |
---|
3107 | dualRowPivot_ = new ClpDualRowDantzig(); |
---|
3108 | dualRowPivot_->setModel(this); |
---|
3109 | #else |
---|
3110 | ClpDualRowPivot * savePivot = NULL; |
---|
3111 | #endif |
---|
3112 | if (!returnCode) { |
---|
3113 | assert (objective_->type()==1); |
---|
3114 | objective_->setType(2); // in case matrix empty |
---|
3115 | returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0); |
---|
3116 | objective_->setType(1); |
---|
3117 | if (!returnCode) { |
---|
3118 | double saveDualBound=dualBound_; |
---|
3119 | dualBound_=CoinMax(dualBound_,1.0e15); |
---|
3120 | swapped=true; |
---|
3121 | double * temp; |
---|
3122 | memcpy(saveLower,lower_,numberTotal*sizeof(double)); |
---|
3123 | temp=saveLower; |
---|
3124 | saveLower=lower_; |
---|
3125 | lower_=temp; |
---|
3126 | //columnLowerWork_ = lower_; |
---|
3127 | //rowLowerWork_ = lower_ + numberColumns_; |
---|
3128 | memcpy(saveUpper,upper_,numberTotal*sizeof(double)); |
---|
3129 | temp=saveUpper; |
---|
3130 | saveUpper=upper_; |
---|
3131 | upper_=temp; |
---|
3132 | //columnUpperWork_ = upper_; |
---|
3133 | //rowUpperWork_ = upper_ + numberColumns_; |
---|
3134 | if (rowScale_) { |
---|
3135 | // scale saved and change arrays |
---|
3136 | double * lowerChange = lower_+numberTotal; |
---|
3137 | double * upperChange = upper_+numberTotal; |
---|
3138 | double * lowerSave = lowerChange+numberTotal; |
---|
3139 | double * upperSave = upperChange+numberTotal; |
---|
3140 | for (int i=0;i<numberColumns_;i++) { |
---|
3141 | double multiplier = inverseColumnScale_[i]; |
---|
3142 | if (lowerSave[i]>-1.0e20) |
---|
3143 | lowerSave[i] *= multiplier; |
---|
3144 | if (upperSave[i]<1.0e20) |
---|
3145 | upperSave[i] *= multiplier; |
---|
3146 | lowerChange[i] *= multiplier; |
---|
3147 | upperChange[i] *= multiplier; |
---|
3148 | } |
---|
3149 | lowerChange += numberColumns_; |
---|
3150 | upperChange += numberColumns_; |
---|
3151 | lowerSave += numberColumns_; |
---|
3152 | upperSave += numberColumns_; |
---|
3153 | for (int i=0;i<numberRows_;i++) { |
---|
3154 | double multiplier = rowScale_[i]; |
---|
3155 | if (lowerSave[i]>-1.0e20) |
---|
3156 | lowerSave[i] *= multiplier; |
---|
3157 | if (upperSave[i]<1.0e20) |
---|
3158 | upperSave[i] *= multiplier; |
---|
3159 | lowerChange[i] *= multiplier; |
---|
3160 | upperChange[i] *= multiplier; |
---|
3161 | } |
---|
3162 | } |
---|
3163 | //double saveEndingTheta = endingTheta; |
---|
3164 | double * saveDuals = NULL; |
---|
3165 | reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data); |
---|
3166 | if (numberPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-4) { |
---|
3167 | // probably can get rid of this if we adjust every change in theta |
---|
3168 | //printf("INFEAS_A %d %g\n",numberPrimalInfeasibilities_, |
---|
3169 | // sumPrimalInfeasibilities_); |
---|
3170 | int pass=100; |
---|
3171 | while(sumPrimalInfeasibilities_) { |
---|
3172 | pass--; |
---|
3173 | if (!pass) |
---|
3174 | break; |
---|
3175 | problemStatus_=-1; |
---|
3176 | for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) { |
---|
3177 | double value=solution_[iSequence]; |
---|
3178 | // remember scaling |
---|
3179 | if (value<lower_[iSequence]-1.0e-9) { |
---|
3180 | lowerCopy[iSequence]+=value-lower_[iSequence]; |
---|
3181 | lower_[iSequence]=value; |
---|
3182 | } else if (value>upper_[iSequence]+1.0e-9) { |
---|
3183 | upperCopy[iSequence]+=value-upper_[iSequence]; |
---|
3184 | upper_[iSequence]=value; |
---|
3185 | } |
---|
3186 | } |
---|
3187 | reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data); |
---|
3188 | } |
---|
3189 | } |
---|
3190 | if (!problemStatus_) { |
---|
3191 | if (nLowerChange||nUpperChange) { |
---|
3192 | #ifndef ALL_DANTZIG |
---|
3193 | // Dantzig |
---|
3194 | savePivot = dualRowPivot_; |
---|
3195 | dualRowPivot_ = new ClpDualRowDantzig(); |
---|
3196 | dualRowPivot_->setModel(this); |
---|
3197 | #endif |
---|
3198 | //for (int i=0;i<numberRows_+numberColumns_;i++) |
---|
3199 | //setFakeBound(i, noFake); |
---|
3200 | // Now do parametrics |
---|
3201 | handler_->message(CLP_PARAMETRICS_STATS, messages_) |
---|
3202 | << startingTheta << objectiveValue() << CoinMessageEol; |
---|
3203 | bool canSkipFactorization=true; |
---|
3204 | while (!returnCode) { |
---|
3205 | paramData.startingTheta=startingTheta; |
---|
3206 | paramData.endingTheta=endingTheta; |
---|
3207 | returnCode = parametricsLoop(paramData, |
---|
3208 | data,canSkipFactorization); |
---|
3209 | startingTheta=paramData.startingTheta; |
---|
3210 | endingTheta=paramData.endingTheta; |
---|
3211 | canSkipFactorization=false; |
---|
3212 | if (!returnCode) { |
---|
3213 | //startingTheta = endingTheta; |
---|
3214 | //endingTheta = saveEndingTheta; |
---|
3215 | handler_->message(CLP_PARAMETRICS_STATS, messages_) |
---|
3216 | << startingTheta << objectiveValue() << CoinMessageEol; |
---|
3217 | if (startingTheta >= endingTheta-primalTolerance_ |
---|
3218 | ||problemStatus_==2) |
---|
3219 | break; |
---|
3220 | } else if (returnCode == -1) { |
---|
3221 | // trouble - do external solve |
---|
3222 | abort(); //needToDoSomething = true; |
---|
3223 | } else if (problemStatus_==1) { |
---|
3224 | // can't move any further |
---|
3225 | handler_->message(CLP_PARAMETRICS_STATS, messages_) |
---|
3226 | << endingTheta << objectiveValue() << CoinMessageEol; |
---|
3227 | problemStatus_=0; |
---|
3228 | } |
---|
3229 | } |
---|
3230 | } |
---|
3231 | dualBound_ = saveDualBound; |
---|
3232 | //reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data); |
---|
3233 | } else { |
---|
3234 | // check if empty |
---|
3235 | //if (!numberRows_||!matrix_->getNumElements()) { |
---|
3236 | // success |
---|
3237 | #ifdef CLP_USER_DRIVEN |
---|
3238 | //theta_ = endingTheta; |
---|
3239 | //eventHandler_->event(ClpEventHandler::theta); |
---|
3240 | #endif |
---|
3241 | //} |
---|
3242 | } |
---|
3243 | } |
---|
3244 | if (problemStatus_==2) { |
---|
3245 | delete [] ray_; |
---|
3246 | ray_ = new double [numberColumns_]; |
---|
3247 | } |
---|
3248 | if (swapped&&lower_) { |
---|
3249 | double * temp=saveLower; |
---|
3250 | saveLower=lower_; |
---|
3251 | lower_=temp; |
---|
3252 | temp=saveUpper; |
---|
3253 | saveUpper=upper_; |
---|
3254 | upper_=temp; |
---|
3255 | } |
---|
3256 | reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0); |
---|
3257 | } |
---|
3258 | if (!scalingFlag_) { |
---|
3259 | memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double)); |
---|
3260 | memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double)); |
---|
3261 | memcpy(rowLower_,lowerCopy+numberColumns_, |
---|
3262 | numberRows_*sizeof(double)); |
---|
3263 | memcpy(rowUpper_,upperCopy+numberColumns_, |
---|
3264 | numberRows_*sizeof(double)); |
---|
3265 | } else { |
---|
3266 | // extra for unscaled |
---|
3267 | double * unscaledCopy; |
---|
3268 | unscaledCopy = lowerCopy + numberTotal; |
---|
3269 | memcpy(columnLower_,unscaledCopy,numberColumns_*sizeof(double)); |
---|
3270 | memcpy(rowLower_,unscaledCopy+numberColumns_, |
---|
3271 | numberRows_*sizeof(double)); |
---|
3272 | unscaledCopy = upperCopy + numberTotal; |
---|
3273 | memcpy(columnUpper_,unscaledCopy,numberColumns_*sizeof(double)); |
---|
3274 | memcpy(rowUpper_,unscaledCopy+numberColumns_, |
---|
3275 | numberRows_*sizeof(double)); |
---|
3276 | } |
---|
3277 | delete [] saveLower; |
---|
3278 | delete [] saveUpper; |
---|
3279 | #ifdef ALL_DANTZIG |
---|
3280 | if (savePivot) { |
---|
3281 | #endif |
---|
3282 | delete dualRowPivot_; |
---|
3283 | dualRowPivot_ = savePivot; |
---|
3284 | #ifdef ALL_DANTZIG |
---|
3285 | } |
---|
3286 | #endif |
---|
3287 | // Restore any saved stuff |
---|
3288 | restoreData(data); |
---|
3289 | perturbation_ = savePerturbation; |
---|
3290 | delete rowArray_[4]; |
---|
3291 | rowArray_[4]=NULL; |
---|
3292 | delete rowArray_[5]; |
---|
3293 | rowArray_[5]=NULL; |
---|
3294 | char line[100]; |
---|
3295 | sprintf(line,"Ending theta %g\n", endingTheta); |
---|
3296 | handler_->message(CLP_GENERAL,messages_) |
---|
3297 | << line << CoinMessageEol; |
---|
3298 | return problemStatus_; |
---|
3299 | } |
---|
3300 | int |
---|
3301 | ClpSimplexOther::parametricsLoop(parametricsData & paramData, |
---|
3302 | ClpDataSave & data,bool canSkipFactorization) |
---|
3303 | { |
---|
3304 | double & startingTheta = paramData.startingTheta; |
---|
3305 | double & endingTheta = paramData.endingTheta; |
---|
3306 | int numberTotal = numberRows_+numberColumns_; |
---|
3307 | // stuff is already at starting |
---|
3308 | const int * lowerList = paramData.lowerList; |
---|
3309 | const int * upperList = paramData.upperList; |
---|
3310 | problemStatus_ = -1; |
---|
3311 | //double saveEndingTheta=endingTheta; |
---|
3312 | |
---|
3313 | // This says whether to restore things etc |
---|
3314 | // startup will have factorized so can skip |
---|
3315 | int factorType = 0; |
---|
3316 | // Start check for cycles |
---|
3317 | progress_.startCheck(); |
---|
3318 | // Say change made on first iteration |
---|
3319 | changeMade_ = 1; |
---|
3320 | /* |
---|
3321 | Status of problem: |
---|
3322 | 0 - optimal |
---|
3323 | 1 - infeasible |
---|
3324 | 2 - unbounded |
---|
3325 | -1 - iterating |
---|
3326 | -2 - factorization wanted |
---|
3327 | -3 - redo checking without factorization |
---|
3328 | -4 - looks infeasible |
---|
3329 | */ |
---|
3330 | while (problemStatus_ < 0) { |
---|
3331 | int iRow, iColumn; |
---|
3332 | // clear |
---|
3333 | for (iRow = 0; iRow < 6; iRow++) { |
---|
3334 | rowArray_[iRow]->clear(); |
---|
3335 | } |
---|
3336 | |
---|
3337 | for (iColumn = 0; iColumn < 2; iColumn++) { |
---|
3338 | columnArray_[iColumn]->clear(); |
---|
3339 | } |
---|
3340 | |
---|
3341 | // give matrix (and model costs and bounds a chance to be |
---|
3342 | // refreshed (normally null) |
---|
3343 | matrix_->refresh(this); |
---|
3344 | // may factorize, checks if problem finished |
---|
3345 | if (!canSkipFactorization) |
---|
3346 | statusOfProblemInParametrics(factorType, data); |
---|
3347 | canSkipFactorization=false; |
---|
3348 | if (numberPrimalInfeasibilities_) { |
---|
3349 | if (largestPrimalError_>1.0e3&&startingTheta>1.0e10) { |
---|
3350 | // treat as success |
---|
3351 | problemStatus_=0; |
---|
3352 | endingTheta=startingTheta; |
---|
3353 | break; |
---|
3354 | } |
---|
3355 | // probably can get rid of this if we adjust every change in theta |
---|
3356 | //printf("INFEAS %d %g\n",numberPrimalInfeasibilities_, |
---|
3357 | // sumPrimalInfeasibilities_); |
---|
3358 | const double * lowerChange = lower_+numberTotal; |
---|
3359 | const double * upperChange = upper_+numberTotal; |
---|
3360 | const double * startLower = lowerChange+numberTotal; |
---|
3361 | const double * startUpper = upperChange+numberTotal; |
---|
3362 | //startingTheta -= 1.0e-7; |
---|
3363 | int nLowerChange = lowerList[-1]; |
---|
3364 | for (int i = 0; i < nLowerChange; i++) { |
---|
3365 | int iSequence = lowerList[i]; |
---|
3366 | lower_[iSequence] = startLower[iSequence] + startingTheta * lowerChange[iSequence]; |
---|
3367 | } |
---|
3368 | int nUpperChange = upperList[-1]; |
---|
3369 | for (int i = 0; i < nUpperChange; i++) { |
---|
3370 | int iSequence = upperList[i]; |
---|
3371 | upper_[iSequence] = startUpper[iSequence] + startingTheta * upperChange[iSequence]; |
---|
3372 | } |
---|
3373 | // adjust rhs in case dual uses |
---|
3374 | memcpy(columnLower_,lower_,numberColumns_*sizeof(double)); |
---|
3375 | memcpy(rowLower_,lower_+numberColumns_,numberRows_*sizeof(double)); |
---|
3376 | memcpy(columnUpper_,upper_,numberColumns_*sizeof(double)); |
---|
3377 | memcpy(rowUpper_,upper_+numberColumns_,numberRows_*sizeof(double)); |
---|
3378 | if (rowScale_) { |
---|
3379 | for (int i=0;i<numberColumns_;i++) { |
---|
3380 | double multiplier = columnScale_[i]; |
---|
3381 | if (columnLower_[i]>-1.0e20) |
---|
3382 | columnLower_[i] *= multiplier; |
---|
3383 | if (columnUpper_[i]<1.0e20) |
---|
3384 | columnUpper_[i] *= multiplier; |
---|
3385 | } |
---|
3386 | for (int i=0;i<numberRows_;i++) { |
---|
3387 | double multiplier = inverseRowScale_[i]; |
---|
3388 | if (rowLower_[i]>-1.0e20) |
---|
3389 | rowLower_[i] *= multiplier; |
---|
3390 | if (rowUpper_[i]<1.0e20) |
---|
3391 | rowUpper_[i] *= multiplier; |
---|
3392 | } |
---|
3393 | } |
---|
3394 | double * saveDuals = NULL; |
---|
3395 | problemStatus_=-1; |
---|
3396 | ClpObjective * saveObjective = objective_; |
---|
3397 | reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data); |
---|
3398 | if (saveObjective!=objective_) { |
---|
3399 | delete objective_; |
---|
3400 | objective_=saveObjective; |
---|
3401 | } |
---|
3402 | int pass=100; |
---|
3403 | double moved=0.0; |
---|
3404 | while(sumPrimalInfeasibilities_) { |
---|
3405 | //printf("INFEAS pass %d %d %g\n",100-pass,numberPrimalInfeasibilities_, |
---|
3406 | // sumPrimalInfeasibilities_); |
---|
3407 | pass--; |
---|
3408 | if (!pass) |
---|
3409 | break; |
---|
3410 | problemStatus_=-1; |
---|
3411 | for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) { |
---|
3412 | double value=solution_[iSequence]; |
---|
3413 | if (value<lower_[iSequence]-1.0e-9) { |
---|
3414 | moved += lower_[iSequence]-value; |
---|
3415 | lower_[iSequence]=value; |
---|
3416 | } else if (value>upper_[iSequence]+1.0e-9) { |
---|
3417 | moved += upper_[iSequence]-value; |
---|
3418 | upper_[iSequence]=value; |
---|
3419 | } |
---|
3420 | } |
---|
3421 | if (!moved) { |
---|
3422 | for (int iSequence=0;iSequence<numberColumns_;iSequence++) { |
---|
3423 | double value=solution_[iSequence]; |
---|
3424 | if (value<lower_[iSequence]-1.0e-9) { |
---|
3425 | moved += lower_[iSequence]-value; |
---|
3426 | lower_[iSequence]=value; |
---|
3427 | } else if (value>upper_[iSequence]+1.0e-9) { |
---|
3428 | moved += upper_[iSequence]-value; |
---|
3429 | upper_[iSequence]=value; |
---|
3430 | } |
---|
3431 | } |
---|
3432 | } |
---|
3433 | assert (moved); |
---|
3434 | reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data); |
---|
3435 | } |
---|
3436 | // adjust |
---|
3437 | //printf("Should adjust - moved %g\n",moved); |
---|
3438 | } |
---|
3439 | // Say good factorization |
---|
3440 | factorType = 1; |
---|
3441 | if (data.sparseThreshold_) { |
---|
3442 | // use default at present |
---|
3443 | factorization_->sparseThreshold(0); |
---|
3444 | factorization_->goSparse(); |
---|
3445 | } |
---|
3446 | |
---|
3447 | // exit if victory declared |
---|
3448 | if (problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7 ) |
---|
3449 | break; |
---|
3450 | |
---|
3451 | // test for maximum iterations |
---|
3452 | if (hitMaximumIterations()) { |
---|
3453 | problemStatus_ = 3; |
---|
3454 | break; |
---|
3455 | } |
---|
3456 | #ifdef CLP_USER_DRIVEN |
---|
3457 | // Check event |
---|
3458 | { |
---|
3459 | int status = eventHandler_->event(ClpEventHandler::endOfFactorization); |
---|
3460 | if (status >= 0) { |
---|
3461 | problemStatus_ = 5; |
---|
3462 | secondaryStatus_ = ClpEventHandler::endOfFactorization; |
---|
3463 | break; |
---|
3464 | } |
---|
3465 | } |
---|
3466 | #endif |
---|
3467 | // Do iterations |
---|
3468 | problemStatus_=-1; |
---|
3469 | whileIterating(paramData, 0.0, |
---|
3470 | NULL); |
---|
3471 | //startingTheta = endingTheta; |
---|
3472 | //endingTheta = saveEndingTheta; |
---|
3473 | } |
---|
3474 | if (!problemStatus_/*||problemStatus_==2*/) { |
---|
3475 | theta_ = endingTheta; |
---|
3476 | #ifdef CLP_USER_DRIVEN |
---|
3477 | { |
---|
3478 | double saveTheta=theta_; |
---|
3479 | theta_ = endingTheta; |
---|
3480 | int status=eventHandler_->event(ClpEventHandler::theta); |
---|
3481 | if (status>=0&&status<10) { |
---|
3482 | endingTheta=theta_; |
---|
3483 | theta_=saveTheta; |
---|
3484 | problemStatus_=-1; |
---|
3485 | } else { |
---|
3486 | if (status>=10) { |
---|
3487 | problemStatus_=status-10; |
---|
3488 | startingTheta=endingTheta; |
---|
3489 | } |
---|
3490 | theta_=saveTheta; |
---|
3491 | } |
---|
3492 | } |
---|
3493 | #endif |
---|
3494 | return 0; |
---|
3495 | } else if (problemStatus_ == 10) { |
---|
3496 | return -1; |
---|
3497 | } else { |
---|
3498 | return problemStatus_; |
---|
3499 | } |
---|
3500 | } |
---|
3501 | /* Checks if finished. Updates status */ |
---|
3502 | void |
---|
3503 | ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData) |
---|
3504 | { |
---|
3505 | if (type == 2) { |
---|
3506 | // trouble - go to recovery |
---|
3507 | problemStatus_ = 10; |
---|
3508 | return; |
---|
3509 | } |
---|
3510 | if (problemStatus_ > -3 || factorization_->pivots()) { |
---|
3511 | // factorize |
---|
3512 | // later on we will need to recover from singularities |
---|
3513 | // also we could skip if first time |
---|
3514 | if (type) { |
---|
3515 | // is factorization okay? |
---|
3516 | if (internalFactorize(1)) { |
---|
3517 | // trouble - go to recovery |
---|
3518 | problemStatus_ = 10; |
---|
3519 | return; |
---|
3520 | } |
---|
3521 | } |
---|
3522 | if (problemStatus_ != -4 || factorization_->pivots() > 10) |
---|
3523 | problemStatus_ = -3; |
---|
3524 | } |
---|
3525 | // at this stage status is -3 or -4 if looks infeasible |
---|
3526 | // get primal and dual solutions |
---|
3527 | gutsOfSolution(NULL, NULL); |
---|
3528 | double realDualInfeasibilities = sumDualInfeasibilities_; |
---|
3529 | // If bad accuracy treat as singular |
---|
3530 | if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) { |
---|
3531 | // trouble - go to recovery |
---|
3532 | problemStatus_ = 10; |
---|
3533 | return; |
---|
3534 | } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) { |
---|
3535 | // Can reduce tolerance |
---|
3536 | double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_); |
---|
3537 | factorization_->pivotTolerance(newTolerance); |
---|
3538 | } |
---|
3539 | // Check if looping |
---|
3540 | int loop; |
---|
3541 | if (type != 2) |
---|
3542 | loop = progress_.looping(); |
---|
3543 | else |
---|
3544 | loop = -1; |
---|
3545 | if (loop >= 0) { |
---|
3546 | problemStatus_ = loop; //exit if in loop |
---|
3547 | if (!problemStatus_) { |
---|
3548 | // declaring victory |
---|
3549 | numberPrimalInfeasibilities_ = 0; |
---|
3550 | sumPrimalInfeasibilities_ = 0.0; |
---|
3551 | } else { |
---|
3552 | problemStatus_ = 10; // instead - try other algorithm |
---|
3553 | } |
---|
3554 | return; |
---|
3555 | } else if (loop < -1) { |
---|
3556 | // something may have changed |
---|
3557 | gutsOfSolution(NULL, NULL); |
---|
3558 | } |
---|
3559 | progressFlag_ = 0; //reset progress flag |
---|
3560 | if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) { |
---|
3561 | handler_->message(CLP_SIMPLEX_STATUS, messages_) |
---|
3562 | << numberIterations_ << objectiveValue(); |
---|
3563 | handler_->printing(sumPrimalInfeasibilities_ > 0.0) |
---|
3564 | << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_; |
---|
3565 | handler_->printing(sumDualInfeasibilities_ > 0.0) |
---|
3566 | << sumDualInfeasibilities_ << numberDualInfeasibilities_; |
---|
3567 | handler_->printing(numberDualInfeasibilitiesWithoutFree_ |
---|
3568 | < numberDualInfeasibilities_) |
---|
3569 | << numberDualInfeasibilitiesWithoutFree_; |
---|
3570 | handler_->message() << CoinMessageEol; |
---|
3571 | } |
---|
3572 | #ifdef CLP_USER_DRIVEN |
---|
3573 | if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-7) { |
---|
3574 | int status=eventHandler_->event(ClpEventHandler::slightlyInfeasible); |
---|
3575 | if (status>=0) { |
---|
3576 | // fix up |
---|
3577 | for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) { |
---|
3578 | double value=solution_[iSequence]; |
---|
3579 | if (value<=lower_[iSequence]-primalTolerance_) { |
---|
3580 | lower_[iSequence]=value; |
---|
3581 | } else if (value>=upper_[iSequence]+primalTolerance_) { |
---|
3582 | upper_[iSequence]=value; |
---|
3583 | } |
---|
3584 | } |
---|
3585 | numberPrimalInfeasibilities_ = 0; |
---|
3586 | sumPrimalInfeasibilities_ = 0.0; |
---|
3587 | } |
---|
3588 | } |
---|
3589 | #endif |
---|
3590 | /* If we are primal feasible and any dual infeasibilities are on |
---|
3591 | free variables then it is better to go to primal */ |
---|
3592 | if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ && |
---|
3593 | numberDualInfeasibilities_) { |
---|
3594 | problemStatus_ = 10; |
---|
3595 | return; |
---|
3596 | } |
---|
3597 | |
---|
3598 | // check optimal |
---|
3599 | // give code benefit of doubt |
---|
3600 | if (sumOfRelaxedDualInfeasibilities_ == 0.0 && |
---|
3601 | sumOfRelaxedPrimalInfeasibilities_ == 0.0) { |
---|
3602 | // say optimal (with these bounds etc) |
---|
3603 | numberDualInfeasibilities_ = 0; |
---|
3604 | sumDualInfeasibilities_ = 0.0; |
---|
3605 | numberPrimalInfeasibilities_ = 0; |
---|
3606 | sumPrimalInfeasibilities_ = 0.0; |
---|
3607 | } |
---|
3608 | if (dualFeasible() || problemStatus_ == -4) { |
---|
3609 | progress_.modifyObjective(objectiveValue_ |
---|
3610 | - sumDualInfeasibilities_ * dualBound_); |
---|
3611 | } |
---|
3612 | if (numberPrimalInfeasibilities_) { |
---|
3613 | if (problemStatus_ == -4 || problemStatus_ == -5) { |
---|
3614 | problemStatus_ = 1; // infeasible |
---|
3615 | } |
---|
3616 | } else if (numberDualInfeasibilities_) { |
---|
3617 | // clean up |
---|
3618 | problemStatus_ = 10; |
---|
3619 | } else { |
---|
3620 | problemStatus_ = 0; |
---|
3621 | } |
---|
3622 | lastGoodIteration_ = numberIterations_; |
---|
3623 | if (problemStatus_ < 0) { |
---|
3624 | sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful |
---|
3625 | if (sumDualInfeasibilities_) |
---|
3626 | numberDualInfeasibilities_ = 1; |
---|
3627 | } |
---|
3628 | // Allow matrices to be sorted etc |
---|
3629 | int fake = -999; // signal sort |
---|
3630 | matrix_->correctSequence(this, fake, fake); |
---|
3631 | } |
---|
3632 | //static double lastThetaX=0.0; |
---|
3633 | /* This has the flow between re-factorizations |
---|
3634 | Reasons to come out: |
---|
3635 | -1 iterations etc |
---|
3636 | -2 inaccuracy |
---|
3637 | -3 slight inaccuracy (and done iterations) |
---|
3638 | +0 looks optimal (might be unbounded - but we will investigate) |
---|
3639 | +1 looks infeasible |
---|
3640 | +3 max iterations |
---|
3641 | +4 accuracy problems |
---|
3642 | */ |
---|
3643 | int |
---|
3644 | ClpSimplexOther::whileIterating(parametricsData & paramData, double /*reportIncrement*/, |
---|
3645 | const double * /*changeObjective*/) |
---|
3646 | { |
---|
3647 | double & startingTheta = paramData.startingTheta; |
---|
3648 | double & endingTheta = paramData.endingTheta; |
---|
3649 | const double * lowerChange = paramData.lowerChange; |
---|
3650 | const double * upperChange = paramData.upperChange; |
---|
3651 | int numberTotal = numberColumns_ + numberRows_; |
---|
3652 | const int * lowerList = paramData.lowerList; |
---|
3653 | const int * upperList = paramData.upperList; |
---|
3654 | //#define CLP_PARAMETRIC_DENSE_ARRAYS 2 |
---|
3655 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
3656 | double * lowerGap = paramData.lowerGap; |
---|
3657 | double * upperGap = paramData.upperGap; |
---|
3658 | double * lowerCoefficient = paramData.lowerCoefficient; |
---|
3659 | double * upperCoefficient = paramData.upperCoefficient; |
---|
3660 | #endif |
---|
3661 | // do basic pointers |
---|
3662 | int * backwardBasic = paramData.backwardBasic; |
---|
3663 | for (int i=0;i<numberTotal;i++) |
---|
3664 | backwardBasic[i]=-1; |
---|
3665 | for (int i=0;i<numberRows_;i++) { |
---|
3666 | int iPivot=pivotVariable_[i]; |
---|
3667 | backwardBasic[iPivot]=i; |
---|
3668 | } |
---|
3669 | { |
---|
3670 | int i; |
---|
3671 | for (i = 0; i < 4; i++) { |
---|
3672 | rowArray_[i]->clear(); |
---|
3673 | } |
---|
3674 | for (i = 0; i < 2; i++) { |
---|
3675 | columnArray_[i]->clear(); |
---|
3676 | } |
---|
3677 | } |
---|
3678 | // if can't trust much and long way from optimal then relax |
---|
3679 | if (largestPrimalError_ > 10.0) |
---|
3680 | factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0)); |
---|
3681 | else |
---|
3682 | factorization_->relaxAccuracyCheck(1.0); |
---|
3683 | // status stays at -1 while iterating, >=0 finished, -2 to invert |
---|
3684 | // status -3 to go to top without an invert |
---|
3685 | int returnCode = -1; |
---|
3686 | double lastTheta = startingTheta; |
---|
3687 | double useTheta = startingTheta; |
---|
3688 | while (problemStatus_ == -1) { |
---|
3689 | double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50); |
---|
3690 | // Get theta for bounds - we know can't crossover |
---|
3691 | int pivotType = nextTheta(1, increaseTheta, paramData, |
---|
3692 | NULL); |
---|
3693 | useTheta += theta_; |
---|
3694 | double change = useTheta - lastTheta; |
---|
3695 | if (paramData.firstIteration) { |
---|
3696 | // redo rhs etc to make as accurate as possible |
---|
3697 | paramData.firstIteration=false; |
---|
3698 | if (change>1.0e-14) { |
---|
3699 | startingTheta=useTheta; |
---|
3700 | lastTheta=startingTheta; |
---|
3701 | change=0.0; |
---|
3702 | // restore rhs |
---|
3703 | const double * saveLower = paramData.lowerChange+2*numberTotal; |
---|
3704 | memcpy(columnLower_,saveLower,numberColumns_*sizeof(double)); |
---|
3705 | memcpy(rowLower_,saveLower+numberColumns_,numberRows_*sizeof(double)); |
---|
3706 | const double * saveUpper = paramData.upperChange+2*numberTotal; |
---|
3707 | memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double)); |
---|
3708 | memcpy(rowUpper_,saveUpper+numberColumns_,numberRows_*sizeof(double)); |
---|
3709 | paramData.startingTheta=startingTheta; |
---|
3710 | computeRhsEtc(paramData); |
---|
3711 | redoInternalArrays(); |
---|
3712 | // Update solution |
---|
3713 | rowArray_[4]->clear(); |
---|
3714 | for (int i=0;i<numberTotal;i++) { |
---|
3715 | if (getStatus(i)==atLowerBound||getStatus(i)==isFixed) |
---|
3716 | solution_[i]=lower_[i]; |
---|
3717 | else if (getStatus(i)==atUpperBound) |
---|
3718 | solution_[i]=upper_[i]; |
---|
3719 | } |
---|
3720 | gutsOfSolution(NULL,NULL); |
---|
3721 | } |
---|
3722 | } |
---|
3723 | if (change>1.0e-14) { |
---|
3724 | int n; |
---|
3725 | n=lowerList[-1]; |
---|
3726 | for (int i=0;i<n;i++) { |
---|
3727 | int iSequence = lowerList[i]; |
---|
3728 | double thisChange = change * lowerChange[iSequence]; |
---|
3729 | double newValue = lower_[iSequence] + thisChange; |
---|
3730 | lower_[iSequence] = newValue; |
---|
3731 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
3732 | if (getStatus(iSequence)==basic) { |
---|
3733 | int iRow=backwardBasic[iSequence]; |
---|
3734 | lowerGap[iRow] -= thisChange; |
---|
3735 | } else if(getStatus(iSequence)==atLowerBound) { |
---|
3736 | solution_[iSequence] = newValue; |
---|
3737 | } |
---|
3738 | #else |
---|
3739 | if(getStatus(iSequence)==atLowerBound) { |
---|
3740 | solution_[iSequence] = newValue; |
---|
3741 | } |
---|
3742 | #endif |
---|
3743 | #if 0 |
---|
3744 | // may have to adjust other bound |
---|
3745 | double otherValue = upper_[iSequence]; |
---|
3746 | if (otherValue-newValue<dualBound_) { |
---|
3747 | //originalBound(iSequence,useTheta,lowerChange,upperChange); |
---|
3748 | //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence); |
---|
3749 | //ClpTraceDebug (fabs(lower_[iSequence]-newValue)<1.0e-5); |
---|
3750 | } |
---|
3751 | #endif |
---|
3752 | } |
---|
3753 | n=upperList[-1]; |
---|
3754 | for (int i=0;i<n;i++) { |
---|
3755 | int iSequence = upperList[i]; |
---|
3756 | double thisChange = change * upperChange[iSequence]; |
---|
3757 | double newValue = upper_[iSequence] + thisChange; |
---|
3758 | upper_[iSequence] = newValue; |
---|
3759 | if(getStatus(iSequence)==atUpperBound|| |
---|
3760 | getStatus(iSequence)==isFixed) { |
---|
3761 | solution_[iSequence] = newValue; |
---|
3762 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
3763 | } else if (getStatus(iSequence)==basic) { |
---|
3764 | int iRow=backwardBasic[iSequence]; |
---|
3765 | upperGap[iRow] += thisChange; |
---|
3766 | #endif |
---|
3767 | } |
---|
3768 | // may have to adjust other bound |
---|
3769 | double otherValue = lower_[iSequence]; |
---|
3770 | if (newValue-otherValue<dualBound_) { |
---|
3771 | //originalBound(iSequence,useTheta,lowerChange,upperChange); |
---|
3772 | //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence); |
---|
3773 | //ClpTraceDebug (fabs(upper_[iSequence]-newValue)<1.0e-5); |
---|
3774 | } |
---|
3775 | } |
---|
3776 | } |
---|
3777 | sequenceIn_=-1; |
---|
3778 | if (pivotType) { |
---|
3779 | if (useTheta>lastTheta+1.0e-9) { |
---|
3780 | handler_->message(CLP_PARAMETRICS_STATS, messages_) |
---|
3781 | << useTheta << objectiveValue() << CoinMessageEol; |
---|
3782 | lastTheta = useTheta; |
---|
3783 | } |
---|
3784 | problemStatus_ = -2; |
---|
3785 | if (!factorization_->pivots()&&pivotRow_<0) |
---|
3786 | problemStatus_=2; |
---|
3787 | #ifdef CLP_USER_DRIVEN |
---|
3788 | { |
---|
3789 | double saveTheta=theta_; |
---|
3790 | theta_ = endingTheta; |
---|
3791 | if (problemStatus_==2&&theta_>paramData.acceptableMaxTheta) |
---|
3792 | theta_=COIN_DBL_MAX; // we have finished |
---|
3793 | int status=eventHandler_->event(ClpEventHandler::theta); |
---|
3794 | if (status>=0&&status<10) { |
---|
3795 | endingTheta=theta_; |
---|
3796 | problemStatus_=-1; |
---|
3797 | continue; |
---|
3798 | } else { |
---|
3799 | if (status>=10) |
---|
3800 | problemStatus_=status-10; |
---|
3801 | if (status<0) |
---|
3802 | startingTheta = useTheta; |
---|
3803 | } |
---|
3804 | theta_=saveTheta; |
---|
3805 | } |
---|
3806 | #else |
---|
3807 | startingTheta = useTheta; |
---|
3808 | #endif |
---|
3809 | return 4; |
---|
3810 | } |
---|
3811 | // choose row to go out |
---|
3812 | //reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1); |
---|
3813 | if (pivotRow_ >= 0) { |
---|
3814 | // we found a pivot row |
---|
3815 | if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) { |
---|
3816 | handler_->message(CLP_SIMPLEX_PIVOTROW, messages_) |
---|
3817 | << pivotRow_ |
---|
3818 | << CoinMessageEol; |
---|
3819 | } |
---|
3820 | // check accuracy of weights |
---|
3821 | dualRowPivot_->checkAccuracy(); |
---|
3822 | // do ratio test for normal iteration |
---|
3823 | double bestPossiblePivot = bestPivot(); |
---|
3824 | if (sequenceIn_ >= 0) { |
---|
3825 | // normal iteration |
---|
3826 | // update the incoming column |
---|
3827 | double btranAlpha = -alpha_ * directionOut_; // for check |
---|
3828 | #ifndef COIN_FAC_NEW |
---|
3829 | unpackPacked(rowArray_[1]); |
---|
3830 | #else |
---|
3831 | unpack(rowArray_[1]); |
---|
3832 | #endif |
---|
3833 | // and update dual weights (can do in parallel - with extra array) |
---|
3834 | rowArray_[2]->clear(); |
---|
3835 | alpha_ = dualRowPivot_->updateWeights(rowArray_[0], |
---|
3836 | rowArray_[2], |
---|
3837 | rowArray_[3], |
---|
3838 | rowArray_[1]); |
---|
3839 | // see if update stable |
---|
3840 | #ifdef CLP_DEBUG |
---|
3841 | if ((handler_->logLevel() & 32)) |
---|
3842 | printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_); |
---|
3843 | #endif |
---|
3844 | double checkValue = 1.0e-7; |
---|
3845 | // if can't trust much and long way from optimal then relax |
---|
3846 | if (largestPrimalError_ > 10.0) |
---|
3847 | checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_); |
---|
3848 | if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || |
---|
3849 | fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) { |
---|
3850 | handler_->message(CLP_DUAL_CHECK, messages_) |
---|
3851 | << btranAlpha |
---|
3852 | << alpha_ |
---|
3853 | << CoinMessageEol; |
---|
3854 | // clear arrays |
---|
3855 | rowArray_[4]->clear(); |
---|
3856 | if (factorization_->pivots()) { |
---|
3857 | dualRowPivot_->unrollWeights(); |
---|
3858 | problemStatus_ = -2; // factorize now |
---|
3859 | rowArray_[0]->clear(); |
---|
3860 | rowArray_[1]->clear(); |
---|
3861 | columnArray_[0]->clear(); |
---|
3862 | returnCode = -2; |
---|
3863 | break; |
---|
3864 | } else { |
---|
3865 | // take on more relaxed criterion |
---|
3866 | double test; |
---|
3867 | if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8) |
---|
3868 | test = 1.0e-1 * fabs(alpha_); |
---|
3869 | else |
---|
3870 | test = 1.0e-4 * (1.0 + fabs(alpha_)); |
---|
3871 | if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || |
---|
3872 | fabs(btranAlpha - alpha_) > test) { |
---|
3873 | dualRowPivot_->unrollWeights(); |
---|
3874 | // need to reject something |
---|
3875 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
3876 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
3877 | << x << sequenceWithin(sequenceOut_) |
---|
3878 | << CoinMessageEol; |
---|
3879 | setFlagged(sequenceOut_); |
---|
3880 | progress_.clearBadTimes(); |
---|
3881 | lastBadIteration_ = numberIterations_; // say be more cautious |
---|
3882 | rowArray_[0]->clear(); |
---|
3883 | rowArray_[1]->clear(); |
---|
3884 | columnArray_[0]->clear(); |
---|
3885 | if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) { |
---|
3886 | //printf("I think should declare infeasible\n"); |
---|
3887 | problemStatus_ = 1; |
---|
3888 | returnCode = 1; |
---|
3889 | break; |
---|
3890 | } |
---|
3891 | continue; |
---|
3892 | } |
---|
3893 | } |
---|
3894 | } |
---|
3895 | // update duals BEFORE replaceColumn so can do updateColumn |
---|
3896 | double objectiveChange = 0.0; |
---|
3897 | // do duals first as variables may flip bounds |
---|
3898 | // rowArray_[0] and columnArray_[0] may have flips |
---|
3899 | // so use rowArray_[3] for work array from here on |
---|
3900 | int nswapped = 0; |
---|
3901 | //rowArray_[0]->cleanAndPackSafe(1.0e-60); |
---|
3902 | //columnArray_[0]->cleanAndPackSafe(1.0e-60); |
---|
3903 | #if CLP_CAN_HAVE_ZERO_OBJ |
---|
3904 | if ((specialOptions_&2097152)==0) { |
---|
3905 | #endif |
---|
3906 | nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], |
---|
3907 | rowArray_[2], theta_, |
---|
3908 | objectiveChange, false); |
---|
3909 | assert (!nswapped); |
---|
3910 | #if CLP_CAN_HAVE_ZERO_OBJ |
---|
3911 | } else { |
---|
3912 | rowArray_[0]->clear(); |
---|
3913 | rowArray_[2]->clear(); |
---|
3914 | columnArray_[0]->clear(); |
---|
3915 | } |
---|
3916 | #endif |
---|
3917 | // which will change basic solution |
---|
3918 | if (nswapped) { |
---|
3919 | abort(); //needs testing |
---|
3920 | factorization_->updateColumn(rowArray_[3], rowArray_[2]); |
---|
3921 | dualRowPivot_->updatePrimalSolution(rowArray_[2], |
---|
3922 | 1.0, objectiveChange); |
---|
3923 | // recompute dualOut_ |
---|
3924 | valueOut_ = solution_[sequenceOut_]; |
---|
3925 | if (directionOut_ < 0) { |
---|
3926 | dualOut_ = valueOut_ - upperOut_; |
---|
3927 | } else { |
---|
3928 | dualOut_ = lowerOut_ - valueOut_; |
---|
3929 | } |
---|
3930 | } |
---|
3931 | // amount primal will move |
---|
3932 | double movement = -dualOut_ * directionOut_ / alpha_; |
---|
3933 | // so objective should increase by fabs(dj)*movement |
---|
3934 | // but we already have objective change - so check will be good |
---|
3935 | if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) { |
---|
3936 | #ifdef CLP_DEBUG |
---|
3937 | if (handler_->logLevel() & 32) |
---|
3938 | printf("movement %g, swap change %g, rest %g * %g\n", |
---|
3939 | objectiveChange + fabs(movement * dualIn_), |
---|
3940 | objectiveChange, movement, dualIn_); |
---|
3941 | #endif |
---|
3942 | assert (objectiveChange + fabs(movement * dualIn_) >= -1.0e-5); |
---|
3943 | if(factorization_->pivots()) { |
---|
3944 | // going backwards - factorize |
---|
3945 | dualRowPivot_->unrollWeights(); |
---|
3946 | problemStatus_ = -2; // factorize now |
---|
3947 | returnCode = -2; |
---|
3948 | break; |
---|
3949 | } |
---|
3950 | } |
---|
3951 | CoinAssert(fabs(dualOut_) < 1.0e50); |
---|
3952 | // if stable replace in basis |
---|
3953 | int updateStatus = factorization_->replaceColumn(this, |
---|
3954 | rowArray_[2], |
---|
3955 | rowArray_[1], |
---|
3956 | pivotRow_, |
---|
3957 | alpha_); |
---|
3958 | // if no pivots, bad update but reasonable alpha - take and invert |
---|
3959 | if (updateStatus == 2 && |
---|
3960 | !factorization_->pivots() && fabs(alpha_) > 1.0e-5) |
---|
3961 | updateStatus = 4; |
---|
3962 | if (updateStatus == 1 || updateStatus == 4) { |
---|
3963 | // slight error |
---|
3964 | if (factorization_->pivots() > 5 || updateStatus == 4) { |
---|
3965 | problemStatus_ = -2; // factorize now |
---|
3966 | returnCode = -3; |
---|
3967 | } |
---|
3968 | } else if (updateStatus == 2) { |
---|
3969 | // major error |
---|
3970 | dualRowPivot_->unrollWeights(); |
---|
3971 | // later we may need to unwind more e.g. fake bounds |
---|
3972 | if (factorization_->pivots()) { |
---|
3973 | problemStatus_ = -2; // factorize now |
---|
3974 | returnCode = -2; |
---|
3975 | break; |
---|
3976 | } else { |
---|
3977 | // need to reject something |
---|
3978 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
3979 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
3980 | << x << sequenceWithin(sequenceOut_) |
---|
3981 | << CoinMessageEol; |
---|
3982 | setFlagged(sequenceOut_); |
---|
3983 | progress_.clearBadTimes(); |
---|
3984 | lastBadIteration_ = numberIterations_; // say be more cautious |
---|
3985 | rowArray_[0]->clear(); |
---|
3986 | rowArray_[1]->clear(); |
---|
3987 | columnArray_[0]->clear(); |
---|
3988 | // make sure dual feasible |
---|
3989 | // look at all rows and columns |
---|
3990 | double objectiveChange = 0.0; |
---|
3991 | reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1], |
---|
3992 | 0.0, objectiveChange, true); |
---|
3993 | continue; |
---|
3994 | } |
---|
3995 | } else if (updateStatus == 3) { |
---|
3996 | // out of memory |
---|
3997 | // increase space if not many iterations |
---|
3998 | if (factorization_->pivots() < |
---|
3999 | 0.5 * factorization_->maximumPivots() && |
---|
4000 | factorization_->pivots() < 200) |
---|
4001 | factorization_->areaFactor( |
---|
4002 | factorization_->areaFactor() * 1.1); |
---|
4003 | problemStatus_ = -2; // factorize now |
---|
4004 | } else if (updateStatus == 5) { |
---|
4005 | problemStatus_ = -2; // factorize now |
---|
4006 | } |
---|
4007 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4008 | int * lowerActive = paramData.lowerActive; |
---|
4009 | int * upperActive = paramData.upperActive; |
---|
4010 | #endif |
---|
4011 | // update change vector |
---|
4012 | { |
---|
4013 | double * work = rowArray_[1]->denseVector(); |
---|
4014 | int number = rowArray_[1]->getNumElements(); |
---|
4015 | int * which = rowArray_[1]->getIndices(); |
---|
4016 | assert (!rowArray_[4]->packedMode()); |
---|
4017 | #ifndef COIN_FAC_NEW |
---|
4018 | assert (rowArray_[1]->packedMode()); |
---|
4019 | #else |
---|
4020 | assert (!rowArray_[1]->packedMode()); |
---|
4021 | #endif |
---|
4022 | double pivotValue = rowArray_[4]->denseVector()[pivotRow_]; |
---|
4023 | double multiplier = -pivotValue/alpha_; |
---|
4024 | double * array=rowArray_[4]->denseVector(); |
---|
4025 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4026 | int lowerN=lowerActive[-1]; |
---|
4027 | int upperN=upperActive[-1]; |
---|
4028 | #endif |
---|
4029 | if (multiplier) { |
---|
4030 | for (int i = 0; i < number; i++) { |
---|
4031 | int iRow = which[i]; |
---|
4032 | #ifndef COIN_FAC_NEW |
---|
4033 | double alpha=multiplier*work[i]; |
---|
4034 | #else |
---|
4035 | double alpha=multiplier*work[iRow]; |
---|
4036 | #endif |
---|
4037 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4038 | double alpha3 = alpha+array[iRow]; |
---|
4039 | int iSequence = pivotVariable_[iRow]; |
---|
4040 | double oldLower = lowerCoefficient[iRow]; |
---|
4041 | double oldUpper = upperCoefficient[iRow]; |
---|
4042 | if (lower_[iSequence]>-1.0e30) { |
---|
4043 | //lowerGap[iRow]=value-lower_[iSequence]; |
---|
4044 | double alpha2 = alpha3 + lowerChange[iSequence]; |
---|
4045 | if (alpha2>1.0e-8) { |
---|
4046 | lowerCoefficient[iRow]=alpha2; |
---|
4047 | if (!oldLower) |
---|
4048 | lowerActive[lowerN++]=iRow; |
---|
4049 | } else { |
---|
4050 | if (oldLower) |
---|
4051 | lowerCoefficient[iRow]=COIN_DBL_MIN; |
---|
4052 | } |
---|
4053 | } else { |
---|
4054 | if (oldLower) |
---|
4055 | lowerCoefficient[iRow]=COIN_DBL_MIN; |
---|
4056 | } |
---|
4057 | if (upper_[iSequence]<1.0e30) { |
---|
4058 | //upperGap[iRow]=-(value-upper_[iSequence]); |
---|
4059 | double alpha2 = -(alpha3+upperChange[iSequence]); |
---|
4060 | if (alpha2>1.0e-8) { |
---|
4061 | upperCoefficient[iRow]=alpha2; |
---|
4062 | if (!oldUpper) |
---|
4063 | upperActive[upperN++]=iRow; |
---|
4064 | } else { |
---|
4065 | if (oldUpper) |
---|
4066 | upperCoefficient[iRow]=COIN_DBL_MIN; |
---|
4067 | } |
---|
4068 | } else { |
---|
4069 | if (oldUpper) |
---|
4070 | upperCoefficient[iRow]=COIN_DBL_MIN; |
---|
4071 | } |
---|
4072 | #endif |
---|
4073 | rowArray_[4]->quickAdd(iRow,alpha); |
---|
4074 | } |
---|
4075 | } |
---|
4076 | pivotValue = array[pivotRow_]; |
---|
4077 | // we want pivot to be -multiplier |
---|
4078 | rowArray_[4]->quickAdd(pivotRow_,-multiplier-pivotValue); |
---|
4079 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4080 | assert (lowerN>=0&&lowerN<=numberRows_); |
---|
4081 | lowerActive[-1]=lowerN; |
---|
4082 | upperActive[-1]=upperN; |
---|
4083 | #endif |
---|
4084 | } |
---|
4085 | // update primal solution |
---|
4086 | #if CLP_CAN_HAVE_ZERO_OBJ |
---|
4087 | if ((specialOptions_&2097152)!=0) |
---|
4088 | theta_=0.0; |
---|
4089 | #endif |
---|
4090 | if (theta_ < 0.0) { |
---|
4091 | #ifdef CLP_DEBUG |
---|
4092 | if (handler_->logLevel() & 32) |
---|
4093 | printf("negative theta %g\n", theta_); |
---|
4094 | #endif |
---|
4095 | theta_ = 0.0; |
---|
4096 | } |
---|
4097 | // do actual flips |
---|
4098 | reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]); |
---|
4099 | //rowArray_[1]->expand(); |
---|
4100 | #ifndef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4101 | dualRowPivot_->updatePrimalSolution(rowArray_[1], |
---|
4102 | movement, |
---|
4103 | objectiveChange); |
---|
4104 | #else |
---|
4105 | // do by hand |
---|
4106 | { |
---|
4107 | double * work = rowArray_[1]->denseVector(); |
---|
4108 | int number = rowArray_[1]->getNumElements(); |
---|
4109 | int * which = rowArray_[1]->getIndices(); |
---|
4110 | int i; |
---|
4111 | if (rowArray_[1]->packedMode()) { |
---|
4112 | for (i = 0; i < number; i++) { |
---|
4113 | int iRow = which[i]; |
---|
4114 | int iSequence = pivotVariable_[iRow]; |
---|
4115 | double value = solution_[iSequence]; |
---|
4116 | double change = movement * work[i]; |
---|
4117 | value -= change; |
---|
4118 | if (lower_[iSequence]>-1.0e30) |
---|
4119 | lowerGap[iRow]=value-lower_[iSequence]; |
---|
4120 | if (upper_[iSequence]<1.0e30) |
---|
4121 | upperGap[iRow]=-(value-upper_[iSequence]); |
---|
4122 | solution_[iSequence] = value; |
---|
4123 | objectiveChange -= change * cost_[iSequence]; |
---|
4124 | work[i] = 0.0; |
---|
4125 | } |
---|
4126 | } else { |
---|
4127 | for (i = 0; i < number; i++) { |
---|
4128 | int iRow = which[i]; |
---|
4129 | int iSequence = pivotVariable_[iRow]; |
---|
4130 | double value = solution_[iSequence]; |
---|
4131 | double change = movement * work[iRow]; |
---|
4132 | value -= change; |
---|
4133 | solution_[iSequence] = value; |
---|
4134 | objectiveChange -= change * cost_[iSequence]; |
---|
4135 | work[iRow] = 0.0; |
---|
4136 | } |
---|
4137 | } |
---|
4138 | rowArray_[1]->setNumElements(0); |
---|
4139 | } |
---|
4140 | #endif |
---|
4141 | // modify dualout |
---|
4142 | dualOut_ /= alpha_; |
---|
4143 | dualOut_ *= -directionOut_; |
---|
4144 | //setStatus(sequenceIn_,basic); |
---|
4145 | dj_[sequenceIn_] = 0.0; |
---|
4146 | //double oldValue = valueIn_; |
---|
4147 | if (directionIn_ == -1) { |
---|
4148 | // as if from upper bound |
---|
4149 | valueIn_ = upperIn_ + dualOut_; |
---|
4150 | } else { |
---|
4151 | // as if from lower bound |
---|
4152 | valueIn_ = lowerIn_ + dualOut_; |
---|
4153 | } |
---|
4154 | objectiveChange = 0.0; |
---|
4155 | #if CLP_CAN_HAVE_ZERO_OBJ |
---|
4156 | if ((specialOptions_&2097152)==0) { |
---|
4157 | #endif |
---|
4158 | for (int i=0;i<numberTotal;i++) |
---|
4159 | objectiveChange += solution_[i]*cost_[i]; |
---|
4160 | objectiveChange -= objectiveValue_; |
---|
4161 | #if CLP_CAN_HAVE_ZERO_OBJ |
---|
4162 | } |
---|
4163 | #endif |
---|
4164 | // outgoing |
---|
4165 | originalBound(sequenceOut_,useTheta,lowerChange,upperChange); |
---|
4166 | lowerOut_=lower_[sequenceOut_]; |
---|
4167 | upperOut_=upper_[sequenceOut_]; |
---|
4168 | // set dj to zero unless values pass |
---|
4169 | if (directionOut_ > 0) { |
---|
4170 | valueOut_ = lowerOut_; |
---|
4171 | dj_[sequenceOut_] = theta_; |
---|
4172 | #if CLP_CAN_HAVE_ZERO_OBJ>1 |
---|
4173 | #ifdef COIN_REUSE_RANDOM |
---|
4174 | if ((specialOptions_&2097152)!=0) { |
---|
4175 | dj_[sequenceOut_] = 1.0e-9*(1.0+CoinDrand48());; |
---|
4176 | } |
---|
4177 | #endif |
---|
4178 | #endif |
---|
4179 | } else { |
---|
4180 | valueOut_ = upperOut_; |
---|
4181 | dj_[sequenceOut_] = -theta_; |
---|
4182 | #if CLP_CAN_HAVE_ZERO_OBJ>1 |
---|
4183 | #ifdef COIN_REUSE_RANDOM |
---|
4184 | if ((specialOptions_&2097152)!=0) { |
---|
4185 | dj_[sequenceOut_] = -1.0e-9*(1.0+CoinDrand48());; |
---|
4186 | } |
---|
4187 | #endif |
---|
4188 | #endif |
---|
4189 | } |
---|
4190 | solution_[sequenceOut_] = valueOut_; |
---|
4191 | int whatNext = housekeeping(objectiveChange); |
---|
4192 | reinterpret_cast<ClpSimplexDual *>(this)->originalBound(sequenceIn_); |
---|
4193 | assert (backwardBasic[sequenceOut_]==pivotRow_); |
---|
4194 | backwardBasic[sequenceOut_]=-1; |
---|
4195 | backwardBasic[sequenceIn_]=pivotRow_; |
---|
4196 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4197 | double value = solution_[sequenceIn_]; |
---|
4198 | double alpha = rowArray_[4]->denseVector()[pivotRow_]; |
---|
4199 | double oldLower = lowerCoefficient[pivotRow_]; |
---|
4200 | double oldUpper = upperCoefficient[pivotRow_]; |
---|
4201 | if (lower_[sequenceIn_]>-1.0e30) { |
---|
4202 | lowerGap[pivotRow_]=value-lower_[sequenceIn_]; |
---|
4203 | double alpha2 = alpha + lowerChange[sequenceIn_]; |
---|
4204 | if (alpha2>1.0e-8) { |
---|
4205 | lowerCoefficient[pivotRow_]=alpha2; |
---|
4206 | if (!oldLower) { |
---|
4207 | int lowerN=lowerActive[-1]; |
---|
4208 | assert (lowerN>=0&&lowerN<numberRows_); |
---|
4209 | lowerActive[lowerN]=pivotRow_; |
---|
4210 | lowerActive[-1]=lowerN+1; |
---|
4211 | } |
---|
4212 | } else { |
---|
4213 | if (oldLower) |
---|
4214 | lowerCoefficient[pivotRow_]=COIN_DBL_MIN; |
---|
4215 | } |
---|
4216 | } else { |
---|
4217 | if (oldLower) |
---|
4218 | lowerCoefficient[pivotRow_]=COIN_DBL_MIN; |
---|
4219 | } |
---|
4220 | if (upper_[sequenceIn_]<1.0e30) { |
---|
4221 | upperGap[pivotRow_]=-(value-upper_[sequenceIn_]); |
---|
4222 | double alpha2 = -(alpha+upperChange[sequenceIn_]); |
---|
4223 | if (alpha2>1.0e-8) { |
---|
4224 | upperCoefficient[pivotRow_]=alpha2; |
---|
4225 | if (!oldUpper) { |
---|
4226 | int upperN=upperActive[-1]; |
---|
4227 | assert (upperN>=0&&upperN<numberRows_); |
---|
4228 | upperActive[upperN]=pivotRow_; |
---|
4229 | upperActive[-1]=upperN+1; |
---|
4230 | } |
---|
4231 | } else { |
---|
4232 | if (oldUpper) |
---|
4233 | upperCoefficient[pivotRow_]=COIN_DBL_MIN; |
---|
4234 | } |
---|
4235 | } else { |
---|
4236 | if (oldUpper) |
---|
4237 | upperCoefficient[pivotRow_]=COIN_DBL_MIN; |
---|
4238 | } |
---|
4239 | #endif |
---|
4240 | { |
---|
4241 | char in[200],out[200]; |
---|
4242 | int iSequence=sequenceIn_; |
---|
4243 | if (iSequence<numberColumns_) { |
---|
4244 | if (lengthNames_) |
---|
4245 | strcpy(in,columnNames_[iSequence].c_str()); |
---|
4246 | else |
---|
4247 | sprintf(in,"C%7.7d",iSequence); |
---|
4248 | } else { |
---|
4249 | iSequence -= numberColumns_; |
---|
4250 | if (lengthNames_) |
---|
4251 | strcpy(in,rowNames_[iSequence].c_str()); |
---|
4252 | else |
---|
4253 | sprintf(in,"R%7.7d",iSequence); |
---|
4254 | } |
---|
4255 | iSequence=sequenceOut_; |
---|
4256 | if (iSequence<numberColumns_) { |
---|
4257 | if (lengthNames_) |
---|
4258 | strcpy(out,columnNames_[iSequence].c_str()); |
---|
4259 | else |
---|
4260 | sprintf(out,"C%7.7d",iSequence); |
---|
4261 | } else { |
---|
4262 | iSequence -= numberColumns_; |
---|
4263 | if (lengthNames_) |
---|
4264 | strcpy(out,rowNames_[iSequence].c_str()); |
---|
4265 | else |
---|
4266 | sprintf(out,"R%7.7d",iSequence); |
---|
4267 | } |
---|
4268 | handler_->message(CLP_PARAMETRICS_STATS2, messages_) |
---|
4269 | << useTheta << objectiveValue() |
---|
4270 | << in << out << CoinMessageEol; |
---|
4271 | } |
---|
4272 | if (useTheta>lastTheta+1.0e-9) { |
---|
4273 | handler_->message(CLP_PARAMETRICS_STATS, messages_) |
---|
4274 | << useTheta << objectiveValue() << CoinMessageEol; |
---|
4275 | lastTheta = useTheta; |
---|
4276 | } |
---|
4277 | // and set bounds correctly |
---|
4278 | originalBound(sequenceIn_,useTheta,lowerChange,upperChange); |
---|
4279 | reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_); |
---|
4280 | if (whatNext == 1) { |
---|
4281 | problemStatus_ = -2; // refactorize |
---|
4282 | } else if (whatNext == 2) { |
---|
4283 | // maximum iterations or equivalent |
---|
4284 | problemStatus_ = 3; |
---|
4285 | returnCode = 3; |
---|
4286 | break; |
---|
4287 | } |
---|
4288 | #ifdef CLP_USER_DRIVEN |
---|
4289 | // Check event |
---|
4290 | { |
---|
4291 | int status = eventHandler_->event(ClpEventHandler::endOfIteration); |
---|
4292 | if (status >= 0) { |
---|
4293 | problemStatus_ = 5; |
---|
4294 | secondaryStatus_ = ClpEventHandler::endOfIteration; |
---|
4295 | returnCode = 4; |
---|
4296 | break; |
---|
4297 | } |
---|
4298 | } |
---|
4299 | #endif |
---|
4300 | } else { |
---|
4301 | // no incoming column is valid |
---|
4302 | #ifdef CLP_USER_DRIVEN |
---|
4303 | rowArray_[0]->clear(); |
---|
4304 | columnArray_[0]->clear(); |
---|
4305 | theta_ = useTheta; |
---|
4306 | lastTheta = useTheta; |
---|
4307 | int action = eventHandler_->event(ClpEventHandler::noTheta); |
---|
4308 | if (action>=0) { |
---|
4309 | endingTheta=theta_; |
---|
4310 | theta_ = 0.0; |
---|
4311 | //adjust [4] from handler - but |
---|
4312 | //rowArray_[4]->clear(); // temp |
---|
4313 | if (action>=0&&action<10) |
---|
4314 | problemStatus_=-1; // carry on |
---|
4315 | else if (action==15) |
---|
4316 | problemStatus_ =5; // say stopped |
---|
4317 | returnCode = 1; |
---|
4318 | if (action==0||action>=10) |
---|
4319 | break; |
---|
4320 | else |
---|
4321 | continue; |
---|
4322 | } else { |
---|
4323 | theta_ = 0.0; |
---|
4324 | } |
---|
4325 | #endif |
---|
4326 | pivotRow_ = -1; |
---|
4327 | #ifdef CLP_DEBUG |
---|
4328 | if (handler_->logLevel() & 32) |
---|
4329 | printf("** no column pivot\n"); |
---|
4330 | #endif |
---|
4331 | if (factorization_->pivots() < 10) { |
---|
4332 | // If we have just factorized and infeasibility reasonable say infeas |
---|
4333 | if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) { |
---|
4334 | if (valueOut_ > upperOut_ + 1.0e-3 || valueOut_ < lowerOut_ - 1.0e-3 |
---|
4335 | || (specialOptions_ & 64) == 0) { |
---|
4336 | // say infeasible |
---|
4337 | problemStatus_ = 1; |
---|
4338 | // unless primal feasible!!!! |
---|
4339 | //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_, |
---|
4340 | // numberDualInfeasibilities_,sumDualInfeasibilities_); |
---|
4341 | if (numberDualInfeasibilities_) |
---|
4342 | problemStatus_ = 10; |
---|
4343 | rowArray_[0]->clear(); |
---|
4344 | columnArray_[0]->clear(); |
---|
4345 | } |
---|
4346 | } |
---|
4347 | // If special option set - put off as long as possible |
---|
4348 | if ((specialOptions_ & 64) == 0) { |
---|
4349 | problemStatus_ = -4; //say looks infeasible |
---|
4350 | } else { |
---|
4351 | // flag |
---|
4352 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
4353 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
4354 | << x << sequenceWithin(sequenceOut_) |
---|
4355 | << CoinMessageEol; |
---|
4356 | setFlagged(sequenceOut_); |
---|
4357 | if (!factorization_->pivots()) { |
---|
4358 | rowArray_[0]->clear(); |
---|
4359 | columnArray_[0]->clear(); |
---|
4360 | continue; |
---|
4361 | } |
---|
4362 | } |
---|
4363 | } |
---|
4364 | rowArray_[0]->clear(); |
---|
4365 | columnArray_[0]->clear(); |
---|
4366 | returnCode = 1; |
---|
4367 | break; |
---|
4368 | } |
---|
4369 | } else { |
---|
4370 | // no pivot row |
---|
4371 | #ifdef CLP_USER_DRIVEN |
---|
4372 | { |
---|
4373 | double saveTheta=theta_; |
---|
4374 | theta_ = endingTheta; |
---|
4375 | int status=eventHandler_->event(ClpEventHandler::theta); |
---|
4376 | if (status>=0&&status<10) { |
---|
4377 | endingTheta=theta_; |
---|
4378 | theta_=saveTheta; |
---|
4379 | continue; |
---|
4380 | } else { |
---|
4381 | theta_=saveTheta; |
---|
4382 | } |
---|
4383 | } |
---|
4384 | #endif |
---|
4385 | #ifdef CLP_DEBUG |
---|
4386 | if (handler_->logLevel() & 32) |
---|
4387 | printf("** no row pivot\n"); |
---|
4388 | #endif |
---|
4389 | int numberPivots = factorization_->pivots(); |
---|
4390 | bool specialCase; |
---|
4391 | int useNumberFake; |
---|
4392 | returnCode = 0; |
---|
4393 | if (numberPivots < 20 && |
---|
4394 | (specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100 |
---|
4395 | && dualBound_ > 1.0e8) { |
---|
4396 | specialCase = true; |
---|
4397 | // as dual bound high - should be okay |
---|
4398 | useNumberFake = 0; |
---|
4399 | } else { |
---|
4400 | specialCase = false; |
---|
4401 | useNumberFake = numberFake_; |
---|
4402 | } |
---|
4403 | if (!numberPivots || specialCase) { |
---|
4404 | // may have crept through - so may be optimal |
---|
4405 | // check any flagged variables |
---|
4406 | int iRow; |
---|
4407 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
4408 | int iPivot = pivotVariable_[iRow]; |
---|
4409 | if (flagged(iPivot)) |
---|
4410 | break; |
---|
4411 | } |
---|
4412 | if (iRow < numberRows_ && numberPivots) { |
---|
4413 | // try factorization |
---|
4414 | returnCode = -2; |
---|
4415 | } |
---|
4416 | |
---|
4417 | if (useNumberFake || numberDualInfeasibilities_) { |
---|
4418 | // may be dual infeasible |
---|
4419 | problemStatus_ = -5; |
---|
4420 | } else { |
---|
4421 | if (iRow < numberRows_) { |
---|
4422 | problemStatus_ = -5; |
---|
4423 | } else { |
---|
4424 | if (numberPivots) { |
---|
4425 | // objective may be wrong |
---|
4426 | objectiveValue_ = innerProduct(cost_, |
---|
4427 | numberColumns_ + numberRows_, |
---|
4428 | solution_); |
---|
4429 | objectiveValue_ += objective_->nonlinearOffset(); |
---|
4430 | objectiveValue_ /= (objectiveScale_ * rhsScale_); |
---|
4431 | if ((specialOptions_ & 16384) == 0) { |
---|
4432 | // and dual_ may be wrong (i.e. for fixed or basic) |
---|
4433 | CoinIndexedVector * arrayVector = rowArray_[1]; |
---|
4434 | arrayVector->clear(); |
---|
4435 | int iRow; |
---|
4436 | double * array = arrayVector->denseVector(); |
---|
4437 | /* Use dual_ instead of array |
---|
4438 | Even though dual_ is only numberRows_ long this is |
---|
4439 | okay as gets permuted to longer rowArray_[2] |
---|
4440 | */ |
---|
4441 | arrayVector->setDenseVector(dual_); |
---|
4442 | int * index = arrayVector->getIndices(); |
---|
4443 | int number = 0; |
---|
4444 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
4445 | int iPivot = pivotVariable_[iRow]; |
---|
4446 | double value = cost_[iPivot]; |
---|
4447 | dual_[iRow] = value; |
---|
4448 | if (value) { |
---|
4449 | index[number++] = iRow; |
---|
4450 | } |
---|
4451 | } |
---|
4452 | arrayVector->setNumElements(number); |
---|
4453 | // Extended duals before "updateTranspose" |
---|
4454 | matrix_->dualExpanded(this, arrayVector, NULL, 0); |
---|
4455 | // Btran basic costs |
---|
4456 | rowArray_[2]->clear(); |
---|
4457 | factorization_->updateColumnTranspose(rowArray_[2], arrayVector); |
---|
4458 | // and return vector |
---|
4459 | arrayVector->setDenseVector(array); |
---|
4460 | } |
---|
4461 | } |
---|
4462 | problemStatus_ = 0; |
---|
4463 | sumPrimalInfeasibilities_ = 0.0; |
---|
4464 | if ((specialOptions_&(1024 + 16384)) != 0) { |
---|
4465 | CoinIndexedVector * arrayVector = rowArray_[1]; |
---|
4466 | arrayVector->clear(); |
---|
4467 | double * rhs = arrayVector->denseVector(); |
---|
4468 | times(1.0, solution_, rhs); |
---|
4469 | bool bad2 = false; |
---|
4470 | int i; |
---|
4471 | for ( i = 0; i < numberRows_; i++) { |
---|
4472 | if (rhs[i] < rowLowerWork_[i] - primalTolerance_ || |
---|
4473 | rhs[i] > rowUpperWork_[i] + primalTolerance_) { |
---|
4474 | bad2 = true; |
---|
4475 | } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) { |
---|
4476 | } |
---|
4477 | rhs[i] = 0.0; |
---|
4478 | } |
---|
4479 | for ( i = 0; i < numberColumns_; i++) { |
---|
4480 | if (solution_[i] < columnLowerWork_[i] - primalTolerance_ || |
---|
4481 | solution_[i] > columnUpperWork_[i] + primalTolerance_) { |
---|
4482 | bad2 = true; |
---|
4483 | } |
---|
4484 | } |
---|
4485 | if (bad2) { |
---|
4486 | problemStatus_ = -3; |
---|
4487 | returnCode = -2; |
---|
4488 | // Force to re-factorize early next time |
---|
4489 | int numberPivots = factorization_->pivots(); |
---|
4490 | forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); |
---|
4491 | } |
---|
4492 | } |
---|
4493 | } |
---|
4494 | } |
---|
4495 | } else { |
---|
4496 | problemStatus_ = -3; |
---|
4497 | returnCode = -2; |
---|
4498 | // Force to re-factorize early next time |
---|
4499 | int numberPivots = factorization_->pivots(); |
---|
4500 | forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); |
---|
4501 | } |
---|
4502 | break; |
---|
4503 | } |
---|
4504 | } |
---|
4505 | startingTheta = lastTheta+theta_; |
---|
4506 | return returnCode; |
---|
4507 | } |
---|
4508 | // Compute new rowLower_ etc (return negative if infeasible - otherwise largest change) |
---|
4509 | double |
---|
4510 | ClpSimplexOther::computeRhsEtc( parametricsData & paramData) |
---|
4511 | { |
---|
4512 | double maxTheta = COIN_DBL_MAX; |
---|
4513 | double largestChange=0.0; |
---|
4514 | double startingTheta = paramData.startingTheta; |
---|
4515 | const double * lowerChange = paramData.lowerChange+ |
---|
4516 | paramData.unscaledChangesOffset; |
---|
4517 | const double * upperChange = paramData.upperChange+ |
---|
4518 | paramData.unscaledChangesOffset; |
---|
4519 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
4520 | double lower = rowLower_[iRow]; |
---|
4521 | double upper = rowUpper_[iRow]; |
---|
4522 | double chgLower = lowerChange[numberColumns_+iRow]; |
---|
4523 | largestChange=CoinMax(largestChange,fabs(chgLower)); |
---|
4524 | double chgUpper = upperChange[numberColumns_+iRow]; |
---|
4525 | largestChange=CoinMax(largestChange,fabs(chgUpper)); |
---|
4526 | if (lower > -1.0e30 && upper < 1.0e30) { |
---|
4527 | if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) { |
---|
4528 | maxTheta = (upper - lower) / (chgLower - chgUpper); |
---|
4529 | } |
---|
4530 | } |
---|
4531 | lower+=startingTheta*chgLower; |
---|
4532 | upper+=startingTheta*chgUpper; |
---|
4533 | #ifndef CLP_USER_DRIVEN |
---|
4534 | if (lower > upper) { |
---|
4535 | maxTheta = -1.0; |
---|
4536 | break; |
---|
4537 | } |
---|
4538 | #endif |
---|
4539 | rowLower_[iRow]=lower; |
---|
4540 | rowUpper_[iRow]=upper; |
---|
4541 | } |
---|
4542 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
4543 | double lower = columnLower_[iColumn]; |
---|
4544 | double upper = columnUpper_[iColumn]; |
---|
4545 | double chgLower = lowerChange[iColumn]; |
---|
4546 | largestChange=CoinMax(largestChange,fabs(chgLower)); |
---|
4547 | double chgUpper = upperChange[iColumn]; |
---|
4548 | largestChange=CoinMax(largestChange,fabs(chgUpper)); |
---|
4549 | if (lower > -1.0e30 && upper < 1.0e30) { |
---|
4550 | if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) { |
---|
4551 | maxTheta = (upper - lower) / (chgLower - chgUpper); |
---|
4552 | } |
---|
4553 | } |
---|
4554 | lower+=startingTheta*chgLower; |
---|
4555 | upper+=startingTheta*chgUpper; |
---|
4556 | #ifndef CLP_USER_DRIVEN |
---|
4557 | if (lower > upper) { |
---|
4558 | maxTheta = -1.0; |
---|
4559 | break; |
---|
4560 | } |
---|
4561 | #endif |
---|
4562 | columnLower_[iColumn]=lower; |
---|
4563 | columnUpper_[iColumn]=upper; |
---|
4564 | } |
---|
4565 | #ifndef CLP_USER_DRIVEN |
---|
4566 | paramData.maxTheta=maxTheta; |
---|
4567 | if (maxTheta<0) |
---|
4568 | largestChange=-1.0; // signal infeasible |
---|
4569 | #else |
---|
4570 | // maxTheta already set |
---|
4571 | /* given largest change element choose acceptable end |
---|
4572 | be safe and make sure difference < 0.1*tolerance */ |
---|
4573 | double acceptableDifference=0.1*primalTolerance_/ |
---|
4574 | CoinMax(largestChange,1.0); |
---|
4575 | paramData.acceptableMaxTheta=maxTheta-acceptableDifference; |
---|
4576 | #endif |
---|
4577 | return largestChange; |
---|
4578 | } |
---|
4579 | // Redo lower_ from rowLower_ etc |
---|
4580 | void |
---|
4581 | ClpSimplexOther::redoInternalArrays() |
---|
4582 | { |
---|
4583 | double * lowerSave = lower_; |
---|
4584 | double * upperSave = upper_; |
---|
4585 | memcpy(lowerSave,columnLower_,numberColumns_*sizeof(double)); |
---|
4586 | memcpy(lowerSave+numberColumns_,rowLower_,numberRows_*sizeof(double)); |
---|
4587 | memcpy(upperSave,columnUpper_,numberColumns_*sizeof(double)); |
---|
4588 | memcpy(upperSave+numberColumns_,rowUpper_,numberRows_*sizeof(double)); |
---|
4589 | if (rowScale_) { |
---|
4590 | // scale arrays |
---|
4591 | for (int i=0;i<numberColumns_;i++) { |
---|
4592 | double multiplier = inverseColumnScale_[i]; |
---|
4593 | if (lowerSave[i]>-1.0e20) |
---|
4594 | lowerSave[i] *= multiplier; |
---|
4595 | if (upperSave[i]<1.0e20) |
---|
4596 | upperSave[i] *= multiplier; |
---|
4597 | } |
---|
4598 | lowerSave += numberColumns_; |
---|
4599 | upperSave += numberColumns_; |
---|
4600 | for (int i=0;i<numberRows_;i++) { |
---|
4601 | double multiplier = rowScale_[i]; |
---|
4602 | if (lowerSave[i]>-1.0e20) |
---|
4603 | lowerSave[i] *= multiplier; |
---|
4604 | if (upperSave[i]<1.0e20) |
---|
4605 | upperSave[i] *= multiplier; |
---|
4606 | } |
---|
4607 | } |
---|
4608 | } |
---|
4609 | #if 0 |
---|
4610 | static int zzzzzz=0; |
---|
4611 | int zzzzzzOther=0; |
---|
4612 | #endif |
---|
4613 | // Finds best possible pivot |
---|
4614 | double |
---|
4615 | ClpSimplexOther::bestPivot(bool justColumns) |
---|
4616 | { |
---|
4617 | // Get good size for pivot |
---|
4618 | // Allow first few iterations to take tiny |
---|
4619 | double acceptablePivot = 1.0e-9; |
---|
4620 | if (numberIterations_ > 100) |
---|
4621 | acceptablePivot = 1.0e-8; |
---|
4622 | if (factorization_->pivots() > 10 || |
---|
4623 | (factorization_->pivots() && sumDualInfeasibilities_)) |
---|
4624 | acceptablePivot = 1.0e-5; // if we have iterated be more strict |
---|
4625 | else if (factorization_->pivots() > 5) |
---|
4626 | acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict |
---|
4627 | else if (factorization_->pivots()) |
---|
4628 | acceptablePivot = 1.0e-8; // relax |
---|
4629 | double bestPossiblePivot = 1.0; |
---|
4630 | // get sign for finding row of tableau |
---|
4631 | // normal iteration |
---|
4632 | // create as packed |
---|
4633 | double direction = directionOut_; |
---|
4634 | #ifndef COIN_FAC_NEW |
---|
4635 | rowArray_[0]->createPacked(1, &pivotRow_, &direction); |
---|
4636 | #else |
---|
4637 | rowArray_[0]->createOneUnpackedElement(pivotRow_, direction); |
---|
4638 | #endif |
---|
4639 | factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]); |
---|
4640 | // put row of tableau in rowArray[0] and columnArray[0] |
---|
4641 | matrix_->transposeTimes(this, -1.0, |
---|
4642 | rowArray_[0], rowArray_[3], columnArray_[0]); |
---|
4643 | sequenceIn_=-1; |
---|
4644 | if (justColumns) |
---|
4645 | rowArray_[0]->clear(); |
---|
4646 | // do ratio test for normal iteration |
---|
4647 | bestPossiblePivot = |
---|
4648 | reinterpret_cast<ClpSimplexDual *> |
---|
4649 | ( this)->dualColumn(rowArray_[0], |
---|
4650 | columnArray_[0], columnArray_[1], |
---|
4651 | rowArray_[3], acceptablePivot, NULL); |
---|
4652 | return bestPossiblePivot; |
---|
4653 | } |
---|
4654 | // Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none) |
---|
4655 | int |
---|
4656 | ClpSimplexOther::nextTheta(int /*type*/, double maxTheta, parametricsData & paramData, |
---|
4657 | const double * /*changeObjective*/) |
---|
4658 | { |
---|
4659 | const double * lowerChange = paramData.lowerChange; |
---|
4660 | const double * upperChange = paramData.upperChange; |
---|
4661 | const int * lowerList = paramData.lowerList; |
---|
4662 | const int * upperList = paramData.upperList; |
---|
4663 | int iSequence; |
---|
4664 | bool toLower = false; |
---|
4665 | //assert (type==1); |
---|
4666 | // may need to decide based on model? |
---|
4667 | bool needFullUpdate = rowArray_[4]->getNumElements()==0; |
---|
4668 | double * array = rowArray_[4]->denseVector(); |
---|
4669 | //rowArray_[4]->checkClean(); |
---|
4670 | const int * row = matrix_->getIndices(); |
---|
4671 | const int * columnLength = matrix_->getVectorLengths(); |
---|
4672 | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
---|
4673 | const double * elementByColumn = matrix_->getElements(); |
---|
4674 | #if 0 |
---|
4675 | double tempArray[5000]; |
---|
4676 | bool checkIt=false; |
---|
4677 | if (factorization_->pivots()&&!needFullUpdate&&sequenceIn_<0) { |
---|
4678 | memcpy(tempArray,array,numberRows_*sizeof(double)); |
---|
4679 | checkIt=true; |
---|
4680 | needFullUpdate=true; |
---|
4681 | } |
---|
4682 | #endif |
---|
4683 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4684 | double * lowerGap = paramData.lowerGap; |
---|
4685 | double * upperGap = paramData.upperGap; |
---|
4686 | double * lowerCoefficient = paramData.lowerCoefficient; |
---|
4687 | double * upperCoefficient = paramData.upperCoefficient; |
---|
4688 | int * lowerActive=paramData.lowerActive; |
---|
4689 | int * upperActive=paramData.upperActive; |
---|
4690 | #endif |
---|
4691 | if (!factorization_->pivots()||needFullUpdate) { |
---|
4692 | //zzzzzz=0; |
---|
4693 | rowArray_[4]->clear(); |
---|
4694 | // get change |
---|
4695 | if (!rowScale_) { |
---|
4696 | int n; |
---|
4697 | n=lowerList[-2]; |
---|
4698 | int i; |
---|
4699 | for (i=0;i<n;i++) { |
---|
4700 | int iSequence = lowerList[i]; |
---|
4701 | assert (iSequence<numberColumns_); |
---|
4702 | if (getColumnStatus(iSequence)==atLowerBound) { |
---|
4703 | double value=lowerChange[iSequence]; |
---|
4704 | for (CoinBigIndex j = columnStart[iSequence]; |
---|
4705 | j < columnStart[iSequence] + columnLength[iSequence]; j++) { |
---|
4706 | rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value); |
---|
4707 | } |
---|
4708 | } |
---|
4709 | } |
---|
4710 | n=lowerList[-1]; |
---|
4711 | const double * change = lowerChange+numberColumns_; |
---|
4712 | for (;i<n;i++) { |
---|
4713 | int iSequence = lowerList[i]-numberColumns_; |
---|
4714 | assert (iSequence>=0); |
---|
4715 | if (getRowStatus(iSequence)==atLowerBound) { |
---|
4716 | double value=change[iSequence]; |
---|
4717 | rowArray_[4]->quickAdd(iSequence, -value); |
---|
4718 | } |
---|
4719 | } |
---|
4720 | n=upperList[-2]; |
---|
4721 | for (i=0;i<n;i++) { |
---|
4722 | int iSequence = upperList[i]; |
---|
4723 | assert (iSequence<numberColumns_); |
---|
4724 | if (getColumnStatus(iSequence)==atUpperBound) { |
---|
4725 | double value=upperChange[iSequence]; |
---|
4726 | for (CoinBigIndex j = columnStart[iSequence]; |
---|
4727 | j < columnStart[iSequence] + columnLength[iSequence]; j++) { |
---|
4728 | rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value); |
---|
4729 | } |
---|
4730 | } |
---|
4731 | } |
---|
4732 | n=upperList[-1]; |
---|
4733 | change = upperChange+numberColumns_; |
---|
4734 | for (;i<n;i++) { |
---|
4735 | int iSequence = upperList[i]-numberColumns_; |
---|
4736 | assert (iSequence>=0); |
---|
4737 | if (getRowStatus(iSequence)==atUpperBound) { |
---|
4738 | double value=change[iSequence]; |
---|
4739 | rowArray_[4]->quickAdd(iSequence, -value); |
---|
4740 | } |
---|
4741 | } |
---|
4742 | } else { |
---|
4743 | int n; |
---|
4744 | n=lowerList[-2]; |
---|
4745 | int i; |
---|
4746 | for (i=0;i<n;i++) { |
---|
4747 | int iSequence = lowerList[i]; |
---|
4748 | assert (iSequence<numberColumns_); |
---|
4749 | if (getColumnStatus(iSequence)==atLowerBound) { |
---|
4750 | double value=lowerChange[iSequence]; |
---|
4751 | // apply scaling |
---|
4752 | double scale = columnScale_[iSequence]; |
---|
4753 | for (CoinBigIndex j = columnStart[iSequence]; |
---|
4754 | j < columnStart[iSequence] + columnLength[iSequence]; j++) { |
---|
4755 | int iRow = row[j]; |
---|
4756 | rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value); |
---|
4757 | } |
---|
4758 | } |
---|
4759 | } |
---|
4760 | n=lowerList[-1]; |
---|
4761 | const double * change = lowerChange+numberColumns_; |
---|
4762 | for (;i<n;i++) { |
---|
4763 | int iSequence = lowerList[i]-numberColumns_; |
---|
4764 | assert (iSequence>=0); |
---|
4765 | if (getRowStatus(iSequence)==atLowerBound) { |
---|
4766 | double value=change[iSequence]; |
---|
4767 | rowArray_[4]->quickAdd(iSequence, -value); |
---|
4768 | } |
---|
4769 | } |
---|
4770 | n=upperList[-2]; |
---|
4771 | for (i=0;i<n;i++) { |
---|
4772 | int iSequence = upperList[i]; |
---|
4773 | assert (iSequence<numberColumns_); |
---|
4774 | if (getColumnStatus(iSequence)==atUpperBound) { |
---|
4775 | double value=upperChange[iSequence]; |
---|
4776 | // apply scaling |
---|
4777 | double scale = columnScale_[iSequence]; |
---|
4778 | for (CoinBigIndex j = columnStart[iSequence]; |
---|
4779 | j < columnStart[iSequence] + columnLength[iSequence]; j++) { |
---|
4780 | int iRow = row[j]; |
---|
4781 | rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value); |
---|
4782 | } |
---|
4783 | } |
---|
4784 | } |
---|
4785 | n=upperList[-1]; |
---|
4786 | change = upperChange+numberColumns_; |
---|
4787 | for (;i<n;i++) { |
---|
4788 | int iSequence = upperList[i]-numberColumns_; |
---|
4789 | assert (iSequence>=0); |
---|
4790 | if (getRowStatus(iSequence)==atUpperBound) { |
---|
4791 | double value=change[iSequence]; |
---|
4792 | rowArray_[4]->quickAdd(iSequence, -value); |
---|
4793 | } |
---|
4794 | } |
---|
4795 | } |
---|
4796 | // ftran it |
---|
4797 | factorization_->updateColumn(rowArray_[0], rowArray_[4]); |
---|
4798 | #if 0 |
---|
4799 | if (checkIt) { |
---|
4800 | for (int i=0;i<numberRows_;i++) { |
---|
4801 | assert (fabs(tempArray[i]-array[i])<1.0e-8); |
---|
4802 | } |
---|
4803 | } |
---|
4804 | #endif |
---|
4805 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4806 | /* later for sparse - keep like CoinIndexedvector |
---|
4807 | and just redo here */ |
---|
4808 | int lowerN=0; |
---|
4809 | int upperN=0; |
---|
4810 | memset(lowerCoefficient,0,numberRows_*sizeof(double)); |
---|
4811 | memset(upperCoefficient,0,numberRows_*sizeof(double)); |
---|
4812 | for (int iRow=0;iRow<numberRows_;iRow++) { |
---|
4813 | iSequence = pivotVariable_[iRow]; |
---|
4814 | double currentSolution = solution_[iSequence]; |
---|
4815 | double alpha = array[iRow]; |
---|
4816 | double thetaCoefficientLower = lowerChange[iSequence] + alpha; |
---|
4817 | double thetaCoefficientUpper = upperChange[iSequence] + alpha; |
---|
4818 | if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) { |
---|
4819 | double currentLower = lower_[iSequence]; |
---|
4820 | ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_); |
---|
4821 | double gap=currentSolution-currentLower; |
---|
4822 | lowerGap[iRow]=gap; |
---|
4823 | lowerCoefficient[iRow]=thetaCoefficientLower; |
---|
4824 | lowerActive[lowerN++]=iRow; |
---|
4825 | //} else { |
---|
4826 | //lowerCoefficient[iRow]=0.0; |
---|
4827 | } |
---|
4828 | if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) { |
---|
4829 | double currentUpper = upper_[iSequence]; |
---|
4830 | ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_); |
---|
4831 | double gap2=-(currentSolution-currentUpper); //positive |
---|
4832 | upperGap[iRow]=gap2; |
---|
4833 | upperCoefficient[iRow]=-thetaCoefficientUpper; |
---|
4834 | upperActive[upperN++]=iRow; |
---|
4835 | } |
---|
4836 | } |
---|
4837 | assert (lowerN>=0&&lowerN<=numberRows_); |
---|
4838 | lowerActive[-1]=lowerN; |
---|
4839 | upperActive[-1]=upperN; |
---|
4840 | #endif |
---|
4841 | } else if (sequenceIn_>=0) { |
---|
4842 | //assert (sequenceIn_>=0); |
---|
4843 | assert (sequenceOut_>=0); |
---|
4844 | assert (sequenceIn_!=sequenceOut_); |
---|
4845 | double change = (directionIn_>0) ? -lowerChange[sequenceIn_] : -upperChange[sequenceIn_]; |
---|
4846 | int needed=0; |
---|
4847 | assert (!rowArray_[5]->getNumElements()); |
---|
4848 | if (change) { |
---|
4849 | if (sequenceIn_<numberColumns_) { |
---|
4850 | if (!rowScale_) { |
---|
4851 | for (CoinBigIndex i = columnStart[sequenceIn_]; |
---|
4852 | i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) { |
---|
4853 | rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change); |
---|
4854 | } |
---|
4855 | } else { |
---|
4856 | // apply scaling |
---|
4857 | double scale = columnScale_[sequenceIn_]; |
---|
4858 | for (CoinBigIndex i = columnStart[sequenceIn_]; |
---|
4859 | i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) { |
---|
4860 | int iRow = row[i]; |
---|
4861 | rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change); |
---|
4862 | } |
---|
4863 | } |
---|
4864 | } else { |
---|
4865 | rowArray_[5]->insert(sequenceIn_-numberColumns_,-change); |
---|
4866 | } |
---|
4867 | needed++; |
---|
4868 | } |
---|
4869 | if (getStatus(sequenceOut_)==atLowerBound) |
---|
4870 | change=lowerChange[sequenceOut_]; |
---|
4871 | else |
---|
4872 | change=upperChange[sequenceOut_]; |
---|
4873 | if (change) { |
---|
4874 | if (sequenceOut_<numberColumns_) { |
---|
4875 | if (!rowScale_) { |
---|
4876 | for (CoinBigIndex i = columnStart[sequenceOut_]; |
---|
4877 | i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) { |
---|
4878 | rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change); |
---|
4879 | } |
---|
4880 | } else { |
---|
4881 | // apply scaling |
---|
4882 | double scale = columnScale_[sequenceOut_]; |
---|
4883 | for (CoinBigIndex i = columnStart[sequenceOut_]; |
---|
4884 | i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) { |
---|
4885 | int iRow = row[i]; |
---|
4886 | rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change); |
---|
4887 | } |
---|
4888 | } |
---|
4889 | } else { |
---|
4890 | rowArray_[5]->quickAdd(sequenceOut_-numberColumns_,-change); |
---|
4891 | } |
---|
4892 | needed++; |
---|
4893 | } |
---|
4894 | //printf("seqin %d seqout %d needed %d\n", |
---|
4895 | // sequenceIn_,sequenceOut_,needed); |
---|
4896 | if (needed) { |
---|
4897 | // ftran it |
---|
4898 | factorization_->updateColumn(rowArray_[0], rowArray_[5]); |
---|
4899 | // add |
---|
4900 | double * array5 = rowArray_[5]->denseVector(); |
---|
4901 | int * index5 = rowArray_[5]->getIndices(); |
---|
4902 | int number5 = rowArray_[5]->getNumElements(); |
---|
4903 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4904 | int lowerN=lowerActive[-1]; |
---|
4905 | int upperN=upperActive[-1]; |
---|
4906 | int nIn4=rowArray_[4]->getNumElements(); |
---|
4907 | int * index4 = rowArray_[4]->getIndices(); |
---|
4908 | #endif |
---|
4909 | for (int i = 0; i < number5; i++) { |
---|
4910 | int iPivot = index5[i]; |
---|
4911 | #ifndef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4912 | rowArray_[4]->quickAdd(iPivot,array5[iPivot]); |
---|
4913 | #else |
---|
4914 | /* later for sparse - modify here */ |
---|
4915 | int iSequence = pivotVariable_[iPivot]; |
---|
4916 | double currentSolution = solution_[iSequence]; |
---|
4917 | double currentAlpha = array[iPivot]; |
---|
4918 | double alpha5 = array5[iPivot]; |
---|
4919 | double alpha = currentAlpha+alpha5; |
---|
4920 | if (currentAlpha) { |
---|
4921 | if (alpha) { |
---|
4922 | array[iPivot] = alpha; |
---|
4923 | } else { |
---|
4924 | array[iPivot] = COIN_DBL_MIN; |
---|
4925 | } |
---|
4926 | } else { |
---|
4927 | index4[nIn4++] = iPivot; |
---|
4928 | array[iPivot] = alpha; |
---|
4929 | } |
---|
4930 | double thetaCoefficientLower = lowerChange[iSequence] + alpha; |
---|
4931 | double thetaCoefficientUpper = upperChange[iSequence] + alpha; |
---|
4932 | double oldLower = lowerCoefficient[iPivot]; |
---|
4933 | double oldUpper = upperCoefficient[iPivot]; |
---|
4934 | if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) { |
---|
4935 | double currentLower = lower_[iSequence]; |
---|
4936 | ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_); |
---|
4937 | double gap=currentSolution-currentLower; |
---|
4938 | lowerGap[iPivot]=gap; |
---|
4939 | lowerCoefficient[iPivot]=thetaCoefficientLower; |
---|
4940 | if (!oldLower) |
---|
4941 | lowerActive[lowerN++]=iPivot; |
---|
4942 | } else { |
---|
4943 | if (oldLower) |
---|
4944 | lowerCoefficient[iPivot]=COIN_DBL_MIN; |
---|
4945 | } |
---|
4946 | if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) { |
---|
4947 | double currentUpper = upper_[iSequence]; |
---|
4948 | ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_); |
---|
4949 | double gap2=-(currentSolution-currentUpper); //positive |
---|
4950 | upperGap[iPivot]=gap2; |
---|
4951 | upperCoefficient[iPivot]=-thetaCoefficientUpper; |
---|
4952 | if (!oldUpper) |
---|
4953 | upperActive[upperN++]=iPivot; |
---|
4954 | } else { |
---|
4955 | if (oldUpper) |
---|
4956 | upperCoefficient[iPivot]=COIN_DBL_MIN; |
---|
4957 | } |
---|
4958 | #endif |
---|
4959 | array5[iPivot]=0.0; |
---|
4960 | } |
---|
4961 | rowArray_[5]->setNumElements(0); |
---|
4962 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
4963 | rowArray_[4]->setNumElements(nIn4); |
---|
4964 | assert (lowerN>=0&&lowerN<=numberRows_); |
---|
4965 | lowerActive[-1]=lowerN; |
---|
4966 | upperActive[-1]=upperN; |
---|
4967 | #endif |
---|
4968 | } |
---|
4969 | } |
---|
4970 | const int * index = rowArray_[4]->getIndices(); |
---|
4971 | int number = rowArray_[4]->getNumElements(); |
---|
4972 | #define TESTXX 0 |
---|
4973 | #ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX |
---|
4974 | int * markDone = reinterpret_cast<int *>(paramData.markDone); |
---|
4975 | int nToZero=(numberRows_+numberColumns_+COIN_ANY_BITS_PER_INT-1)>>COIN_ANY_SHIFT_PER_INT; |
---|
4976 | memset(markDone,0,nToZero*sizeof(int)); |
---|
4977 | const int * backwardBasic = paramData.backwardBasic; |
---|
4978 | #endif |
---|
4979 | // first ones with alpha |
---|
4980 | double theta1=maxTheta; |
---|
4981 | int pivotRow1=-1; |
---|
4982 | #ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX |
---|
4983 | int pivotRow2=-1; |
---|
4984 | double theta2=maxTheta; |
---|
4985 | #endif |
---|
4986 | #ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX |
---|
4987 | for (int i=0;i<number;i++) { |
---|
4988 | int iPivot=index[i]; |
---|
4989 | iSequence = pivotVariable_[iPivot]; |
---|
4990 | //assert(!markDone[iSequence]); |
---|
4991 | int word = iSequence >> COIN_ANY_SHIFT_PER_INT; |
---|
4992 | int bit = iSequence & COIN_ANY_MASK_PER_INT; |
---|
4993 | markDone[word] |= ( 1 << bit ); |
---|
4994 | // solution value will be sol - theta*alpha |
---|
4995 | // bounds will be bounds + change *theta |
---|
4996 | double currentSolution = solution_[iSequence]; |
---|
4997 | double alpha = array[iPivot]; |
---|
4998 | double thetaCoefficientLower = lowerChange[iSequence] + alpha; |
---|
4999 | double thetaCoefficientUpper = upperChange[iSequence] + alpha; |
---|
5000 | if (thetaCoefficientLower > 1.0e-8) { |
---|
5001 | double currentLower = lower_[iSequence]; |
---|
5002 | ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_); |
---|
5003 | assert (currentSolution >= currentLower - 100.0*primalTolerance_); |
---|
5004 | double gap=currentSolution-currentLower; |
---|
5005 | if (thetaCoefficientLower*theta1>gap) { |
---|
5006 | theta1 = gap/thetaCoefficientLower; |
---|
5007 | //toLower=true; |
---|
5008 | pivotRow1=iPivot; |
---|
5009 | } |
---|
5010 | } |
---|
5011 | if (thetaCoefficientUpper < -1.0e-8) { |
---|
5012 | double currentUpper = upper_[iSequence]; |
---|
5013 | ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_); |
---|
5014 | assert (currentSolution <= currentUpper + 100.0*primalTolerance_); |
---|
5015 | double gap2=currentSolution-currentUpper; //negative |
---|
5016 | if (thetaCoefficientUpper*theta2<gap2) { |
---|
5017 | theta2 = gap2/thetaCoefficientUpper; |
---|
5018 | //toLower=false; |
---|
5019 | pivotRow2=iPivot; |
---|
5020 | } |
---|
5021 | } |
---|
5022 | } |
---|
5023 | // now others |
---|
5024 | int nLook=lowerList[-1]; |
---|
5025 | for (int i=0;i<nLook;i++) { |
---|
5026 | int iSequence = lowerList[i]; |
---|
5027 | int word = iSequence >> COIN_ANY_SHIFT_PER_INT; |
---|
5028 | int bit = iSequence & COIN_ANY_MASK_PER_INT; |
---|
5029 | if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) { |
---|
5030 | double currentSolution = solution_[iSequence]; |
---|
5031 | double currentLower = lower_[iSequence]; |
---|
5032 | ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_); |
---|
5033 | double thetaCoefficient = lowerChange[iSequence]; |
---|
5034 | if (thetaCoefficient > 0.0) { |
---|
5035 | double gap=currentSolution-currentLower; |
---|
5036 | if (thetaCoefficient*theta1>gap) { |
---|
5037 | theta1 = gap/thetaCoefficient; |
---|
5038 | //toLower=true; |
---|
5039 | pivotRow1 = backwardBasic[iSequence]; |
---|
5040 | } |
---|
5041 | } |
---|
5042 | } |
---|
5043 | } |
---|
5044 | nLook=upperList[-1]; |
---|
5045 | for (int i=0;i<nLook;i++) { |
---|
5046 | int iSequence = upperList[i]; |
---|
5047 | int word = iSequence >> COIN_ANY_SHIFT_PER_INT; |
---|
5048 | int bit = iSequence & COIN_ANY_MASK_PER_INT; |
---|
5049 | if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) { |
---|
5050 | double currentSolution = solution_[iSequence]; |
---|
5051 | double currentUpper = upper_[iSequence]; |
---|
5052 | ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_); |
---|
5053 | double thetaCoefficient = upperChange[iSequence]; |
---|
5054 | if (thetaCoefficient < 0) { |
---|
5055 | double gap=currentSolution-currentUpper; //negative |
---|
5056 | if (thetaCoefficient*theta2<gap) { |
---|
5057 | theta2 = gap/thetaCoefficient; |
---|
5058 | //toLower=false; |
---|
5059 | pivotRow2 = backwardBasic[iSequence]; |
---|
5060 | } |
---|
5061 | } |
---|
5062 | } |
---|
5063 | } |
---|
5064 | if (theta2<theta1) { |
---|
5065 | theta_=theta2; |
---|
5066 | toLower=false; |
---|
5067 | pivotRow_=pivotRow2; |
---|
5068 | } else { |
---|
5069 | theta_=theta1; |
---|
5070 | toLower=true; |
---|
5071 | pivotRow_=pivotRow1; |
---|
5072 | } |
---|
5073 | #if 0 //TESTXX |
---|
5074 | #ifdef CLP_PARAMETRIC_DENSE_ARRAYS |
---|
5075 | { |
---|
5076 | double * checkArray = new double[numberRows_]; |
---|
5077 | memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double)); |
---|
5078 | int lowerN=lowerActive[-1]; |
---|
5079 | for (int i=0;i<lowerN;i++) { |
---|
5080 | int iRow=lowerActive[i]; |
---|
5081 | int iSequence = pivotVariable_[iRow]; |
---|
5082 | double alpha = array[iRow]; |
---|
5083 | double thetaCoefficient = lowerChange[iSequence] + alpha; |
---|
5084 | if (thetaCoefficient > 1.0e-8&&lower_[iSequence]>-1.0e30) { |
---|
5085 | assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5); |
---|
5086 | if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) { |
---|
5087 | abort(); |
---|
5088 | } |
---|
5089 | } else { |
---|
5090 | assert (fabs(checkArray[iRow])<1.0e-12); |
---|
5091 | if (fabs(checkArray[iRow])>1.0e-12) { |
---|
5092 | abort(); |
---|
5093 | } |
---|
5094 | } |
---|
5095 | checkArray[iRow]=0.0; |
---|
5096 | } |
---|
5097 | for (int i=0;i<numberRows_;i++) { |
---|
5098 | assert (!checkArray[i]); |
---|
5099 | if (checkArray[i]) |
---|
5100 | abort(); |
---|
5101 | } |
---|
5102 | memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double)); |
---|
5103 | int upperN=upperActive[-1]; |
---|
5104 | for (int i=0;i<upperN;i++) { |
---|
5105 | int iRow=upperActive[i]; |
---|
5106 | int iSequence = pivotVariable_[iRow]; |
---|
5107 | double alpha = array[iRow]; |
---|
5108 | double thetaCoefficient = -(upperChange[iSequence] + alpha); |
---|
5109 | if (thetaCoefficient > 1.0e-8&&upper_[iSequence]<1.0e30) { |
---|
5110 | assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5); |
---|
5111 | if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) { |
---|
5112 | abort(); |
---|
5113 | } |
---|
5114 | } else { |
---|
5115 | assert (fabs(checkArray[iRow])<1.0e-12); |
---|
5116 | if (fabs(checkArray[iRow])>1.0e-12) { |
---|
5117 | abort(); |
---|
5118 | } |
---|
5119 | } |
---|
5120 | checkArray[iRow]=0.0; |
---|
5121 | } |
---|
5122 | for (int i=0;i<numberRows_;i++) { |
---|
5123 | assert (!checkArray[i]); |
---|
5124 | if (checkArray[i]) |
---|
5125 | abort(); |
---|
5126 | } |
---|
5127 | delete [] checkArray; |
---|
5128 | } |
---|
5129 | double theta3=maxTheta; |
---|
5130 | int pivotRow3=-1; |
---|
5131 | int lowerN=lowerActive[-1]; |
---|
5132 | for (int i=0;i<lowerN;i++) { |
---|
5133 | int iRow=lowerActive[i]; |
---|
5134 | double lowerC = lowerCoefficient[iRow]; |
---|
5135 | double gap=lowerGap[iRow]; |
---|
5136 | if (toLower&&iRow==pivotRow_) { |
---|
5137 | assert (lowerC*theta3>gap-1.0e-8); |
---|
5138 | if (lowerC*theta3<gap-1.0e-8) |
---|
5139 | abort(); |
---|
5140 | } |
---|
5141 | if (lowerC*theta3>gap&&lowerC!=COIN_DBL_MIN) { |
---|
5142 | theta3 = gap/lowerC; |
---|
5143 | pivotRow3=iRow; |
---|
5144 | } |
---|
5145 | } |
---|
5146 | int pivotRow4=pivotRow3; |
---|
5147 | double theta4=theta3; |
---|
5148 | int upperN=upperActive[-1]; |
---|
5149 | for (int i=0;i<upperN;i++) { |
---|
5150 | int iRow=upperActive[i]; |
---|
5151 | double upperC = upperCoefficient[iRow]; |
---|
5152 | double gap=upperGap[iRow]; |
---|
5153 | if (!toLower&&iRow==pivotRow_) { |
---|
5154 | assert (upperC*theta3>gap-1.0e-8); |
---|
5155 | if (upperC*theta3<gap-1.0e-8) |
---|
5156 | abort(); |
---|
5157 | } |
---|
5158 | if (upperC*theta4>gap&&upperC!=COIN_DBL_MIN) { |
---|
5159 | theta4 = gap/upperC; |
---|
5160 | pivotRow4=iRow; |
---|
5161 | } |
---|
5162 | } |
---|
5163 | bool toLower3; |
---|
5164 | if (theta4<theta3) { |
---|
5165 | theta3=theta4; |
---|
5166 | toLower3=false; |
---|
5167 | pivotRow3=pivotRow4; |
---|
5168 | } else { |
---|
5169 | toLower3=true; |
---|
5170 | } |
---|
5171 | if (fabs(theta3-theta_)>1.0e-8) |
---|
5172 | abort(); |
---|
5173 | if (toLower!=toLower3||pivotRow_!=pivotRow3) { |
---|
5174 | printf("bad piv - good %d %g %s, bad %d %g %s\n",pivotRow_,theta_,toLower ? "toLower" : "toUpper", |
---|
5175 | pivotRow3,theta3,toLower3 ? "toLower" : "toUpper"); |
---|
5176 | //zzzzzz++; |
---|
5177 | if (true/*zzzzzz>zzzzzzOther*/) { |
---|
5178 | printf("Swapping\n"); |
---|
5179 | pivotRow_=pivotRow3; |
---|
5180 | theta_=theta3; |
---|
5181 | toLower=toLower3; |
---|
5182 | } |
---|
5183 | } |
---|
5184 | #endif |
---|
5185 | #endif |
---|
5186 | #else |
---|
5187 | #if 0 //CLP_PARAMETRIC_DENSE_ARRAYS==2 |
---|
5188 | { |
---|
5189 | double * checkArray = new double[numberRows_]; |
---|
5190 | memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double)); |
---|
5191 | int lowerN=lowerActive[-1]; |
---|
5192 | for (int i=0;i<lowerN;i++) { |
---|
5193 | int iRow=lowerActive[i]; |
---|
5194 | checkArray[iRow]=0.0; |
---|
5195 | } |
---|
5196 | for (int i=0;i<numberRows_;i++) { |
---|
5197 | assert (!checkArray[i]); |
---|
5198 | if (checkArray[i]) |
---|
5199 | abort(); |
---|
5200 | } |
---|
5201 | memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double)); |
---|
5202 | int upperN=upperActive[-1]; |
---|
5203 | for (int i=0;i<upperN;i++) { |
---|
5204 | int iRow=upperActive[i]; |
---|
5205 | checkArray[iRow]=0.0; |
---|
5206 | } |
---|
5207 | for (int i=0;i<numberRows_;i++) { |
---|
5208 | assert (!checkArray[i]); |
---|
5209 | if (checkArray[i]) |
---|
5210 | abort(); |
---|
5211 | } |
---|
5212 | delete [] checkArray; |
---|
5213 | } |
---|
5214 | #endif |
---|
5215 | int lowerN=lowerActive[-1]; |
---|
5216 | for (int i=0;i<lowerN;i++) { |
---|
5217 | int iRow=lowerActive[i]; |
---|
5218 | double lowerC = lowerCoefficient[iRow]; |
---|
5219 | double gap=lowerGap[iRow]; |
---|
5220 | if (lowerC*theta1>gap&&lowerC!=COIN_DBL_MIN) { |
---|
5221 | theta1 = gap/lowerC; |
---|
5222 | pivotRow1=iRow; |
---|
5223 | } |
---|
5224 | } |
---|
5225 | pivotRow_=pivotRow1; |
---|
5226 | theta_=theta1; |
---|
5227 | int upperN=upperActive[-1]; |
---|
5228 | for (int i=0;i<upperN;i++) { |
---|
5229 | int iRow=upperActive |
---|