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