1 | /* $Id: ClpSimplexDual.cpp 1785 2011-08-25 10:17:49Z 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 | #if CAN_HAVE_ZERO_OBJ>1 |
---|
1226 | if ((specialOptions_&2097152)!=0) |
---|
1227 | theta_=0.0; |
---|
1228 | #endif |
---|
1229 | } else { |
---|
1230 | // Make sure direction plausible |
---|
1231 | CoinAssert (upperOut_ < 1.0e50 || lowerOut_ > -1.0e50); |
---|
1232 | // If in integer cleanup do direction using duals |
---|
1233 | // may be wrong way round |
---|
1234 | if(ifValuesPass == 2) { |
---|
1235 | if (dual_[pivotRow_] > 0.0) { |
---|
1236 | // this will give a -1 in pivot row (as slacks are -1.0) |
---|
1237 | directionOut_ = 1; |
---|
1238 | } else { |
---|
1239 | directionOut_ = -1; |
---|
1240 | } |
---|
1241 | } |
---|
1242 | if (directionOut_ < 0 && fabs(valueOut_ - upperOut_) > dualBound_ + primalTolerance_) { |
---|
1243 | if (fabs(valueOut_ - upperOut_) > fabs(valueOut_ - lowerOut_)) |
---|
1244 | directionOut_ = 1; |
---|
1245 | } else if (directionOut_ > 0 && fabs(valueOut_ - lowerOut_) > dualBound_ + primalTolerance_) { |
---|
1246 | if (fabs(valueOut_ - upperOut_) < fabs(valueOut_ - lowerOut_)) |
---|
1247 | directionOut_ = -1; |
---|
1248 | } |
---|
1249 | double direction = directionOut_; |
---|
1250 | rowArray_[0]->createPacked(1, &pivotRow_, &direction); |
---|
1251 | factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]); |
---|
1252 | // put row of tableau in rowArray[0] and columnArray[0] |
---|
1253 | if (!scaledMatrix_) { |
---|
1254 | matrix_->transposeTimes(this, -1.0, |
---|
1255 | rowArray_[0], rowArray_[3], columnArray_[0]); |
---|
1256 | } else { |
---|
1257 | double * saveR = rowScale_; |
---|
1258 | double * saveC = columnScale_; |
---|
1259 | rowScale_ = NULL; |
---|
1260 | columnScale_ = NULL; |
---|
1261 | scaledMatrix_->transposeTimes(this, -1.0, |
---|
1262 | rowArray_[0], rowArray_[3], columnArray_[0]); |
---|
1263 | rowScale_ = saveR; |
---|
1264 | columnScale_ = saveC; |
---|
1265 | } |
---|
1266 | acceptablePivot *= 10.0; |
---|
1267 | // do ratio test |
---|
1268 | if (ifValuesPass == 1) { |
---|
1269 | checkPossibleValuesMove(rowArray_[0], columnArray_[0], |
---|
1270 | acceptablePivot); |
---|
1271 | } else { |
---|
1272 | checkPossibleCleanup(rowArray_[0], columnArray_[0], |
---|
1273 | acceptablePivot); |
---|
1274 | if (sequenceIn_ < 0) { |
---|
1275 | rowArray_[0]->clear(); |
---|
1276 | columnArray_[0]->clear(); |
---|
1277 | continue; // can't do anything |
---|
1278 | } |
---|
1279 | } |
---|
1280 | |
---|
1281 | // recompute true dualOut_ |
---|
1282 | if (directionOut_ < 0) { |
---|
1283 | dualOut_ = valueOut_ - upperOut_; |
---|
1284 | } else { |
---|
1285 | dualOut_ = lowerOut_ - valueOut_; |
---|
1286 | } |
---|
1287 | // check what happened if was values pass |
---|
1288 | // may want to move part way i.e. movement |
---|
1289 | bool normalIteration = (sequenceIn_ != sequenceOut_); |
---|
1290 | |
---|
1291 | clearPivoted(sequenceOut_); // make sure won't be done again |
---|
1292 | // see if end of values pass |
---|
1293 | if (!numberCandidates) { |
---|
1294 | int iRow; |
---|
1295 | delete [] candidateList; |
---|
1296 | delete [] givenDuals; |
---|
1297 | candidate = -2; // -2 signals end |
---|
1298 | givenDuals = NULL; |
---|
1299 | candidateList = NULL; |
---|
1300 | ifValuesPass = 1; |
---|
1301 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1302 | int iPivot = pivotVariable_[iRow]; |
---|
1303 | //assert (fabs(dj_[iPivot]),1.0e-5); |
---|
1304 | clearPivoted(iPivot); |
---|
1305 | } |
---|
1306 | } |
---|
1307 | if (!normalIteration) { |
---|
1308 | //rowArray_[0]->cleanAndPackSafe(1.0e-60); |
---|
1309 | //columnArray_[0]->cleanAndPackSafe(1.0e-60); |
---|
1310 | updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_); |
---|
1311 | if (candidate == -2) |
---|
1312 | problemStatus_ = -2; |
---|
1313 | continue; // skip rest of iteration |
---|
1314 | } else { |
---|
1315 | // recompute dualOut_ |
---|
1316 | if (directionOut_ < 0) { |
---|
1317 | dualOut_ = valueOut_ - upperOut_; |
---|
1318 | } else { |
---|
1319 | dualOut_ = lowerOut_ - valueOut_; |
---|
1320 | } |
---|
1321 | } |
---|
1322 | } |
---|
1323 | if (sequenceIn_ >= 0) { |
---|
1324 | // normal iteration |
---|
1325 | // update the incoming column |
---|
1326 | double btranAlpha = -alpha_ * directionOut_; // for check |
---|
1327 | unpackPacked(rowArray_[1]); |
---|
1328 | // moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]); |
---|
1329 | // and update dual weights (can do in parallel - with extra array) |
---|
1330 | alpha_ = dualRowPivot_->updateWeights(rowArray_[0], |
---|
1331 | rowArray_[2], |
---|
1332 | rowArray_[3], |
---|
1333 | rowArray_[1]); |
---|
1334 | // see if update stable |
---|
1335 | #ifdef CLP_DEBUG |
---|
1336 | if ((handler_->logLevel() & 32)) |
---|
1337 | printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_); |
---|
1338 | #endif |
---|
1339 | double checkValue = 1.0e-7; |
---|
1340 | // if can't trust much and long way from optimal then relax |
---|
1341 | if (largestPrimalError_ > 10.0) |
---|
1342 | checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_); |
---|
1343 | if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || |
---|
1344 | fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) { |
---|
1345 | handler_->message(CLP_DUAL_CHECK, messages_) |
---|
1346 | << btranAlpha |
---|
1347 | << alpha_ |
---|
1348 | << CoinMessageEol; |
---|
1349 | if (factorization_->pivots()) { |
---|
1350 | dualRowPivot_->unrollWeights(); |
---|
1351 | problemStatus_ = -2; // factorize now |
---|
1352 | rowArray_[0]->clear(); |
---|
1353 | rowArray_[1]->clear(); |
---|
1354 | columnArray_[0]->clear(); |
---|
1355 | returnCode = -2; |
---|
1356 | break; |
---|
1357 | } else { |
---|
1358 | // take on more relaxed criterion |
---|
1359 | double test; |
---|
1360 | if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8) |
---|
1361 | test = 1.0e-1 * fabs(alpha_); |
---|
1362 | else |
---|
1363 | test = 1.0e-4 * (1.0 + fabs(alpha_)); |
---|
1364 | if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || |
---|
1365 | fabs(btranAlpha - alpha_) > test) { |
---|
1366 | dualRowPivot_->unrollWeights(); |
---|
1367 | // need to reject something |
---|
1368 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
1369 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
1370 | << x << sequenceWithin(sequenceOut_) |
---|
1371 | << CoinMessageEol; |
---|
1372 | #ifdef COIN_DEVELOP |
---|
1373 | printf("flag a %g %g\n", btranAlpha, alpha_); |
---|
1374 | #endif |
---|
1375 | //#define FEB_TRY |
---|
1376 | #if 1 //def FEB_TRY |
---|
1377 | // Make safer? |
---|
1378 | factorization_->saferTolerances (-0.99, -1.03); |
---|
1379 | #endif |
---|
1380 | setFlagged(sequenceOut_); |
---|
1381 | progress_.clearBadTimes(); |
---|
1382 | lastBadIteration_ = numberIterations_; // say be more cautious |
---|
1383 | rowArray_[0]->clear(); |
---|
1384 | rowArray_[1]->clear(); |
---|
1385 | columnArray_[0]->clear(); |
---|
1386 | if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) { |
---|
1387 | //printf("I think should declare infeasible\n"); |
---|
1388 | problemStatus_ = 1; |
---|
1389 | returnCode = 1; |
---|
1390 | break; |
---|
1391 | } |
---|
1392 | continue; |
---|
1393 | } |
---|
1394 | } |
---|
1395 | } |
---|
1396 | // update duals BEFORE replaceColumn so can do updateColumn |
---|
1397 | double objectiveChange = 0.0; |
---|
1398 | // do duals first as variables may flip bounds |
---|
1399 | // rowArray_[0] and columnArray_[0] may have flips |
---|
1400 | // so use rowArray_[3] for work array from here on |
---|
1401 | int nswapped = 0; |
---|
1402 | //rowArray_[0]->cleanAndPackSafe(1.0e-60); |
---|
1403 | //columnArray_[0]->cleanAndPackSafe(1.0e-60); |
---|
1404 | if (candidate == -1) { |
---|
1405 | #if CLP_CAN_HAVE_ZERO_OBJ>1 |
---|
1406 | if ((specialOptions_&2097152)==0) { |
---|
1407 | #endif |
---|
1408 | // make sure incoming doesn't count |
---|
1409 | Status saveStatus = getStatus(sequenceIn_); |
---|
1410 | setStatus(sequenceIn_, basic); |
---|
1411 | nswapped = updateDualsInDual(rowArray_[0], columnArray_[0], |
---|
1412 | rowArray_[2], theta_, |
---|
1413 | objectiveChange, false); |
---|
1414 | setStatus(sequenceIn_, saveStatus); |
---|
1415 | #if CLP_CAN_HAVE_ZERO_OBJ>1 |
---|
1416 | } else { |
---|
1417 | rowArray_[0]->clear(); |
---|
1418 | rowArray_[2]->clear(); |
---|
1419 | columnArray_[0]->clear(); |
---|
1420 | } |
---|
1421 | #endif |
---|
1422 | } else { |
---|
1423 | updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_); |
---|
1424 | } |
---|
1425 | double oldDualOut = dualOut_; |
---|
1426 | // which will change basic solution |
---|
1427 | if (nswapped) { |
---|
1428 | if (rowArray_[2]->getNumElements()) { |
---|
1429 | factorization_->updateColumn(rowArray_[3], rowArray_[2]); |
---|
1430 | dualRowPivot_->updatePrimalSolution(rowArray_[2], |
---|
1431 | 1.0, objectiveChange); |
---|
1432 | } |
---|
1433 | // recompute dualOut_ |
---|
1434 | valueOut_ = solution_[sequenceOut_]; |
---|
1435 | if (directionOut_ < 0) { |
---|
1436 | dualOut_ = valueOut_ - upperOut_; |
---|
1437 | } else { |
---|
1438 | dualOut_ = lowerOut_ - valueOut_; |
---|
1439 | } |
---|
1440 | #if 0 |
---|
1441 | if (dualOut_ < 0.0) { |
---|
1442 | #ifdef CLP_DEBUG |
---|
1443 | if (handler_->logLevel() & 32) { |
---|
1444 | printf(" dualOut_ %g %g save %g\n", dualOut_, averagePrimalInfeasibility, saveDualOut); |
---|
1445 | printf("values %g %g %g %g %g %g %g\n", lowerOut_, valueOut_, upperOut_, |
---|
1446 | objectiveChange,); |
---|
1447 | } |
---|
1448 | #endif |
---|
1449 | if (upperOut_ == lowerOut_) |
---|
1450 | dualOut_ = 0.0; |
---|
1451 | } |
---|
1452 | if(dualOut_ < -CoinMax(1.0e-12 * averagePrimalInfeasibility, 1.0e-8) |
---|
1453 | && factorization_->pivots() > 100 && |
---|
1454 | getStatus(sequenceIn_) != isFree) { |
---|
1455 | // going backwards - factorize |
---|
1456 | dualRowPivot_->unrollWeights(); |
---|
1457 | problemStatus_ = -2; // factorize now |
---|
1458 | returnCode = -2; |
---|
1459 | break; |
---|
1460 | } |
---|
1461 | #endif |
---|
1462 | } |
---|
1463 | // amount primal will move |
---|
1464 | double movement = -dualOut_ * directionOut_ / alpha_; |
---|
1465 | double movementOld = oldDualOut * directionOut_ / alpha_; |
---|
1466 | // so objective should increase by fabs(dj)*movement |
---|
1467 | // but we already have objective change - so check will be good |
---|
1468 | if (objectiveChange + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(objectiveValue_))) { |
---|
1469 | #ifdef CLP_DEBUG |
---|
1470 | if (handler_->logLevel() & 32) |
---|
1471 | printf("movement %g, swap change %g, rest %g * %g\n", |
---|
1472 | objectiveChange + fabs(movement * dualIn_), |
---|
1473 | objectiveChange, movement, dualIn_); |
---|
1474 | #endif |
---|
1475 | if(factorization_->pivots()) { |
---|
1476 | // going backwards - factorize |
---|
1477 | dualRowPivot_->unrollWeights(); |
---|
1478 | problemStatus_ = -2; // factorize now |
---|
1479 | returnCode = -2; |
---|
1480 | break; |
---|
1481 | } |
---|
1482 | } |
---|
1483 | // if stable replace in basis |
---|
1484 | int updateStatus = factorization_->replaceColumn(this, |
---|
1485 | rowArray_[2], |
---|
1486 | rowArray_[1], |
---|
1487 | pivotRow_, |
---|
1488 | alpha_, |
---|
1489 | (moreSpecialOptions_ & 16) != 0, |
---|
1490 | acceptablePivot); |
---|
1491 | // If looks like bad pivot - refactorize |
---|
1492 | if (fabs(dualOut_) > 1.0e50) |
---|
1493 | updateStatus = 2; |
---|
1494 | // if no pivots, bad update but reasonable alpha - take and invert |
---|
1495 | if (updateStatus == 2 && |
---|
1496 | !factorization_->pivots() && fabs(alpha_) > 1.0e-5) |
---|
1497 | updateStatus = 4; |
---|
1498 | if (updateStatus == 1 || updateStatus == 4) { |
---|
1499 | // slight error |
---|
1500 | if (factorization_->pivots() > 5 || updateStatus == 4) { |
---|
1501 | problemStatus_ = -2; // factorize now |
---|
1502 | returnCode = -3; |
---|
1503 | } |
---|
1504 | } else if (updateStatus == 2) { |
---|
1505 | // major error |
---|
1506 | dualRowPivot_->unrollWeights(); |
---|
1507 | // later we may need to unwind more e.g. fake bounds |
---|
1508 | if (factorization_->pivots() && |
---|
1509 | ((moreSpecialOptions_ & 16) == 0 || factorization_->pivots() > 4)) { |
---|
1510 | problemStatus_ = -2; // factorize now |
---|
1511 | returnCode = -2; |
---|
1512 | moreSpecialOptions_ |= 16; |
---|
1513 | break; |
---|
1514 | } else { |
---|
1515 | // need to reject something |
---|
1516 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
1517 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
1518 | << x << sequenceWithin(sequenceOut_) |
---|
1519 | << CoinMessageEol; |
---|
1520 | #ifdef COIN_DEVELOP |
---|
1521 | printf("flag b %g\n", alpha_); |
---|
1522 | #endif |
---|
1523 | setFlagged(sequenceOut_); |
---|
1524 | progress_.clearBadTimes(); |
---|
1525 | lastBadIteration_ = numberIterations_; // say be more cautious |
---|
1526 | rowArray_[0]->clear(); |
---|
1527 | rowArray_[1]->clear(); |
---|
1528 | columnArray_[0]->clear(); |
---|
1529 | // make sure dual feasible |
---|
1530 | // look at all rows and columns |
---|
1531 | double objectiveChange = 0.0; |
---|
1532 | updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1], |
---|
1533 | 0.0, objectiveChange, true); |
---|
1534 | rowArray_[1]->clear(); |
---|
1535 | columnArray_[0]->clear(); |
---|
1536 | continue; |
---|
1537 | } |
---|
1538 | } else if (updateStatus == 3) { |
---|
1539 | // out of memory |
---|
1540 | // increase space if not many iterations |
---|
1541 | if (factorization_->pivots() < |
---|
1542 | 0.5 * factorization_->maximumPivots() && |
---|
1543 | factorization_->pivots() < 200) |
---|
1544 | factorization_->areaFactor( |
---|
1545 | factorization_->areaFactor() * 1.1); |
---|
1546 | problemStatus_ = -2; // factorize now |
---|
1547 | } else if (updateStatus == 5) { |
---|
1548 | problemStatus_ = -2; // factorize now |
---|
1549 | } |
---|
1550 | // update primal solution |
---|
1551 | if (theta_ < 0.0 && candidate == -1) { |
---|
1552 | #ifdef CLP_DEBUG |
---|
1553 | if (handler_->logLevel() & 32) |
---|
1554 | printf("negative theta %g\n", theta_); |
---|
1555 | #endif |
---|
1556 | theta_ = 0.0; |
---|
1557 | } |
---|
1558 | // do actual flips |
---|
1559 | flipBounds(rowArray_[0], columnArray_[0]); |
---|
1560 | //rowArray_[1]->expand(); |
---|
1561 | dualRowPivot_->updatePrimalSolution(rowArray_[1], |
---|
1562 | movement, |
---|
1563 | objectiveChange); |
---|
1564 | #ifdef CLP_DEBUG |
---|
1565 | double oldobj = objectiveValue_; |
---|
1566 | #endif |
---|
1567 | // modify dualout |
---|
1568 | dualOut_ /= alpha_; |
---|
1569 | dualOut_ *= -directionOut_; |
---|
1570 | //setStatus(sequenceIn_,basic); |
---|
1571 | dj_[sequenceIn_] = 0.0; |
---|
1572 | double oldValue = valueIn_; |
---|
1573 | if (directionIn_ == -1) { |
---|
1574 | // as if from upper bound |
---|
1575 | valueIn_ = upperIn_ + dualOut_; |
---|
1576 | } else { |
---|
1577 | // as if from lower bound |
---|
1578 | valueIn_ = lowerIn_ + dualOut_; |
---|
1579 | } |
---|
1580 | objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue); |
---|
1581 | // outgoing |
---|
1582 | // set dj to zero unless values pass |
---|
1583 | if (directionOut_ > 0) { |
---|
1584 | valueOut_ = lowerOut_; |
---|
1585 | if (candidate == -1) |
---|
1586 | dj_[sequenceOut_] = theta_; |
---|
1587 | } else { |
---|
1588 | valueOut_ = upperOut_; |
---|
1589 | if (candidate == -1) |
---|
1590 | dj_[sequenceOut_] = -theta_; |
---|
1591 | } |
---|
1592 | solution_[sequenceOut_] = valueOut_; |
---|
1593 | int whatNext = housekeeping(objectiveChange); |
---|
1594 | #ifdef CLP_REPORT_PROGRESS |
---|
1595 | if (ixxxxxx > ixxyyyy - 5) { |
---|
1596 | handler_->setLogLevel(63); |
---|
1597 | int nTotal = numberColumns_ + numberRows_; |
---|
1598 | double oldObj = 0.0; |
---|
1599 | double newObj = 0.0; |
---|
1600 | for (int i = 0; i < nTotal; i++) { |
---|
1601 | if (savePSol[i]) |
---|
1602 | oldObj += savePSol[i] * saveCost[i]; |
---|
1603 | if (solution_[i]) |
---|
1604 | newObj += solution_[i] * cost_[i]; |
---|
1605 | bool printIt = false; |
---|
1606 | if (cost_[i] != saveCost[i]) |
---|
1607 | printIt = true; |
---|
1608 | if (status_[i] != saveStat[i]) |
---|
1609 | printIt = true; |
---|
1610 | if (printIt) |
---|
1611 | printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n", |
---|
1612 | i, saveStat[i], saveCost[i], savePSol[i], |
---|
1613 | status_[i], cost_[i], solution_[i]); |
---|
1614 | // difference |
---|
1615 | savePSol[i] = solution_[i] - savePSol[i]; |
---|
1616 | } |
---|
1617 | printf("pivots %d, old obj %g new %g\n", |
---|
1618 | factorization_->pivots(), |
---|
1619 | oldObj, newObj); |
---|
1620 | memset(saveDj, 0, numberRows_ * sizeof(double)); |
---|
1621 | times(1.0, savePSol, saveDj); |
---|
1622 | double largest = 1.0e-6; |
---|
1623 | int k = -1; |
---|
1624 | for (int i = 0; i < numberRows_; i++) { |
---|
1625 | saveDj[i] -= savePSol[i+numberColumns_]; |
---|
1626 | if (fabs(saveDj[i]) > largest) { |
---|
1627 | largest = fabs(saveDj[i]); |
---|
1628 | k = i; |
---|
1629 | } |
---|
1630 | } |
---|
1631 | if (k >= 0) |
---|
1632 | printf("Not null %d %g\n", k, largest); |
---|
1633 | } |
---|
1634 | #endif |
---|
1635 | #ifdef VUB |
---|
1636 | { |
---|
1637 | if ((sequenceIn_ < numberColumns_ && vub[sequenceIn_] >= 0) || toVub[sequenceIn_] >= 0 || |
---|
1638 | (sequenceOut_ < numberColumns_ && vub[sequenceOut_] >= 0) || toVub[sequenceOut_] >= 0) { |
---|
1639 | int inSequence = sequenceIn_; |
---|
1640 | int inVub = -1; |
---|
1641 | if (sequenceIn_ < numberColumns_) |
---|
1642 | inVub = vub[sequenceIn_]; |
---|
1643 | int inBack = toVub[inSequence]; |
---|
1644 | int inSlack = -1; |
---|
1645 | if (inSequence >= numberColumns_ && inBack >= 0) { |
---|
1646 | inSlack = inSequence - numberColumns_; |
---|
1647 | inSequence = inBack; |
---|
1648 | inBack = toVub[inSequence]; |
---|
1649 | } |
---|
1650 | if (inVub >= 0) |
---|
1651 | printf("Vub %d in ", inSequence); |
---|
1652 | if (inBack >= 0 && inSlack < 0) |
---|
1653 | printf("%d (descendent of %d) in ", inSequence, inBack); |
---|
1654 | if (inSlack >= 0) |
---|
1655 | printf("slack for row %d -> %d (descendent of %d) in ", inSlack, inSequence, inBack); |
---|
1656 | int outSequence = sequenceOut_; |
---|
1657 | int outVub = -1; |
---|
1658 | if (sequenceOut_ < numberColumns_) |
---|
1659 | outVub = vub[sequenceOut_]; |
---|
1660 | int outBack = toVub[outSequence]; |
---|
1661 | int outSlack = -1; |
---|
1662 | if (outSequence >= numberColumns_ && outBack >= 0) { |
---|
1663 | outSlack = outSequence - numberColumns_; |
---|
1664 | outSequence = outBack; |
---|
1665 | outBack = toVub[outSequence]; |
---|
1666 | } |
---|
1667 | if (outVub >= 0) |
---|
1668 | printf("Vub %d out ", outSequence); |
---|
1669 | if (outBack >= 0 && outSlack < 0) |
---|
1670 | printf("%d (descendent of %d) out ", outSequence, outBack); |
---|
1671 | if (outSlack >= 0) |
---|
1672 | printf("slack for row %d -> %d (descendent of %d) out ", outSlack, outSequence, outBack); |
---|
1673 | printf("\n"); |
---|
1674 | } |
---|
1675 | } |
---|
1676 | #endif |
---|
1677 | #if 0 |
---|
1678 | if (numberIterations_ > 206033) |
---|
1679 | handler_->setLogLevel(63); |
---|
1680 | if (numberIterations_ > 210567) |
---|
1681 | exit(77); |
---|
1682 | #endif |
---|
1683 | if (!givenDuals && ifValuesPass && ifValuesPass != 2) { |
---|
1684 | handler_->message(CLP_END_VALUES_PASS, messages_) |
---|
1685 | << numberIterations_; |
---|
1686 | whatNext = 1; |
---|
1687 | } |
---|
1688 | #ifdef CHECK_ACCURACY |
---|
1689 | if (whatNext) { |
---|
1690 | CoinMemcpyN(solution_, (numberRows_ + numberColumns_), zzzzzz); |
---|
1691 | } |
---|
1692 | #endif |
---|
1693 | //if (numberIterations_==1890) |
---|
1694 | //whatNext=1; |
---|
1695 | //if (numberIterations_>2000) |
---|
1696 | //exit(77); |
---|
1697 | // and set bounds correctly |
---|
1698 | originalBound(sequenceIn_); |
---|
1699 | changeBound(sequenceOut_); |
---|
1700 | #ifdef CLP_DEBUG |
---|
1701 | if (objectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16)) |
---|
1702 | printf("obj backwards %g %g\n", objectiveValue_, oldobj); |
---|
1703 | #endif |
---|
1704 | #if 0 |
---|
1705 | { |
---|
1706 | for (int i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
1707 | FakeBound bound = getFakeBound(i); |
---|
1708 | if (bound == ClpSimplexDual::upperFake) { |
---|
1709 | assert (upper_[i] < 1.0e20); |
---|
1710 | } else if (bound == ClpSimplexDual::lowerFake) { |
---|
1711 | assert (lower_[i] > -1.0e20); |
---|
1712 | } else if (bound == ClpSimplexDual::bothFake) { |
---|
1713 | assert (upper_[i] < 1.0e20); |
---|
1714 | assert (lower_[i] > -1.0e20); |
---|
1715 | } |
---|
1716 | } |
---|
1717 | } |
---|
1718 | #endif |
---|
1719 | if (whatNext == 1 || candidate == -2) { |
---|
1720 | problemStatus_ = -2; // refactorize |
---|
1721 | } else if (whatNext == 2) { |
---|
1722 | // maximum iterations or equivalent |
---|
1723 | problemStatus_ = 3; |
---|
1724 | returnCode = 3; |
---|
1725 | break; |
---|
1726 | } |
---|
1727 | // Check event |
---|
1728 | { |
---|
1729 | int status = eventHandler_->event(ClpEventHandler::endOfIteration); |
---|
1730 | if (status >= 0) { |
---|
1731 | problemStatus_ = 5; |
---|
1732 | secondaryStatus_ = ClpEventHandler::endOfIteration; |
---|
1733 | returnCode = 4; |
---|
1734 | break; |
---|
1735 | } |
---|
1736 | } |
---|
1737 | } else { |
---|
1738 | #ifdef CLP_INVESTIGATE_SERIAL |
---|
1739 | z_thinks = 1; |
---|
1740 | #endif |
---|
1741 | // no incoming column is valid |
---|
1742 | pivotRow_ = -1; |
---|
1743 | #ifdef CLP_DEBUG |
---|
1744 | if (handler_->logLevel() & 32) |
---|
1745 | printf("** no column pivot\n"); |
---|
1746 | #endif |
---|
1747 | if (factorization_->pivots() < 2 && acceptablePivot_ <= 1.0e-8) { |
---|
1748 | //&&goodAccuracy()) { |
---|
1749 | // If not in branch and bound etc save ray |
---|
1750 | delete [] ray_; |
---|
1751 | if ((specialOptions_&(1024 | 4096)) == 0 || (specialOptions_ & 32) != 0) { |
---|
1752 | // create ray anyway |
---|
1753 | ray_ = new double [ numberRows_]; |
---|
1754 | rowArray_[0]->expand(); // in case packed |
---|
1755 | CoinMemcpyN(rowArray_[0]->denseVector(), numberRows_, ray_); |
---|
1756 | } else { |
---|
1757 | ray_ = NULL; |
---|
1758 | } |
---|
1759 | // If we have just factorized and infeasibility reasonable say infeas |
---|
1760 | double dualTest = ((specialOptions_ & 4096) != 0) ? 1.0e8 : 1.0e13; |
---|
1761 | if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > dualTest) { |
---|
1762 | double testValue = 1.0e-4; |
---|
1763 | if (!factorization_->pivots() && numberPrimalInfeasibilities_ == 1) |
---|
1764 | testValue = 1.0e-6; |
---|
1765 | if (valueOut_ > upperOut_ + testValue || valueOut_ < lowerOut_ - testValue |
---|
1766 | || (specialOptions_ & 64) == 0) { |
---|
1767 | // say infeasible |
---|
1768 | problemStatus_ = 1; |
---|
1769 | // unless primal feasible!!!! |
---|
1770 | //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_, |
---|
1771 | // numberDualInfeasibilities_,sumDualInfeasibilities_); |
---|
1772 | //#define TEST_CLP_NODE |
---|
1773 | #ifndef TEST_CLP_NODE |
---|
1774 | // Should be correct - but ... |
---|
1775 | int numberFake = numberAtFakeBound(); |
---|
1776 | double sumPrimal = (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_; |
---|
1777 | if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-5 || |
---|
1778 | (sumPrimal < 1.0e5 && (specialOptions_ & 1024) != 0 && factorization_->pivots())) { |
---|
1779 | if (sumPrimal > 50.0 && factorization_->pivots() > 2) { |
---|
1780 | problemStatus_ = -4; |
---|
1781 | #ifdef COIN_DEVELOP |
---|
1782 | printf("status to -4 at %d - primalinf %g pivots %d\n", |
---|
1783 | __LINE__, sumPrimalInfeasibilities_, |
---|
1784 | factorization_->pivots()); |
---|
1785 | #endif |
---|
1786 | } else { |
---|
1787 | problemStatus_ = 10; |
---|
1788 | #if COIN_DEVELOP>1 |
---|
1789 | 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", |
---|
1790 | __LINE__, numberPrimalInfeasibilities_, |
---|
1791 | sumPrimalInfeasibilities_, |
---|
1792 | numberDualInfeasibilities_, sumDualInfeasibilities_, |
---|
1793 | numberFake_, dualBound_, factorization_->pivots(), |
---|
1794 | (specialOptions_ & 1024) != 0 ? 1 : 0, |
---|
1795 | (specialOptions_ & 2048) != 0 ? 1 : 0, |
---|
1796 | (specialOptions_ & 4096) != 0 ? 1 : 0, |
---|
1797 | (specialOptions_ & 8192) != 0 ? 1 : 0, |
---|
1798 | (specialOptions_ & 16384) != 0 ? 1 : 0 |
---|
1799 | ); |
---|
1800 | #endif |
---|
1801 | // Get rid of objective |
---|
1802 | if ((specialOptions_ & 16384) == 0) |
---|
1803 | objective_ = new ClpLinearObjective(NULL, numberColumns_); |
---|
1804 | } |
---|
1805 | } |
---|
1806 | #else |
---|
1807 | if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-6) { |
---|
1808 | #ifdef COIN_DEVELOP |
---|
1809 | printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n", |
---|
1810 | __LINE__, numberPrimalInfeasibilities_, |
---|
1811 | sumPrimalInfeasibilities_, |
---|
1812 | numberDualInfeasibilities_, sumDualInfeasibilities_, |
---|
1813 | numberFake_, dualBound_, factorization_->pivots()); |
---|
1814 | #endif |
---|
1815 | if ((specialOptions_ & 1024) != 0 && factorization_->pivots()) { |
---|
1816 | problemStatus_ = 10; |
---|
1817 | #if COIN_DEVELOP>1 |
---|
1818 | printf("returning at %d\n", __LINE__); |
---|
1819 | #endif |
---|
1820 | // Get rid of objective |
---|
1821 | if ((specialOptions_ & 16384) == 0) |
---|
1822 | objective_ = new ClpLinearObjective(NULL, numberColumns_); |
---|
1823 | } |
---|
1824 | } |
---|
1825 | #endif |
---|
1826 | rowArray_[0]->clear(); |
---|
1827 | columnArray_[0]->clear(); |
---|
1828 | returnCode = 1; |
---|
1829 | break; |
---|
1830 | } |
---|
1831 | } |
---|
1832 | // If special option set - put off as long as possible |
---|
1833 | if ((specialOptions_ & 64) == 0 || (moreSpecialOptions_ & 64) != 0) { |
---|
1834 | if (factorization_->pivots() == 0) |
---|
1835 | problemStatus_ = -4; //say looks infeasible |
---|
1836 | } else { |
---|
1837 | // flag |
---|
1838 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
1839 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
1840 | << x << sequenceWithin(sequenceOut_) |
---|
1841 | << CoinMessageEol; |
---|
1842 | #ifdef COIN_DEVELOP |
---|
1843 | printf("flag c\n"); |
---|
1844 | #endif |
---|
1845 | setFlagged(sequenceOut_); |
---|
1846 | if (!factorization_->pivots()) { |
---|
1847 | rowArray_[0]->clear(); |
---|
1848 | columnArray_[0]->clear(); |
---|
1849 | continue; |
---|
1850 | } |
---|
1851 | } |
---|
1852 | } |
---|
1853 | if (factorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8) |
---|
1854 | acceptablePivot_ = 1.0e-8; |
---|
1855 | rowArray_[0]->clear(); |
---|
1856 | columnArray_[0]->clear(); |
---|
1857 | returnCode = 1; |
---|
1858 | break; |
---|
1859 | } |
---|
1860 | } else { |
---|
1861 | #ifdef CLP_INVESTIGATE_SERIAL |
---|
1862 | z_thinks = 0; |
---|
1863 | #endif |
---|
1864 | // no pivot row |
---|
1865 | #ifdef CLP_DEBUG |
---|
1866 | if (handler_->logLevel() & 32) |
---|
1867 | printf("** no row pivot\n"); |
---|
1868 | #endif |
---|
1869 | // If in branch and bound try and get rid of fixed variables |
---|
1870 | if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) { |
---|
1871 | assert (!candidateList); |
---|
1872 | candidateList = new int[numberRows_]; |
---|
1873 | int iRow; |
---|
1874 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1875 | int iPivot = pivotVariable_[iRow]; |
---|
1876 | if (flagged(iPivot) || !pivoted(iPivot)) |
---|
1877 | continue; |
---|
1878 | assert (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]); |
---|
1879 | candidateList[numberCandidates++] = iRow; |
---|
1880 | } |
---|
1881 | // and set first candidate |
---|
1882 | if (!numberCandidates) { |
---|
1883 | delete [] candidateList; |
---|
1884 | candidateList = NULL; |
---|
1885 | int iRow; |
---|
1886 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1887 | int iPivot = pivotVariable_[iRow]; |
---|
1888 | clearPivoted(iPivot); |
---|
1889 | } |
---|
1890 | } else { |
---|
1891 | ifValuesPass = 2; |
---|
1892 | continue; |
---|
1893 | } |
---|
1894 | } |
---|
1895 | int numberPivots = factorization_->pivots(); |
---|
1896 | bool specialCase; |
---|
1897 | int useNumberFake; |
---|
1898 | returnCode = 0; |
---|
1899 | if (numberPivots <= CoinMax(dontFactorizePivots_, 20) && |
---|
1900 | (specialOptions_ & 2048) != 0 && (true || !numberChanged_ || perturbation_ == 101) |
---|
1901 | && dualBound_ >= 1.0e8) { |
---|
1902 | specialCase = true; |
---|
1903 | // as dual bound high - should be okay |
---|
1904 | useNumberFake = 0; |
---|
1905 | } else { |
---|
1906 | specialCase = false; |
---|
1907 | useNumberFake = numberFake_; |
---|
1908 | } |
---|
1909 | if (!numberPivots || specialCase) { |
---|
1910 | // may have crept through - so may be optimal |
---|
1911 | // check any flagged variables |
---|
1912 | int iRow; |
---|
1913 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1914 | int iPivot = pivotVariable_[iRow]; |
---|
1915 | if (flagged(iPivot)) |
---|
1916 | break; |
---|
1917 | } |
---|
1918 | if (iRow < numberRows_ && numberPivots) { |
---|
1919 | // try factorization |
---|
1920 | returnCode = -2; |
---|
1921 | } |
---|
1922 | |
---|
1923 | if (useNumberFake || numberDualInfeasibilities_) { |
---|
1924 | // may be dual infeasible |
---|
1925 | if ((specialOptions_ & 1024) == 0) |
---|
1926 | problemStatus_ = -5; |
---|
1927 | else if (!useNumberFake && numberPrimalInfeasibilities_ |
---|
1928 | && !numberPivots) |
---|
1929 | problemStatus_ = 1; |
---|
1930 | } else { |
---|
1931 | if (iRow < numberRows_) { |
---|
1932 | #ifdef COIN_DEVELOP |
---|
1933 | std::cout << "Flagged variables at end - infeasible?" << std::endl; |
---|
1934 | printf("Probably infeasible - pivot was %g\n", alpha_); |
---|
1935 | #endif |
---|
1936 | //if (fabs(alpha_)<1.0e-4) { |
---|
1937 | //problemStatus_=1; |
---|
1938 | //} else { |
---|
1939 | #ifdef CLP_DEBUG |
---|
1940 | abort(); |
---|
1941 | #endif |
---|
1942 | //} |
---|
1943 | problemStatus_ = -5; |
---|
1944 | } else { |
---|
1945 | problemStatus_ = 0; |
---|
1946 | #ifndef CLP_CHECK_NUMBER_PIVOTS |
---|
1947 | #define CLP_CHECK_NUMBER_PIVOTS 10 |
---|
1948 | #endif |
---|
1949 | #if CLP_CHECK_NUMBER_PIVOTS < 20 |
---|
1950 | if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) { |
---|
1951 | #ifndef NDEBUG_CLP |
---|
1952 | int nTotal = numberRows_ + numberColumns_; |
---|
1953 | double * comp = CoinCopyOfArray(solution_, nTotal); |
---|
1954 | #endif |
---|
1955 | computePrimals(rowActivityWork_, columnActivityWork_); |
---|
1956 | #ifndef NDEBUG_CLP |
---|
1957 | double largest = 1.0e-5; |
---|
1958 | int bad = -1; |
---|
1959 | for (int i = 0; i < nTotal; i++) { |
---|
1960 | double value = solution_[i]; |
---|
1961 | double larger = CoinMax(fabs(value), fabs(comp[i])); |
---|
1962 | double tol = 1.0e-5 + 1.0e-5 * larger; |
---|
1963 | double diff = fabs(value - comp[i]); |
---|
1964 | if (diff - tol > largest) { |
---|
1965 | bad = i; |
---|
1966 | largest = diff - tol; |
---|
1967 | } |
---|
1968 | } |
---|
1969 | if (bad >= 0) |
---|
1970 | COIN_DETAIL_PRINT(printf("bad %d old %g new %g\n", bad, comp[bad], solution_[bad])); |
---|
1971 | #endif |
---|
1972 | checkPrimalSolution(rowActivityWork_, columnActivityWork_); |
---|
1973 | if (numberPrimalInfeasibilities_) { |
---|
1974 | #ifdef CLP_INVESTIGATE |
---|
1975 | printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_, |
---|
1976 | sumPrimalInfeasibilities_); |
---|
1977 | #endif |
---|
1978 | problemStatus_ = -1; |
---|
1979 | returnCode = -2; |
---|
1980 | } |
---|
1981 | #ifndef NDEBUG_CLP |
---|
1982 | memcpy(solution_, comp, nTotal * sizeof(double)); |
---|
1983 | delete [] comp; |
---|
1984 | #endif |
---|
1985 | } |
---|
1986 | #endif |
---|
1987 | if (!problemStatus_) { |
---|
1988 | // make it look OK |
---|
1989 | numberPrimalInfeasibilities_ = 0; |
---|
1990 | sumPrimalInfeasibilities_ = 0.0; |
---|
1991 | numberDualInfeasibilities_ = 0; |
---|
1992 | sumDualInfeasibilities_ = 0.0; |
---|
1993 | // May be perturbed |
---|
1994 | if (perturbation_ == 101 || numberChanged_) { |
---|
1995 | numberChanged_ = 0; // Number of variables with changed costs |
---|
1996 | perturbation_ = 102; // stop any perturbations |
---|
1997 | //double changeCost; |
---|
1998 | //changeBounds(1,NULL,changeCost); |
---|
1999 | createRim4(false); |
---|
2000 | // make sure duals are current |
---|
2001 | computeDuals(givenDuals); |
---|
2002 | checkDualSolution(); |
---|
2003 | progress_.modifyObjective(-COIN_DBL_MAX); |
---|
2004 | if (numberDualInfeasibilities_) { |
---|
2005 | problemStatus_ = 10; // was -3; |
---|
2006 | } else { |
---|
2007 | computeObjectiveValue(true); |
---|
2008 | } |
---|
2009 | } else if (numberPivots) { |
---|
2010 | computeObjectiveValue(true); |
---|
2011 | } |
---|
2012 | if (numberPivots < -1000) { |
---|
2013 | // objective may be wrong |
---|
2014 | objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_); |
---|
2015 | objectiveValue_ += objective_->nonlinearOffset(); |
---|
2016 | objectiveValue_ /= (objectiveScale_ * rhsScale_); |
---|
2017 | if ((specialOptions_ & 16384) == 0) { |
---|
2018 | // and dual_ may be wrong (i.e. for fixed or basic) |
---|
2019 | CoinIndexedVector * arrayVector = rowArray_[1]; |
---|
2020 | arrayVector->clear(); |
---|
2021 | int iRow; |
---|
2022 | double * array = arrayVector->denseVector(); |
---|
2023 | /* Use dual_ instead of array |
---|
2024 | Even though dual_ is only numberRows_ long this is |
---|
2025 | okay as gets permuted to longer rowArray_[2] |
---|
2026 | */ |
---|
2027 | arrayVector->setDenseVector(dual_); |
---|
2028 | int * index = arrayVector->getIndices(); |
---|
2029 | int number = 0; |
---|
2030 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2031 | int iPivot = pivotVariable_[iRow]; |
---|
2032 | double value = cost_[iPivot]; |
---|
2033 | dual_[iRow] = value; |
---|
2034 | if (value) { |
---|
2035 | index[number++] = iRow; |
---|
2036 | } |
---|
2037 | } |
---|
2038 | arrayVector->setNumElements(number); |
---|
2039 | // Extended duals before "updateTranspose" |
---|
2040 | matrix_->dualExpanded(this, arrayVector, NULL, 0); |
---|
2041 | // Btran basic costs |
---|
2042 | rowArray_[2]->clear(); |
---|
2043 | factorization_->updateColumnTranspose(rowArray_[2], arrayVector); |
---|
2044 | // and return vector |
---|
2045 | arrayVector->setDenseVector(array); |
---|
2046 | } |
---|
2047 | } |
---|
2048 | sumPrimalInfeasibilities_ = 0.0; |
---|
2049 | } |
---|
2050 | if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_) { |
---|
2051 | CoinIndexedVector * arrayVector = rowArray_[1]; |
---|
2052 | arrayVector->clear(); |
---|
2053 | double * rhs = arrayVector->denseVector(); |
---|
2054 | times(1.0, solution_, rhs); |
---|
2055 | #ifdef CHECK_ACCURACY |
---|
2056 | bool bad = false; |
---|
2057 | #endif |
---|
2058 | bool bad2 = false; |
---|
2059 | int i; |
---|
2060 | for ( i = 0; i < numberRows_; i++) { |
---|
2061 | if (rhs[i] < rowLowerWork_[i] - primalTolerance_ || |
---|
2062 | rhs[i] > rowUpperWork_[i] + primalTolerance_) { |
---|
2063 | bad2 = true; |
---|
2064 | #ifdef CHECK_ACCURACY |
---|
2065 | printf("row %d out of bounds %g, %g correct %g bad %g\n", i, |
---|
2066 | rowLowerWork_[i], rowUpperWork_[i], |
---|
2067 | rhs[i], rowActivityWork_[i]); |
---|
2068 | #endif |
---|
2069 | } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) { |
---|
2070 | #ifdef CHECK_ACCURACY |
---|
2071 | bad = true; |
---|
2072 | printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]); |
---|
2073 | #endif |
---|
2074 | } |
---|
2075 | rhs[i] = 0.0; |
---|
2076 | } |
---|
2077 | for ( i = 0; i < numberColumns_; i++) { |
---|
2078 | if (solution_[i] < columnLowerWork_[i] - primalTolerance_ || |
---|
2079 | solution_[i] > columnUpperWork_[i] + primalTolerance_) { |
---|
2080 | bad2 = true; |
---|
2081 | #ifdef CHECK_ACCURACY |
---|
2082 | printf("column %d out of bounds %g, %g correct %g bad %g\n", i, |
---|
2083 | columnLowerWork_[i], columnUpperWork_[i], |
---|
2084 | solution_[i], columnActivityWork_[i]); |
---|
2085 | #endif |
---|
2086 | } |
---|
2087 | } |
---|
2088 | if (bad2) { |
---|
2089 | problemStatus_ = -3; |
---|
2090 | returnCode = -2; |
---|
2091 | // Force to re-factorize early next time |
---|
2092 | int numberPivots = factorization_->pivots(); |
---|
2093 | forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); |
---|
2094 | } |
---|
2095 | } |
---|
2096 | } |
---|
2097 | } |
---|
2098 | } else { |
---|
2099 | problemStatus_ = -3; |
---|
2100 | returnCode = -2; |
---|
2101 | // Force to re-factorize early next time |
---|
2102 | int numberPivots = factorization_->pivots(); |
---|
2103 | forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); |
---|
2104 | } |
---|
2105 | break; |
---|
2106 | } |
---|
2107 | } |
---|
2108 | if (givenDuals) { |
---|
2109 | CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals); |
---|
2110 | // get rid of any values pass array |
---|
2111 | delete [] candidateList; |
---|
2112 | } |
---|
2113 | delete [] dubiousWeights; |
---|
2114 | #ifdef CLP_REPORT_PROGRESS |
---|
2115 | if (ixxxxxx > ixxyyyy - 5) { |
---|
2116 | int nTotal = numberColumns_ + numberRows_; |
---|
2117 | double oldObj = 0.0; |
---|
2118 | double newObj = 0.0; |
---|
2119 | for (int i = 0; i < nTotal; i++) { |
---|
2120 | if (savePSol[i]) |
---|
2121 | oldObj += savePSol[i] * saveCost[i]; |
---|
2122 | if (solution_[i]) |
---|
2123 | newObj += solution_[i] * cost_[i]; |
---|
2124 | bool printIt = false; |
---|
2125 | if (cost_[i] != saveCost[i]) |
---|
2126 | printIt = true; |
---|
2127 | if (status_[i] != saveStat[i]) |
---|
2128 | printIt = true; |
---|
2129 | if (printIt) |
---|
2130 | printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n", |
---|
2131 | i, saveStat[i], saveCost[i], savePSol[i], |
---|
2132 | status_[i], cost_[i], solution_[i]); |
---|
2133 | // difference |
---|
2134 | savePSol[i] = solution_[i] - savePSol[i]; |
---|
2135 | } |
---|
2136 | printf("exit pivots %d, old obj %g new %g\n", |
---|
2137 | factorization_->pivots(), |
---|
2138 | oldObj, newObj); |
---|
2139 | memset(saveDj, 0, numberRows_ * sizeof(double)); |
---|
2140 | times(1.0, savePSol, saveDj); |
---|
2141 | double largest = 1.0e-6; |
---|
2142 | int k = -1; |
---|
2143 | for (int i = 0; i < numberRows_; i++) { |
---|
2144 | saveDj[i] -= savePSol[i+numberColumns_]; |
---|
2145 | if (fabs(saveDj[i]) > largest) { |
---|
2146 | largest = fabs(saveDj[i]); |
---|
2147 | k = i; |
---|
2148 | } |
---|
2149 | } |
---|
2150 | if (k >= 0) |
---|
2151 | printf("Not null %d %g\n", k, largest); |
---|
2152 | } |
---|
2153 | delete [] savePSol ; |
---|
2154 | delete [] saveDj ; |
---|
2155 | delete [] saveCost ; |
---|
2156 | delete [] saveStat ; |
---|
2157 | #endif |
---|
2158 | return returnCode; |
---|
2159 | } |
---|
2160 | /* The duals are updated by the given arrays. |
---|
2161 | Returns number of infeasibilities. |
---|
2162 | rowArray and columnarray will have flipped |
---|
2163 | The output vector has movement (row length array) */ |
---|
2164 | int |
---|
2165 | ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray, |
---|
2166 | CoinIndexedVector * columnArray, |
---|
2167 | CoinIndexedVector * outputArray, |
---|
2168 | double theta, |
---|
2169 | double & objectiveChange, |
---|
2170 | bool fullRecompute) |
---|
2171 | { |
---|
2172 | |
---|
2173 | outputArray->clear(); |
---|
2174 | |
---|
2175 | |
---|
2176 | int numberInfeasibilities = 0; |
---|
2177 | int numberRowInfeasibilities = 0; |
---|
2178 | |
---|
2179 | // get a tolerance |
---|
2180 | double tolerance = dualTolerance_; |
---|
2181 | // we can't really trust infeasibilities if there is dual error |
---|
2182 | double error = CoinMin(1.0e-2, largestDualError_); |
---|
2183 | // allow tolerance at least slightly bigger than standard |
---|
2184 | tolerance = tolerance + error; |
---|
2185 | |
---|
2186 | double changeObj = 0.0; |
---|
2187 | |
---|
2188 | // Coding is very similar but we can save a bit by splitting |
---|
2189 | // Do rows |
---|
2190 | if (!fullRecompute) { |
---|
2191 | int i; |
---|
2192 | double * COIN_RESTRICT reducedCost = djRegion(0); |
---|
2193 | const double * COIN_RESTRICT lower = lowerRegion(0); |
---|
2194 | const double * COIN_RESTRICT upper = upperRegion(0); |
---|
2195 | const double * COIN_RESTRICT cost = costRegion(0); |
---|
2196 | double * COIN_RESTRICT work; |
---|
2197 | int number; |
---|
2198 | int * COIN_RESTRICT which; |
---|
2199 | const unsigned char * COIN_RESTRICT statusArray = status_ + numberColumns_; |
---|
2200 | assert(rowArray->packedMode()); |
---|
2201 | work = rowArray->denseVector(); |
---|
2202 | number = rowArray->getNumElements(); |
---|
2203 | which = rowArray->getIndices(); |
---|
2204 | double multiplier[] = { -1.0, 1.0}; |
---|
2205 | for (i = 0; i < number; i++) { |
---|
2206 | int iSequence = which[i]; |
---|
2207 | double alphaI = work[i]; |
---|
2208 | work[i] = 0.0; |
---|
2209 | int iStatus = (statusArray[iSequence] & 3) - 1; |
---|
2210 | if (iStatus) { |
---|
2211 | double value = reducedCost[iSequence] - theta * alphaI; |
---|
2212 | reducedCost[iSequence] = value; |
---|
2213 | double mult = multiplier[iStatus-1]; |
---|
2214 | value *= mult; |
---|
2215 | if (value < -tolerance) { |
---|
2216 | // flipping bounds |
---|
2217 | double movement = mult * (lower[iSequence] - upper[iSequence]); |
---|
2218 | which[numberInfeasibilities++] = iSequence; |
---|
2219 | #ifndef NDEBUG |
---|
2220 | if (fabs(movement) >= 1.0e30) |
---|
2221 | resetFakeBounds(-1000 - iSequence); |
---|
2222 | #endif |
---|
2223 | #ifdef CLP_DEBUG |
---|
2224 | if ((handler_->logLevel() & 32)) |
---|
2225 | printf("%d %d, new dj %g, alpha %g, movement %g\n", |
---|
2226 | 0, iSequence, value, alphaI, movement); |
---|
2227 | #endif |
---|
2228 | changeObj -= movement * cost[iSequence]; |
---|
2229 | outputArray->quickAdd(iSequence, movement); |
---|
2230 | } |
---|
2231 | } |
---|
2232 | } |
---|
2233 | // Do columns |
---|
2234 | reducedCost = djRegion(1); |
---|
2235 | lower = lowerRegion(1); |
---|
2236 | upper = upperRegion(1); |
---|
2237 | cost = costRegion(1); |
---|
2238 | // set number of infeasibilities in row array |
---|
2239 | numberRowInfeasibilities = numberInfeasibilities; |
---|
2240 | rowArray->setNumElements(numberInfeasibilities); |
---|
2241 | numberInfeasibilities = 0; |
---|
2242 | work = columnArray->denseVector(); |
---|
2243 | number = columnArray->getNumElements(); |
---|
2244 | which = columnArray->getIndices(); |
---|
2245 | if ((moreSpecialOptions_ & 8) != 0) { |
---|
2246 | const unsigned char * COIN_RESTRICT statusArray = status_; |
---|
2247 | for (i = 0; i < number; i++) { |
---|
2248 | int iSequence = which[i]; |
---|
2249 | double alphaI = work[i]; |
---|
2250 | work[i] = 0.0; |
---|
2251 | |
---|
2252 | int iStatus = (statusArray[iSequence] & 3) - 1; |
---|
2253 | if (iStatus) { |
---|
2254 | double value = reducedCost[iSequence] - theta * alphaI; |
---|
2255 | reducedCost[iSequence] = value; |
---|
2256 | double mult = multiplier[iStatus-1]; |
---|
2257 | value *= mult; |
---|
2258 | if (value < -tolerance) { |
---|
2259 | // flipping bounds |
---|
2260 | double movement = mult * (upper[iSequence] - lower[iSequence]); |
---|
2261 | which[numberInfeasibilities++] = iSequence; |
---|
2262 | #ifndef NDEBUG |
---|
2263 | if (fabs(movement) >= 1.0e30) |
---|
2264 | resetFakeBounds(-1000 - iSequence); |
---|
2265 | #endif |
---|
2266 | #ifdef CLP_DEBUG |
---|
2267 | if ((handler_->logLevel() & 32)) |
---|
2268 | printf("%d %d, new dj %g, alpha %g, movement %g\n", |
---|
2269 | 1, iSequence, value, alphaI, movement); |
---|
2270 | #endif |
---|
2271 | changeObj += movement * cost[iSequence]; |
---|
2272 | matrix_->add(this, outputArray, iSequence, movement); |
---|
2273 | } |
---|
2274 | } |
---|
2275 | } |
---|
2276 | } else { |
---|
2277 | for (i = 0; i < number; i++) { |
---|
2278 | int iSequence = which[i]; |
---|
2279 | double alphaI = work[i]; |
---|
2280 | work[i] = 0.0; |
---|
2281 | |
---|
2282 | Status status = getStatus(iSequence); |
---|
2283 | if (status == atLowerBound) { |
---|
2284 | double value = reducedCost[iSequence] - theta * alphaI; |
---|
2285 | reducedCost[iSequence] = value; |
---|
2286 | double movement = 0.0; |
---|
2287 | |
---|
2288 | if (value < -tolerance) { |
---|
2289 | // to upper bound |
---|
2290 | which[numberInfeasibilities++] = iSequence; |
---|
2291 | movement = upper[iSequence] - lower[iSequence]; |
---|
2292 | #ifndef NDEBUG |
---|
2293 | if (fabs(movement) >= 1.0e30) |
---|
2294 | resetFakeBounds(-1000 - iSequence); |
---|
2295 | #endif |
---|
2296 | #ifdef CLP_DEBUG |
---|
2297 | if ((handler_->logLevel() & 32)) |
---|
2298 | printf("%d %d, new dj %g, alpha %g, movement %g\n", |
---|
2299 | 1, iSequence, value, alphaI, movement); |
---|
2300 | #endif |
---|
2301 | changeObj += movement * cost[iSequence]; |
---|
2302 | matrix_->add(this, outputArray, iSequence, movement); |
---|
2303 | } |
---|
2304 | } else if (status == atUpperBound) { |
---|
2305 | double value = reducedCost[iSequence] - theta * alphaI; |
---|
2306 | reducedCost[iSequence] = value; |
---|
2307 | double movement = 0.0; |
---|
2308 | |
---|
2309 | if (value > tolerance) { |
---|
2310 | // to lower bound (if swap) |
---|
2311 | which[numberInfeasibilities++] = iSequence; |
---|
2312 | movement = lower[iSequence] - upper[iSequence]; |
---|
2313 | #ifndef NDEBUG |
---|
2314 | if (fabs(movement) >= 1.0e30) |
---|
2315 | resetFakeBounds(-1000 - iSequence); |
---|
2316 | #endif |
---|
2317 | #ifdef CLP_DEBUG |
---|
2318 | if ((handler_->logLevel() & 32)) |
---|
2319 | printf("%d %d, new dj %g, alpha %g, movement %g\n", |
---|
2320 | 1, iSequence, value, alphaI, movement); |
---|
2321 | #endif |
---|
2322 | changeObj += movement * cost[iSequence]; |
---|
2323 | matrix_->add(this, outputArray, iSequence, movement); |
---|
2324 | } |
---|
2325 | } else if (status == isFree) { |
---|
2326 | double value = reducedCost[iSequence] - theta * alphaI; |
---|
2327 | reducedCost[iSequence] = value; |
---|
2328 | } |
---|
2329 | } |
---|
2330 | } |
---|
2331 | } else { |
---|
2332 | double * COIN_RESTRICT solution = solutionRegion(0); |
---|
2333 | double * COIN_RESTRICT reducedCost = djRegion(0); |
---|
2334 | const double * COIN_RESTRICT lower = lowerRegion(0); |
---|
2335 | const double * COIN_RESTRICT upper = upperRegion(0); |
---|
2336 | const double * COIN_RESTRICT cost = costRegion(0); |
---|
2337 | int * COIN_RESTRICT which; |
---|
2338 | which = rowArray->getIndices(); |
---|
2339 | int iSequence; |
---|
2340 | for (iSequence = 0; iSequence < numberRows_; iSequence++) { |
---|
2341 | double value = reducedCost[iSequence]; |
---|
2342 | |
---|
2343 | Status status = getStatus(iSequence + numberColumns_); |
---|
2344 | // more likely to be at upper bound ? |
---|
2345 | if (status == atUpperBound) { |
---|
2346 | double movement = 0.0; |
---|
2347 | //#define NO_SWAP7 |
---|
2348 | if (value > tolerance) { |
---|
2349 | // to lower bound (if swap) |
---|
2350 | // put back alpha |
---|
2351 | which[numberInfeasibilities++] = iSequence; |
---|
2352 | movement = lower[iSequence] - upper[iSequence]; |
---|
2353 | changeObj += movement * cost[iSequence]; |
---|
2354 | outputArray->quickAdd(iSequence, -movement); |
---|
2355 | #ifndef NDEBUG |
---|
2356 | if (fabs(movement) >= 1.0e30) |
---|
2357 | resetFakeBounds(-1000 - iSequence); |
---|
2358 | #endif |
---|
2359 | #ifndef NO_SWAP7 |
---|
2360 | } else if (value > -tolerance) { |
---|
2361 | // at correct bound but may swap |
---|
2362 | FakeBound bound = getFakeBound(iSequence + numberColumns_); |
---|
2363 | if (bound == ClpSimplexDual::upperFake) { |
---|
2364 | movement = lower[iSequence] - upper[iSequence]; |
---|
2365 | #ifndef NDEBUG |
---|
2366 | if (fabs(movement) >= 1.0e30) |
---|
2367 | resetFakeBounds(-1000 - iSequence); |
---|
2368 | #endif |
---|
2369 | setStatus(iSequence + numberColumns_, atLowerBound); |
---|
2370 | solution[iSequence] = lower[iSequence]; |
---|
2371 | changeObj += movement * cost[iSequence]; |
---|
2372 | //numberFake_--; |
---|
2373 | //setFakeBound(iSequence+numberColumns_,noFake); |
---|
2374 | } |
---|
2375 | #endif |
---|
2376 | } |
---|
2377 | } else if (status == atLowerBound) { |
---|
2378 | double movement = 0.0; |
---|
2379 | |
---|
2380 | if (value < -tolerance) { |
---|
2381 | // to upper bound |
---|
2382 | // put back alpha |
---|
2383 | which[numberInfeasibilities++] = iSequence; |
---|
2384 | movement = upper[iSequence] - lower[iSequence]; |
---|
2385 | #ifndef NDEBUG |
---|
2386 | if (fabs(movement) >= 1.0e30) |
---|
2387 | resetFakeBounds(-1000 - iSequence); |
---|
2388 | #endif |
---|
2389 | changeObj += movement * cost[iSequence]; |
---|
2390 | outputArray->quickAdd(iSequence, -movement); |
---|
2391 | #ifndef NO_SWAP7 |
---|
2392 | } else if (value < tolerance) { |
---|
2393 | // at correct bound but may swap |
---|
2394 | FakeBound bound = getFakeBound(iSequence + numberColumns_); |
---|
2395 | if (bound == ClpSimplexDual::lowerFake) { |
---|
2396 | movement = upper[iSequence] - lower[iSequence]; |
---|
2397 | #ifndef NDEBUG |
---|
2398 | if (fabs(movement) >= 1.0e30) |
---|
2399 | resetFakeBounds(-1000 - iSequence); |
---|
2400 | #endif |
---|
2401 | setStatus(iSequence + numberColumns_, atUpperBound); |
---|
2402 | solution[iSequence] = upper[iSequence]; |
---|
2403 | changeObj += movement * cost[iSequence]; |
---|
2404 | //numberFake_--; |
---|
2405 | //setFakeBound(iSequence+numberColumns_,noFake); |
---|
2406 | } |
---|
2407 | #endif |
---|
2408 | } |
---|
2409 | } |
---|
2410 | } |
---|
2411 | // Do columns |
---|
2412 | solution = solutionRegion(1); |
---|
2413 | reducedCost = djRegion(1); |
---|
2414 | lower = lowerRegion(1); |
---|
2415 | upper = upperRegion(1); |
---|
2416 | cost = costRegion(1); |
---|
2417 | // set number of infeasibilities in row array |
---|
2418 | numberRowInfeasibilities = numberInfeasibilities; |
---|
2419 | rowArray->setNumElements(numberInfeasibilities); |
---|
2420 | numberInfeasibilities = 0; |
---|
2421 | which = columnArray->getIndices(); |
---|
2422 | for (iSequence = 0; iSequence < numberColumns_; iSequence++) { |
---|
2423 | double value = reducedCost[iSequence]; |
---|
2424 | |
---|
2425 | Status status = getStatus(iSequence); |
---|
2426 | if (status == atLowerBound) { |
---|
2427 | double movement = 0.0; |
---|
2428 | |
---|
2429 | if (value < -tolerance) { |
---|
2430 | // to upper bound |
---|
2431 | // put back alpha |
---|
2432 | which[numberInfeasibilities++] = iSequence; |
---|
2433 | movement = upper[iSequence] - lower[iSequence]; |
---|
2434 | #ifndef NDEBUG |
---|
2435 | if (fabs(movement) >= 1.0e30) |
---|
2436 | resetFakeBounds(-1000 - iSequence); |
---|
2437 | #endif |
---|
2438 | changeObj += movement * cost[iSequence]; |
---|
2439 | matrix_->add(this, outputArray, iSequence, movement); |
---|
2440 | #ifndef NO_SWAP7 |
---|
2441 | } else if (value < tolerance) { |
---|
2442 | // at correct bound but may swap |
---|
2443 | FakeBound bound = getFakeBound(iSequence); |
---|
2444 | if (bound == ClpSimplexDual::lowerFake) { |
---|
2445 | movement = upper[iSequence] - lower[iSequence]; |
---|
2446 | #ifndef NDEBUG |
---|
2447 | if (fabs(movement) >= 1.0e30) |
---|
2448 | resetFakeBounds(-1000 - iSequence); |
---|
2449 | #endif |
---|
2450 | setStatus(iSequence, atUpperBound); |
---|
2451 | solution[iSequence] = upper[iSequence]; |
---|
2452 | changeObj += movement * cost[iSequence]; |
---|
2453 | //numberFake_--; |
---|
2454 | //setFakeBound(iSequence,noFake); |
---|
2455 | } |
---|
2456 | #endif |
---|
2457 | } |
---|
2458 | } else if (status == atUpperBound) { |
---|
2459 | double movement = 0.0; |
---|
2460 | |
---|
2461 | if (value > tolerance) { |
---|
2462 | // to lower bound (if swap) |
---|
2463 | // put back alpha |
---|
2464 | which[numberInfeasibilities++] = iSequence; |
---|
2465 | movement = lower[iSequence] - upper[iSequence]; |
---|
2466 | #ifndef NDEBUG |
---|
2467 | if (fabs(movement) >= 1.0e30) |
---|
2468 | resetFakeBounds(-1000 - iSequence); |
---|
2469 | #endif |
---|
2470 | changeObj += movement * cost[iSequence]; |
---|
2471 | matrix_->add(this, outputArray, iSequence, movement); |
---|
2472 | #ifndef NO_SWAP7 |
---|
2473 | } else if (value > -tolerance) { |
---|
2474 | // at correct bound but may swap |
---|
2475 | FakeBound bound = getFakeBound(iSequence); |
---|
2476 | if (bound == ClpSimplexDual::upperFake) { |
---|
2477 | movement = lower[iSequence] - upper[iSequence]; |
---|
2478 | #ifndef NDEBUG |
---|
2479 | if (fabs(movement) >= 1.0e30) |
---|
2480 | resetFakeBounds(-1000 - iSequence); |
---|
2481 | #endif |
---|
2482 | setStatus(iSequence, atLowerBound); |
---|
2483 | solution[iSequence] = lower[iSequence]; |
---|
2484 | changeObj += movement * cost[iSequence]; |
---|
2485 | //numberFake_--; |
---|
2486 | //setFakeBound(iSequence,noFake); |
---|
2487 | } |
---|
2488 | #endif |
---|
2489 | } |
---|
2490 | } |
---|
2491 | } |
---|
2492 | } |
---|
2493 | |
---|
2494 | #ifdef CLP_DEBUG |
---|
2495 | if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0) |
---|
2496 | printf("%d fake after full update\n", numberFake_); |
---|
2497 | #endif |
---|
2498 | // set number of infeasibilities |
---|
2499 | columnArray->setNumElements(numberInfeasibilities); |
---|
2500 | numberInfeasibilities += numberRowInfeasibilities; |
---|
2501 | if (fullRecompute) { |
---|
2502 | // do actual flips |
---|
2503 | flipBounds(rowArray, columnArray); |
---|
2504 | } |
---|
2505 | objectiveChange += changeObj; |
---|
2506 | return numberInfeasibilities; |
---|
2507 | } |
---|
2508 | void |
---|
2509 | ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray, |
---|
2510 | CoinIndexedVector * columnArray, |
---|
2511 | double theta) |
---|
2512 | { |
---|
2513 | |
---|
2514 | // use a tighter tolerance except for all being okay |
---|
2515 | double tolerance = dualTolerance_; |
---|
2516 | |
---|
2517 | // Coding is very similar but we can save a bit by splitting |
---|
2518 | // Do rows |
---|
2519 | { |
---|
2520 | int i; |
---|
2521 | double * reducedCost = djRegion(0); |
---|
2522 | double * work; |
---|
2523 | int number; |
---|
2524 | int * which; |
---|
2525 | work = rowArray->denseVector(); |
---|
2526 | number = rowArray->getNumElements(); |
---|
2527 | which = rowArray->getIndices(); |
---|
2528 | for (i = 0; i < number; i++) { |
---|
2529 | int iSequence = which[i]; |
---|
2530 | double alphaI = work[i]; |
---|
2531 | double value = reducedCost[iSequence] - theta * alphaI; |
---|
2532 | work[i] = 0.0; |
---|
2533 | reducedCost[iSequence] = value; |
---|
2534 | |
---|
2535 | Status status = getStatus(iSequence + numberColumns_); |
---|
2536 | // more likely to be at upper bound ? |
---|
2537 | if (status == atUpperBound) { |
---|
2538 | |
---|
2539 | if (value > tolerance) |
---|
2540 | reducedCost[iSequence] = 0.0; |
---|
2541 | } else if (status == atLowerBound) { |
---|
2542 | |
---|
2543 | if (value < -tolerance) { |
---|
2544 | reducedCost[iSequence] = 0.0; |
---|
2545 | } |
---|
2546 | } |
---|
2547 | } |
---|
2548 | } |
---|
2549 | rowArray->setNumElements(0); |
---|
2550 | |
---|
2551 | // Do columns |
---|
2552 | { |
---|
2553 | int i; |
---|
2554 | double * reducedCost = djRegion(1); |
---|
2555 | double * work; |
---|
2556 | int number; |
---|
2557 | int * which; |
---|
2558 | work = columnArray->denseVector(); |
---|
2559 | number = columnArray->getNumElements(); |
---|
2560 | which = columnArray->getIndices(); |
---|
2561 | |
---|
2562 | for (i = 0; i < number; i++) { |
---|
2563 | int iSequence = which[i]; |
---|
2564 | double alphaI = work[i]; |
---|
2565 | double value = reducedCost[iSequence] - theta * alphaI; |
---|
2566 | work[i] = 0.0; |
---|
2567 | reducedCost[iSequence] = value; |
---|
2568 | |
---|
2569 | Status status = getStatus(iSequence); |
---|
2570 | if (status == atLowerBound) { |
---|
2571 | if (value < -tolerance) |
---|
2572 | reducedCost[iSequence] = 0.0; |
---|
2573 | } else if (status == atUpperBound) { |
---|
2574 | if (value > tolerance) |
---|
2575 | reducedCost[iSequence] = 0.0; |
---|
2576 | } |
---|
2577 | } |
---|
2578 | } |
---|
2579 | columnArray->setNumElements(0); |
---|
2580 | } |
---|
2581 | /* |
---|
2582 | Chooses dual pivot row |
---|
2583 | Would be faster with separate region to scan |
---|
2584 | and will have this (with square of infeasibility) when steepest |
---|
2585 | For easy problems we can just choose one of the first rows we look at |
---|
2586 | */ |
---|
2587 | void |
---|
2588 | ClpSimplexDual::dualRow(int alreadyChosen) |
---|
2589 | { |
---|
2590 | // get pivot row using whichever method it is |
---|
2591 | int chosenRow = -1; |
---|
2592 | #ifdef FORCE_FOLLOW |
---|
2593 | bool forceThis = false; |
---|
2594 | if (!fpFollow && strlen(forceFile)) { |
---|
2595 | fpFollow = fopen(forceFile, "r"); |
---|
2596 | assert (fpFollow); |
---|
2597 | } |
---|
2598 | if (fpFollow) { |
---|
2599 | if (numberIterations_ <= force_iteration) { |
---|
2600 | // read to next Clp0102 |
---|
2601 | char temp[300]; |
---|
2602 | while (fgets(temp, 250, fpFollow)) { |
---|
2603 | if (strncmp(temp, "Clp0102", 7)) |
---|
2604 | continue; |
---|
2605 | char cin, cout; |
---|
2606 | sscanf(temp + 9, "%d%*f%*s%*c%c%d%*s%*c%c%d", |
---|
2607 | &force_iteration, &cin, &force_in, &cout, &force_out); |
---|
2608 | if (cin == 'R') |
---|
2609 | force_in += numberColumns_; |
---|
2610 | if (cout == 'R') |
---|
2611 | force_out += numberColumns_; |
---|
2612 | forceThis = true; |
---|
2613 | assert (numberIterations_ == force_iteration - 1); |
---|
2614 | printf("Iteration %d will force %d out and %d in\n", |
---|
2615 | force_iteration, force_out, force_in); |
---|
2616 | alreadyChosen = force_out; |
---|
2617 | break; |
---|
2618 | } |
---|
2619 | } else { |
---|
2620 | // use old |
---|
2621 | forceThis = true; |
---|
2622 | } |
---|
2623 | if (!forceThis) { |
---|
2624 | fclose(fpFollow); |
---|
2625 | fpFollow = NULL; |
---|
2626 | forceFile = ""; |
---|
2627 | } |
---|
2628 | } |
---|
2629 | #endif |
---|
2630 | double freeAlpha = 0.0; |
---|
2631 | if (alreadyChosen < 0) { |
---|
2632 | // first see if any free variables and put them in basis |
---|
2633 | int nextFree = nextSuperBasic(); |
---|
2634 | //nextFree=-1; //off |
---|
2635 | if (nextFree >= 0) { |
---|
2636 | // unpack vector and find a good pivot |
---|
2637 | unpack(rowArray_[1], nextFree); |
---|
2638 | factorization_->updateColumn(rowArray_[2], rowArray_[1]); |
---|
2639 | |
---|
2640 | double * work = rowArray_[1]->denseVector(); |
---|
2641 | int number = rowArray_[1]->getNumElements(); |
---|
2642 | int * which = rowArray_[1]->getIndices(); |
---|
2643 | double bestFeasibleAlpha = 0.0; |
---|
2644 | int bestFeasibleRow = -1; |
---|
2645 | double bestInfeasibleAlpha = 0.0; |
---|
2646 | int bestInfeasibleRow = -1; |
---|
2647 | int i; |
---|
2648 | |
---|
2649 | for (i = 0; i < number; i++) { |
---|
2650 | int iRow = which[i]; |
---|
2651 | double alpha = fabs(work[iRow]); |
---|
2652 | if (alpha > 1.0e-3) { |
---|
2653 | int iSequence = pivotVariable_[iRow]; |
---|
2654 | double value = solution_[iSequence]; |
---|
2655 | double lower = lower_[iSequence]; |
---|
2656 | double upper = upper_[iSequence]; |
---|
2657 | double infeasibility = 0.0; |
---|
2658 | if (value > upper) |
---|
2659 | infeasibility = value - upper; |
---|
2660 | else if (value < lower) |
---|
2661 | infeasibility = lower - value; |
---|
2662 | if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) { |
---|
2663 | if (!flagged(iSequence)) { |
---|
2664 | bestInfeasibleAlpha = infeasibility * alpha; |
---|
2665 | bestInfeasibleRow = iRow; |
---|
2666 | } |
---|
2667 | } |
---|
2668 | if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) { |
---|
2669 | bestFeasibleAlpha = alpha; |
---|
2670 | bestFeasibleRow = iRow; |
---|
2671 | } |
---|
2672 | } |
---|
2673 | } |
---|
2674 | if (bestInfeasibleRow >= 0) |
---|
2675 | chosenRow = bestInfeasibleRow; |
---|
2676 | else if (bestFeasibleAlpha > 1.0e-2) |
---|
2677 | chosenRow = bestFeasibleRow; |
---|
2678 | if (chosenRow >= 0) { |
---|
2679 | pivotRow_ = chosenRow; |
---|
2680 | freeAlpha = work[chosenRow]; |
---|
2681 | } |
---|
2682 | rowArray_[1]->clear(); |
---|
2683 | } |
---|
2684 | } else { |
---|
2685 | // in values pass |
---|
2686 | chosenRow = alreadyChosen; |
---|
2687 | #ifdef FORCE_FOLLOW |
---|
2688 | if(forceThis) { |
---|
2689 | alreadyChosen = -1; |
---|
2690 | chosenRow = -1; |
---|
2691 | for (int i = 0; i < numberRows_; i++) { |
---|
2692 | if (pivotVariable_[i] == force_out) { |
---|
2693 | chosenRow = i; |
---|
2694 | break; |
---|
2695 | } |
---|
2696 | } |
---|
2697 | assert (chosenRow >= 0); |
---|
2698 | } |
---|
2699 | #endif |
---|
2700 | pivotRow_ = chosenRow; |
---|
2701 | } |
---|
2702 | if (chosenRow < 0) |
---|
2703 | pivotRow_ = dualRowPivot_->pivotRow(); |
---|
2704 | |
---|
2705 | if (pivotRow_ >= 0) { |
---|
2706 | sequenceOut_ = pivotVariable_[pivotRow_]; |
---|
2707 | valueOut_ = solution_[sequenceOut_]; |
---|
2708 | lowerOut_ = lower_[sequenceOut_]; |
---|
2709 | upperOut_ = upper_[sequenceOut_]; |
---|
2710 | if (alreadyChosen < 0) { |
---|
2711 | // if we have problems we could try other way and hope we get a |
---|
2712 | // zero pivot? |
---|
2713 | if (valueOut_ > upperOut_) { |
---|
2714 | directionOut_ = -1; |
---|
2715 | dualOut_ = valueOut_ - upperOut_; |
---|
2716 | } else if (valueOut_ < lowerOut_) { |
---|
2717 | directionOut_ = 1; |
---|
2718 | dualOut_ = lowerOut_ - valueOut_; |
---|
2719 | } else { |
---|
2720 | #if 1 |
---|
2721 | // odd (could be free) - it's feasible - go to nearest |
---|
2722 | if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) { |
---|
2723 | directionOut_ = 1; |
---|
2724 | dualOut_ = lowerOut_ - valueOut_; |
---|
2725 | } else { |
---|
2726 | directionOut_ = -1; |
---|
2727 | dualOut_ = valueOut_ - upperOut_; |
---|
2728 | } |
---|
2729 | #else |
---|
2730 | // odd (could be free) - it's feasible - improve obj |
---|
2731 | printf("direction from alpha of %g is %d\n", |
---|
2732 | freeAlpha, freeAlpha > 0.0 ? 1 : -1); |
---|
2733 | if (valueOut_ - lowerOut_ > 1.0e20) |
---|
2734 | freeAlpha = 1.0; |
---|
2735 | else if(upperOut_ - valueOut_ > 1.0e20) |
---|
2736 | freeAlpha = -1.0; |
---|
2737 | //if (valueOut_-lowerOut_<upperOut_-valueOut_) { |
---|
2738 | if (freeAlpha < 0.0) { |
---|
2739 | directionOut_ = 1; |
---|
2740 | dualOut_ = lowerOut_ - valueOut_; |
---|
2741 | } else { |
---|
2742 | directionOut_ = -1; |
---|
2743 | dualOut_ = valueOut_ - upperOut_; |
---|
2744 | } |
---|
2745 | printf("direction taken %d - bounds %g %g %g\n", |
---|
2746 | directionOut_, lowerOut_, valueOut_, upperOut_); |
---|
2747 | #endif |
---|
2748 | } |
---|
2749 | #ifdef CLP_DEBUG |
---|
2750 | assert(dualOut_ >= 0.0); |
---|
2751 | #endif |
---|
2752 | } else { |
---|
2753 | // in values pass so just use sign of dj |
---|
2754 | // We don't want to go through any barriers so set dualOut low |
---|
2755 | // free variables will never be here |
---|
2756 | dualOut_ = 1.0e-6; |
---|
2757 | if (dj_[sequenceOut_] > 0.0) { |
---|
2758 | // this will give a -1 in pivot row (as slacks are -1.0) |
---|
2759 | directionOut_ = 1; |
---|
2760 | } else { |
---|
2761 | directionOut_ = -1; |
---|
2762 | } |
---|
2763 | } |
---|
2764 | } |
---|
2765 | return ; |
---|
2766 | } |
---|
2767 | // Checks if any fake bounds active - if so returns number and modifies |
---|
2768 | // dualBound_ and everything. |
---|
2769 | // Free variables will be left as free |
---|
2770 | // Returns number of bounds changed if >=0 |
---|
2771 | // Returns -1 if not initialize and no effect |
---|
2772 | // Fills in changeVector which can be used to see if unbounded |
---|
2773 | // and cost of change vector |
---|
2774 | int |
---|
2775 | ClpSimplexDual::changeBounds(int initialize, |
---|
2776 | CoinIndexedVector * outputArray, |
---|
2777 | double & changeCost) |
---|
2778 | { |
---|
2779 | numberFake_ = 0; |
---|
2780 | if (!initialize) { |
---|
2781 | int numberInfeasibilities; |
---|
2782 | double newBound; |
---|
2783 | newBound = 5.0 * dualBound_; |
---|
2784 | numberInfeasibilities = 0; |
---|
2785 | changeCost = 0.0; |
---|
2786 | // put back original bounds and then check |
---|
2787 | createRim1(false); |
---|
2788 | int iSequence; |
---|
2789 | // bounds will get bigger - just look at ones at bounds |
---|
2790 | for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { |
---|
2791 | double lowerValue = lower_[iSequence]; |
---|
2792 | double upperValue = upper_[iSequence]; |
---|
2793 | double value = solution_[iSequence]; |
---|
2794 | setFakeBound(iSequence, ClpSimplexDual::noFake); |
---|
2795 | switch(getStatus(iSequence)) { |
---|
2796 | |
---|
2797 | case basic: |
---|
2798 | case ClpSimplex::isFixed: |
---|
2799 | break; |
---|
2800 | case isFree: |
---|
2801 | case superBasic: |
---|
2802 | break; |
---|
2803 | case atUpperBound: |
---|
2804 | if (fabs(value - upperValue) > primalTolerance_) |
---|
2805 | numberInfeasibilities++; |
---|
2806 | break; |
---|
2807 | case atLowerBound: |
---|
2808 | if (fabs(value - lowerValue) > primalTolerance_) |
---|
2809 | numberInfeasibilities++; |
---|
2810 | break; |
---|
2811 | } |
---|
2812 | } |
---|
2813 | // If dual infeasible then carry on |
---|
2814 | if (numberInfeasibilities) { |
---|
2815 | handler_->message(CLP_DUAL_CHECKB, messages_) |
---|
2816 | << newBound |
---|
2817 | << CoinMessageEol; |
---|
2818 | int iSequence; |
---|
2819 | for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { |
---|
2820 | double lowerValue = lower_[iSequence]; |
---|
2821 | double upperValue = upper_[iSequence]; |
---|
2822 | double newLowerValue; |
---|
2823 | double newUpperValue; |
---|
2824 | Status status = getStatus(iSequence); |
---|
2825 | if (status == atUpperBound || |
---|
2826 | status == atLowerBound) { |
---|
2827 | double value = solution_[iSequence]; |
---|
2828 | if (value - lowerValue <= upperValue - value) { |
---|
2829 | newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound); |
---|
2830 | newUpperValue = CoinMin(upperValue, newLowerValue + newBound); |
---|
2831 | } else { |
---|
2832 | newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound); |
---|
2833 | newLowerValue = CoinMax(lowerValue, newUpperValue - newBound); |
---|
2834 | } |
---|
2835 | lower_[iSequence] = newLowerValue; |
---|
2836 | upper_[iSequence] = newUpperValue; |
---|
2837 | if (newLowerValue > lowerValue) { |
---|
2838 | if (newUpperValue < upperValue) { |
---|
2839 | setFakeBound(iSequence, ClpSimplexDual::bothFake); |
---|
2840 | #ifdef CLP_INVESTIGATE |
---|
2841 | abort(); // No idea what should happen here - I have never got here |
---|
2842 | #endif |
---|
2843 | numberFake_++; |
---|
2844 | } else { |
---|
2845 | setFakeBound(iSequence, ClpSimplexDual::lowerFake); |
---|
2846 | numberFake_++; |
---|
2847 | } |
---|
2848 | } else { |
---|
2849 | if (newUpperValue < upperValue) { |
---|
2850 | setFakeBound(iSequence, ClpSimplexDual::upperFake); |
---|
2851 | numberFake_++; |
---|
2852 | } |
---|
2853 | } |
---|
2854 | if (status == atUpperBound) |
---|
2855 | solution_[iSequence] = newUpperValue; |
---|
2856 | else |
---|
2857 | solution_[iSequence] = newLowerValue; |
---|
2858 | double movement = solution_[iSequence] - value; |
---|
2859 | if (movement && outputArray) { |
---|
2860 | if (iSequence >= numberColumns_) { |
---|
2861 | outputArray->quickAdd(iSequence, -movement); |
---|
2862 | changeCost += movement * cost_[iSequence]; |
---|
2863 | } else { |
---|
2864 | matrix_->add(this, outputArray, iSequence, movement); |
---|
2865 | changeCost += movement * cost_[iSequence]; |
---|
2866 | } |
---|
2867 | } |
---|
2868 | } |
---|
2869 | } |
---|
2870 | dualBound_ = newBound; |
---|
2871 | } else { |
---|
2872 | numberInfeasibilities = -1; |
---|
2873 | } |
---|
2874 | return numberInfeasibilities; |
---|
2875 | } else if (initialize == 1 || initialize == 3) { |
---|
2876 | int iSequence; |
---|
2877 | if (initialize == 3) { |
---|
2878 | for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { |
---|
2879 | setFakeBound(iSequence, ClpSimplexDual::noFake); |
---|
2880 | } |
---|
2881 | } |
---|
2882 | double testBound = 0.999999 * dualBound_; |
---|
2883 | for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { |
---|
2884 | Status status = getStatus(iSequence); |
---|
2885 | if (status == atUpperBound || |
---|
2886 | status == atLowerBound) { |
---|
2887 | double lowerValue = lower_[iSequence]; |
---|
2888 | double upperValue = upper_[iSequence]; |
---|
2889 | double value = solution_[iSequence]; |
---|
2890 | if (lowerValue > -largeValue_ || upperValue < largeValue_) { |
---|
2891 | if (true || lowerValue - value > -0.5 * dualBound_ || |
---|
2892 | upperValue - value < 0.5 * dualBound_) { |
---|
2893 | if (fabs(lowerValue - value) <= fabs(upperValue - value)) { |
---|
2894 | if (upperValue > lowerValue + testBound) { |
---|
2895 | if (getFakeBound(iSequence) == ClpSimplexDual::noFake) |
---|
2896 | numberFake_++; |
---|
2897 | upper_[iSequence] = lowerValue + dualBound_; |
---|
2898 | setFakeBound(iSequence, ClpSimplexDual::upperFake); |
---|
2899 | } |
---|
2900 | } else { |
---|
2901 | if (lowerValue < upperValue - testBound) { |
---|
2902 | if (getFakeBound(iSequence) == ClpSimplexDual::noFake) |
---|
2903 | numberFake_++; |
---|
2904 | lower_[iSequence] = upperValue - dualBound_; |
---|
2905 | setFakeBound(iSequence, ClpSimplexDual::lowerFake); |
---|
2906 | } |
---|
2907 | } |
---|
2908 | } else { |
---|
2909 | if (getFakeBound(iSequence) == ClpSimplexDual::noFake) |
---|
2910 | numberFake_++; |
---|
2911 | lower_[iSequence] = -0.5 * dualBound_; |
---|
2912 | upper_[iSequence] = 0.5 * dualBound_; |
---|
2913 | setFakeBound(iSequence, ClpSimplexDual::bothFake); |
---|
2914 | abort(); |
---|
2915 | } |
---|
2916 | if (status == atUpperBound) |
---|
2917 | solution_[iSequence] = upper_[iSequence]; |
---|
2918 | else |
---|
2919 | solution_[iSequence] = lower_[iSequence]; |
---|
2920 | } else { |
---|
2921 | // set non basic free variables to fake bounds |
---|
2922 | // I don't think we should ever get here |
---|
2923 | CoinAssert(!("should not be here")); |
---|
2924 | lower_[iSequence] = -0.5 * dualBound_; |
---|
2925 | upper_[iSequence] = 0.5 * dualBound_; |
---|
2926 | setFakeBound(iSequence, ClpSimplexDual::bothFake); |
---|
2927 | numberFake_++; |
---|
2928 | setStatus(iSequence, atUpperBound); |
---|
2929 | solution_[iSequence] = 0.5 * dualBound_; |
---|
2930 | } |
---|
2931 | } else if (status == basic) { |
---|
2932 | // make sure not at fake bound and bounds correct |
---|
2933 | setFakeBound(iSequence, ClpSimplexDual::noFake); |
---|
2934 | double gap = upper_[iSequence] - lower_[iSequence]; |
---|
2935 | if (gap > 0.5 * dualBound_ && gap < 2.0 * dualBound_) { |
---|
2936 | if (iSequence < numberColumns_) { |
---|
2937 | if (columnScale_) { |
---|
2938 | double multiplier = rhsScale_ * inverseColumnScale_[iSequence]; |
---|
2939 | // lower |
---|
2940 | double value = columnLower_[iSequence]; |
---|
2941 | if (value > -1.0e30) { |
---|
2942 | value *= multiplier; |
---|
2943 | } |
---|
2944 | lower_[iSequence] = value; |
---|
2945 | // upper |
---|
2946 | value = columnUpper_[iSequence]; |
---|
2947 | if (value < 1.0e30) { |
---|
2948 | value *= multiplier; |
---|
2949 | } |
---|
2950 | upper_[iSequence] = value; |
---|
2951 | } else { |
---|
2952 | lower_[iSequence] = columnLower_[iSequence];; |
---|
2953 | upper_[iSequence] = columnUpper_[iSequence];; |
---|
2954 | } |
---|
2955 | } else { |
---|
2956 | int iRow = iSequence - numberColumns_; |
---|
2957 | if (rowScale_) { |
---|
2958 | // lower |
---|
2959 | double multiplier = rhsScale_ * rowScale_[iRow]; |
---|
2960 | double value = rowLower_[iRow]; |
---|
2961 | if (value > -1.0e30) { |
---|
2962 | value *= multiplier; |
---|
2963 | } |
---|
2964 | lower_[iSequence] = value; |
---|
2965 | // upper |
---|
2966 | value = rowUpper_[iRow]; |
---|
2967 | if (value < 1.0e30) { |
---|
2968 | value *= multiplier; |
---|
2969 | } |
---|
2970 | upper_[iSequence] = value; |
---|
2971 | } else { |
---|
2972 | lower_[iSequence] = rowLower_[iRow];; |
---|
2973 | upper_[iSequence] = rowUpper_[iRow];; |
---|
2974 | } |
---|
2975 | } |
---|
2976 | } |
---|
2977 | } |
---|
2978 | } |
---|
2979 | |
---|
2980 | return 1; |
---|
2981 | } else { |
---|
2982 | // just reset changed ones |
---|
2983 | if (columnScale_) { |
---|
2984 | int iSequence; |
---|
2985 | for (iSequence = 0; iSequence < numberColumns_; iSequence++) { |
---|
2986 | FakeBound fakeStatus = getFakeBound(iSequence); |
---|
2987 | if (fakeStatus != noFake) { |
---|
2988 | if ((static_cast<int> (fakeStatus) & 1) != 0) { |
---|
2989 | // lower |
---|
2990 | double value = columnLower_[iSequence]; |
---|
2991 | if (value > -1.0e30) { |
---|
2992 | double multiplier = rhsScale_ * inverseColumnScale_[iSequence]; |
---|
2993 | value *= multiplier; |
---|
2994 | } |
---|
2995 | columnLowerWork_[iSequence] = value; |
---|
2996 | } |
---|
2997 | if ((static_cast<int> (fakeStatus) & 2) != 0) { |
---|
2998 | // upper |
---|
2999 | double value = columnUpper_[iSequence]; |
---|
3000 | if (value < 1.0e30) { |
---|
3001 | double multiplier = rhsScale_ * inverseColumnScale_[iSequence]; |
---|
3002 | value *= multiplier; |
---|
3003 | } |
---|
3004 | columnUpperWork_[iSequence] = value; |
---|
3005 | } |
---|
3006 | } |
---|
3007 | } |
---|
3008 | for (iSequence = 0; iSequence < numberRows_; iSequence++) { |
---|
3009 | FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_); |
---|
3010 | if (fakeStatus != noFake) { |
---|
3011 | if ((static_cast<int> (fakeStatus) & 1) != 0) { |
---|
3012 | // lower |
---|
3013 | double value = rowLower_[iSequence]; |
---|
3014 | if (value > -1.0e30) { |
---|
3015 | double multiplier = rhsScale_ * rowScale_[iSequence]; |
---|
3016 | value *= multiplier; |
---|
3017 | } |
---|
3018 | rowLowerWork_[iSequence] = value; |
---|
3019 | } |
---|
3020 | if ((static_cast<int> (fakeStatus) & 2) != 0) { |
---|
3021 | // upper |
---|
3022 | double value = rowUpper_[iSequence]; |
---|
3023 | if (value < 1.0e30) { |
---|
3024 | double multiplier = rhsScale_ * rowScale_[iSequence]; |
---|
3025 | value *= multiplier; |
---|
3026 | } |
---|
3027 | rowUpperWork_[iSequence] = value; |
---|
3028 | } |
---|
3029 | } |
---|
3030 | } |
---|
3031 | } else { |
---|
3032 | int iSequence; |
---|
3033 | for (iSequence = 0; iSequence < numberColumns_; iSequence++) { |
---|
3034 | FakeBound fakeStatus = getFakeBound(iSequence); |
---|
3035 | if ((static_cast<int> (fakeStatus) & 1) != 0) { |
---|
3036 | // lower |
---|
3037 | columnLowerWork_[iSequence] = columnLower_[iSequence]; |
---|
3038 | } |
---|
3039 | if ((static_cast<int> (fakeStatus) & 2) != 0) { |
---|
3040 | // upper |
---|
3041 | columnUpperWork_[iSequence] = columnUpper_[iSequence]; |
---|
3042 | } |
---|
3043 | } |
---|
3044 | for (iSequence = 0; iSequence < numberRows_; iSequence++) { |
---|
3045 | FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_); |
---|
3046 | if ((static_cast<int> (fakeStatus) & 1) != 0) { |
---|
3047 | // lower |
---|
3048 | rowLowerWork_[iSequence] = rowLower_[iSequence]; |
---|
3049 | } |
---|
3050 | if ((static_cast<int> (fakeStatus) & 2) != 0) { |
---|
3051 | // upper |
---|
3052 | rowUpperWork_[iSequence] = rowUpper_[iSequence]; |
---|
3053 | } |
---|
3054 | } |
---|
3055 | } |
---|
3056 | return 0; |
---|
3057 | } |
---|
3058 | } |
---|
3059 | int |
---|
3060 | ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray, |
---|
3061 | const CoinIndexedVector * columnArray, |
---|
3062 | CoinIndexedVector * spareArray, |
---|
3063 | double acceptablePivot, |
---|
3064 | double & upperReturn, double &bestReturn, double & badFree) |
---|
3065 | { |
---|
3066 | // do first pass to get possibles |
---|
3067 | double * spare = spareArray->denseVector(); |
---|
3068 | int * index = spareArray->getIndices(); |
---|
3069 | const double * work; |
---|
3070 | int number; |
---|
3071 | const int * which; |
---|
3072 | const double * reducedCost; |
---|
3073 | // We can also see if infeasible or pivoting on free |
---|
3074 | double tentativeTheta = 1.0e15; |
---|
3075 | double upperTheta = 1.0e31; |
---|
3076 | double freePivot = acceptablePivot; |
---|
3077 | double bestPossible = 0.0; |
---|
3078 | int numberRemaining = 0; |
---|
3079 | int i; |
---|
3080 | badFree = 0.0; |
---|
3081 | if ((moreSpecialOptions_ & 8) != 0) { |
---|
3082 | // No free or super basic |
---|
3083 | double multiplier[] = { -1.0, 1.0}; |
---|
3084 | double dualT = - dualTolerance_; |
---|
3085 | for (int iSection = 0; iSection < 2; iSection++) { |
---|
3086 | |
---|
3087 | int addSequence; |
---|
3088 | unsigned char * statusArray; |
---|
3089 | if (!iSection) { |
---|
3090 | work = rowArray->denseVector(); |
---|
3091 | number = rowArray->getNumElements(); |
---|
3092 | which = rowArray->getIndices(); |
---|
3093 | reducedCost = rowReducedCost_; |
---|
3094 | addSequence = numberColumns_; |
---|
3095 | statusArray = status_ + numberColumns_; |
---|
3096 | } else { |
---|
3097 | work = columnArray->denseVector(); |
---|
3098 | number = columnArray->getNumElements(); |
---|
3099 | which = columnArray->getIndices(); |
---|
3100 | reducedCost = reducedCostWork_; |
---|
3101 | addSequence = 0; |
---|
3102 | statusArray = status_; |
---|
3103 | } |
---|
3104 | |
---|
3105 | for (i = 0; i < number; i++) { |
---|
3106 | int iSequence = which[i]; |
---|
3107 | double alpha; |
---|
3108 | double oldValue; |
---|
3109 | double value; |
---|
3110 | |
---|
3111 | assert (getStatus(iSequence + addSequence) != isFree |
---|
3112 | && getStatus(iSequence + addSequence) != superBasic); |
---|
3113 | int iStatus = (statusArray[iSequence] & 3) - 1; |
---|
3114 | if (iStatus) { |
---|
3115 | double mult = multiplier[iStatus-1]; |
---|
3116 | alpha = work[i] * mult; |
---|
3117 | if (alpha > 0.0) { |
---|
3118 | oldValue = reducedCost[iSequence] * mult; |
---|
3119 | value = oldValue - tentativeTheta * alpha; |
---|
3120 | if (value < dualT) { |
---|
3121 | bestPossible = CoinMax(bestPossible, alpha); |
---|
3122 | value = oldValue - upperTheta * alpha; |
---|
3123 | if (value < dualT && alpha >= acceptablePivot) { |
---|
3124 | upperTheta = (oldValue - dualT) / alpha; |
---|
3125 | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
---|
3126 | } |
---|
3127 | // add to list |
---|
3128 | spare[numberRemaining] = alpha * mult; |
---|
3129 | index[numberRemaining++] = iSequence + addSequence; |
---|
3130 | } |
---|
3131 | } |
---|
3132 | } |
---|
3133 | } |
---|
3134 | } |
---|
3135 | } else { |
---|
3136 | // some free or super basic |
---|
3137 | for (int iSection = 0; iSection < 2; iSection++) { |
---|
3138 | |
---|
3139 | int addSequence; |
---|
3140 | |
---|
3141 | if (!iSection) { |
---|
3142 | work = rowArray->denseVector(); |
---|
3143 | number = rowArray->getNumElements(); |
---|
3144 | which = rowArray->getIndices(); |
---|
3145 | reducedCost = rowReducedCost_; |
---|
3146 | addSequence = numberColumns_; |
---|
3147 | } else { |
---|
3148 | work = columnArray->denseVector(); |
---|
3149 | number = columnArray->getNumElements(); |
---|
3150 | which = columnArray->getIndices(); |
---|
3151 | reducedCost = reducedCostWork_; |
---|
3152 | addSequence = 0; |
---|
3153 | } |
---|
3154 | |
---|
3155 | for (i = 0; i < number; i++) { |
---|
3156 | int iSequence = which[i]; |
---|
3157 | double alpha; |
---|
3158 | double oldValue; |
---|
3159 | double value; |
---|
3160 | bool keep; |
---|
3161 | |
---|
3162 | switch(getStatus(iSequence + addSequence)) { |
---|
3163 | |
---|
3164 | case basic: |
---|
3165 | case ClpSimplex::isFixed: |
---|
3166 | break; |
---|
3167 | case isFree: |
---|
3168 | case superBasic: |
---|
3169 | alpha = work[i]; |
---|
3170 | bestPossible = CoinMax(bestPossible, fabs(alpha)); |
---|
3171 | oldValue = reducedCost[iSequence]; |
---|
3172 | // If free has to be very large - should come in via dualRow |
---|
3173 | //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3) |
---|
3174 | //break; |
---|
3175 | if (oldValue > dualTolerance_) { |
---|
3176 | keep = true; |
---|
3177 | } else if (oldValue < -dualTolerance_) { |
---|
3178 | keep = true; |
---|
3179 | } else { |
---|
3180 | if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) { |
---|
3181 | keep = true; |
---|
3182 | } else { |
---|
3183 | keep = false; |
---|
3184 | badFree = CoinMax(badFree, fabs(alpha)); |
---|
3185 | } |
---|
3186 | } |
---|
3187 | if (keep) { |
---|
3188 | // free - choose largest |
---|
3189 | if (fabs(alpha) > freePivot) { |
---|
3190 | freePivot = fabs(alpha); |
---|
3191 | sequenceIn_ = iSequence + addSequence; |
---|
3192 | theta_ = oldValue / alpha; |
---|
3193 | alpha_ = alpha; |
---|
3194 | } |
---|
3195 | } |
---|
3196 | break; |
---|
3197 | case atUpperBound: |
---|
3198 | alpha = work[i]; |
---|
3199 | oldValue = reducedCost[iSequence]; |
---|
3200 | value = oldValue - tentativeTheta * alpha; |
---|
3201 | //assert (oldValue<=dualTolerance_*1.0001); |
---|
3202 | if (value > dualTolerance_) { |
---|
3203 | bestPossible = CoinMax(bestPossible, -alpha); |
---|
3204 | value = oldValue - upperTheta * alpha; |
---|
3205 | if (value > dualTolerance_ && -alpha >= acceptablePivot) { |
---|
3206 | upperTheta = (oldValue - dualTolerance_) / alpha; |
---|
3207 | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
---|
3208 | } |
---|
3209 | // add to list |
---|
3210 | spare[numberRemaining] = alpha; |
---|
3211 | index[numberRemaining++] = iSequence + addSequence; |
---|
3212 | } |
---|
3213 | break; |
---|
3214 | case atLowerBound: |
---|
3215 | alpha = work[i]; |
---|
3216 | oldValue = reducedCost[iSequence]; |
---|
3217 | value = oldValue - tentativeTheta * alpha; |
---|
3218 | //assert (oldValue>=-dualTolerance_*1.0001); |
---|
3219 | if (value < -dualTolerance_) { |
---|
3220 | bestPossible = CoinMax(bestPossible, alpha); |
---|
3221 | value = oldValue - upperTheta * alpha; |
---|
3222 | if (value < -dualTolerance_ && alpha >= acceptablePivot) { |
---|
3223 | upperTheta = (oldValue + dualTolerance_) / alpha; |
---|
3224 | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
---|
3225 | } |
---|
3226 | // add to list |
---|
3227 | spare[numberRemaining] = alpha; |
---|
3228 | index[numberRemaining++] = iSequence + addSequence; |
---|
3229 | } |
---|
3230 | break; |
---|
3231 | } |
---|
3232 | } |
---|
3233 | } |
---|
3234 | } |
---|
3235 | upperReturn = upperTheta; |
---|
3236 | bestReturn = bestPossible; |
---|
3237 | return numberRemaining; |
---|
3238 | } |
---|
3239 | /* |
---|
3240 | Row array has row part of pivot row (as duals so sign may be switched) |
---|
3241 | Column array has column part. |
---|
3242 | This chooses pivot column. |
---|
3243 | Spare array will be needed when we start getting clever. |
---|
3244 | We will check for basic so spare array will never overflow. |
---|
3245 | If necessary will modify costs |
---|
3246 | */ |
---|
3247 | double |
---|
3248 | ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray, |
---|
3249 | CoinIndexedVector * columnArray, |
---|
3250 | CoinIndexedVector * spareArray, |
---|
3251 | CoinIndexedVector * spareArray2, |
---|
3252 | double acceptablePivot, |
---|
3253 | CoinBigIndex * /*dubiousWeights*/) |
---|
3254 | { |
---|
3255 | int numberPossiblySwapped = 0; |
---|
3256 | int numberRemaining = 0; |
---|
3257 | |
---|
3258 | double totalThru = 0.0; // for when variables flip |
---|
3259 | //double saveAcceptable=acceptablePivot; |
---|
3260 | //acceptablePivot=1.0e-9; |
---|
3261 | |
---|
3262 | double bestEverPivot = acceptablePivot; |
---|
3263 | int lastSequence = -1; |
---|
3264 | double lastPivot = 0.0; |
---|
3265 | double upperTheta; |
---|
3266 | double newTolerance = dualTolerance_; |
---|
3267 | //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance]; |
---|
3268 | // will we need to increase tolerance |
---|
3269 | bool thisIncrease = false; |
---|
3270 | // If we think we need to modify costs (not if something from broad sweep) |
---|
3271 | bool modifyCosts = false; |
---|
3272 | // Increase in objective due to swapping bounds (may be negative) |
---|
3273 | double increaseInObjective = 0.0; |
---|
3274 | |
---|
3275 | // use spareArrays to put ones looked at in |
---|
3276 | // we are going to flip flop between |
---|
3277 | int iFlip = 0; |
---|
3278 | // Possible list of pivots |
---|
3279 | int interesting[2]; |
---|
3280 | // where possible swapped ones are |
---|
3281 | int swapped[2]; |
---|
3282 | // for zeroing out arrays after |
---|
3283 | int marker[2][2]; |
---|
3284 | // pivot elements |
---|
3285 | double * array[2], * spare, * spare2; |
---|
3286 | // indices |
---|
3287 | int * indices[2], * index, * index2; |
---|
3288 | spareArray2->clear(); |
---|
3289 | array[0] = spareArray->denseVector(); |
---|
3290 | indices[0] = spareArray->getIndices(); |
---|
3291 | spare = array[0]; |
---|
3292 | index = indices[0]; |
---|
3293 | array[1] = spareArray2->denseVector(); |
---|
3294 | indices[1] = spareArray2->getIndices(); |
---|
3295 | int i; |
---|
3296 | |
---|
3297 | // initialize lists |
---|
3298 | for (i = 0; i < 2; i++) { |
---|
3299 | interesting[i] = 0; |
---|
3300 | swapped[i] = numberColumns_; |
---|
3301 | marker[i][0] = 0; |
---|
3302 | marker[i][1] = numberColumns_; |
---|
3303 | } |
---|
3304 | /* |
---|
3305 | First we get a list of possible pivots. We can also see if the |
---|
3306 | problem looks infeasible or whether we want to pivot in free variable. |
---|
3307 | This may make objective go backwards but can only happen a finite |
---|
3308 | number of times and I do want free variables basic. |
---|
3309 | |
---|
3310 | Then we flip back and forth. At the start of each iteration |
---|
3311 | interesting[iFlip] should have possible candidates and swapped[iFlip] |
---|
3312 | will have pivots if we decide to take a previous pivot. |
---|
3313 | At end of each iteration interesting[1-iFlip] should have |
---|
3314 | candidates if we go through this theta and swapped[1-iFlip] |
---|
3315 | pivots if we don't go through. |
---|
3316 | |
---|
3317 | At first we increase theta and see what happens. We start |
---|
3318 | theta at a reasonable guess. If in right area then we do bit by bit. |
---|
3319 | |
---|
3320 | */ |
---|
3321 | |
---|
3322 | // do first pass to get possibles |
---|
3323 | upperTheta = 1.0e31; |
---|
3324 | double bestPossible = 0.0; |
---|
3325 | double badFree = 0.0; |
---|
3326 | alpha_ = 0.0; |
---|
3327 | if (spareIntArray_[0] >= 0) { |
---|
3328 | numberRemaining = dualColumn0(rowArray, columnArray, spareArray, |
---|
3329 | acceptablePivot, upperTheta, bestPossible, badFree); |
---|
3330 | } else { |
---|
3331 | // already done |
---|
3332 | numberRemaining = spareArray->getNumElements(); |
---|
3333 | spareArray->setNumElements(0); |
---|
3334 | upperTheta = spareDoubleArray_[0]; |
---|
3335 | bestPossible = spareDoubleArray_[1]; |
---|
3336 | if (spareIntArray_[0] == -1) { |
---|
3337 | theta_ = spareDoubleArray_[2]; |
---|
3338 | alpha_ = spareDoubleArray_[3]; |
---|
3339 | sequenceIn_ = spareIntArray_[1]; |
---|
3340 | } else { |
---|
3341 | #if 0 |
---|
3342 | int n = numberRemaining; |
---|
3343 | double u = upperTheta; |
---|
3344 | double b = bestPossible; |
---|
3345 | upperTheta = 1.0e31; |
---|
3346 | bestPossible = 0.0; |
---|
3347 | numberRemaining = dualColumn0(rowArray, columnArray, spareArray, |
---|
3348 | acceptablePivot, upperTheta, bestPossible, badFree); |
---|
3349 | assert (n == numberRemaining); |
---|
3350 | assert (fabs(b - bestPossible) < 1.0e-7); |
---|
3351 | assert (fabs(u - upperTheta) < 1.0e-7); |
---|
3352 | #endif |
---|
3353 | } |
---|
3354 | } |
---|
3355 | // switch off |
---|
3356 | spareIntArray_[0] = 0; |
---|
3357 | // We can also see if infeasible or pivoting on free |
---|
3358 | double tentativeTheta = 1.0e25; |
---|
3359 | interesting[0] = numberRemaining; |
---|
3360 | marker[0][0] = numberRemaining; |
---|
3361 | |
---|
3362 | if (!numberRemaining && sequenceIn_ < 0) |
---|
3363 | return 0.0; // Looks infeasible |
---|
3364 | |
---|
3365 | // If sum of bad small pivots too much |
---|
3366 | #define MORE_CAREFUL |
---|
3367 | #ifdef MORE_CAREFUL |
---|
3368 | bool badSumPivots = false; |
---|
3369 | #endif |
---|
3370 | if (sequenceIn_ >= 0) { |
---|
3371 | // free variable - always choose |
---|
3372 | } else { |
---|
3373 | |
---|
3374 | theta_ = 1.0e50; |
---|
3375 | // now flip flop between spare arrays until reasonable theta |
---|
3376 | tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7); |
---|
3377 | |
---|
3378 | // loops increasing tentative theta until can't go through |
---|
3379 | |
---|
3380 | while (tentativeTheta < 1.0e22) { |
---|
3381 | double thruThis = 0.0; |
---|
3382 | |
---|
3383 | double bestPivot = acceptablePivot; |
---|
3384 | int bestSequence = -1; |
---|
3385 | |
---|
3386 | numberPossiblySwapped = numberColumns_; |
---|
3387 | numberRemaining = 0; |
---|
3388 | |
---|
3389 | upperTheta = 1.0e50; |
---|
3390 | |
---|
3391 | spare = array[iFlip]; |
---|
3392 | index = indices[iFlip]; |
---|
3393 | spare2 = array[1-iFlip]; |
---|
3394 | index2 = indices[1-iFlip]; |
---|
3395 | |
---|
3396 | // try 3 different ways |
---|
3397 | // 1 bias increase by ones with slightly wrong djs |
---|
3398 | // 2 bias by all |
---|
3399 | // 3 bias by all - tolerance |
---|
3400 | #define TRYBIAS 3 |
---|
3401 | |
---|
3402 | |
---|
3403 | double increaseInThis = 0.0; //objective increase in this loop |
---|
3404 | |
---|
3405 | for (i = 0; i < interesting[iFlip]; i++) { |
---|
3406 | int iSequence = index[i]; |
---|
3407 | double alpha = spare[i]; |
---|
3408 | double oldValue = dj_[iSequence]; |
---|
3409 | double value = oldValue - tentativeTheta * alpha; |
---|
3410 | |
---|
3411 | if (alpha < 0.0) { |
---|
3412 | //at upper bound |
---|
3413 | if (value > newTolerance) { |
---|
3414 | double range = upper_[iSequence] - lower_[iSequence]; |
---|
3415 | thruThis -= range * alpha; |
---|
3416 | #if TRYBIAS==1 |
---|
3417 | if (oldValue > 0.0) |
---|
3418 | increaseInThis -= oldValue * range; |
---|
3419 | #elif TRYBIAS==2 |
---|
3420 | increaseInThis -= oldValue * range; |
---|
3421 | #else |
---|
3422 | increaseInThis -= (oldValue + dualTolerance_) * range; |
---|
3423 | #endif |
---|
3424 | // goes on swapped list (also means candidates if too many) |
---|
3425 | spare2[--numberPossiblySwapped] = alpha; |
---|
3426 | index2[numberPossiblySwapped] = iSequence; |
---|
3427 | if (fabs(alpha) > bestPivot) { |
---|
3428 | bestPivot = fabs(alpha); |
---|
3429 | bestSequence = numberPossiblySwapped; |
---|
3430 | } |
---|
3431 | } else { |
---|
3432 | value = oldValue - upperTheta * alpha; |
---|
3433 | if (value > newTolerance && -alpha >= acceptablePivot) |
---|
3434 | upperTheta = (oldValue - newTolerance) / alpha; |
---|
3435 | spare2[numberRemaining] = alpha; |
---|
3436 | index2[numberRemaining++] = iSequence; |
---|
3437 | } |
---|
3438 | } else { |
---|
3439 | // at lower bound |
---|
3440 | if (value < -newTolerance) { |
---|
3441 | double range = upper_[iSequence] - lower_[iSequence]; |
---|
3442 | thruThis += range * alpha; |
---|
3443 | //?? is this correct - and should we look at good ones |
---|
3444 | #if TRYBIAS==1 |
---|
3445 | if (oldValue < 0.0) |
---|
3446 | increaseInThis += oldValue * range; |
---|
3447 | #elif TRYBIAS==2 |
---|
3448 | increaseInThis += oldValue * range; |
---|
3449 | #else |
---|
3450 | increaseInThis += (oldValue - dualTolerance_) * range; |
---|
3451 | #endif |
---|
3452 | // goes on swapped list (also means candidates if too many) |
---|
3453 | spare2[--numberPossiblySwapped] = alpha; |
---|
3454 | index2[numberPossiblySwapped] = iSequence; |
---|
3455 | if (fabs(alpha) > bestPivot) { |
---|
3456 | bestPivot = fabs(alpha); |
---|
3457 | bestSequence = numberPossiblySwapped; |
---|
3458 | } |
---|
3459 | } else { |
---|
3460 | value = oldValue - upperTheta * alpha; |
---|
3461 | if (value < -newTolerance && alpha >= acceptablePivot) |
---|
3462 | upperTheta = (oldValue + newTolerance) / alpha; |
---|
3463 | spare2[numberRemaining] = alpha; |
---|
3464 | index2[numberRemaining++] = iSequence; |
---|
3465 | } |
---|
3466 | } |
---|
3467 | } |
---|
3468 | swapped[1-iFlip] = numberPossiblySwapped; |
---|
3469 | interesting[1-iFlip] = numberRemaining; |
---|
3470 | marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining); |
---|
3471 | marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped); |
---|
3472 | |
---|
3473 | if (totalThru + thruThis >= fabs(dualOut_) || |
---|
3474 | increaseInObjective + increaseInThis < 0.0) { |
---|
3475 | // We should be pivoting in this batch |
---|
3476 | // so compress down to this lot |
---|
3477 | numberRemaining = 0; |
---|
3478 | for (i = numberColumns_ - 1; i >= swapped[1-iFlip]; i--) { |
---|
3479 | spare[numberRemaining] = spare2[i]; |
---|
3480 | index[numberRemaining++] = index2[i]; |
---|
3481 | } |
---|
3482 | interesting[iFlip] = numberRemaining; |
---|
3483 | int iTry; |
---|
3484 | #define MAXTRY 100 |
---|
3485 | // first get ratio with tolerance |
---|
3486 | for (iTry = 0; iTry < MAXTRY; iTry++) { |
---|
3487 | |
---|
3488 | upperTheta = 1.0e50; |
---|
3489 | numberPossiblySwapped = numberColumns_; |
---|
3490 | numberRemaining = 0; |
---|
3491 | |
---|
3492 | increaseInThis = 0.0; //objective increase in this loop |
---|
3493 | |
---|
3494 | thruThis = 0.0; |
---|
3495 | |
---|
3496 | spare = array[iFlip]; |
---|
3497 | index = indices[iFlip]; |
---|
3498 | spare2 = array[1-iFlip]; |
---|
3499 | index2 = indices[1-iFlip]; |
---|
3500 | for (i = 0; i < interesting[iFlip]; i++) { |
---|
3501 | int iSequence = index[i]; |
---|
3502 | double alpha = spare[i]; |
---|
3503 | double oldValue = dj_[iSequence]; |
---|
3504 | double value = oldValue - upperTheta * alpha; |
---|
3505 | |
---|
3506 | if (alpha < 0.0) { |
---|
3507 | //at upper bound |
---|
3508 | if (value > newTolerance) { |
---|
3509 | if (-alpha >= acceptablePivot) { |
---|
3510 | upperTheta = (oldValue - newTolerance) / alpha; |
---|
3511 | } |
---|
3512 | } |
---|
3513 | } else { |
---|
3514 | // at lower bound |
---|
3515 | if (value < -newTolerance) { |
---|
3516 | if (alpha >= acceptablePivot) { |
---|
3517 | upperTheta = (oldValue + newTolerance) / alpha; |
---|
3518 | } |
---|
3519 | } |
---|
3520 | } |
---|
3521 | } |
---|
3522 | bestPivot = acceptablePivot; |
---|
3523 | sequenceIn_ = -1; |
---|
3524 | #ifdef DUBIOUS_WEIGHTS |
---|
3525 | double bestWeight = COIN_DBL_MAX; |
---|
3526 | #endif |
---|
3527 | double largestPivot = acceptablePivot; |
---|
3528 | // now choose largest and sum all ones which will go through |
---|
3529 | //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]); |
---|
3530 | // Sum of bad small pivots |
---|
3531 | #ifdef MORE_CAREFUL |
---|
3532 | double sumBadPivots = 0.0; |
---|
3533 | badSumPivots = false; |
---|
3534 | #endif |
---|
3535 | // Make sure upperTheta will work (-O2 and above gives problems) |
---|
3536 | upperTheta *= 1.0000000001; |
---|
3537 | for (i = 0; i < interesting[iFlip]; i++) { |
---|
3538 | int iSequence = index[i]; |
---|
3539 | double alpha = spare[i]; |
---|
3540 | double value = dj_[iSequence] - upperTheta * alpha; |
---|
3541 | double badDj = 0.0; |
---|
3542 | |
---|
3543 | bool addToSwapped = false; |
---|
3544 | |
---|
3545 | if (alpha < 0.0) { |
---|
3546 | //at upper bound |
---|
3547 | if (value >= 0.0) { |
---|
3548 | addToSwapped = true; |
---|
3549 | #if TRYBIAS==1 |
---|
3550 | badDj = -CoinMax(dj_[iSequence], 0.0); |
---|
3551 | #elif TRYBIAS==2 |
---|
3552 | badDj = -dj_[iSequence]; |
---|
3553 | #else |
---|
3554 | badDj = -dj_[iSequence] - dualTolerance_; |
---|
3555 | #endif |
---|
3556 | } |
---|
3557 | } else { |
---|
3558 | // at lower bound |
---|
3559 | if (value <= 0.0) { |
---|
3560 | addToSwapped = true; |
---|
3561 | #if TRYBIAS==1 |
---|
3562 | badDj = CoinMin(dj_[iSequence], 0.0); |
---|
3563 | #elif TRYBIAS==2 |
---|
3564 | badDj = dj_[iSequence]; |
---|
3565 | #else |
---|
3566 | badDj = dj_[iSequence] - dualTolerance_; |
---|
3567 | #endif |
---|
3568 | } |
---|
3569 | } |
---|
3570 | if (!addToSwapped) { |
---|
3571 | // add to list of remaining |
---|
3572 | spare2[numberRemaining] = alpha; |
---|
3573 | index2[numberRemaining++] = iSequence; |
---|
3574 | } else { |
---|
3575 | // add to list of swapped |
---|
3576 | spare2[--numberPossiblySwapped] = alpha; |
---|
3577 | index2[numberPossiblySwapped] = iSequence; |
---|
3578 | // select if largest pivot |
---|
3579 | bool take = false; |
---|
3580 | double absAlpha = fabs(alpha); |
---|
3581 | #ifdef DUBIOUS_WEIGHTS |
---|
3582 | // User could do anything to break ties here |
---|
3583 | double weight; |
---|
3584 | if (dubiousWeights) |
---|
3585 | weight = dubiousWeights[iSequence]; |
---|
3586 | else |
---|
3587 | weight = 1.0; |
---|
3588 | weight += randomNumberGenerator_.randomDouble() * 1.0e-2; |
---|
3589 | if (absAlpha > 2.0 * bestPivot) { |
---|
3590 | take = true; |
---|
3591 | } else if (absAlpha > largestPivot) { |
---|
3592 | // could multiply absAlpha and weight |
---|
3593 | if (weight * bestPivot < bestWeight * absAlpha) |
---|
3594 | take = true; |
---|
3595 | } |
---|
3596 | #else |
---|
3597 | if (absAlpha > bestPivot) |
---|
3598 | take = true; |
---|
3599 | #endif |
---|
3600 | #ifdef MORE_CAREFUL |
---|
3601 | if (absAlpha < acceptablePivot && upperTheta < 1.0e20) { |
---|
3602 | if (alpha < 0.0) { |
---|
3603 | //at upper bound |
---|
3604 | if (value > dualTolerance_) { |
---|
3605 | double gap = upper_[iSequence] - lower_[iSequence]; |
---|
3606 | if (gap < 1.0e20) |
---|
3607 | sumBadPivots += value * gap; |
---|
3608 | else |
---|
3609 | sumBadPivots += 1.0e20; |
---|
3610 | //printf("bad %d alpha %g dj at upper %g\n", |
---|
3611 | // iSequence,alpha,value); |
---|
3612 | } |
---|
3613 | } else { |
---|
3614 | //at lower bound |
---|
3615 | if (value < -dualTolerance_) { |
---|
3616 | double gap = upper_[iSequence] - lower_[iSequence]; |
---|
3617 | if (gap < 1.0e20) |
---|
3618 | sumBadPivots -= value * gap; |
---|
3619 | else |
---|
3620 | sumBadPivots += 1.0e20; |
---|
3621 | //printf("bad %d alpha %g dj at lower %g\n", |
---|
3622 | // iSequence,alpha,value); |
---|
3623 | } |
---|
3624 | } |
---|
3625 | } |
---|
3626 | #endif |
---|
3627 | #ifdef FORCE_FOLLOW |
---|
3628 | if (iSequence == force_in) { |
---|
3629 | printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot); |
---|
3630 | take = true; |
---|
3631 | } |
---|
3632 | #endif |
---|
3633 | if (take) { |
---|
3634 | sequenceIn_ = numberPossiblySwapped; |
---|
3635 | bestPivot = absAlpha; |
---|
3636 | theta_ = dj_[iSequence] / alpha; |
---|
3637 | largestPivot = CoinMax(largestPivot, 0.5 * bestPivot); |
---|
3638 | #ifdef DUBIOUS_WEIGHTS |
---|
3639 | bestWeight = weight; |
---|
3640 | #endif |
---|
3641 | //printf(" taken seq %d alpha %g weight %d\n", |
---|
3642 | // iSequence,absAlpha,dubiousWeights[iSequence]); |
---|
3643 | } else { |
---|
3644 | //printf(" not taken seq %d alpha %g weight %d\n", |
---|
3645 | // iSequence,absAlpha,dubiousWeights[iSequence]); |
---|
3646 | } |
---|
3647 | double range = upper_[iSequence] - lower_[iSequence]; |
---|
3648 | thruThis += range * fabs(alpha); |
---|
3649 | increaseInThis += badDj * range; |
---|
3650 | } |
---|
3651 | } |
---|
3652 | marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining); |
---|
3653 | marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped); |
---|
3654 | #ifdef MORE_CAREFUL |
---|
3655 | // If we have done pivots and things look bad set alpha_ 0.0 to force factorization |
---|
3656 | if (sumBadPivots > 1.0e4) { |
---|
3657 | if (handler_->logLevel() > 1) |
---|
3658 | printf("maybe forcing re-factorization - sum %g %d pivots\n", sumBadPivots, |
---|
3659 | factorization_->pivots()); |
---|
3660 | if(factorization_->pivots() > 3) { |
---|
3661 | badSumPivots = true; |
---|
3662 | break; |
---|
3663 | } |
---|
3664 | } |
---|
3665 | #endif |
---|
3666 | swapped[1-iFlip] = numberPossiblySwapped; |
---|
3667 | interesting[1-iFlip] = numberRemaining; |
---|
3668 | // If we stop now this will be increase in objective (I think) |
---|
3669 | double increase = (fabs(dualOut_) - totalThru) * theta_; |
---|
3670 | increase += increaseInObjective; |
---|
3671 | if (theta_ < 0.0) |
---|
3672 | thruThis += fabs(dualOut_); // force using this one |
---|
3673 | if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) { |
---|
3674 | // back |
---|
3675 | // We may need to be more careful - we could do by |
---|
3676 | // switch so we always do fine grained? |
---|
3677 | bestPivot = 0.0; |
---|
3678 | } else { |
---|
3679 | // add in |
---|
3680 | totalThru += thruThis; |
---|
3681 | increaseInObjective += increaseInThis; |
---|
3682 | } |
---|
3683 | if (bestPivot < 0.1 * bestEverPivot && |
---|
3684 | bestEverPivot > 1.0e-6 && |
---|
3685 | (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) { |
---|
3686 | // back to previous one |
---|
3687 | sequenceIn_ = lastSequence; |
---|
3688 | // swap regions |
---|
3689 | iFlip = 1 - iFlip; |
---|
3690 | break; |
---|
3691 | } else if (sequenceIn_ == -1 && upperTheta > largeValue_) { |
---|
3692 | if (lastPivot > acceptablePivot) { |
---|
3693 | // back to previous one |
---|
3694 | sequenceIn_ = lastSequence; |
---|
3695 | // swap regions |
---|
3696 | iFlip = 1 - iFlip; |
---|
3697 | } else { |
---|
3698 | // can only get here if all pivots too small |
---|
3699 | } |
---|
3700 | break; |
---|
3701 | } else if (totalThru >= fabs(dualOut_)) { |
---|
3702 | modifyCosts = true; // fine grain - we can modify costs |
---|
3703 | break; // no point trying another loop |
---|
3704 | } else { |
---|
3705 | lastSequence = sequenceIn_; |
---|
3706 | if (bestPivot > bestEverPivot) |
---|
3707 | bestEverPivot = bestPivot; |
---|
3708 | iFlip = 1 - iFlip; |
---|
3709 | modifyCosts = true; // fine grain - we can modify costs |
---|
3710 | } |
---|
3711 | } |
---|
3712 | if (iTry == MAXTRY) |
---|
3713 | iFlip = 1 - iFlip; // flip back |
---|
3714 | break; |
---|
3715 | } else { |
---|
3716 | // skip this lot |
---|
3717 | if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) { |
---|
3718 | bestEverPivot = bestPivot; |
---|
3719 | lastSequence = bestSequence; |
---|
3720 | } else { |
---|
3721 | // keep old swapped |
---|
3722 | CoinMemcpyN(array[iFlip] + swapped[iFlip], |
---|
3723 | numberColumns_ - swapped[iFlip], array[1-iFlip] + swapped[iFlip]); |
---|
3724 | CoinMemcpyN(indices[iFlip] + swapped[iFlip], |
---|
3725 | numberColumns_ - swapped[iFlip], indices[1-iFlip] + swapped[iFlip]); |
---|
3726 | marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], swapped[iFlip]); |
---|
3727 | swapped[1-iFlip] = swapped[iFlip]; |
---|
3728 | } |
---|
3729 | increaseInObjective += increaseInThis; |
---|
3730 | iFlip = 1 - iFlip; // swap regions |
---|
3731 | tentativeTheta = 2.0 * upperTheta; |
---|
3732 | totalThru += thruThis; |
---|
3733 | } |
---|
3734 | } |
---|
3735 | |
---|
3736 | // can get here without sequenceIn_ set but with lastSequence |
---|
3737 | if (sequenceIn_ < 0 && lastSequence >= 0) { |
---|
3738 | // back to previous one |
---|
3739 | sequenceIn_ = lastSequence; |
---|
3740 | // swap regions |
---|
3741 | iFlip = 1 - iFlip; |
---|
3742 | } |
---|
3743 | |
---|
3744 | #define MINIMUMTHETA 1.0e-18 |
---|
3745 | // Movement should be minimum for anti-degeneracy - unless |
---|
3746 | // fixed variable out |
---|
3747 | double minimumTheta; |
---|
3748 | if (upperOut_ > lowerOut_) |
---|
3749 | minimumTheta = MINIMUMTHETA; |
---|
3750 | else |
---|
3751 | minimumTheta = 0.0; |
---|
3752 | if (sequenceIn_ >= 0) { |
---|
3753 | // at this stage sequenceIn_ is just pointer into index array |
---|
3754 | // flip just so we can use iFlip |
---|
3755 | iFlip = 1 - iFlip; |
---|
3756 | spare = array[iFlip]; |
---|
3757 | index = indices[iFlip]; |
---|
3758 | double oldValue; |
---|
3759 | alpha_ = spare[sequenceIn_]; |
---|
3760 | sequenceIn_ = indices[iFlip][sequenceIn_]; |
---|
3761 | oldValue = dj_[sequenceIn_]; |
---|
3762 | theta_ = CoinMax(oldValue / alpha_, 0.0); |
---|
3763 | if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) { |
---|
3764 | // can't pivot to zero |
---|
3765 | #if 0 |
---|
3766 | if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) { |
---|
3767 | theta_ = minimumTheta; |
---|
3768 | } else if (oldValue - minimumTheta*alpha_ >= -newTolerance) { |
---|
3769 | theta_ = minimumTheta; |
---|
3770 | thisIncrease = true; |
---|
3771 | } else { |
---|
3772 | theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0); |
---|
3773 | thisIncrease = true; |
---|
3774 | } |
---|
3775 | #else |
---|
3776 | theta_ = minimumTheta; |
---|
3777 | #endif |
---|
3778 | } |
---|
3779 | // may need to adjust costs so all dual feasible AND pivoted is exactly 0 |
---|
3780 | //int costOffset = numberRows_+numberColumns_; |
---|
3781 | if (modifyCosts && !badSumPivots) { |
---|
3782 | int i; |
---|
3783 | for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) { |
---|
3784 | int iSequence = index[i]; |
---|
3785 | double alpha = spare[i]; |
---|
3786 | double value = dj_[iSequence] - theta_ * alpha; |
---|
3787 | |
---|
3788 | // can't be free here |
---|
3789 | |
---|
3790 | if (alpha < 0.0) { |
---|
3791 | //at upper bound |
---|
3792 | if (value > dualTolerance_) { |
---|
3793 | thisIncrease = true; |
---|
3794 | #if CLP_CAN_HAVE_ZERO_OBJ<2 |
---|
3795 | #define MODIFYCOST 2 |
---|
3796 | #endif |
---|
3797 | #if MODIFYCOST |
---|
3798 | // modify cost to hit new tolerance |
---|
3799 | double modification = alpha * theta_ - dj_[iSequence] |
---|
3800 | + newTolerance; |
---|
3801 | if ((specialOptions_&(2048 + 4096 + 16384)) != 0) { |
---|
3802 | if ((specialOptions_ & 16384) != 0) { |
---|
3803 | if (fabs(modification) < 1.0e-8) |
---|
3804 | modification = 0.0; |
---|
3805 | } else if ((specialOptions_ & 2048) != 0) { |
---|
3806 | if (fabs(modification) < 1.0e-10) |
---|
3807 | modification = 0.0; |
---|
3808 | } else { |
---|
3809 | if (fabs(modification) < 1.0e-12) |
---|
3810 | modification = 0.0; |
---|
3811 | } |
---|
3812 | } |
---|
3813 | dj_[iSequence] += modification; |
---|
3814 | cost_[iSequence] += modification; |
---|
3815 | if (modification) |
---|
3816 | numberChanged_ ++; // Say changed costs |
---|
3817 | //cost_[iSequence+costOffset] += modification; // save change |
---|
3818 | #endif |
---|
3819 | } |
---|
3820 | } else { |
---|
3821 | // at lower bound |
---|
3822 | if (-value > dualTolerance_) { |
---|
3823 | thisIncrease = true; |
---|
3824 | #if MODIFYCOST |
---|
3825 | // modify cost to hit new tolerance |
---|
3826 | double modification = alpha * theta_ - dj_[iSequence] |
---|
3827 | - newTolerance; |
---|
3828 | //modification = CoinMax(modification,-dualTolerance_); |
---|
3829 | //assert (fabs(modification)<1.0e-7); |
---|
3830 | if ((specialOptions_&(2048 + 4096)) != 0) { |
---|
3831 | if ((specialOptions_ & 2048) != 0) { |
---|
3832 | if (fabs(modification) < 1.0e-10) |
---|
3833 | modification = 0.0; |
---|
3834 | } else { |
---|
3835 | if (fabs(modification) < 1.0e-12) |
---|
3836 | modification = 0.0; |
---|
3837 | } |
---|
3838 | } |
---|
3839 | dj_[iSequence] += modification; |
---|
3840 | cost_[iSequence] += modification; |
---|
3841 | if (modification) |
---|
3842 | numberChanged_ ++; // Say changed costs |
---|
3843 | //cost_[iSequence+costOffset] += modification; // save change |
---|
3844 | #endif |
---|
3845 | } |
---|
3846 | } |
---|
3847 | } |
---|
3848 | } |
---|
3849 | } |
---|
3850 | } |
---|
3851 | |
---|
3852 | #ifdef MORE_CAREFUL |
---|
3853 | // If we have done pivots and things look bad set alpha_ 0.0 to force factorization |
---|
3854 | if ((badSumPivots || |
---|
3855 | fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) { |
---|
3856 | if (handler_->logLevel() > 1) |
---|
3857 | printf("forcing re-factorization\n"); |
---|
3858 | sequenceIn_ = -1; |
---|
3859 | } |
---|
3860 | #endif |
---|
3861 | if (sequenceIn_ >= 0) { |
---|
3862 | lowerIn_ = lower_[sequenceIn_]; |
---|
3863 | upperIn_ = upper_[sequenceIn_]; |
---|
3864 | valueIn_ = solution_[sequenceIn_]; |
---|
3865 | dualIn_ = dj_[sequenceIn_]; |
---|
3866 | |
---|
3867 | if (numberTimesOptimal_) { |
---|
3868 | // can we adjust cost back closer to original |
---|
3869 | //*** add coding |
---|
3870 | } |
---|
3871 | #if MODIFYCOST>1 |
---|
3872 | // modify cost to hit zero exactly |
---|
3873 | // so (dualIn_+modification)==theta_*alpha_ |
---|
3874 | double modification = theta_ * alpha_ - dualIn_; |
---|
3875 | // But should not move objective too much ?? |
---|
3876 | #define DONT_MOVE_OBJECTIVE |
---|
3877 | #ifdef DONT_MOVE_OBJECTIVE |
---|
3878 | double moveObjective = fabs(modification * solution_[sequenceIn_]); |
---|
3879 | double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3); |
---|
3880 | if (moveObjective > smallMove) { |
---|
3881 | if (handler_->logLevel() > 1) |
---|
3882 | printf("would move objective by %g - original mod %g sol value %g\n", moveObjective, |
---|
3883 | modification, solution_[sequenceIn_]); |
---|
3884 | modification *= smallMove / moveObjective; |
---|
3885 | } |
---|
3886 | #endif |
---|
3887 | if (badSumPivots) |
---|
3888 | modification = 0.0; |
---|
3889 | if ((specialOptions_&(2048 + 4096)) != 0) { |
---|
3890 | if ((specialOptions_ & 16384) != 0) { |
---|
3891 | // in fast dual |
---|
3892 | if (fabs(modification) < 1.0e-7) |
---|
3893 | modification = 0.0; |
---|
3894 | } else if ((specialOptions_ & 2048) != 0) { |
---|
3895 | if (fabs(modification) < 1.0e-10) |
---|
3896 | modification = 0.0; |
---|
3897 | } else { |
---|
3898 | if (fabs(modification) < 1.0e-12) |
---|
3899 | modification = 0.0; |
---|
3900 | } |
---|
3901 | } |
---|
3902 | dualIn_ += modification; |
---|
3903 | dj_[sequenceIn_] = dualIn_; |
---|
3904 | cost_[sequenceIn_] += modification; |
---|
3905 | if (modification) |
---|
3906 | numberChanged_ ++; // Say changed costs |
---|
3907 | //int costOffset = numberRows_+numberColumns_; |
---|
3908 | //cost_[sequenceIn_+costOffset] += modification; // save change |
---|
3909 | //assert (fabs(modification)<1.0e-6); |
---|
3910 | #ifdef CLP_DEBUG |
---|
3911 | if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15) |
---|
3912 | printf("exact %d new cost %g, change %g\n", sequenceIn_, |
---|
3913 | cost_[sequenceIn_], modification); |
---|
3914 | #endif |
---|
3915 | #endif |
---|
3916 | |
---|
3917 | if (alpha_ < 0.0) { |
---|
3918 | // as if from upper bound |
---|
3919 | directionIn_ = -1; |
---|
3920 | upperIn_ = valueIn_; |
---|
3921 | } else { |
---|
3922 | // as if from lower bound |
---|
3923 | directionIn_ = 1; |
---|
3924 | lowerIn_ = valueIn_; |
---|
3925 | } |
---|
3926 | } else { |
---|
3927 | // no pivot |
---|
3928 | bestPossible = 0.0; |
---|
3929 | alpha_ = 0.0; |
---|
3930 | } |
---|
3931 | //if (thisIncrease) |
---|
3932 | //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance]; |
---|
3933 | |
---|
3934 | // clear arrays |
---|
3935 | |
---|
3936 | for (i = 0; i < 2; i++) { |
---|
3937 | CoinZeroN(array[i], marker[i][0]); |
---|
3938 | CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]); |
---|
3939 | } |
---|
3940 | return bestPossible; |
---|
3941 | } |
---|
3942 | #ifdef CLP_ALL_ONE_FILE |
---|
3943 | #undef MAXTRY |
---|
3944 | #endif |
---|
3945 | /* Checks if tentative optimal actually means unbounded |
---|
3946 | Returns -3 if not, 2 if is unbounded */ |
---|
3947 | int |
---|
3948 | ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray, |
---|
3949 | CoinIndexedVector * spare, |
---|
3950 | double changeCost) |
---|
3951 | { |
---|
3952 | int status = 2; // say unbounded |
---|
3953 | factorization_->updateColumn(spare, ray); |
---|
3954 | // get reduced cost |
---|
3955 | int i; |
---|
3956 | int number = ray->getNumElements(); |
---|
3957 | int * index = ray->getIndices(); |
---|
3958 | double * array = ray->denseVector(); |
---|
3959 | for (i = 0; i < number; i++) { |
---|
3960 | int iRow = index[i]; |
---|
3961 | int iPivot = pivotVariable_[iRow]; |
---|
3962 | changeCost -= cost(iPivot) * array[iRow]; |
---|
3963 | } |
---|
3964 | double way; |
---|
3965 | if (changeCost > 0.0) { |
---|
3966 | //try going down |
---|
3967 | way = 1.0; |
---|
3968 | } else if (changeCost < 0.0) { |
---|
3969 | //try going up |
---|
3970 | way = -1.0; |
---|
3971 | } else { |
---|
3972 | #ifdef CLP_DEBUG |
---|
3973 | printf("can't decide on up or down\n"); |
---|
3974 | #endif |
---|
3975 | way = 0.0; |
---|
3976 | status = -3; |
---|
3977 | } |
---|
3978 | double movement = 1.0e10 * way; // some largish number |
---|
3979 | double zeroTolerance = 1.0e-14 * dualBound_; |
---|
3980 | for (i = 0; i < number; i++) { |
---|
3981 | int iRow = index[i]; |
---|
3982 | int iPivot = pivotVariable_[iRow]; |
---|
3983 | double arrayValue = array[iRow]; |
---|
3984 | if (fabs(arrayValue) < zeroTolerance) |
---|
3985 | arrayValue = 0.0; |
---|
3986 | double newValue = solution(iPivot) + movement * arrayValue; |
---|
3987 | if (newValue > upper(iPivot) + primalTolerance_ || |
---|
3988 | newValue < lower(iPivot) - primalTolerance_) |
---|
3989 | status = -3; // not unbounded |
---|
3990 | } |
---|
3991 | if (status == 2) { |
---|
3992 | // create ray |
---|
3993 | delete [] ray_; |
---|
3994 | ray_ = new double [numberColumns_]; |
---|
3995 | CoinZeroN(ray_, numberColumns_); |
---|
3996 | for (i = 0; i < number; i++) { |
---|
3997 | int iRow = index[i]; |
---|
3998 | int iPivot = pivotVariable_[iRow]; |
---|
3999 | double arrayValue = array[iRow]; |
---|
4000 | if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance) |
---|
4001 | ray_[iPivot] = way * array[iRow]; |
---|
4002 | } |
---|
4003 | } |
---|
4004 | ray->clear(); |
---|
4005 | return status; |
---|
4006 | } |
---|
4007 | //static int count_alpha=0; |
---|
4008 | /* Checks if finished. Updates status */ |
---|
4009 | void |
---|
4010 | ClpSimplexDual::statusOfProblemInDual(int & lastCleaned, int type, |
---|
4011 | double * givenDuals, ClpDataSave & saveData, |
---|
4012 | int ifValuesPass) |
---|
4013 | { |
---|
4014 | #ifdef CLP_INVESTIGATE_SERIAL |
---|
4015 | if (z_thinks > 0 && z_thinks < 2) |
---|
4016 | z_thinks += 2; |
---|
4017 | #endif |
---|
4018 | bool arraysNotCreated = (type==0); |
---|
4019 | // If lots of iterations then adjust costs if large ones |
---|
4020 | if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) { |
---|
4021 | double largest = 0.0; |
---|
4022 | for (int i = 0; i < numberRows_; i++) { |
---|
4023 | int iColumn = pivotVariable_[i]; |
---|
4024 | largest = CoinMax(largest, fabs(cost_[iColumn])); |
---|
4025 | } |
---|
4026 | if (largest > 1.0e6) { |
---|
4027 | objectiveScale_ = 1.0e6 / largest; |
---|
4028 | for (int i = 0; i < numberRows_ + numberColumns_; i++) |
---|
4029 | cost_[i] *= objectiveScale_; |
---|
4030 | } |
---|
4031 | } |
---|
4032 | int numberPivots = factorization_->pivots(); |
---|
4033 | double realDualInfeasibilities = 0.0; |
---|
4034 | if (type == 2) { |
---|
4035 | if (alphaAccuracy_ != -1.0) |
---|
4036 | alphaAccuracy_ = -2.0; |
---|
4037 | // trouble - restore solution |
---|
4038 | CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_); |
---|
4039 | CoinMemcpyN(savedSolution_ + numberColumns_ , |
---|
4040 | numberRows_, rowActivityWork_); |
---|
4041 | CoinMemcpyN(savedSolution_ , |
---|
4042 | numberColumns_, columnActivityWork_); |
---|
4043 | // restore extra stuff |
---|
4044 | int dummy; |
---|
4045 | matrix_->generalExpanded(this, 6, dummy); |
---|
4046 | forceFactorization_ = 1; // a bit drastic but .. |
---|
4047 | changeMade_++; // say something changed |
---|
4048 | // get correct bounds on all variables |
---|
4049 | resetFakeBounds(0); |
---|
4050 | } |
---|
4051 | int tentativeStatus = problemStatus_; |
---|
4052 | double changeCost; |
---|
4053 | bool unflagVariables = true; |
---|
4054 | bool weightsSaved = false; |
---|
4055 | bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_; |
---|
4056 | int dontFactorizePivots = dontFactorizePivots_; |
---|
4057 | if (type == 3) { |
---|
4058 | type = 1; |
---|
4059 | dontFactorizePivots = 1; |
---|
4060 | } |
---|
4061 | if (alphaAccuracy_ < 0.0 || !numberPivots || alphaAccuracy_ > 1.0e4 || numberPivots > 20) { |
---|
4062 | if (problemStatus_ > -3 || numberPivots > dontFactorizePivots) { |
---|
4063 | // factorize |
---|
4064 | // later on we will need to recover from singularities |
---|
4065 | // also we could skip if first time |
---|
4066 | // save dual weights |
---|
4067 | dualRowPivot_->saveWeights(this, 1); |
---|
4068 | weightsSaved = true; |
---|
4069 | if (type) { |
---|
4070 | // is factorization okay? |
---|
4071 | if (internalFactorize(1)) { |
---|
4072 | // no - restore previous basis |
---|
4073 | unflagVariables = false; |
---|
4074 | assert (type == 1); |
---|
4075 | changeMade_++; // say something changed |
---|
4076 | // Keep any flagged variables |
---|
4077 | int i; |
---|
4078 | for (i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
4079 | if (flagged(i)) |
---|
4080 | saveStatus_[i] |= 64; //say flagged |
---|
4081 | } |
---|
4082 | CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_); |
---|
4083 | CoinMemcpyN(savedSolution_ + numberColumns_ , |
---|
4084 | numberRows_, rowActivityWork_); |
---|
4085 | CoinMemcpyN(savedSolution_ , |
---|
4086 | numberColumns_, columnActivityWork_); |
---|
4087 | // restore extra stuff |
---|
4088 | int dummy; |
---|
4089 | matrix_->generalExpanded(this, 6, dummy); |
---|
4090 | // get correct bounds on all variables |
---|
4091 | resetFakeBounds(1); |
---|
4092 | // need to reject something |
---|
4093 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
4094 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
4095 | << x << sequenceWithin(sequenceOut_) |
---|
4096 | << CoinMessageEol; |
---|
4097 | #ifdef COIN_DEVELOP |
---|
4098 | printf("flag d\n"); |
---|
4099 | #endif |
---|
4100 | setFlagged(sequenceOut_); |
---|
4101 | progress_.clearBadTimes(); |
---|
4102 | |
---|
4103 | // Go to safe |
---|
4104 | factorization_->pivotTolerance(0.99); |
---|
4105 | forceFactorization_ = 1; // a bit drastic but .. |
---|
4106 | type = 2; |
---|
4107 | //assert (internalFactorize(1)==0); |
---|
4108 | if (internalFactorize(1)) { |
---|
4109 | CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_); |
---|
4110 | CoinMemcpyN(savedSolution_ + numberColumns_ , |
---|
4111 | numberRows_, rowActivityWork_); |
---|
4112 | CoinMemcpyN(savedSolution_ , |
---|
4113 | numberColumns_, columnActivityWork_); |
---|
4114 | // restore extra stuff |
---|
4115 | int dummy; |
---|
4116 | matrix_->generalExpanded(this, 6, dummy); |
---|
4117 | // debug |
---|
4118 | int returnCode = internalFactorize(1); |
---|
4119 | while (returnCode) { |
---|
4120 | // ouch |
---|
4121 | // switch off dense |
---|
4122 | int saveDense = factorization_->denseThreshold(); |
---|
4123 | factorization_->setDenseThreshold(0); |
---|
4124 | // Go to safe |
---|
4125 | factorization_->pivotTolerance(0.99); |
---|
4126 | // make sure will do safe factorization |
---|
4127 | pivotVariable_[0] = -1; |
---|
4128 | returnCode = internalFactorize(2); |
---|
4129 | factorization_->setDenseThreshold(saveDense); |
---|
4130 | } |
---|
4131 | // get correct bounds on all variables |
---|
4132 | resetFakeBounds(1); |
---|
4133 | } |
---|
4134 | } |
---|
4135 | } |
---|
4136 | if (problemStatus_ != -4 || numberPivots > 10) |
---|
4137 | problemStatus_ = -3; |
---|
4138 | } |
---|
4139 | } else { |
---|
4140 | //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_); |
---|
4141 | //count_alpha++; |
---|
4142 | //if ((count_alpha%5000)==0) |
---|
4143 | //printf("count alpha %d\n",count_alpha); |
---|
4144 | } |
---|
4145 | // at this stage status is -3 or -4 if looks infeasible |
---|
4146 | // get primal and dual solutions |
---|
4147 | #if 0 |
---|
4148 | { |
---|
4149 | int numberTotal = numberRows_ + numberColumns_; |
---|
4150 | double * saveSol = CoinCopyOfArray(solution_, numberTotal); |
---|
4151 | double * saveDj = CoinCopyOfArray(dj_, numberTotal); |
---|
4152 | double tolerance = type ? 1.0e-4 : 1.0e-8; |
---|
4153 | // always if values pass |
---|
4154 | double saveObj = objectiveValue_; |
---|
4155 | double sumPrimal = sumPrimalInfeasibilities_; |
---|
4156 | int numberPrimal = numberPrimalInfeasibilities_; |
---|
4157 | double sumDual = sumDualInfeasibilities_; |
---|
4158 | int numberDual = numberDualInfeasibilities_; |
---|
4159 | gutsOfSolution(givenDuals, NULL); |
---|
4160 | int j; |
---|
4161 | double largestPrimal = tolerance; |
---|
4162 | int iPrimal = -1; |
---|
4163 | for (j = 0; j < numberTotal; j++) { |
---|
4164 | double difference = solution_[j] - saveSol[j]; |
---|
4165 | if (fabs(difference) > largestPrimal) { |
---|
4166 | iPrimal = j; |
---|
4167 | largestPrimal = fabs(difference); |
---|
4168 | } |
---|
4169 | } |
---|
4170 | double largestDual = tolerance; |
---|
4171 | int iDual = -1; |
---|
4172 | for (j = 0; j < numberTotal; j++) { |
---|
4173 | double difference = dj_[j] - saveDj[j]; |
---|
4174 | if (fabs(difference) > largestDual && upper_[j] > lower_[j]) { |
---|
4175 | iDual = j; |
---|
4176 | largestDual = fabs(difference); |
---|
4177 | } |
---|
4178 | } |
---|
4179 | if (!type) { |
---|
4180 | if (fabs(saveObj - objectiveValue_) > 1.0e-5 || |
---|
4181 | numberPrimal != numberPrimalInfeasibilities_ || numberPrimal != 1 || |
---|
4182 | fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 || iPrimal >= 0 || |
---|
4183 | numberDual != numberDualInfeasibilities_ || numberDual != 0 || |
---|
4184 | fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 || iDual >= 0) |
---|
4185 | 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", |
---|
4186 | type, numberIterations_, numberPivots, |
---|
4187 | numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_, |
---|
4188 | largestPrimal, iPrimal, |
---|
4189 | numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_, |
---|
4190 | largestDual, iDual, |
---|
4191 | saveObj, objectiveValue_); |
---|
4192 | } else { |
---|
4193 | if (fabs(saveObj - objectiveValue_) > 1.0e-5 || |
---|
4194 | numberPrimalInfeasibilities_ || iPrimal >= 0 || |
---|
4195 | numberDualInfeasibilities_ || iDual >= 0) |
---|
4196 | 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", |
---|
4197 | type, numberIterations_, numberPivots, |
---|
4198 | numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_, |
---|
4199 | largestPrimal, iPrimal, |
---|
4200 | numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_, |
---|
4201 | largestDual, iDual, |
---|
4202 | saveObj, objectiveValue_); |
---|
4203 | } |
---|
4204 | delete [] saveSol; |
---|
4205 | delete [] saveDj; |
---|
4206 | } |
---|
4207 | #else |
---|
4208 | if (type || ifValuesPass) |
---|
4209 | gutsOfSolution(givenDuals, NULL); |
---|
4210 | #endif |
---|
4211 | // If bad accuracy treat as singular |
---|
4212 | if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) { |
---|
4213 | // restore previous basis |
---|
4214 | unflagVariables = false; |
---|
4215 | changeMade_++; // say something changed |
---|
4216 | // Keep any flagged variables |
---|
4217 | int i; |
---|
4218 | for (i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
4219 | if (flagged(i)) |
---|
4220 | saveStatus_[i] |= 64; //say flagged |
---|
4221 | } |
---|
4222 | CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_); |
---|
4223 | CoinMemcpyN(savedSolution_ + numberColumns_ , |
---|
4224 | numberRows_, rowActivityWork_); |
---|
4225 | CoinMemcpyN(savedSolution_ , |
---|
4226 | numberColumns_, columnActivityWork_); |
---|
4227 | // restore extra stuff |
---|
4228 | int dummy; |
---|
4229 | matrix_->generalExpanded(this, 6, dummy); |
---|
4230 | // get correct bounds on all variables |
---|
4231 | resetFakeBounds(1); |
---|
4232 | // need to reject something |
---|
4233 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
4234 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
4235 | << x << sequenceWithin(sequenceOut_) |
---|
4236 | << CoinMessageEol; |
---|
4237 | #ifdef COIN_DEVELOP |
---|
4238 | printf("flag e\n"); |
---|
4239 | #endif |
---|
4240 | setFlagged(sequenceOut_); |
---|
4241 | progress_.clearBadTimes(); |
---|
4242 | |
---|
4243 | // Go to safer |
---|
4244 | double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99); |
---|
4245 | factorization_->pivotTolerance(newTolerance); |
---|
4246 | forceFactorization_ = 1; // a bit drastic but .. |
---|
4247 | if (alphaAccuracy_ != -1.0) |
---|
4248 | alphaAccuracy_ = -2.0; |
---|
4249 | type = 2; |
---|
4250 | //assert (internalFactorize(1)==0); |
---|
4251 | if (internalFactorize(1)) { |
---|
4252 | CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_); |
---|
4253 | CoinMemcpyN(savedSolution_ + numberColumns_ , |
---|
4254 | numberRows_, rowActivityWork_); |
---|
4255 | CoinMemcpyN(savedSolution_ , |
---|
4256 | numberColumns_, columnActivityWork_); |
---|
4257 | // restore extra stuff |
---|
4258 | int dummy; |
---|
4259 | matrix_->generalExpanded(this, 6, dummy); |
---|
4260 | // debug |
---|
4261 | int returnCode = internalFactorize(1); |
---|
4262 | while (returnCode) { |
---|
4263 | // ouch |
---|
4264 | // switch off dense |
---|
4265 | int saveDense = factorization_->denseThreshold(); |
---|
4266 | factorization_->setDenseThreshold(0); |
---|
4267 | // Go to safe |
---|
4268 | factorization_->pivotTolerance(0.99); |
---|
4269 | // make sure will do safe factorization |
---|
4270 | pivotVariable_[0] = -1; |
---|
4271 | returnCode = internalFactorize(2); |
---|
4272 | factorization_->setDenseThreshold(saveDense); |
---|
4273 | } |
---|
4274 | // get correct bounds on all variables |
---|
4275 | resetFakeBounds(1); |
---|
4276 | } |
---|
4277 | // get primal and dual solutions |
---|
4278 | gutsOfSolution(givenDuals, NULL); |
---|
4279 | } else if (goodAccuracy()) { |
---|
4280 | // Can reduce tolerance |
---|
4281 | double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_); |
---|
4282 | factorization_->pivotTolerance(newTolerance); |
---|
4283 | } |
---|
4284 | bestObjectiveValue_ = CoinMax(bestObjectiveValue_, |
---|
4285 | objectiveValue_ - bestPossibleImprovement_); |
---|
4286 | bool reallyBadProblems = false; |
---|
4287 | // Double check infeasibility if no action |
---|
4288 | if (progress_.lastIterationNumber(0) == numberIterations_) { |
---|
4289 | if (dualRowPivot_->looksOptimal()) { |
---|
4290 | numberPrimalInfeasibilities_ = 0; |
---|
4291 | sumPrimalInfeasibilities_ = 0.0; |
---|
4292 | } |
---|
4293 | #if 1 |
---|
4294 | } else { |
---|
4295 | double thisObj = objectiveValue_ - bestPossibleImprovement_; |
---|
4296 | #ifdef CLP_INVESTIGATE |
---|
4297 | assert (bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100); |
---|
4298 | if (bestPossibleImprovement_) |
---|
4299 | printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_, |
---|
4300 | thisObj); |
---|
4301 | #endif |
---|
4302 | double lastObj = progress_.lastObjective(0); |
---|
4303 | #ifndef NDEBUG |
---|
4304 | #ifdef COIN_DEVELOP |
---|
4305 | resetFakeBounds(-1); |
---|
4306 | #endif |
---|
4307 | #endif |
---|
4308 | #ifdef CLP_REPORT_PROGRESS |
---|
4309 | ixxxxxx++; |
---|
4310 | if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) { |
---|
4311 | char temp[20]; |
---|
4312 | sprintf(temp, "sol%d.out", ixxxxxx); |
---|
4313 | printf("sol%d.out\n", ixxxxxx); |
---|
4314 | FILE * fp = fopen(temp, "w"); |
---|
4315 | int nTotal = numberRows_ + numberColumns_; |
---|
4316 | for (int i = 0; i < nTotal; i++) |
---|
4317 | fprintf(fp, "%d %d %g %g %g %g %g\n", |
---|
4318 | i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]); |
---|
4319 | fclose(fp); |
---|
4320 | } |
---|
4321 | #endif |
---|
4322 | if(!ifValuesPass && firstFree_ < 0) { |
---|
4323 | double testTol = 5.0e-3; |
---|
4324 | if (progress_.timesFlagged() > 10) { |
---|
4325 | testTol *= pow(2.0, progress_.timesFlagged() - 8); |
---|
4326 | } else if (progress_.timesFlagged() > 5) { |
---|
4327 | testTol *= 5.0; |
---|
4328 | } |
---|
4329 | if (lastObj > thisObj + |
---|
4330 | testTol*(fabs(thisObj) + fabs(lastObj)) + testTol) { |
---|
4331 | int maxFactor = factorization_->maximumPivots(); |
---|
4332 | if ((specialOptions_ & 1048576) == 0) { |
---|
4333 | if (progress_.timesFlagged() > 10) |
---|
4334 | progress_.incrementReallyBadTimes(); |
---|
4335 | if (maxFactor > 10 - 9) { |
---|
4336 | #ifdef COIN_DEVELOP |
---|
4337 | printf("lastobj %g thisobj %g\n", lastObj, thisObj); |
---|
4338 | #endif |
---|
4339 | //if (forceFactorization_<0) |
---|
4340 | //forceFactorization_= maxFactor; |
---|
4341 | //forceFactorization_ = CoinMax(1,(forceFactorization_>>1)); |
---|
4342 | if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 && |
---|
4343 | largestPrimalError_ < 1.0e2) { |
---|
4344 | // Just save costs |
---|
4345 | // save extra copy of cost_ |
---|
4346 | int nTotal = numberRows_ + numberColumns_; |
---|
4347 | double * temp = new double [2*nTotal]; |
---|
4348 | memcpy(temp, cost_, nTotal * sizeof(double)); |
---|
4349 | memcpy(temp + nTotal, cost_, nTotal * sizeof(double)); |
---|
4350 | delete [] cost_; |
---|
4351 | cost_ = temp; |
---|
4352 | objectiveWork_ = cost_; |
---|
4353 | rowObjectiveWork_ = cost_ + numberColumns_; |
---|
4354 | progressFlag_ |= 4; |
---|
4355 | } else { |
---|
4356 | forceFactorization_ = 1; |
---|
4357 | #ifdef COIN_DEVELOP |
---|
4358 | printf("Reducing factorization frequency - bad backwards\n"); |
---|
4359 | #endif |
---|
4360 | #if 1 |
---|
4361 | unflagVariables = false; |
---|
4362 | changeMade_++; // say something changed |
---|
4363 | int nTotal = numberRows_ + numberColumns_; |
---|
4364 | CoinMemcpyN(saveStatus_, nTotal, status_); |
---|
4365 | CoinMemcpyN(savedSolution_ + numberColumns_ , |
---|
4366 | numberRows_, rowActivityWork_); |
---|
4367 | CoinMemcpyN(savedSolution_ , |
---|
4368 | numberColumns_, columnActivityWork_); |
---|
4369 | if ((progressFlag_ & 4) == 0) { |
---|
4370 | // save extra copy of cost_ |
---|
4371 | double * temp = new double [2*nTotal]; |
---|
4372 | memcpy(temp, cost_, nTotal * sizeof(double)); |
---|
4373 | memcpy(temp + nTotal, cost_, nTotal * sizeof(double)); |
---|
4374 | delete [] cost_; |
---|
4375 | cost_ = temp; |
---|
4376 | objectiveWork_ = cost_; |
---|
4377 | rowObjectiveWork_ = cost_ + numberColumns_; |
---|
4378 | progressFlag_ |= 4; |
---|
4379 | } else { |
---|
4380 | memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double)); |
---|
4381 | } |
---|
4382 | // restore extra stuff |
---|
4383 | int dummy; |
---|
4384 | matrix_->generalExpanded(this, 6, dummy); |
---|
4385 | double pivotTolerance = factorization_->pivotTolerance(); |
---|
4386 | if(pivotTolerance < 0.2) |
---|
4387 | factorization_->pivotTolerance(0.2); |
---|
4388 | else if(progress_.timesFlagged() > 2) |
---|
4389 | factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99)); |
---|
4390 | if (alphaAccuracy_ != -1.0) |
---|
4391 | alphaAccuracy_ = -2.0; |
---|
4392 | if (internalFactorize(1)) { |
---|
4393 | CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_); |
---|
4394 | CoinMemcpyN(savedSolution_ + numberColumns_ , |
---|
4395 | numberRows_, rowActivityWork_); |
---|
4396 | CoinMemcpyN(savedSolution_ , |
---|
4397 | numberColumns_, columnActivityWork_); |
---|
4398 | // restore extra stuff |
---|
4399 | int dummy; |
---|
4400 | matrix_->generalExpanded(this, 6, dummy); |
---|
4401 | // debug |
---|
4402 | int returnCode = internalFactorize(1); |
---|
4403 | while (returnCode) { |
---|
4404 | // ouch |
---|
4405 | // switch off dense |
---|
4406 | int saveDense = factorization_->denseThreshold(); |
---|
4407 | factorization_->setDenseThreshold(0); |
---|
4408 | // Go to safe |
---|
4409 | factorization_->pivotTolerance(0.99); |
---|
4410 | // make sure will do safe factorization |
---|
4411 | pivotVariable_[0] = -1; |
---|
4412 | returnCode = internalFactorize(2); |
---|
4413 | factorization_->setDenseThreshold(saveDense); |
---|
4414 | } |
---|
4415 | } |
---|
4416 | resetFakeBounds(0); |
---|
4417 | type = 2; // so will restore weights |
---|
4418 | // get primal and dual solutions |
---|
4419 | gutsOfSolution(givenDuals, NULL); |
---|
4420 | if (numberPivots < 2) { |
---|
4421 | // need to reject something |
---|
4422 | char x = isColumn(sequenceOut_) ? 'C' : 'R'; |
---|
4423 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
4424 | << x << sequenceWithin(sequenceOut_) |
---|
4425 | << CoinMessageEol; |
---|
4426 | #ifdef COIN_DEVELOP |
---|
4427 | printf("flag d\n"); |
---|
4428 | #endif |
---|
4429 | setFlagged(sequenceOut_); |
---|
4430 | progress_.clearBadTimes(); |
---|
4431 | progress_.incrementTimesFlagged(); |
---|
4432 | } |
---|
4433 | if (numberPivots < 10) |
---|
4434 | reallyBadProblems = true; |
---|
4435 | #ifdef COIN_DEVELOP |
---|
4436 | printf("obj now %g\n", objectiveValue_); |
---|
4437 | #endif |
---|
4438 | progress_.modifyObjective(objectiveValue_ |
---|
4439 | - bestPossibleImprovement_); |
---|
4440 | #endif |
---|
4441 | } |
---|
4442 | } |
---|
4443 | } else { |
---|
4444 | // in fast dual give up |
---|
4445 | #ifdef COIN_DEVELOP |
---|
4446 | printf("In fast dual?\n"); |
---|
4447 | #endif |
---|
4448 | problemStatus_ = 3; |
---|
4449 | } |
---|
4450 | } else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) { |
---|
4451 | numberTimesOptimal_ = 0; |
---|
4452 | } |
---|
4453 | } |
---|
4454 | #endif |
---|
4455 | } |
---|
4456 | // Up tolerance if looks a bit odd |
---|
4457 | if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) { |
---|
4458 | if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) { |
---|
4459 | int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1); |
---|
4460 | if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) { |
---|
4461 | if (factorization_->pivotTolerance() < 0.9) { |
---|
4462 | // up tolerance |
---|
4463 | factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91)); |
---|
4464 | //printf("tol now %g\n",factorization_->pivotTolerance()); |
---|
4465 | progress_.clearIterationNumbers(); |
---|
4466 | } |
---|
4467 | } |
---|
4468 | } |
---|
4469 | } |
---|
4470 | // Check if looping |
---|
4471 | int loop; |
---|
4472 | if (!givenDuals && type != 2) |
---|
4473 | loop = progress_.looping(); |
---|
4474 | else |
---|
4475 | loop = -1; |
---|
4476 | if (progress_.reallyBadTimes() > 10) { |
---|
4477 | problemStatus_ = 10; // instead - try other algorithm |
---|
4478 | #if COIN_DEVELOP>2 |
---|
4479 | printf("returning at %d\n", __LINE__); |
---|
4480 | #endif |
---|
4481 | } |
---|
4482 | int situationChanged = 0; |
---|
4483 | if (loop >= 0) { |
---|
4484 | problemStatus_ = loop; //exit if in loop |
---|
4485 | if (!problemStatus_) { |
---|
4486 | // declaring victory |
---|
4487 | numberPrimalInfeasibilities_ = 0; |
---|
4488 | sumPrimalInfeasibilities_ = 0.0; |
---|
4489 | } else { |
---|
4490 | problemStatus_ = 10; // instead - try other algorithm |
---|
4491 | #if COIN_DEVELOP>2 |
---|
4492 | printf("returning at %d\n", __LINE__); |
---|
4493 | #endif |
---|
4494 | } |
---|
4495 | return; |
---|
4496 | } else if (loop < -1) { |
---|
4497 | // something may have changed |
---|
4498 | gutsOfSolution(NULL, NULL); |
---|
4499 | situationChanged = 1; |
---|
4500 | } |
---|
4501 | // really for free variables in |
---|
4502 | if((progressFlag_ & 2) != 0) { |
---|
4503 | situationChanged = 2; |
---|
4504 | } |
---|
4505 | progressFlag_ &= (~3); //reset progress flag |
---|
4506 | if ((progressFlag_ & 4) != 0) { |
---|
4507 | // save copy of cost_ |
---|
4508 | int nTotal = numberRows_ + numberColumns_; |
---|
4509 | memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double)); |
---|
4510 | } |
---|
4511 | /*if (!numberIterations_&&sumDualInfeasibilities_) |
---|
4512 | printf("OBJ %g sumPinf %g sumDinf %g\n", |
---|
4513 | objectiveValue(),sumPrimalInfeasibilities_, |
---|
4514 | sumDualInfeasibilities_);*/ |
---|
4515 | // mark as having gone optimal if looks like it |
---|
4516 | if (!numberPrimalInfeasibilities_&& |
---|
4517 | !numberDualInfeasibilities_) |
---|
4518 | progressFlag_ |= 8; |
---|
4519 | if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) { |
---|
4520 | handler_->message(CLP_SIMPLEX_STATUS, messages_) |
---|
4521 | << numberIterations_ << objectiveValue(); |
---|
4522 | handler_->printing(sumPrimalInfeasibilities_ > 0.0) |
---|
4523 | << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_; |
---|
4524 | handler_->printing(sumDualInfeasibilities_ > 0.0) |
---|
4525 | << sumDualInfeasibilities_ << numberDualInfeasibilities_; |
---|
4526 | handler_->printing(numberDualInfeasibilitiesWithoutFree_ |
---|
4527 | < numberDualInfeasibilities_) |
---|
4528 | << numberDualInfeasibilitiesWithoutFree_; |
---|
4529 | handler_->message() << CoinMessageEol; |
---|
4530 | } |
---|
4531 | #if 0 |
---|
4532 | printf("IT %d %g %g(%d) %g(%d)\n", |
---|
4533 | numberIterations_, objectiveValue(), |
---|
4534 | sumPrimalInfeasibilities_, numberPrimalInfeasibilities_, |
---|
4535 | sumDualInfeasibilities_, numberDualInfeasibilities_); |
---|
4536 | #endif |
---|
4537 | double approximateObjective = objectiveValue_; |
---|
4538 | #ifdef CLP_REPORT_PROGRESS |
---|
4539 | if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) { |
---|
4540 | char temp[20]; |
---|
4541 | sprintf(temp, "x_sol%d.out", ixxxxxx); |
---|
4542 | FILE * fp = fopen(temp, "w"); |
---|
4543 | int nTotal = numberRows_ + numberColumns_; |
---|
4544 | for (int i = 0; i < nTotal; i++) |
---|
4545 | fprintf(fp, "%d %d %g %g %g %g %g\n", |
---|
4546 | i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]); |
---|
4547 | fclose(fp); |
---|
4548 | if (ixxxxxx == ixxyyyy) |
---|
4549 | exit(6); |
---|
4550 | } |
---|
4551 | #endif |
---|
4552 | realDualInfeasibilities = sumDualInfeasibilities_; |
---|
4553 | double saveTolerance = dualTolerance_; |
---|
4554 | // If we need to carry on cleaning variables |
---|
4555 | if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) { |
---|
4556 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
4557 | int iPivot = pivotVariable_[iRow]; |
---|
4558 | if (!flagged(iPivot) && pivoted(iPivot)) { |
---|
4559 | // carry on |
---|
4560 | numberPrimalInfeasibilities_ = -1; |
---|
4561 | sumOfRelaxedPrimalInfeasibilities_ = 1.0; |
---|
4562 | sumPrimalInfeasibilities_ = 1.0; |
---|
4563 | break; |
---|
4564 | } |
---|
4565 | } |
---|
4566 | } |
---|
4567 | /* If we are primal feasible and any dual infeasibilities are on |
---|
4568 | free variables then it is better to go to primal */ |
---|
4569 | if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ && |
---|
4570 | numberDualInfeasibilities_) |
---|
4571 | problemStatus_ = 10; |
---|
4572 | // dual bound coming in |
---|
4573 | //double saveDualBound = dualBound_; |
---|
4574 | bool needCleanFake = false; |
---|
4575 | while (problemStatus_ <= -3) { |
---|
4576 | int cleanDuals = 0; |
---|
4577 | if (situationChanged != 0) |
---|
4578 | cleanDuals = 1; |
---|
4579 | int numberChangedBounds = 0; |
---|
4580 | int doOriginalTolerance = 0; |
---|
4581 | if ( lastCleaned == numberIterations_) |
---|
4582 | doOriginalTolerance = 1; |
---|
4583 | // check optimal |
---|
4584 | // give code benefit of doubt |
---|
4585 | if (sumOfRelaxedDualInfeasibilities_ == 0.0 && |
---|
4586 | sumOfRelaxedPrimalInfeasibilities_ == 0.0) { |
---|
4587 | // say optimal (with these bounds etc) |
---|
4588 | numberDualInfeasibilities_ = 0; |
---|
4589 | sumDualInfeasibilities_ = 0.0; |
---|
4590 | numberPrimalInfeasibilities_ = 0; |
---|
4591 | sumPrimalInfeasibilities_ = 0.0; |
---|
4592 | } |
---|
4593 | //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) { |
---|
4594 | if (dualFeasible() || problemStatus_ == -4) { |
---|
4595 | progress_.modifyObjective(objectiveValue_ |
---|
4596 | - bestPossibleImprovement_); |
---|
4597 | #ifdef COIN_DEVELOP |
---|
4598 | if (sumDualInfeasibilities_ || bestPossibleImprovement_) |
---|
4599 | printf("improve %g dualinf %g -> %g\n", |
---|
4600 | bestPossibleImprovement_, sumDualInfeasibilities_, |
---|
4601 | sumDualInfeasibilities_ * dualBound_); |
---|
4602 | #endif |
---|
4603 | // see if cutoff reached |
---|
4604 | double limit = 0.0; |
---|
4605 | getDblParam(ClpDualObjectiveLimit, limit); |
---|
4606 | #if 0 |
---|
4607 | if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ > |
---|
4608 | limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) { |
---|
4609 | //looks infeasible on objective |
---|
4610 | if (perturbation_ == 101) { |
---|
4611 | cleanDuals = 1; |
---|
4612 | // Save costs |
---|
4613 | int numberTotal = numberRows_ + numberColumns_; |
---|
4614 | double * saveCost = CoinCopyOfArray(cost_, numberTotal); |
---|
4615 | // make sure fake bounds are back |
---|
4616 | changeBounds(1, NULL, changeCost); |
---|
4617 | createRim4(false); |
---|
4618 | // make sure duals are current |
---|
4619 | computeDuals(givenDuals); |
---|
4620 | checkDualSolution(); |
---|
4621 | if(objectiveValue()*optimizationDirection_ > |
---|
4622 | limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) { |
---|
4623 | perturbation_ = 102; // stop any perturbations |
---|
4624 | printf("cutoff test succeeded\n"); |
---|
4625 | } else { |
---|
4626 | printf("cutoff test failed\n"); |
---|
4627 | // put back |
---|
4628 | memcpy(cost_, saveCost, numberTotal * sizeof(double)); |
---|
4629 | // make sure duals are current |
---|
4630 | computeDuals(givenDuals); |
---|
4631 | checkDualSolution(); |
---|
4632 | progress_.modifyObjective(-COIN_DBL_MAX); |
---|
4633 | problemStatus_ = -1; |
---|
4634 | } |
---|
4635 | delete [] saveCost; |
---|
4636 | } |
---|
4637 | } |
---|
4638 | #endif |
---|
4639 | if (primalFeasible() && !givenDuals) { |
---|
4640 | // may be optimal - or may be bounds are wrong |
---|
4641 | handler_->message(CLP_DUAL_BOUNDS, messages_) |
---|
4642 | << dualBound_ |
---|
4643 | << CoinMessageEol; |
---|
4644 | // save solution in case unbounded |
---|
4645 | double * saveColumnSolution = NULL; |
---|
4646 | double * saveRowSolution = NULL; |
---|
4647 | bool inCbc = (specialOptions_ & (0x01000000 | 16384)) != 0; |
---|
4648 | if (!inCbc) { |
---|
4649 | saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_); |
---|
4650 | saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_); |
---|
4651 | } |
---|
4652 | numberChangedBounds = changeBounds(0, rowArray_[3], changeCost); |
---|
4653 | if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) { |
---|
4654 | //looks optimal - do we need to reset tolerance |
---|
4655 | if (perturbation_ == 101) { |
---|
4656 | perturbation_ = 102; // stop any perturbations |
---|
4657 | cleanDuals = 1; |
---|
4658 | // make sure fake bounds are back |
---|
4659 | //computeObjectiveValue(); |
---|
4660 | changeBounds(1, NULL, changeCost); |
---|
4661 | //computeObjectiveValue(); |
---|
4662 | createRim4(false); |
---|
4663 | // make sure duals are current |
---|
4664 | computeDuals(givenDuals); |
---|
4665 | checkDualSolution(); |
---|
4666 | progress_.modifyObjective(-COIN_DBL_MAX); |
---|
4667 | #define DUAL_TRY_FASTER |
---|
4668 | #ifdef DUAL_TRY_FASTER |
---|
4669 | if (numberDualInfeasibilities_) { |
---|
4670 | #endif |
---|
4671 | numberChanged_ = 1; // force something to happen |
---|
4672 | lastCleaned = numberIterations_ - 1; |
---|
4673 | #ifdef DUAL_TRY_FASTER |
---|
4674 | } else { |
---|
4675 | //double value = objectiveValue_; |
---|
4676 | computeObjectiveValue(true); |
---|
4677 | //printf("old %g new %g\n",value,objectiveValue_); |
---|
4678 | //numberChanged_=1; |
---|
4679 | } |
---|
4680 | #endif |
---|
4681 | } |
---|
4682 | if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 && |
---|
4683 | (numberChanged_ || (specialOptions_ & 4096) == 0)) { |
---|
4684 | #if CLP_CAN_HAVE_ZERO_OBJ |
---|
4685 | if ((specialOptions_&2097152)==0) { |
---|
4686 | #endif |
---|
4687 | doOriginalTolerance = 2; |
---|
4688 | numberTimesOptimal_++; |
---|
4689 | changeMade_++; // say something changed |
---|
4690 | if (numberTimesOptimal_ == 1) { |
---|
4691 | dualTolerance_ = dblParam_[ClpDualTolerance]; |
---|
4692 | } else { |
---|
4693 | if (numberTimesOptimal_ == 2) { |
---|
4694 | // better to have small tolerance even if slower |
---|
4695 | factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15)); |
---|
4696 | } |
---|
4697 | dualTolerance_ = dblParam_[ClpDualTolerance]; |
---|
4698 | dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1); |
---|
4699 | } |
---|
4700 | cleanDuals = 2; // If nothing changed optimal else primal |
---|
4701 | #if CLP_CAN_HAVE_ZERO_OBJ |
---|
4702 | } else { |
---|
4703 | // no cost - skip checks |
---|
4704 | problemStatus_=0; |
---|
4705 | } |
---|
4706 | #endif |
---|
4707 | } else { |
---|
4708 | problemStatus_ = 0; // optimal |
---|
4709 | if (lastCleaned < numberIterations_ && numberChanged_) { |
---|
4710 | handler_->message(CLP_SIMPLEX_GIVINGUP, messages_) |
---|
4711 | << CoinMessageEol; |
---|
4712 | } |
---|
4713 | } |
---|
4714 | } else { |
---|
4715 | cleanDuals = 1; |
---|
4716 | if (doOriginalTolerance == 1) { |
---|
4717 | // check unbounded |
---|
4718 | // find a variable with bad dj |
---|
4719 | int iSequence; |
---|
4720 | int iChosen = -1; |
---|
4721 | if (!inCbc) { |
---|
4722 | double largest = 100.0 * primalTolerance_; |
---|
4723 | for (iSequence = 0; iSequence < numberRows_ + numberColumns_; |
---|
4724 | iSequence++) { |
---|
4725 | double djValue = dj_[iSequence]; |
---|
4726 | double originalLo = originalLower(iSequence); |
---|
4727 | double originalUp = originalUpper(iSequence); |
---|
4728 | if (fabs(djValue) > fabs(largest)) { |
---|
4729 | if (getStatus(iSequence) != basic) { |
---|
4730 | if (djValue > 0 && originalLo < -1.0e20) { |
---|
4731 | if (djValue > fabs(largest)) { |
---|
4732 | largest = djValue; |
---|
4733 | iChosen = iSequence; |
---|
4734 | } |
---|
4735 | } else if (djValue < 0 && originalUp > 1.0e20) { |
---|
4736 | if (-djValue > fabs(largest)) { |
---|
4737 | largest = djValue; |
---|
4738 | iChosen = iSequence; |
---|
4739 | } |
---|
4740 | } |
---|
4741 | } |
---|
4742 | } |
---|
4743 | } |
---|
4744 | } |
---|
4745 | if (iChosen >= 0) { |
---|
4746 | int iSave = sequenceIn_; |
---|
4747 | sequenceIn_ = iChosen; |
---|
4748 | unpack(rowArray_[1]); |
---|
4749 | sequenceIn_ = iSave; |
---|
4750 | // if dual infeasibilities then must be free vector so add in dual |
---|
4751 | if (numberDualInfeasibilities_) { |
---|
4752 | if (fabs(changeCost) > 1.0e-5) |
---|
4753 | COIN_DETAIL_PRINT(printf("Odd free/unbounded combo\n")); |
---|
4754 | changeCost += cost_[iChosen]; |
---|
4755 | } |
---|
4756 | problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0], |
---|
4757 | changeCost); |
---|
4758 | rowArray_[1]->clear(); |
---|
4759 | } else { |
---|
4760 | problemStatus_ = -3; |
---|
4761 | } |
---|
4762 | if (problemStatus_ == 2 && perturbation_ == 101) { |
---|
4763 | perturbation_ = 102; // stop any perturbations |
---|
4764 | cleanDuals = 1; |
---|
4765 | createRim4(false); |
---|
4766 | progress_.modifyObjective(-COIN_DBL_MAX); |
---|
4767 | problemStatus_ = -1; |
---|
4768 | } |
---|
4769 | if (problemStatus_ == 2) { |
---|
4770 | // it is unbounded - restore solution |
---|
4771 | // but first add in changes to non-basic |
---|
4772 | int iColumn; |
---|
4773 | double * original = columnArray_[0]->denseVector(); |
---|
4774 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
4775 | if(getColumnStatus(iColumn) != basic) |
---|
4776 | ray_[iColumn] += |
---|
4777 | saveColumnSolution[iColumn] - original[iColumn]; |
---|
4778 | columnActivityWork_[iColumn] = original[iColumn]; |
---|
4779 | } |
---|
4780 | CoinMemcpyN(saveRowSolution, numberRows_, |
---|
4781 | rowActivityWork_); |
---|
4782 | } |
---|
4783 | } else { |
---|
4784 | doOriginalTolerance = 2; |
---|
4785 | rowArray_[0]->clear(); |
---|
4786 | } |
---|
4787 | } |
---|
4788 | delete [] saveColumnSolution; |
---|
4789 | delete [] saveRowSolution; |
---|
4790 | } |
---|
4791 | if (problemStatus_ == -4 || problemStatus_ == -5) { |
---|
4792 | // may be infeasible - or may be bounds are wrong |
---|
4793 | numberChangedBounds = changeBounds(0, NULL, changeCost); |
---|
4794 | needCleanFake = true; |
---|
4795 | /* Should this be here as makes no difference to being feasible. |
---|
4796 | But seems to make a difference to run times. */ |
---|
4797 | if (perturbation_ == 101 && 0) { |
---|
4798 | perturbation_ = 102; // stop any perturbations |
---|
4799 | cleanDuals = 1; |
---|
4800 | numberChangedBounds = 1; |
---|
4801 | // make sure fake bounds are back |
---|
4802 | changeBounds(1, NULL, changeCost); |
---|
4803 | needCleanFake = true; |
---|
4804 | createRim4(false); |
---|
4805 | progress_.modifyObjective(-COIN_DBL_MAX); |
---|
4806 | } |
---|
4807 | if ((numberChangedBounds <= 0 || dualBound_ > 1.0e20 || |
---|
4808 | (largestPrimalError_ > 1.0 && dualBound_ > 1.0e17)) && |
---|
4809 | (numberPivots < 4 || sumPrimalInfeasibilities_ > 1.0e-6)) { |
---|
4810 | problemStatus_ = 1; // infeasible |
---|
4811 | if (perturbation_ == 101) { |
---|
4812 | perturbation_ = 102; // stop any perturbations |
---|
4813 | //cleanDuals=1; |
---|
4814 | //numberChangedBounds=1; |
---|
4815 | //createRim4(false); |
---|
4816 | } |
---|
4817 | } else { |
---|
4818 | problemStatus_ = -1; //iterate |
---|
4819 | cleanDuals = 1; |
---|
4820 | if (numberChangedBounds <= 0) |
---|
4821 | doOriginalTolerance = 2; |
---|
4822 | // and delete ray which has been created |
---|
4823 | delete [] ray_; |
---|
4824 | ray_ = NULL; |
---|
4825 | } |
---|
4826 | |
---|
4827 | } |
---|
4828 | } else { |
---|
4829 | cleanDuals = 1; |
---|
4830 | } |
---|
4831 | if (problemStatus_ < 0) { |
---|
4832 | if (doOriginalTolerance == 2) { |
---|
4833 | // put back original tolerance |
---|
4834 | lastCleaned = numberIterations_; |
---|
4835 | numberChanged_ = 0; // Number of variables with changed costs |
---|
4836 | handler_->message(CLP_DUAL_ORIGINAL, messages_) |
---|
4837 | << CoinMessageEol; |
---|
4838 | perturbation_ = 102; // stop any perturbations |
---|
4839 | #if 0 |
---|
4840 | double * xcost = new double[numberRows_+numberColumns_]; |
---|
4841 | double * xlower = new double[numberRows_+numberColumns_]; |
---|
4842 | double * xupper = new double[numberRows_+numberColumns_]; |
---|
4843 | double * xdj = new double[numberRows_+numberColumns_]; |
---|
4844 | double * xsolution = new double[numberRows_+numberColumns_]; |
---|
4845 | CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost); |
---|
4846 | CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower); |
---|
4847 | CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper); |
---|
4848 | CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj); |
---|
4849 | CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution); |
---|
4850 | #endif |
---|
4851 | createRim4(false); |
---|
4852 | progress_.modifyObjective(-COIN_DBL_MAX); |
---|
4853 | // make sure duals are current |
---|
4854 | computeDuals(givenDuals); |
---|
4855 | checkDualSolution(); |
---|
4856 | #if 0 |
---|
4857 | int i; |
---|
4858 | for (i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
4859 | if (cost_[i] != xcost[i]) |
---|
4860 | printf("** %d old cost %g new %g sol %g\n", |
---|
4861 | i, xcost[i], cost_[i], solution_[i]); |
---|
4862 | if (lower_[i] != xlower[i]) |
---|
4863 | printf("** %d old lower %g new %g sol %g\n", |
---|
4864 | i, xlower[i], lower_[i], solution_[i]); |
---|
4865 | if (upper_[i] != xupper[i]) |
---|
4866 | printf("** %d old upper %g new %g sol %g\n", |
---|
4867 | i, xupper[i], upper_[i], solution_[i]); |
---|
4868 | if (dj_[i] != xdj[i]) |
---|
4869 | printf("** %d old dj %g new %g sol %g\n", |
---|
4870 | i, xdj[i], dj_[i], solution_[i]); |
---|
4871 | if (solution_[i] != xsolution[i]) |
---|
4872 | printf("** %d old solution %g new %g sol %g\n", |
---|
4873 | i, xsolution[i], solution_[i], solution_[i]); |
---|
4874 | } |
---|
4875 | //delete [] xcost; |
---|
4876 | //delete [] xupper; |
---|
4877 | //delete [] xlower; |
---|
4878 | //delete [] xdj; |
---|
4879 | //delete [] xsolution; |
---|
4880 | #endif |
---|
4881 | // put back bounds as they were if was optimal |
---|
4882 | if (doOriginalTolerance == 2 && cleanDuals != 2) { |
---|
4883 | changeMade_++; // say something changed |
---|
4884 | /* We may have already changed some bounds in this function |
---|
4885 | so save numberFake_ and add in. |
---|
4886 | |
---|
4887 | Worst that can happen is that we waste a bit of time - but it must be finite. |
---|
4888 | */ |
---|
4889 | //int saveNumberFake = numberFake_; |
---|
4890 | //resetFakeBounds(-1); |
---|
4891 | changeBounds(3, NULL, changeCost); |
---|
4892 | needCleanFake = true; |
---|
4893 | //numberFake_ += saveNumberFake; |
---|
4894 | //resetFakeBounds(-1); |
---|
4895 | cleanDuals = 2; |
---|
4896 | //cleanDuals=1; |
---|
4897 | } |
---|
4898 | #if 0 |
---|
4899 | //int i; |
---|
4900 | for (i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
4901 | if (cost_[i] != xcost[i]) |
---|
4902 | printf("** %d old cost %g new %g sol %g\n", |
---|
4903 | i, xcost[i], cost_[i], solution_[i]); |
---|
4904 | if (lower_[i] != xlower[i]) |
---|
4905 | printf("** %d old lower %g new %g sol %g\n", |
---|
4906 | i, xlower[i], lower_[i], solution_[i]); |
---|
4907 | if (upper_[i] != xupper[i]) |
---|
4908 | printf("** %d old upper %g new %g sol %g\n", |
---|
4909 | i, xupper[i], upper_[i], solution_[i]); |
---|
4910 | if (dj_[i] != xdj[i]) |
---|
4911 | printf("** %d old dj %g new %g sol %g\n", |
---|
4912 | i, xdj[i], dj_[i], solution_[i]); |
---|
4913 | if (solution_[i] != xsolution[i]) |
---|
4914 | printf("** %d old solution %g new %g sol %g\n", |
---|
4915 | i, xsolution[i], solution_[i], solution_[i]); |
---|
4916 | } |
---|
4917 | delete [] xcost; |
---|
4918 | delete [] xupper; |
---|
4919 | delete [] xlower; |
---|
4920 | delete [] xdj; |
---|
4921 | delete [] xsolution; |
---|
4922 | #endif |
---|
4923 | } |
---|
4924 | if (cleanDuals == 1 || (cleanDuals == 2 && !numberDualInfeasibilities_)) { |
---|
4925 | // make sure dual feasible |
---|
4926 | // look at all rows and columns |
---|
4927 | rowArray_[0]->clear(); |
---|
4928 | columnArray_[0]->clear(); |
---|
4929 | double objectiveChange = 0.0; |
---|
4930 | double savePrimalInfeasibilities = sumPrimalInfeasibilities_; |
---|
4931 | if (!numberIterations_) { |
---|
4932 | int nTotal = numberRows_ + numberColumns_; |
---|
4933 | if (arraysNotCreated) { |
---|
4934 | // create save arrays |
---|
4935 | delete [] saveStatus_; |
---|
4936 | delete [] savedSolution_; |
---|
4937 | saveStatus_ = new unsigned char [nTotal]; |
---|
4938 | savedSolution_ = new double [nTotal]; |
---|
4939 | arraysNotCreated = false; |
---|
4940 | } |
---|
4941 | // save arrays |
---|
4942 | CoinMemcpyN(status_, nTotal, saveStatus_); |
---|
4943 | CoinMemcpyN(rowActivityWork_, |
---|
4944 | numberRows_, savedSolution_ + numberColumns_); |
---|
4945 | CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_); |
---|
4946 | } |
---|
4947 | #if 0 |
---|
4948 | double * xcost = new double[numberRows_+numberColumns_]; |
---|
4949 | double * xlower = new double[numberRows_+numberColumns_]; |
---|
4950 | double * xupper = new double[numberRows_+numberColumns_]; |
---|
4951 | double * xdj = new double[numberRows_+numberColumns_]; |
---|
4952 | double * xsolution = new double[numberRows_+numberColumns_]; |
---|
4953 | CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost); |
---|
4954 | CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower); |
---|
4955 | CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper); |
---|
4956 | CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj); |
---|
4957 | CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution); |
---|
4958 | #endif |
---|
4959 | if (givenDuals) |
---|
4960 | dualTolerance_ = 1.0e50; |
---|
4961 | #if CLP_CAN_HAVE_ZERO_OBJ>1 |
---|
4962 | if ((specialOptions_&2097152)==0) { |
---|
4963 | #endif |
---|
4964 | updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1], |
---|
4965 | 0.0, objectiveChange, true); |
---|
4966 | #if CLP_CAN_HAVE_ZERO_OBJ>1 |
---|
4967 | } else { |
---|
4968 | rowArray_[0]->clear(); |
---|
4969 | rowArray_[1]->clear(); |
---|
4970 | columnArray_[0]->clear(); |
---|
4971 | } |
---|
4972 | #endif |
---|
4973 | dualTolerance_ = saveTolerance; |
---|
4974 | #if 0 |
---|
4975 | int i; |
---|
4976 | for (i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
4977 | if (cost_[i] != xcost[i]) |
---|
4978 | printf("** %d old cost %g new %g sol %g\n", |
---|
4979 | i, xcost[i], cost_[i], solution_[i]); |
---|
4980 | if (lower_[i] != xlower[i]) |
---|
4981 | printf("** %d old lower %g new %g sol %g\n", |
---|
4982 | i, xlower[i], lower_[i], solution_[i]); |
---|
4983 | if (upper_[i] != xupper[i]) |
---|
4984 | printf("** %d old upper %g new %g sol %g\n", |
---|
4985 | i, xupper[i], upper_[i], solution_[i]); |
---|
4986 | if (dj_[i] != xdj[i]) |
---|
4987 | printf("** %d old dj %g new %g sol %g\n", |
---|
4988 | i, xdj[i], dj_[i], solution_[i]); |
---|
4989 | if (solution_[i] != xsolution[i]) |
---|
4990 | printf("** %d old solution %g new %g sol %g\n", |
---|
4991 | i, xsolution[i], solution_[i], solution_[i]); |
---|
4992 | } |
---|
4993 | delete [] xcost; |
---|
4994 | delete [] xupper; |
---|
4995 | delete [] xlower; |
---|
4996 | delete [] xdj; |
---|
4997 | delete [] xsolution; |
---|
4998 | #endif |
---|
4999 | // for now - recompute all |
---|
5000 | gutsOfSolution(NULL, NULL); |
---|
5001 | if (givenDuals) |
---|
5002 | dualTolerance_ = 1.0e50; |
---|
5003 | #if CLP_CAN_HAVE_ZERO_OBJ>1 |
---|
5004 | if ((specialOptions_&2097152)==0) { |
---|
5005 | #endif |
---|
5006 | updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1], |
---|
5007 | 0.0, objectiveChange, true); |
---|
5008 | #if CLP_CAN_HAVE_ZERO_OBJ>1 |
---|
5009 | } else { |
---|
5010 | rowArray_[0]->clear(); |
---|
5011 | rowArray_[1]->clear(); |
---|
5012 | columnArray_[0]->clear(); |
---|
5013 | } |
---|
5014 | #endif |
---|
5015 | dualTolerance_ = saveTolerance; |
---|
5016 | if (!numberIterations_ && sumPrimalInfeasibilities_ > |
---|
5017 | 1.0e5*(savePrimalInfeasibilities+1.0e3) && |
---|
5018 | (moreSpecialOptions_ & 256) == 0) { |
---|
5019 | // Use primal |
---|
5020 | int nTotal = numberRows_ + numberColumns_; |
---|
5021 | CoinMemcpyN(saveStatus_, nTotal, status_); |
---|
5022 | CoinMemcpyN(savedSolution_ + numberColumns_ , |
---|
5023 | numberRows_, rowActivityWork_); |
---|
5024 | CoinMemcpyN(savedSolution_ , |
---|
5025 | numberColumns_, columnActivityWork_); |
---|
5026 | problemStatus_ = 10; |
---|
5027 | situationChanged = 0; |
---|
5028 | } |
---|
5029 | //assert(numberDualInfeasibilitiesWithoutFree_==0); |
---|
5030 | if (numberDualInfeasibilities_) { |
---|
5031 | if ((numberPrimalInfeasibilities_ || numberPivots) |
---|
5032 | && problemStatus_!=10) { |
---|
5033 | problemStatus_ = -1; // carry on as normal |
---|
5034 | } else { |
---|
5035 | problemStatus_ = 10; // try primal |
---|
5036 | #if COIN_DEVELOP>1 |
---|
5037 | printf("returning at %d\n", __LINE__); |
---|
5038 | #endif |
---|
5039 | } |
---|
5040 | } else if (situationChanged == 2) { |
---|
5041 | problemStatus_ = -1; // carry on as normal |
---|
5042 | // need to reset bounds |
---|
5043 | changeBounds(3, NULL, changeCost); |
---|
5044 | } |
---|
5045 | situationChanged = 0; |
---|
5046 | } else { |
---|
5047 | // iterate |
---|
5048 | if (cleanDuals != 2) { |
---|
5049 | problemStatus_ = -1; |
---|
5050 | } else { |
---|
5051 | problemStatus_ = 10; // try primal |
---|
5052 | #if COIN_DEVELOP>2 |
---|
5053 | printf("returning at %d\n", __LINE__); |
---|
5054 | #endif |
---|
5055 | } |
---|
5056 | } |
---|
5057 | } |
---|
5058 | } |
---|
5059 | // unflag all variables (we may want to wait a bit?) |
---|
5060 | if ((tentativeStatus != -2 && tentativeStatus != -1) && unflagVariables) { |
---|
5061 | int iRow; |
---|
5062 | int numberFlagged = 0; |
---|
5063 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
5064 | int iPivot = pivotVariable_[iRow]; |
---|
5065 | if (flagged(iPivot)) { |
---|
5066 | numberFlagged++; |
---|
5067 | clearFlagged(iPivot); |
---|
5068 | } |
---|
5069 | } |
---|
5070 | #ifdef COIN_DEVELOP |
---|
5071 | if (numberFlagged) { |
---|
5072 | printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus, |
---|
5073 | problemStatus_, numberPrimalInfeasibilities_, |
---|
5074 | numberTimesOptimal_); |
---|
5075 | } |
---|
5076 | #endif |
---|
5077 | unflagVariables = numberFlagged > 0; |
---|
5078 | if (numberFlagged && !numberPivots) { |
---|
5079 | /* looks like trouble as we have not done any iterations. |
---|
5080 | Try changing pivot tolerance then give it a few goes and give up */ |
---|
5081 | if (factorization_->pivotTolerance() < 0.9) { |
---|
5082 | factorization_->pivotTolerance(0.99); |
---|
5083 | problemStatus_ = -1; |
---|
5084 | } else if (numberTimesOptimal_ < 3) { |
---|
5085 | numberTimesOptimal_++; |
---|
5086 | problemStatus_ = -1; |
---|
5087 | } else { |
---|
5088 | unflagVariables = false; |
---|
5089 | //secondaryStatus_ = 1; // and say probably infeasible |
---|
5090 | if ((moreSpecialOptions_ & 256) == 0) { |
---|
5091 | // try primal |
---|
5092 | problemStatus_ = 10; |
---|
5093 | } else { |
---|
5094 | // almost certainly infeasible |
---|
5095 | problemStatus_ = 1; |
---|
5096 | } |
---|
5097 | #if COIN_DEVELOP>1 |
---|
5098 | printf("returning at %d\n", __LINE__); |
---|
5099 | #endif |
---|
5100 | } |
---|
5101 | } |
---|
5102 | } |
---|
5103 | if (problemStatus_ < 0) { |
---|
5104 | if (needCleanFake) { |
---|
5105 | double dummyChangeCost = 0.0; |
---|
5106 | changeBounds(3, NULL, dummyChangeCost); |
---|
5107 | } |
---|
5108 | #if 0 |
---|
5109 | if (objectiveValue_ < lastObjectiveValue_ - 1.0e-8 * |
---|
5110 | CoinMax(fabs(objectivevalue_), fabs(lastObjectiveValue_))) { |
---|
5111 | } else { |
---|
5112 | lastObjectiveValue_ = objectiveValue_; |
---|
5113 | } |
---|
5114 | #endif |
---|
5115 | if (type == 0 || type == 1) { |
---|
5116 | if (!type && arraysNotCreated) { |
---|
5117 | // create save arrays |
---|
5118 | delete [] saveStatus_; |
---|
5119 | delete [] savedSolution_; |
---|
5120 | saveStatus_ = new unsigned char [numberRows_+numberColumns_]; |
---|
5121 | savedSolution_ = new double [numberRows_+numberColumns_]; |
---|
5122 | } |
---|
5123 | // save arrays |
---|
5124 | CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_); |
---|
5125 | CoinMemcpyN(rowActivityWork_, |
---|
5126 | numberRows_, savedSolution_ + numberColumns_); |
---|
5127 | CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_); |
---|
5128 | // save extra stuff |
---|
5129 | int dummy; |
---|
5130 | matrix_->generalExpanded(this, 5, dummy); |
---|
5131 | } |
---|
5132 | if (weightsSaved) { |
---|
5133 | // restore weights (if saved) - also recompute infeasibility list |
---|
5134 | if (!reallyBadProblems && (largestPrimalError_ < 100.0 || numberPivots > 10)) { |
---|
5135 | if (tentativeStatus > -3) |
---|
5136 | dualRowPivot_->saveWeights(this, (type < 2) ? 2 : 4); |
---|
5137 | else |
---|
5138 | dualRowPivot_->saveWeights(this, 3); |
---|
5139 | } else { |
---|
5140 | // reset weights or scale back |
---|
5141 | dualRowPivot_->saveWeights(this, 6); |
---|
5142 | } |
---|
5143 | } else if (weightsSaved2 && numberPrimalInfeasibilities_) { |
---|
5144 | dualRowPivot_->saveWeights(this, 3); |
---|
5145 | } |
---|
5146 | } |
---|
5147 | // see if cutoff reached |
---|
5148 | double limit = 0.0; |
---|
5149 | getDblParam(ClpDualObjectiveLimit, limit); |
---|
5150 | #if 0 |
---|
5151 | if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ > |
---|
5152 | limit + 100.0) { |
---|
5153 | printf("lim %g obj %g %g - wo perturb %g sum dual %g\n", |
---|
5154 | limit, objectiveValue_, objectiveValue(), computeInternalObjectiveValue(), sumDualInfeasibilities_); |
---|
5155 | } |
---|
5156 | #endif |
---|
5157 | if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ > |
---|
5158 | limit && !numberAtFakeBound()) { |
---|
5159 | bool looksInfeasible = !numberDualInfeasibilities_; |
---|
5160 | if (objectiveValue()*optimizationDirection_ > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 && |
---|
5161 | sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000) |
---|
5162 | looksInfeasible = true; |
---|
5163 | if (looksInfeasible) { |
---|
5164 | // Even if not perturbed internal costs may have changed |
---|
5165 | // be careful |
---|
5166 | if (true || numberIterations_) { |
---|
5167 | if(computeInternalObjectiveValue() > limit) { |
---|
5168 | problemStatus_ = 1; |
---|
5169 | secondaryStatus_ = 1; // and say was on cutoff |
---|
5170 | } |
---|
5171 | } else { |
---|
5172 | problemStatus_ = 1; |
---|
5173 | secondaryStatus_ = 1; // and say was on cutoff |
---|
5174 | } |
---|
5175 | } |
---|
5176 | } |
---|
5177 | // If we are in trouble and in branch and bound give up |
---|
5178 | if ((specialOptions_ & 1024) != 0) { |
---|
5179 | int looksBad = 0; |
---|
5180 | if (largestPrimalError_ * largestDualError_ > 1.0e2) { |
---|
5181 | looksBad = 1; |
---|
5182 | } else if (largestPrimalError_ > 1.0e-2 |
---|
5183 | && objectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) { |
---|
5184 | looksBad = 2; |
---|
5185 | } |
---|
5186 | if (looksBad) { |
---|
5187 | if (factorization_->pivotTolerance() < 0.9) { |
---|
5188 | // up tolerance |
---|
5189 | factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91)); |
---|
5190 | } else if (numberIterations_ > 10000) { |
---|
5191 | if (handler_->logLevel() > 2) |
---|
5192 | printf("bad dual - saying infeasible %d\n", looksBad); |
---|
5193 | problemStatus_ = 1; |
---|
5194 | secondaryStatus_ = 1; // and say was on cutoff |
---|
5195 | } else if (largestPrimalError_ > 1.0e5) { |
---|
5196 | { |
---|
5197 | int iBigB = -1; |
---|
5198 | double bigB = 0.0; |
---|
5199 | int iBigN = -1; |
---|
5200 | double bigN = 0.0; |
---|
5201 | for (int i = 0; i < numberRows_ + numberColumns_; i++) { |
---|
5202 | double value = fabs(solution_[i]); |
---|
5203 | if (getStatus(i) == basic) { |
---|
5204 | if (value > bigB) { |
---|
5205 | bigB = value; |
---|
5206 | iBigB = i; |
---|
5207 | } |
---|
5208 | } else { |
---|
5209 | if (value > bigN) { |
---|
5210 | bigN = value; |
---|
5211 | iBigN = i; |
---|
5212 | } |
---|
5213 | } |
---|
5214 | } |
---|
5215 | #ifdef CLP_INVESTIGATE |
---|
5216 | if (bigB > 1.0e8 || bigN > 1.0e8) { |
---|
5217 | if (handler_->logLevel() > 0) |
---|
5218 | printf("it %d - basic %d %g, nonbasic %d %g\n", |
---|
5219 | numberIterations_, iBigB, bigB, iBigN, bigN); |
---|
5220 | } |
---|
5221 | #endif |
---|
5222 | } |
---|
5223 | #if COIN_DEVELOP!=2 |
---|
5224 | if (handler_->logLevel() > 2) |
---|
5225 | #endif |
---|
5226 | printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_); |
---|
5227 | allSlackBasis(true); |
---|
5228 | problemStatus_ = 10; |
---|
5229 | } |
---|
5230 | } |
---|
5231 | } |
---|
5232 | if (problemStatus_ < 0 && !changeMade_) { |
---|
5233 | problemStatus_ = 4; // unknown |
---|
5234 | } |
---|
5235 | lastGoodIteration_ = numberIterations_; |
---|
5236 | if (numberIterations_ > lastBadIteration_ + 100) |
---|
5237 | moreSpecialOptions_ &= ~16; // clear check accuracy flag |
---|
5238 | if (problemStatus_ < 0) { |
---|
5239 | sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful |
---|
5240 | if (sumDualInfeasibilities_) |
---|
5241 | numberDualInfeasibilities_ = 1; |
---|
5242 | } |
---|
5243 | #ifdef CLP_REPORT_PROGRESS |
---|
5244 | if (ixxxxxx > ixxyyyy - 3) { |
---|
5245 | printf("objectiveValue_ %g\n", objectiveValue_); |
---|
5246 | handler_->setLogLevel(63); |
---|
5247 | int nTotal = numberColumns_ + numberRows_; |
---|
5248 | double newObj = 0.0; |
---|
5249 | for (int i = 0; i < nTotal; i++) { |
---|
5250 | if (solution_[i]) |
---|
5251 | newObj += solution_[i] * cost_[i]; |
---|
5252 | } |
---|
5253 | printf("xxx obj %g\n", newObj); |
---|
5254 | // for now - recompute all |
---|
5255 | gutsOfSolution(NULL, NULL); |
---|
5256 | newObj = 0.0; |
---|
5257 | for (int i =< |
---|