1/* $Id: ClpPredictorCorrector.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2003, 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 Implements crude primal dual predictor corrector algorithm
8
9 */
10//#define SOME_DEBUG
11
12#include "CoinPragma.hpp"
13#include <math.h>
14
15#include "CoinHelperFunctions.hpp"
16#include "ClpPredictorCorrector.hpp"
17#include "CoinPackedMatrix.hpp"
18#include "ClpMessage.hpp"
19#include "ClpCholeskyBase.hpp"
20#include "ClpHelperFunctions.hpp"
21#include "ClpQuadraticObjective.hpp"
22#include <cfloat>
23#include <cassert>
24#include <string>
25#include <cstdio>
26#include <iostream>
27#if 0
28static int yyyyyy = 0;
29void ClpPredictorCorrector::saveSolution(std::string fileName)
30{
31 FILE * fp = fopen(fileName.c_str(), "wb");
32 if (fp) {
33 int numberRows = numberRows_;
34 int numberColumns = numberColumns_;
35 fwrite(&numberRows, sizeof(int), 1, fp);
36 fwrite(&numberColumns, sizeof(int), 1, fp);
37 CoinWorkDouble dsave[20];
38 memset(dsave, 0, sizeof(dsave));
39 fwrite(dsave, sizeof(CoinWorkDouble), 20, fp);
40 int msave[20];
41 memset(msave, 0, sizeof(msave));
42 msave[0] = numberIterations_;
43 fwrite(msave, sizeof(int), 20, fp);
44 fwrite(dual_, sizeof(CoinWorkDouble), numberRows, fp);
45 fwrite(errorRegion_, sizeof(CoinWorkDouble), numberRows, fp);
46 fwrite(rhsFixRegion_, sizeof(CoinWorkDouble), numberRows, fp);
47 fwrite(solution_, sizeof(CoinWorkDouble), numberColumns, fp);
48 fwrite(solution_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
49 fwrite(diagonal_, sizeof(CoinWorkDouble), numberColumns, fp);
50 fwrite(diagonal_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
51 fwrite(wVec_, sizeof(CoinWorkDouble), numberColumns, fp);
52 fwrite(wVec_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
53 fwrite(zVec_, sizeof(CoinWorkDouble), numberColumns, fp);
54 fwrite(zVec_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
55 fwrite(upperSlack_, sizeof(CoinWorkDouble), numberColumns, fp);
56 fwrite(upperSlack_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
57 fwrite(lowerSlack_, sizeof(CoinWorkDouble), numberColumns, fp);
58 fwrite(lowerSlack_ + numberColumns, sizeof(CoinWorkDouble), numberRows, fp);
59 fclose(fp);
60 } else {
61 std::cout << "Unable to open file " << fileName << std::endl;
62 }
63}
64#endif
65// Could change on CLP_LONG_CHOLESKY or COIN_LONG_WORK?
66static CoinWorkDouble eScale = 1.0e27;
67static CoinWorkDouble eBaseCaution = 1.0e-12;
68static CoinWorkDouble eBase = 1.0e-12;
69static CoinWorkDouble eRatio = 1.0e40;
70static CoinWorkDouble eRatioCaution = 1.0e25;
71static CoinWorkDouble eDiagonal = 1.0e25;
72static CoinWorkDouble eDiagonalCaution = 1.0e18;
73static CoinWorkDouble eExtra = 1.0e-12;
74
75// main function
76
77int ClpPredictorCorrector::solve ( )
78{
79 problemStatus_ = -1;
80 algorithm_ = 1;
81 //create all regions
82 if (!createWorkingData()) {
83 problemStatus_ = 4;
84 return 2;
85 }
86#if COIN_LONG_WORK
87 // reallocate some regions
88 double * dualSave = dual_;
89 dual_ = reinterpret_cast<double *>(new CoinWorkDouble[numberRows_]);
90 double * reducedCostSave = reducedCost_;
91 reducedCost_ = reinterpret_cast<double *>(new CoinWorkDouble[numberColumns_]);
92#endif
93 //diagonalPerturbation_=1.0e-25;
94 ClpMatrixBase * saveMatrix = NULL;
95 // If quadratic then make copy so we can actually scale or normalize
96#ifndef NO_RTTI
97 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
98#else
99 ClpQuadraticObjective * quadraticObj = NULL;
100 if (objective_->type() == 2)
101 quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
102#endif
103 /* If modeSwitch is :
104 0 - normal
105 1 - bit switch off centering
106 2 - bit always do type 2
107 4 - accept corrector nearly always
108 */
109 int modeSwitch = 0;
110 //if (quadraticObj)
111 //modeSwitch |= 1; // switch off centring for now
112 //if (quadraticObj)
113 //modeSwitch |=4;
114 ClpObjective * saveObjective = NULL;
115 if (quadraticObj) {
116 // check KKT is on
117 if (!cholesky_->kkt()) {
118 //No!
119 handler_->message(CLP_BARRIER_KKT, messages_)
120 << CoinMessageEol;
121 return -1;
122 }
123 saveObjective = objective_;
124 // We are going to make matrix full rather half
125 objective_ = new ClpQuadraticObjective(*quadraticObj, 1);
126 }
127 bool allowIncreasingGap = (modeSwitch & 4) != 0;
128 // If scaled then really scale matrix
129 if (scalingFlag_ > 0 && rowScale_) {
130 saveMatrix = matrix_;
131 matrix_ = matrix_->scaledColumnCopy(this);
132 }
133 //initializeFeasible(); - this just set fixed flag
134 smallestInfeasibility_ = COIN_DBL_MAX;
135 int i;
136 for (i = 0; i < LENGTH_HISTORY; i++)
137 historyInfeasibility_[i] = COIN_DBL_MAX;
138
139 //bool firstTime=true;
140 //firstFactorization(true);
141 int returnCode = cholesky_->order(this);
142 if (returnCode || cholesky_->symbolic()) {
143 COIN_DETAIL_PRINT(printf("Error return from symbolic - probably not enough memory\n"));
144 problemStatus_ = 4;
145 //delete all temporary regions
146 deleteWorkingData();
147 if (saveMatrix) {
148 // restore normal copy
149 delete matrix_;
150 matrix_ = saveMatrix;
151 }
152 // Restore quadratic objective if necessary
153 if (saveObjective) {
154 delete objective_;
155 objective_ = saveObjective;
156 }
157 return -1;
158 }
159 mu_ = 1.0e10;
160 diagonalScaleFactor_ = 1.0;
161 //set iterations
162 numberIterations_ = -1;
163 int numberTotal = numberRows_ + numberColumns_;
164 //initialize solution here
165 if(createSolution() < 0) {
166 COIN_DETAIL_PRINT(printf("Not enough memory\n"));
167 problemStatus_ = 4;
168 //delete all temporary regions
169 deleteWorkingData();
170 if (saveMatrix) {
171 // restore normal copy
172 delete matrix_;
173 matrix_ = saveMatrix;
174 }
175 return -1;
176 }
177 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
178 // Could try centering steps without any original step i.e. just center
179 //firstFactorization(false);
180 CoinZeroN(dualArray, numberRows_);
181 multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, errorRegion_, 0.0);
182 matrix_->times(1.0, solution_, errorRegion_);
183 maximumRHSError_ = maximumAbsElement(errorRegion_, numberRows_);
184 maximumBoundInfeasibility_ = maximumRHSError_;
185 //CoinWorkDouble maximumDualError_=COIN_DBL_MAX;
186 //initialize
187 actualDualStep_ = 0.0;
188 actualPrimalStep_ = 0.0;
189 gonePrimalFeasible_ = false;
190 goneDualFeasible_ = false;
191 //bool hadGoodSolution=false;
192 diagonalNorm_ = solutionNorm_;
193 mu_ = solutionNorm_;
194 int numberFixed = updateSolution(-COIN_DBL_MAX);
195 int numberFixedTotal = numberFixed;
196 //int numberRows_DroppedBefore=0;
197 //CoinWorkDouble extra=eExtra;
198 //CoinWorkDouble maximumPerturbation=COIN_DBL_MAX;
199 //constants for infeas interior point
200 const CoinWorkDouble beta2 = 0.99995;
201 const CoinWorkDouble tau = 0.00002;
202 CoinWorkDouble lastComplementarityGap = COIN_DBL_MAX * 1.0e-20;
203 CoinWorkDouble lastStep = 1.0;
204 // use to see if to take affine
205 CoinWorkDouble checkGap = COIN_DBL_MAX;
206 int lastGoodIteration = 0;
207 CoinWorkDouble bestObjectiveGap = COIN_DBL_MAX;
208 CoinWorkDouble bestObjective = COIN_DBL_MAX;
209 int bestKilled = -1;
210 int saveIteration = -1;
211 int saveIteration2 = -1;
212 bool sloppyOptimal = false;
213 CoinWorkDouble * savePi = NULL;
214 CoinWorkDouble * savePrimal = NULL;
215 CoinWorkDouble * savePi2 = NULL;
216 CoinWorkDouble * savePrimal2 = NULL;
217 // Extra regions for centering
218 CoinWorkDouble * saveX = new CoinWorkDouble[numberTotal];
219 CoinWorkDouble * saveY = new CoinWorkDouble[numberRows_];
220 CoinWorkDouble * saveZ = new CoinWorkDouble[numberTotal];
221 CoinWorkDouble * saveW = new CoinWorkDouble[numberTotal];
222 CoinWorkDouble * saveSL = new CoinWorkDouble[numberTotal];
223 CoinWorkDouble * saveSU = new CoinWorkDouble[numberTotal];
224 // Save smallest mu used in primal dual moves
225 CoinWorkDouble smallestPrimalDualMu = COIN_DBL_MAX;
226 CoinWorkDouble objScale = optimizationDirection_ /
227 (rhsScale_ * objectiveScale_);
228 while (problemStatus_ < 0) {
229 //#define FULL_DEBUG
230#ifdef FULL_DEBUG
231 {
232 int i;
233 printf("row pi artvec rhsfx\n");
234 for (i = 0; i < numberRows_; i++) {
235 printf("%d %g %g %g\n", i, dual_[i], errorRegion_[i], rhsFixRegion_[i]);
236 }
237 printf(" col dsol ddiag dwvec dzvec dbdslu dbdsll\n");
238 for (i = 0; i < numberColumns_ + numberRows_; i++) {
239 printf(" %d %g %g %g %g %g %g\n", i, solution_[i], diagonal_[i], wVec_[i],
240 zVec_[i], upperSlack_[i], lowerSlack_[i]);
241 }
242 }
243#endif
244 complementarityGap_ = complementarityGap(numberComplementarityPairs_,
245 numberComplementarityItems_, 0);
246 handler_->message(CLP_BARRIER_ITERATION, messages_)
247 << numberIterations_
248 << static_cast<double>(primalObjective_ * objScale - dblParam_[ClpObjOffset])
249 << static_cast<double>(dualObjective_ * objScale - dblParam_[ClpObjOffset])
250 << static_cast<double>(complementarityGap_)
251 << numberFixedTotal
252 << cholesky_->rank()
253 << CoinMessageEol;
254#if 0
255 if (numberIterations_ == -1) {
256 saveSolution("xxx.sav");
257 if (yyyyyy)
258 exit(99);
259 }
260#endif
261 // move up history
262 for (i = 1; i < LENGTH_HISTORY; i++)
263 historyInfeasibility_[i-1] = historyInfeasibility_[i];
264 historyInfeasibility_[LENGTH_HISTORY-1] = complementarityGap_;
265 // switch off saved if changes
266 //if (saveIteration+10<numberIterations_&&
267 //complementarityGap_*2.0<historyInfeasibility_[0])
268 //saveIteration=-1;
269 lastStep = CoinMin(actualPrimalStep_, actualDualStep_);
270 CoinWorkDouble goodGapChange;
271 //#define KEEP_GOING_IF_FIXED 5
272#ifndef KEEP_GOING_IF_FIXED
273#define KEEP_GOING_IF_FIXED 10000
274#endif
275 if (!sloppyOptimal) {
276 goodGapChange = 0.93;
277 } else {
278 goodGapChange = 0.7;
279 if (numberFixed > KEEP_GOING_IF_FIXED)
280 goodGapChange = 0.99; // make more likely to carry on
281 }
282 CoinWorkDouble gapO;
283 CoinWorkDouble lastGood = bestObjectiveGap;
284 if (gonePrimalFeasible_ && goneDualFeasible_) {
285 CoinWorkDouble largestObjective;
286 if (CoinAbs(primalObjective_) > CoinAbs(dualObjective_)) {
287 largestObjective = CoinAbs(primalObjective_);
288 } else {
289 largestObjective = CoinAbs(dualObjective_);
290 }
291 if (largestObjective < 1.0) {
292 largestObjective = 1.0;
293 }
294 gapO = CoinAbs(primalObjective_ - dualObjective_) / largestObjective;
295 handler_->message(CLP_BARRIER_OBJECTIVE_GAP, messages_)
296 << static_cast<double>(gapO)
297 << CoinMessageEol;
298 //start saving best
299 bool saveIt = false;
300 if (gapO < bestObjectiveGap) {
301 bestObjectiveGap = gapO;
302#ifndef SAVE_ON_OBJ
303 saveIt = true;
304#endif
305 }
306 if (primalObjective_ < bestObjective) {
307 bestObjective = primalObjective_;
308#ifdef SAVE_ON_OBJ
309 saveIt = true;
310#endif
311 }
312 if (numberFixedTotal > bestKilled) {
313 bestKilled = numberFixedTotal;
314#if KEEP_GOING_IF_FIXED<10
315 saveIt = true;
316#endif
317 }
318 if (saveIt) {
319#if KEEP_GOING_IF_FIXED<10
320 COIN_DETAIL_PRINT(printf("saving\n"));
321#endif
322 saveIteration = numberIterations_;
323 if (!savePi) {
324 savePi = new CoinWorkDouble[numberRows_];
325 savePrimal = new CoinWorkDouble [numberTotal];
326 }
327 CoinMemcpyN(dualArray, numberRows_, savePi);
328 CoinMemcpyN(solution_, numberTotal, savePrimal);
329 } else if(gapO > 2.0 * bestObjectiveGap) {
330 //maybe be more sophisticated e.g. re-initialize having
331 //fixed variables and dropped rows
332 //std::cout <<" gap increasing "<<std::endl;
333 }
334 //std::cout <<"could stop"<<std::endl;
335 //gapO=0.0;
336 if (CoinAbs(primalObjective_ - dualObjective_) < dualTolerance()) {
337 gapO = 0.0;
338 }
339 } else {
340 gapO = COIN_DBL_MAX;
341 if (saveIteration >= 0) {
342 handler_->message(CLP_BARRIER_GONE_INFEASIBLE, messages_)
343 << CoinMessageEol;
344 CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
345 // save alternate
346 if (numberFixedTotal > bestKilled &&
347 maximumBoundInfeasibility_ < 1.0e-6 &&
348 scaledRHSError < 1.0e-2) {
349 bestKilled = numberFixedTotal;
350#if KEEP_GOING_IF_FIXED<10
351 COIN_DETAIL_PRINT(printf("saving alternate\n"));
352#endif
353 saveIteration2 = numberIterations_;
354 if (!savePi2) {
355 savePi2 = new CoinWorkDouble[numberRows_];
356 savePrimal2 = new CoinWorkDouble [numberTotal];
357 }
358 CoinMemcpyN(dualArray, numberRows_, savePi2);
359 CoinMemcpyN(solution_, numberTotal, savePrimal2);
360 }
361 if (sloppyOptimal) {
362 // vaguely optimal
363 if (maximumBoundInfeasibility_ > 1.0e-2 ||
364 scaledRHSError > 1.0e-2 ||
365 maximumDualError_ > objectiveNorm_ * 1.0e-2) {
366 handler_->message(CLP_BARRIER_EXIT2, messages_)
367 << saveIteration
368 << CoinMessageEol;
369 problemStatus_ = 0; // benefit of doubt
370 break;
371 }
372 } else {
373 // not close to optimal but check if getting bad
374 CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
375 if ((maximumBoundInfeasibility_ > 1.0e-1 ||
376 scaledRHSError > 1.0e-1 ||
377 maximumDualError_ > objectiveNorm_ * 1.0e-1)
378 && (numberIterations_ > 50
379 && complementarityGap_ > 0.9 * historyInfeasibility_[0])) {
380 handler_->message(CLP_BARRIER_EXIT2, messages_)
381 << saveIteration
382 << CoinMessageEol;
383 break;
384 }
385 if (complementarityGap_ > 0.95 * checkGap && bestObjectiveGap < 1.0e-3 &&
386 (numberIterations_ > saveIteration + 5 || numberIterations_ > 100)) {
387 handler_->message(CLP_BARRIER_EXIT2, messages_)
388 << saveIteration
389 << CoinMessageEol;
390 break;
391 }
392 }
393 }
394 if (complementarityGap_ > 0.5 * checkGap && primalObjective_ >
395 bestObjective + 1.0e-9 &&
396 (numberIterations_ > saveIteration + 5 || numberIterations_ > 100)) {
397 handler_->message(CLP_BARRIER_EXIT2, messages_)
398 << saveIteration
399 << CoinMessageEol;
400 break;
401 }
402 }
403 if ((gapO < 1.0e-6 || (gapO < 1.0e-4 && complementarityGap_ < 0.1)) && !sloppyOptimal) {
404 sloppyOptimal = true;
405 handler_->message(CLP_BARRIER_CLOSE_TO_OPTIMAL, messages_)
406 << numberIterations_ << static_cast<double>(complementarityGap_)
407 << CoinMessageEol;
408 }
409 int numberBack = quadraticObj ? 10 : 5;
410 //tryJustPredictor=true;
411 //printf("trying just predictor\n");
412 //}
413 if (complementarityGap_ >= 1.05 * lastComplementarityGap) {
414 handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
415 << static_cast<double>(complementarityGap_) << "increasing"
416 << CoinMessageEol;
417 if (saveIteration >= 0 && sloppyOptimal) {
418 handler_->message(CLP_BARRIER_EXIT2, messages_)
419 << saveIteration
420 << CoinMessageEol;
421 break;
422 } else if (numberIterations_ - lastGoodIteration >= numberBack &&
423 complementarityGap_ < 1.0e-6) {
424 break; // not doing very well - give up
425 }
426 } else if (complementarityGap_ < goodGapChange * lastComplementarityGap) {
427 lastGoodIteration = numberIterations_;
428 lastComplementarityGap = complementarityGap_;
429 } else if (numberIterations_ - lastGoodIteration >= numberBack &&
430 complementarityGap_ < 1.0e-3) {
431 handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
432 << static_cast<double>(complementarityGap_) << "not decreasing"
433 << CoinMessageEol;
434 if (gapO > 0.75 * lastGood && numberFixed < KEEP_GOING_IF_FIXED) {
435 break;
436 }
437 } else if (numberIterations_ - lastGoodIteration >= 2 &&
438 complementarityGap_ < 1.0e-6) {
439 handler_->message(CLP_BARRIER_COMPLEMENTARITY, messages_)
440 << static_cast<double>(complementarityGap_) << "not decreasing"
441 << CoinMessageEol;
442 break;
443 }
444 if (numberIterations_ > maximumBarrierIterations_ || hitMaximumIterations()) {
445 handler_->message(CLP_BARRIER_STOPPING, messages_)
446 << CoinMessageEol;
447 problemStatus_ = 3;
448 onStopped(); // set secondary status
449 break;
450 }
451 if (gapO < targetGap_) {
452 problemStatus_ = 0;
453 handler_->message(CLP_BARRIER_EXIT, messages_)
454 << " "
455 << CoinMessageEol;
456 break;//finished
457 }
458 if (complementarityGap_ < 1.0e-12) {
459 problemStatus_ = 0;
460 handler_->message(CLP_BARRIER_EXIT, messages_)
461 << "- small complementarity gap"
462 << CoinMessageEol;
463 break;//finished
464 }
465 if (complementarityGap_ < 1.0e-10 && gapO < 1.0e-10) {
466 problemStatus_ = 0;
467 handler_->message(CLP_BARRIER_EXIT, messages_)
468 << "- objective gap and complementarity gap both small"
469 << CoinMessageEol;
470 break;//finished
471 }
472 if (gapO < 1.0e-9) {
473 CoinWorkDouble value = gapO * complementarityGap_;
474 value *= actualPrimalStep_;
475 value *= actualDualStep_;
476 //std::cout<<value<<std::endl;
477 if (value < 1.0e-17 && numberIterations_ > lastGoodIteration) {
478 problemStatus_ = 0;
479 handler_->message(CLP_BARRIER_EXIT, messages_)
480 << "- objective gap and complementarity gap both smallish and small steps"
481 << CoinMessageEol;
482 break;//finished
483 }
484 }
485 CoinWorkDouble nextGap = COIN_DBL_MAX;
486 int nextNumber = 0;
487 int nextNumberItems = 0;
488 worstDirectionAccuracy_ = 0.0;
489 int newDropped = 0;
490 //Predictor step
491 //prepare for cholesky. Set up scaled diagonal in deltaX
492 // ** for efficiency may be better if scale factor known before
493 CoinWorkDouble norm2 = 0.0;
494 CoinWorkDouble maximumValue;
495 getNorms(diagonal_, numberTotal, maximumValue, norm2);
496 diagonalNorm_ = CoinSqrt(norm2 / numberComplementarityPairs_);
497 diagonalScaleFactor_ = 1.0;
498 CoinWorkDouble maximumAllowable = eScale;
499 //scale so largest is less than allowable ? could do better
500 CoinWorkDouble factor = 0.5;
501 while (maximumValue > maximumAllowable) {
502 diagonalScaleFactor_ *= factor;
503 maximumValue *= factor;
504 } /* endwhile */
505 if (diagonalScaleFactor_ != 1.0) {
506 handler_->message(CLP_BARRIER_SCALING, messages_)
507 << "diagonal" << static_cast<double>(diagonalScaleFactor_)
508 << CoinMessageEol;
509 diagonalNorm_ *= diagonalScaleFactor_;
510 }
511 multiplyAdd(NULL, numberTotal, 0.0, diagonal_,
512 diagonalScaleFactor_);
513 int * rowsDroppedThisTime = new int [numberRows_];
514 newDropped = cholesky_->factorize(diagonal_, rowsDroppedThisTime);
515 if (newDropped) {
516 if (newDropped == -1) {
517 COIN_DETAIL_PRINT(printf("Out of memory\n"));
518 problemStatus_ = 4;
519 //delete all temporary regions
520 deleteWorkingData();
521 if (saveMatrix) {
522 // restore normal copy
523 delete matrix_;
524 matrix_ = saveMatrix;
525 }
526 return -1;
527 } else {
528#ifndef NDEBUG
529 //int newDropped2=cholesky_->factorize(diagonal_,rowsDroppedThisTime);
530 //assert(!newDropped2);
531#endif
532 if (newDropped < 0 && 0) {
533 //replace dropped
534 newDropped = -newDropped;
535 //off 1 to allow for reset all
536 newDropped--;
537 //set all bits false
538 cholesky_->resetRowsDropped();
539 }
540 }
541 }
542 delete [] rowsDroppedThisTime;
543 if (cholesky_->status()) {
544 std::cout << "bad cholesky?" << std::endl;
545 abort();
546 }
547 int phase = 0; // predictor, corrector , primal dual
548 CoinWorkDouble directionAccuracy = 0.0;
549 bool doCorrector = true;
550 bool goodMove = true;
551 //set up for affine direction
552 setupForSolve(phase);
553 if ((modeSwitch & 2) == 0) {
554 directionAccuracy = findDirectionVector(phase);
555 if (directionAccuracy > worstDirectionAccuracy_) {
556 worstDirectionAccuracy_ = directionAccuracy;
557 }
558 if (saveIteration > 0 && directionAccuracy > 1.0) {
559 handler_->message(CLP_BARRIER_EXIT2, messages_)
560 << saveIteration
561 << CoinMessageEol;
562 break;
563 }
564 findStepLength(phase);
565 nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
566 debugMove(0, actualPrimalStep_, actualDualStep_);
567 debugMove(0, 1.0e-2, 1.0e-2);
568 }
569 CoinWorkDouble affineGap = nextGap;
570 int bestPhase = 0;
571 CoinWorkDouble bestNextGap = nextGap;
572 // ?
573 bestNextGap = CoinMax(nextGap, 0.8 * complementarityGap_);
574 if (quadraticObj)
575 bestNextGap = CoinMax(nextGap, 0.99 * complementarityGap_);
576 if (complementarityGap_ > 1.0e-4 * numberComplementarityPairs_) {
577 //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
578 CoinWorkDouble part1 = nextGap / numberComplementarityPairs_;
579 part1 = nextGap / numberComplementarityItems_;
580 CoinWorkDouble part2 = nextGap / complementarityGap_;
581 mu_ = part1 * part2 * part2;
582#if 0
583 CoinWorkDouble papermu = complementarityGap_ / numberComplementarityPairs_;
584 CoinWorkDouble affmu = nextGap / nextNumber;
585 CoinWorkDouble sigma = pow(affmu / papermu, 3);
586 printf("mu %g, papermu %g, affmu %g, sigma %g sigmamu %g\n",
587 mu_, papermu, affmu, sigma, sigma * papermu);
588#endif
589 //printf("paper mu %g\n",(nextGap*nextGap*nextGap)/(complementarityGap_*complementarityGap_*
590 // (CoinWorkDouble) numberComplementarityPairs_));
591 } else {
592 CoinWorkDouble phi;
593 if (numberComplementarityPairs_ <= 5000) {
594 phi = pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_), 2.0);
595 } else {
596 phi = pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_), 1.5);
597 if (phi < 500.0 * 500.0) {
598 phi = 500.0 * 500.0;
599 }
600 }
601 mu_ = complementarityGap_ / phi;
602 }
603 //save information
604 CoinWorkDouble product = affineProduct();
605#if 0
606 //can we do corrector step?
607 CoinWorkDouble xx = complementarityGap_ * (beta2 - tau) + product;
608 if (xx > 0.0) {
609 CoinWorkDouble saveMu = mu_;
610 CoinWorkDouble mu2 = numberComplementarityPairs_;
611 mu2 = xx / mu2;
612 if (mu2 > mu_) {
613 //std::cout<<" could increase to "<<mu2<<std::endl;
614 //was mu2=mu2*0.25;
615 mu2 = mu2 * 0.99;
616 if (mu2 < mu_) {
617 mu_ = mu2;
618 //std::cout<<" changing mu to "<<mu_<<std::endl;
619 } else {
620 //std::cout<<std::endl;
621 }
622 } else {
623 //std::cout<<" should decrease to "<<mu2<<std::endl;
624 mu_ = 0.5 * mu2;
625 //std::cout<<" changing mu to "<<mu_<<std::endl;
626 }
627 handler_->message(CLP_BARRIER_MU, messages_)
628 << saveMu << mu_
629 << CoinMessageEol;
630 } else {
631 //std::cout<<" bad by any standards"<<std::endl;
632 }
633#endif
634 if (complementarityGap_*(beta2 - tau) + product - mu_ * numberComplementarityPairs_ < 0.0 && 0) {
635#ifdef SOME_DEBUG
636 printf("failed 1 product %.18g mu %.18g - %.18g < 0.0, nextGap %.18g\n", product, mu_,
637 complementarityGap_*(beta2 - tau) + product - mu_ * numberComplementarityPairs_,
638 nextGap);
639#endif
640 doCorrector = false;
641 if (nextGap > 0.9 * complementarityGap_ || 1) {
642 goodMove = false;
643 bestNextGap = COIN_DBL_MAX;
644 }
645 //CoinWorkDouble floatNumber = 2.0*numberComplementarityPairs_;
646 //floatNumber = 1.0*numberComplementarityItems_;
647 //mu_=nextGap/floatNumber;
648 handler_->message(CLP_BARRIER_INFO, messages_)
649 << "no corrector step"
650 << CoinMessageEol;
651 } else {
652 phase = 1;
653 }
654 // If bad gap - try standard primal dual
655 if (nextGap > complementarityGap_ * 1.001)
656 goodMove = false;
657 if ((modeSwitch & 2) != 0)
658 goodMove = false;
659 if (goodMove && doCorrector) {
660 CoinMemcpyN(deltaX_, numberTotal, saveX);
661 CoinMemcpyN(deltaY_, numberRows_, saveY);
662 CoinMemcpyN(deltaZ_, numberTotal, saveZ);
663 CoinMemcpyN(deltaW_, numberTotal, saveW);
664 CoinMemcpyN(deltaSL_, numberTotal, saveSL);
665 CoinMemcpyN(deltaSU_, numberTotal, saveSU);
666#ifdef HALVE
667 CoinWorkDouble savePrimalStep = actualPrimalStep_;
668 CoinWorkDouble saveDualStep = actualDualStep_;
669 CoinWorkDouble saveMu = mu_;
670#endif
671 //set up for next step
672 setupForSolve(phase);
673 CoinWorkDouble directionAccuracy2 = findDirectionVector(phase);
674 if (directionAccuracy2 > worstDirectionAccuracy_) {
675 worstDirectionAccuracy_ = directionAccuracy2;
676 }
677 CoinWorkDouble testValue = 1.0e2 * directionAccuracy;
678 if (1.0e2 * projectionTolerance_ > testValue) {
679 testValue = 1.0e2 * projectionTolerance_;
680 }
681 if (primalTolerance() > testValue) {
682 testValue = primalTolerance();
683 }
684 if (maximumRHSError_ > testValue) {
685 testValue = maximumRHSError_;
686 }
687 if (directionAccuracy2 > testValue && numberIterations_ >= -77) {
688 goodMove = false;
689#ifdef SOME_DEBUG
690 printf("accuracy %g phase 1 failed, test value %g\n",
691 directionAccuracy2, testValue);
692#endif
693 }
694 if (goodMove) {
695 phase = 1;
696 CoinWorkDouble norm = findStepLength(phase);
697 nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
698 debugMove(1, actualPrimalStep_, actualDualStep_);
699 //debugMove(1,1.0e-7,1.0e-7);
700 goodMove = checkGoodMove(true, bestNextGap, allowIncreasingGap);
701 if (norm < 0)
702 goodMove = false;
703 if (!goodMove) {
704#ifdef SOME_DEBUG
705 printf("checkGoodMove failed\n");
706#endif
707 }
708 }
709#ifdef HALVE
710 int nHalve = 0;
711 // relax test
712 bestNextGap = CoinMax(bestNextGap, 0.9 * complementarityGap_);
713 while (!goodMove) {
714 mu_ = saveMu;
715 actualPrimalStep_ = savePrimalStep;
716 actualDualStep_ = saveDualStep;
717 int i;
718 //printf("halve %d\n",nHalve);
719 nHalve++;
720 const CoinWorkDouble lambda = 0.5;
721 for (i = 0; i < numberRows_; i++)
722 deltaY_[i] = lambda * deltaY_[i] + (1.0 - lambda) * saveY[i];
723 for (i = 0; i < numberTotal; i++) {
724 deltaX_[i] = lambda * deltaX_[i] + (1.0 - lambda) * saveX[i];
725 deltaZ_[i] = lambda * deltaZ_[i] + (1.0 - lambda) * saveZ[i];
726 deltaW_[i] = lambda * deltaW_[i] + (1.0 - lambda) * saveW[i];
727 deltaSL_[i] = lambda * deltaSL_[i] + (1.0 - lambda) * saveSL[i];
728 deltaSU_[i] = lambda * deltaSU_[i] + (1.0 - lambda) * saveSU[i];
729 }
730 CoinMemcpyN(saveX, numberTotal, deltaX_);
731 CoinMemcpyN(saveY, numberRows_, deltaY_);
732 CoinMemcpyN(saveZ, numberTotal, deltaZ_);
733 CoinMemcpyN(saveW, numberTotal, deltaW_);
734 CoinMemcpyN(saveSL, numberTotal, deltaSL_);
735 CoinMemcpyN(saveSU, numberTotal, deltaSU_);
736 findStepLength(1);
737 nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
738 goodMove = checkGoodMove(true, bestNextGap, allowIncreasingGap);
739 if (nHalve > 10)
740 break;
741 //assert (goodMove);
742 }
743 if (nHalve && handler_->logLevel() > 2)
744 printf("halved %d times\n", nHalve);
745#endif
746 }
747 //bestPhase=-1;
748 //goodMove=false;
749 if (!goodMove) {
750 // Just primal dual step
751 CoinWorkDouble floatNumber;
752 floatNumber = 2.0 * numberComplementarityPairs_;
753 //floatNumber = numberComplementarityItems_;
754 CoinWorkDouble saveMu = mu_; // use one from predictor corrector
755 mu_ = complementarityGap_ / floatNumber;
756 // If going well try small mu
757 mu_ *= CoinSqrt((1.0 - lastStep) / (1.0 + 10.0 * lastStep));
758 CoinWorkDouble mu1 = mu_;
759 CoinWorkDouble phi;
760 if (numberComplementarityPairs_ <= 500) {
761 phi = pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_), 2.0);
762 } else {
763 phi = pow(static_cast<CoinWorkDouble> (numberComplementarityPairs_), 1.5);
764 if (phi < 500.0 * 500.0) {
765 phi = 500.0 * 500.0;
766 }
767 }
768 mu_ = complementarityGap_ / phi;
769 //printf("pd mu %g, alternate %g, smallest %g\n",
770 // mu_,mu1,smallestPrimalDualMu);
771 mu_ = CoinSqrt(mu_ * mu1);
772 mu_ = mu1;
773 if ((numberIterations_ & 1) == 0 || numberIterations_ < 10)
774 mu_ = saveMu;
775 //mu_=CoinMin(smallestPrimalDualMu*0.95,mu_);
776 smallestPrimalDualMu = mu_;
777 // Try simpler
778 floatNumber = numberComplementarityItems_;
779 mu_ = 0.5 * complementarityGap_ / floatNumber;
780 //if ((modeSwitch&2)==0) {
781 //if ((numberIterations_&1)==0)
782 // mu_ *= 0.5;
783 //} else {
784 //mu_ *= 0.8;
785 //}
786 //set up for next step
787 setupForSolve(2);
788 findDirectionVector(2);
789 CoinWorkDouble norm = findStepLength(2);
790 // just for debug
791 bestNextGap = complementarityGap_ * 1.0005;
792 //bestNextGap=COIN_DBL_MAX;
793 nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
794 debugMove(2, actualPrimalStep_, actualDualStep_);
795 //debugMove(2,1.0e-7,1.0e-7);
796 checkGoodMove(false, bestNextGap, allowIncreasingGap);
797 if ((nextGap > 0.9 * complementarityGap_ && bestPhase == 0 && affineGap < nextGap
798 && (numberIterations_ > 80 || (numberIterations_ > 20 && quadraticObj))) || norm < 0.0) {
799 // Back to affine
800 phase = 0;
801 setupForSolve(phase);
802 directionAccuracy = findDirectionVector(phase);
803 findStepLength(phase);
804 nextGap = complementarityGap(nextNumber, nextNumberItems, 1);
805 bestNextGap = complementarityGap_;
806 //checkGoodMove(false, bestNextGap,allowIncreasingGap);
807 }
808 }
809 if (numberIterations_ == 0)
810 smallestPrimalDualMu = mu_;
811 if (!goodMove)
812 mu_ = nextGap / (static_cast<CoinWorkDouble> (nextNumber) * 1.1);
813 //if (quadraticObj)
814 //goodMove=true;
815 //goodMove=false; //TEMP
816 // Do centering steps
817 int numberTries = 0;
818 CoinWorkDouble nextCenterGap = 0.0;
819 int numberGoodTries = 0;
820#ifdef COIN_DETAIL
821 CoinWorkDouble originalDualStep = actualDualStep_;
822 CoinWorkDouble originalPrimalStep = actualPrimalStep_;
823#endif
824 if (actualDualStep_ > 0.9 && actualPrimalStep_ > 0.9)
825 goodMove = false; // don't bother
826 if ((modeSwitch & 1) != 0)
827 goodMove = false;
828 while (goodMove && numberTries < 5) {
829 goodMove = false;
830 numberTries++;
831 CoinMemcpyN(deltaX_, numberTotal, saveX);
832 CoinMemcpyN(deltaY_, numberRows_, saveY);
833 CoinMemcpyN(deltaZ_, numberTotal, saveZ);
834 CoinMemcpyN(deltaW_, numberTotal, saveW);
835 CoinWorkDouble savePrimalStep = actualPrimalStep_;
836 CoinWorkDouble saveDualStep = actualDualStep_;
837 CoinWorkDouble saveMu = mu_;
838 setupForSolve(3);
839 findDirectionVector(3);
840 findStepLength(3);
841 debugMove(3, actualPrimalStep_, actualDualStep_);
842 //debugMove(3,1.0e-7,1.0e-7);
843 CoinWorkDouble xGap = complementarityGap(nextNumber, nextNumberItems, 3);
844 // If one small then that's the one that counts
845 CoinWorkDouble checkDual = saveDualStep;
846 CoinWorkDouble checkPrimal = savePrimalStep;
847 if (checkDual > 5.0 * checkPrimal) {
848 checkDual = 2.0 * checkPrimal;
849 } else if (checkPrimal > 5.0 * checkDual) {
850 checkPrimal = 2.0 * checkDual;
851 }
852 if (actualPrimalStep_ < checkPrimal ||
853 actualDualStep_ < checkDual ||
854 (xGap > nextGap && xGap > 0.9 * complementarityGap_)) {
855 //if (actualPrimalStep_<=checkPrimal||
856 //actualDualStep_<=checkDual) {
857#ifdef SOME_DEBUG
858 printf("PP rejected gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n", xGap,
859 actualPrimalStep_, actualDualStep_, nextGap, savePrimalStep, saveDualStep);
860#endif
861 mu_ = saveMu;
862 actualPrimalStep_ = savePrimalStep;
863 actualDualStep_ = saveDualStep;
864 CoinMemcpyN(saveX, numberTotal, deltaX_);
865 CoinMemcpyN(saveY, numberRows_, deltaY_);
866 CoinMemcpyN(saveZ, numberTotal, deltaZ_);
867 CoinMemcpyN(saveW, numberTotal, deltaW_);
868 } else {
869#ifdef SOME_DEBUG
870 printf("PPphase 3 gap %.18g, steps %.18g %.18g, 2 gap %.18g, steps %.18g %.18g\n", xGap,
871 actualPrimalStep_, actualDualStep_, nextGap, savePrimalStep, saveDualStep);
872#endif
873 numberGoodTries++;
874 nextCenterGap = xGap;
875 // See if big enough change
876 if (actualPrimalStep_ < 1.01 * checkPrimal ||
877 actualDualStep_ < 1.01 * checkDual) {
878 // stop now
879 } else {
880 // carry on
881 goodMove = true;
882 }
883 }
884 }
885 if (numberGoodTries && handler_->logLevel() > 1) {
886 COIN_DETAIL_PRINT(printf("%d centering steps moved from (gap %.18g, dual %.18g, primal %.18g) to (gap %.18g, dual %.18g, primal %.18g)\n",
887 numberGoodTries, static_cast<double>(nextGap), static_cast<double>(originalDualStep),
888 static_cast<double>(originalPrimalStep),
889 static_cast<double>(nextCenterGap), static_cast<double>(actualDualStep_),
890 static_cast<double>(actualPrimalStep_)));
891 }
892 // save last gap
893 checkGap = complementarityGap_;
894 numberFixed = updateSolution(nextGap);
895 numberFixedTotal += numberFixed;
896 } /* endwhile */
897 delete [] saveX;
898 delete [] saveY;
899 delete [] saveZ;
900 delete [] saveW;
901 delete [] saveSL;
902 delete [] saveSU;
903 if (savePi) {
904 if (numberIterations_ - saveIteration > 20 &&
905 numberIterations_ - saveIteration2 < 5) {
906#if KEEP_GOING_IF_FIXED<10
907 std::cout << "Restoring2 from iteration " << saveIteration2 << std::endl;
908#endif
909 CoinMemcpyN(savePi2, numberRows_, dualArray);
910 CoinMemcpyN(savePrimal2, numberTotal, solution_);
911 } else {
912#if KEEP_GOING_IF_FIXED<10
913 std::cout << "Restoring from iteration " << saveIteration << std::endl;
914#endif
915 CoinMemcpyN(savePi, numberRows_, dualArray);
916 CoinMemcpyN(savePrimal, numberTotal, solution_);
917 }
918 delete [] savePi;
919 delete [] savePrimal;
920 }
921 delete [] savePi2;
922 delete [] savePrimal2;
923 //recompute slacks
924 // Split out solution
925 CoinZeroN(rowActivity_, numberRows_);
926 CoinMemcpyN(solution_, numberColumns_, columnActivity_);
927 matrix_->times(1.0, columnActivity_, rowActivity_);
928 //unscale objective
929 multiplyAdd(NULL, numberTotal, 0.0, cost_, scaleFactor_);
930 multiplyAdd(NULL, numberRows_, 0, dualArray, scaleFactor_);
931 checkSolution();
932 //CoinMemcpyN(reducedCost_,numberColumns_,dj_);
933 // If quadratic use last solution
934 // Restore quadratic objective if necessary
935 if (saveObjective) {
936 delete objective_;
937 objective_ = saveObjective;
938 objectiveValue_ = 0.5 * (primalObjective_ + dualObjective_);
939 }
940 handler_->message(CLP_BARRIER_END, messages_)
941 << static_cast<double>(sumPrimalInfeasibilities_)
942 << static_cast<double>(sumDualInfeasibilities_)
943 << static_cast<double>(complementarityGap_)
944 << static_cast<double>(objectiveValue())
945 << CoinMessageEol;
946 //#ifdef SOME_DEBUG
947 if (handler_->logLevel() > 1)
948 COIN_DETAIL_PRINT(printf("ENDRUN status %d after %d iterations\n", problemStatus_, numberIterations_));
949 //#endif
950 //std::cout<<"Absolute primal infeasibility at end "<<sumPrimalInfeasibilities_<<std::endl;
951 //std::cout<<"Absolute dual infeasibility at end "<<sumDualInfeasibilities_<<std::endl;
952 //std::cout<<"Absolute complementarity at end "<<complementarityGap_<<std::endl;
953 //std::cout<<"Primal objective "<<objectiveValue()<<std::endl;
954 //std::cout<<"maximum complementarity "<<worstComplementarity_<<std::endl;
955#if COIN_LONG_WORK
956 // put back dual
957 delete [] dual_;
958 delete [] reducedCost_;
959 dual_ = dualSave;
960 reducedCost_ = reducedCostSave;
961#endif
962 //delete all temporary regions
963 deleteWorkingData();
964#if KEEP_GOING_IF_FIXED<10
965#if 0 //ndef NDEBUG
966 {
967 static int kk = 0;
968 char name[20];
969 sprintf(name, "save.sol.%d", kk);
970 kk++;
971 printf("saving to file %s\n", name);
972 FILE * fp = fopen(name, "wb");
973 int numberWritten;
974 numberWritten = fwrite(&numberColumns_, sizeof(int), 1, fp);
975 assert (numberWritten == 1);
976 numberWritten = fwrite(columnActivity_, sizeof(double), numberColumns_, fp);
977 assert (numberWritten == numberColumns_);
978 fclose(fp);
979 }
980#endif
981#endif
982 if (saveMatrix) {
983 // restore normal copy
984 delete matrix_;
985 matrix_ = saveMatrix;
986 }
987 return problemStatus_;
988}
989// findStepLength.
990//phase - 0 predictor
991// 1 corrector
992// 2 primal dual
993CoinWorkDouble ClpPredictorCorrector::findStepLength( int phase)
994{
995 CoinWorkDouble directionNorm = 0.0;
996 CoinWorkDouble maximumPrimalStep = COIN_DBL_MAX * 1.0e-20;
997 CoinWorkDouble maximumDualStep = COIN_DBL_MAX;
998 int numberTotal = numberRows_ + numberColumns_;
999 CoinWorkDouble tolerance = 1.0e-12;
1000 int chosenPrimalSequence = -1;
1001 int chosenDualSequence = -1;
1002 bool lowPrimal = false;
1003 bool lowDual = false;
1004 // If done many iterations then allow to hit boundary
1005 CoinWorkDouble hitTolerance;
1006 //printf("objective norm %g\n",objectiveNorm_);
1007 if (numberIterations_ < 80 || !gonePrimalFeasible_)
1008 hitTolerance = COIN_DBL_MAX;
1009 else
1010 hitTolerance = CoinMax(1.0e3, 1.0e-3 * objectiveNorm_);
1011 int iColumn;
1012 //printf("dual value %g\n",dual_[0]);
1013 //printf(" X dX lS dlS uS dUs dj Z dZ t dT\n");
1014 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1015 if (!flagged(iColumn)) {
1016 CoinWorkDouble directionElement = deltaX_[iColumn];
1017 if (directionNorm < CoinAbs(directionElement)) {
1018 directionNorm = CoinAbs(directionElement);
1019 }
1020 if (lowerBound(iColumn)) {
1021 CoinWorkDouble delta = - deltaSL_[iColumn];
1022 CoinWorkDouble z1 = deltaZ_[iColumn];
1023 CoinWorkDouble newZ = zVec_[iColumn] + z1;
1024 if (zVec_[iColumn] > tolerance) {
1025 if (zVec_[iColumn] < -z1 * maximumDualStep) {
1026 maximumDualStep = -zVec_[iColumn] / z1;
1027 chosenDualSequence = iColumn;
1028 lowDual = true;
1029 }
1030 }
1031 if (lowerSlack_[iColumn] < maximumPrimalStep * delta) {
1032 CoinWorkDouble newStep = lowerSlack_[iColumn] / delta;
1033 if (newStep > 0.2 || newZ < hitTolerance || delta > 1.0e3 || delta <= 1.0e-6 || dj_[iColumn] < hitTolerance) {
1034 maximumPrimalStep = newStep;
1035 chosenPrimalSequence = iColumn;
1036 lowPrimal = true;
1037 } else {
1038 //printf("small %d delta %g newZ %g step %g\n",iColumn,delta,newZ,newStep);
1039 }
1040 }
1041 }
1042 if (upperBound(iColumn)) {
1043 CoinWorkDouble delta = - deltaSU_[iColumn];
1044 CoinWorkDouble w1 = deltaW_[iColumn];
1045 CoinWorkDouble newT = wVec_[iColumn] + w1;
1046 if (wVec_[iColumn] > tolerance) {
1047 if (wVec_[iColumn] < -w1 * maximumDualStep) {
1048 maximumDualStep = -wVec_[iColumn] / w1;
1049 chosenDualSequence = iColumn;
1050 lowDual = false;
1051 }
1052 }
1053 if (upperSlack_[iColumn] < maximumPrimalStep * delta) {
1054 CoinWorkDouble newStep = upperSlack_[iColumn] / delta;
1055 if (newStep > 0.2 || newT < hitTolerance || delta > 1.0e3 || delta <= 1.0e-6 || dj_[iColumn] > -hitTolerance) {
1056 maximumPrimalStep = newStep;
1057 chosenPrimalSequence = iColumn;
1058 lowPrimal = false;
1059 } else {
1060 //printf("small %d delta %g newT %g step %g\n",iColumn,delta,newT,newStep);
1061 }
1062 }
1063 }
1064 }
1065 }
1066#ifdef SOME_DEBUG
1067 printf("new step - phase %d, norm %.18g, dual step %.18g, primal step %.18g\n",
1068 phase, directionNorm, maximumDualStep, maximumPrimalStep);
1069 if (lowDual)
1070 printf("ld %d %g %g => %g (dj %g,sol %g) ",
1071 chosenDualSequence, zVec_[chosenDualSequence],
1072 deltaZ_[chosenDualSequence], zVec_[chosenDualSequence] +
1073 maximumDualStep * deltaZ_[chosenDualSequence], dj_[chosenDualSequence],
1074 solution_[chosenDualSequence]);
1075 else
1076 printf("ud %d %g %g => %g (dj %g,sol %g) ",
1077 chosenDualSequence, wVec_[chosenDualSequence],
1078 deltaW_[chosenDualSequence], wVec_[chosenDualSequence] +
1079 maximumDualStep * deltaW_[chosenDualSequence], dj_[chosenDualSequence],
1080 solution_[chosenDualSequence]);
1081 if (lowPrimal)
1082 printf("lp %d %g %g => %g (dj %g,sol %g)\n",
1083 chosenPrimalSequence, lowerSlack_[chosenPrimalSequence],
1084 deltaSL_[chosenPrimalSequence], lowerSlack_[chosenPrimalSequence] +
1085 maximumPrimalStep * deltaSL_[chosenPrimalSequence],
1086 dj_[chosenPrimalSequence], solution_[chosenPrimalSequence]);
1087 else
1088 printf("up %d %g %g => %g (dj %g,sol %g)\n",
1089 chosenPrimalSequence, upperSlack_[chosenPrimalSequence],
1090 deltaSU_[chosenPrimalSequence], upperSlack_[chosenPrimalSequence] +
1091 maximumPrimalStep * deltaSU_[chosenPrimalSequence],
1092 dj_[chosenPrimalSequence], solution_[chosenPrimalSequence]);
1093#endif
1094 actualPrimalStep_ = stepLength_ * maximumPrimalStep;
1095 if (phase >= 0 && actualPrimalStep_ > 1.0) {
1096 actualPrimalStep_ = 1.0;
1097 }
1098 actualDualStep_ = stepLength_ * maximumDualStep;
1099 if (phase >= 0 && actualDualStep_ > 1.0) {
1100 actualDualStep_ = 1.0;
1101 }
1102 // See if quadratic objective
1103#ifndef NO_RTTI
1104 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
1105#else
1106 ClpQuadraticObjective * quadraticObj = NULL;
1107 if (objective_->type() == 2)
1108 quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
1109#endif
1110 if (quadraticObj) {
1111 // Use smaller unless very small
1112 CoinWorkDouble smallerStep = CoinMin(actualDualStep_, actualPrimalStep_);
1113 if (smallerStep > 0.0001) {
1114 actualDualStep_ = smallerStep;
1115 actualPrimalStep_ = smallerStep;
1116 }
1117 }
1118#define OFFQ
1119#ifndef OFFQ
1120 if (quadraticObj) {
1121 // Don't bother if phase 0 or 3 or large gap
1122 //if ((phase==1||phase==2||phase==0)&&maximumDualError_>0.1*complementarityGap_
1123 //&&smallerStep>0.001) {
1124 if ((phase == 1 || phase == 2 || phase == 0 || phase == 3)) {
1125 // minimize complementarity + norm*dual inf ? primal inf
1126 // at first - just check better - if not
1127 // Complementarity gap will be a*change*change + b*change +c
1128 CoinWorkDouble a = 0.0;
1129 CoinWorkDouble b = 0.0;
1130 CoinWorkDouble c = 0.0;
1131 /* SQUARE of dual infeasibility will be:
1132 square of dj - ......
1133 */
1134 CoinWorkDouble aq = 0.0;
1135 CoinWorkDouble bq = 0.0;
1136 CoinWorkDouble cq = 0.0;
1137 CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
1138 CoinWorkDouble * linearDjChange = new CoinWorkDouble[numberTotal];
1139 CoinZeroN(linearDjChange, numberColumns_);
1140 multiplyAdd(deltaY_, numberRows_, 1.0, linearDjChange + numberColumns_, 0.0);
1141 matrix_->transposeTimes(-1.0, deltaY_, linearDjChange);
1142 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1143 const int * columnQuadratic = quadratic->getIndices();
1144 const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1145 const int * columnQuadraticLength = quadratic->getVectorLengths();
1146 CoinWorkDouble * quadraticElement = quadratic->getMutableElements();
1147 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1148 CoinWorkDouble oldPrimal = solution_[iColumn];
1149 if (!flagged(iColumn)) {
1150 if (lowerBound(iColumn)) {
1151 CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
1152 c += lowerSlack_[iColumn] * zVec_[iColumn];
1153 b += lowerSlack_[iColumn] * deltaZ_[iColumn] + zVec_[iColumn] * change;
1154 a += deltaZ_[iColumn] * change;
1155 }
1156 if (upperBound(iColumn)) {
1157 CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
1158 c += upperSlack_[iColumn] * wVec_[iColumn];
1159 b += upperSlack_[iColumn] * deltaW_[iColumn] + wVec_[iColumn] * change;
1160 a += deltaW_[iColumn] * change;
1161 }
1162 // new djs are dj_ + change*value
1163 CoinWorkDouble djChange = linearDjChange[iColumn];
1164 if (iColumn < numberColumns_) {
1165 for (CoinBigIndex j = columnQuadraticStart[iColumn];
1166 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1167 int jColumn = columnQuadratic[j];
1168 CoinWorkDouble changeJ = deltaX_[jColumn];
1169 CoinWorkDouble elementValue = quadraticElement[j];
1170 djChange += changeJ * elementValue;
1171 }
1172 }
1173 CoinWorkDouble gammaTerm = gamma2;
1174 if (primalR_) {
1175 gammaTerm += primalR_[iColumn];
1176 }
1177 djChange += gammaTerm;
1178 // and dual infeasibility
1179 CoinWorkDouble oldInf = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn] +
1180 gammaTerm * solution_[iColumn];
1181 CoinWorkDouble changeInf = djChange - deltaZ_[iColumn] + deltaW_[iColumn];
1182 cq += oldInf * oldInf;
1183 bq += 2.0 * oldInf * changeInf;
1184 aq += changeInf * changeInf;
1185 } else {
1186 // fixed
1187 if (lowerBound(iColumn)) {
1188 c += lowerSlack_[iColumn] * zVec_[iColumn];
1189 }
1190 if (upperBound(iColumn)) {
1191 c += upperSlack_[iColumn] * wVec_[iColumn];
1192 }
1193 // new djs are dj_ + change*value
1194 CoinWorkDouble djChange = linearDjChange[iColumn];
1195 if (iColumn < numberColumns_) {
1196 for (CoinBigIndex j = columnQuadraticStart[iColumn];
1197 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1198 int jColumn = columnQuadratic[j];
1199 CoinWorkDouble changeJ = deltaX_[jColumn];
1200 CoinWorkDouble elementValue = quadraticElement[j];
1201 djChange += changeJ * elementValue;
1202 }
1203 }
1204 CoinWorkDouble gammaTerm = gamma2;
1205 if (primalR_) {
1206 gammaTerm += primalR_[iColumn];
1207 }
1208 djChange += gammaTerm;
1209 // and dual infeasibility
1210 CoinWorkDouble oldInf = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn] +
1211 gammaTerm * solution_[iColumn];
1212 CoinWorkDouble changeInf = djChange - deltaZ_[iColumn] + deltaW_[iColumn];
1213 cq += oldInf * oldInf;
1214 bq += 2.0 * oldInf * changeInf;
1215 aq += changeInf * changeInf;
1216 }
1217 }
1218 delete [] linearDjChange;
1219 // ? We want to minimize complementarityGap + solutionNorm_*square of inf ??
1220 // maybe use inf and do line search
1221 // To check see if matches at current step
1222 CoinWorkDouble step = actualPrimalStep_;
1223 //Current gap + solutionNorm_ * CoinSqrt (sum square inf)
1224 CoinWorkDouble multiplier = solutionNorm_;
1225 multiplier *= 0.01;
1226 multiplier = 1.0;
1227 CoinWorkDouble currentInf = multiplier * CoinSqrt(cq);
1228 CoinWorkDouble nextInf = multiplier * CoinSqrt(CoinMax(cq + step * bq + step * step * aq, 0.0));
1229 CoinWorkDouble allowedIncrease = 1.4;
1230#ifdef SOME_DEBUG
1231 printf("lin %g %g %g -> %g\n", a, b, c,
1232 c + b * step + a * step * step);
1233 printf("quad %g %g %g -> %g\n", aq, bq, cq,
1234 cq + bq * step + aq * step * step);
1235 debugMove(7, step, step);
1236 printf ("current dualInf %g, with step of %g is %g\n",
1237 currentInf, step, nextInf);
1238#endif
1239 if (b > -1.0e-6) {
1240 if (phase != 0)
1241 directionNorm = -1.0;
1242 }
1243 if ((phase == 1 || phase == 2 || phase == 0 || phase == 3) && nextInf > 0.1 * complementarityGap_ &&
1244 nextInf > currentInf * allowedIncrease) {
1245 //cq = CoinMax(cq,10.0);
1246 // convert to (x+q)*(x+q) = w
1247 CoinWorkDouble q = bq / (1.0 * aq);
1248 CoinWorkDouble w = CoinMax(q * q + (cq / aq) * (allowedIncrease - 1.0), 0.0);
1249 w = CoinSqrt(w);
1250 CoinWorkDouble stepX = w - q;
1251 step = stepX;
1252 nextInf =
1253 multiplier * CoinSqrt(CoinMax(cq + step * bq + step * step * aq, 0.0));
1254#ifdef SOME_DEBUG
1255 printf ("with step of %g dualInf is %g\n",
1256 step, nextInf);
1257#endif
1258 actualDualStep_ = CoinMin(step, actualDualStep_);
1259 actualPrimalStep_ = CoinMin(step, actualPrimalStep_);
1260 }
1261 }
1262 } else {
1263 // probably pointless as linear
1264 // minimize complementarity
1265 // Complementarity gap will be a*change*change + b*change +c
1266 CoinWorkDouble a = 0.0;
1267 CoinWorkDouble b = 0.0;
1268 CoinWorkDouble c = 0.0;
1269 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1270 CoinWorkDouble oldPrimal = solution_[iColumn];
1271 if (!flagged(iColumn)) {
1272 if (lowerBound(iColumn)) {
1273 CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
1274 c += lowerSlack_[iColumn] * zVec_[iColumn];
1275 b += lowerSlack_[iColumn] * deltaZ_[iColumn] + zVec_[iColumn] * change;
1276 a += deltaZ_[iColumn] * change;
1277 }
1278 if (upperBound(iColumn)) {
1279 CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
1280 c += upperSlack_[iColumn] * wVec_[iColumn];
1281 b += upperSlack_[iColumn] * deltaW_[iColumn] + wVec_[iColumn] * change;
1282 a += deltaW_[iColumn] * change;
1283 }
1284 } else {
1285 // fixed
1286 if (lowerBound(iColumn)) {
1287 c += lowerSlack_[iColumn] * zVec_[iColumn];
1288 }
1289 if (upperBound(iColumn)) {
1290 c += upperSlack_[iColumn] * wVec_[iColumn];
1291 }
1292 }
1293 }
1294 // ? We want to minimize complementarityGap;
1295 // maybe use inf and do line search
1296 // To check see if matches at current step
1297 CoinWorkDouble step = CoinMin(actualPrimalStep_, actualDualStep_);
1298 CoinWorkDouble next = c + b * step + a * step * step;
1299#ifdef SOME_DEBUG
1300 printf("lin %g %g %g -> %g\n", a, b, c,
1301 c + b * step + a * step * step);
1302 debugMove(7, step, step);
1303#endif
1304 if (b > -1.0e-6) {
1305 if (phase == 0) {
1306#ifdef SOME_DEBUG
1307 printf("*** odd phase 0 direction\n");
1308#endif
1309 } else {
1310 directionNorm = -1.0;
1311 }
1312 }
1313 // and with ratio
1314 a = 0.0;
1315 b = 0.0;
1316 CoinWorkDouble ratio = actualDualStep_ / actualPrimalStep_;
1317 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1318 CoinWorkDouble oldPrimal = solution_[iColumn];
1319 if (!flagged(iColumn)) {
1320 if (lowerBound(iColumn)) {
1321 CoinWorkDouble change = oldPrimal + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
1322 b += lowerSlack_[iColumn] * deltaZ_[iColumn] * ratio + zVec_[iColumn] * change;
1323 a += deltaZ_[iColumn] * change * ratio;
1324 }
1325 if (upperBound(iColumn)) {
1326 CoinWorkDouble change = upper_[iColumn] - oldPrimal - deltaX_[iColumn] - upperSlack_[iColumn];
1327 b += upperSlack_[iColumn] * deltaW_[iColumn] * ratio + wVec_[iColumn] * change;
1328 a += deltaW_[iColumn] * change * ratio;
1329 }
1330 }
1331 }
1332 // ? We want to minimize complementarityGap;
1333 // maybe use inf and do line search
1334 // To check see if matches at current step
1335 step = actualPrimalStep_;
1336 CoinWorkDouble next2 = c + b * step + a * step * step;
1337 if (next2 > next) {
1338 actualPrimalStep_ = CoinMin(actualPrimalStep_, actualDualStep_);
1339 actualDualStep_ = actualPrimalStep_;
1340 }
1341#ifdef SOME_DEBUG
1342 printf("linb %g %g %g -> %g\n", a, b, c,
1343 c + b * step + a * step * step);
1344 debugMove(7, actualPrimalStep_, actualDualStep_);
1345#endif
1346 if (b > -1.0e-6) {
1347 if (phase == 0) {
1348#ifdef SOME_DEBUG
1349 printf("*** odd phase 0 direction\n");
1350#endif
1351 } else {
1352 directionNorm = -1.0;
1353 }
1354 }
1355 }
1356#else
1357 //actualPrimalStep_ =0.5*actualDualStep_;
1358#endif
1359#ifdef FULL_DEBUG
1360 if (phase == 3) {
1361 CoinWorkDouble minBeta = 0.1 * mu_;
1362 CoinWorkDouble maxBeta = 10.0 * mu_;
1363 for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
1364 if (!flagged(iColumn)) {
1365 if (lowerBound(iColumn)) {
1366 CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
1367 CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
1368 CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
1369 CoinWorkDouble gapProduct = dualValue * primalValue;
1370 if (delta2Z_[iColumn] < minBeta || delta2Z_[iColumn] > maxBeta)
1371 printf("3lower %d primal %g, dual %g, gap %g, old gap %g\n",
1372 iColumn, primalValue, dualValue, gapProduct, delta2Z_[iColumn]);
1373 }
1374 if (upperBound(iColumn)) {
1375 CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
1376 CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
1377 CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
1378 CoinWorkDouble gapProduct = dualValue * primalValue;
1379 if (delta2W_[iColumn] < minBeta || delta2W_[iColumn] > maxBeta)
1380 printf("3upper %d primal %g, dual %g, gap %g, old gap %g\n",
1381 iColumn, primalValue, dualValue, gapProduct, delta2W_[iColumn]);
1382 }
1383 }
1384 }
1385 }
1386#endif
1387#ifdef SOME_DEBUG_not
1388 {
1389 CoinWorkDouble largestL = 0.0;
1390 CoinWorkDouble smallestL = COIN_DBL_MAX;
1391 CoinWorkDouble largestU = 0.0;
1392 CoinWorkDouble smallestU = COIN_DBL_MAX;
1393 CoinWorkDouble sumL = 0.0;
1394 CoinWorkDouble sumU = 0.0;
1395 int nL = 0;
1396 int nU = 0;
1397 for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
1398 if (!flagged(iColumn)) {
1399 if (lowerBound(iColumn)) {
1400 CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
1401 CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
1402 CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
1403 CoinWorkDouble gapProduct = dualValue * primalValue;
1404 largestL = CoinMax(largestL, gapProduct);
1405 smallestL = CoinMin(smallestL, gapProduct);
1406 nL++;
1407 sumL += gapProduct;
1408 }
1409 if (upperBound(iColumn)) {
1410 CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
1411 CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
1412 CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
1413 CoinWorkDouble gapProduct = dualValue * primalValue;
1414 largestU = CoinMax(largestU, gapProduct);
1415 smallestU = CoinMin(smallestU, gapProduct);
1416 nU++;
1417 sumU += gapProduct;
1418 }
1419 }
1420 }
1421 CoinWorkDouble mu = (sumL + sumU) / (static_cast<CoinWorkDouble> (nL + nU));
1422
1423 CoinWorkDouble minBeta = 0.1 * mu;
1424 CoinWorkDouble maxBeta = 10.0 * mu;
1425 int nBL = 0;
1426 int nAL = 0;
1427 int nBU = 0;
1428 int nAU = 0;
1429 for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
1430 if (!flagged(iColumn)) {
1431 if (lowerBound(iColumn)) {
1432 CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
1433 CoinWorkDouble dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
1434 CoinWorkDouble primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
1435 CoinWorkDouble gapProduct = dualValue * primalValue;
1436 if (gapProduct < minBeta)
1437 nBL++;
1438 else if (gapProduct > maxBeta)
1439 nAL++;
1440 //if (gapProduct<0.1*minBeta)
1441 //printf("Lsmall one %d dual %g primal %g\n",iColumn,
1442 // dualValue,primalValue);
1443 }
1444 if (upperBound(iColumn)) {
1445 CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
1446 CoinWorkDouble dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
1447 CoinWorkDouble primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
1448 CoinWorkDouble gapProduct = dualValue * primalValue;
1449 if (gapProduct < minBeta)
1450 nBU++;
1451 else if (gapProduct > maxBeta)
1452 nAU++;
1453 //if (gapProduct<0.1*minBeta)
1454 //printf("Usmall one %d dual %g primal %g\n",iColumn,
1455 // dualValue,primalValue);
1456 }
1457 }
1458 }
1459 printf("phase %d new mu %.18g new gap %.18g\n", phase, mu, sumL + sumU);
1460 printf(" %d lower, smallest %.18g, %d below - largest %.18g, %d above\n",
1461 nL, smallestL, nBL, largestL, nAL);
1462 printf(" %d upper, smallest %.18g, %d below - largest %.18g, %d above\n",
1463 nU, smallestU, nBU, largestU, nAU);
1464 }
1465#endif
1466 return directionNorm;
1467}
1468/* Does solve. region1 is for deltaX (columns+rows), region2 for deltaPi (rows) */
1469void
1470ClpPredictorCorrector::solveSystem(CoinWorkDouble * region1, CoinWorkDouble * region2,
1471 const CoinWorkDouble * region1In, const CoinWorkDouble * region2In,
1472 const CoinWorkDouble * saveRegion1, const CoinWorkDouble * saveRegion2,
1473 bool gentleRefine)
1474{
1475 int iRow;
1476 int numberTotal = numberRows_ + numberColumns_;
1477 if (region2In) {
1478 // normal
1479 for (iRow = 0; iRow < numberRows_; iRow++)
1480 region2[iRow] = region2In[iRow];
1481 } else {
1482 // initial solution - (diagonal is 1 or 0)
1483 CoinZeroN(region2, numberRows_);
1484 }
1485 int iColumn;
1486 if (cholesky_->type() < 20) {
1487 // not KKT
1488 for (iColumn = 0; iColumn < numberTotal; iColumn++)
1489 region1[iColumn] = region1In[iColumn] * diagonal_[iColumn];
1490 multiplyAdd(region1 + numberColumns_, numberRows_, -1.0, region2, 1.0);
1491 matrix_->times(1.0, region1, region2);
1492 CoinWorkDouble maximumRHS = maximumAbsElement(region2, numberRows_);
1493 CoinWorkDouble scale = 1.0;
1494 CoinWorkDouble unscale = 1.0;
1495 if (maximumRHS > 1.0e-30) {
1496 if (maximumRHS <= 0.5) {
1497 CoinWorkDouble factor = 2.0;
1498 while (maximumRHS <= 0.5) {
1499 maximumRHS *= factor;
1500 scale *= factor;
1501 } /* endwhile */
1502 } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
1503 CoinWorkDouble factor = 0.5;
1504 while (maximumRHS >= 2.0) {
1505 maximumRHS *= factor;
1506 scale *= factor;
1507 } /* endwhile */
1508 }
1509 unscale = diagonalScaleFactor_ / scale;
1510 } else {
1511 //effectively zero
1512 scale = 0.0;
1513 unscale = 0.0;
1514 }
1515 multiplyAdd(NULL, numberRows_, 0.0, region2, scale);
1516 cholesky_->solve(region2);
1517 multiplyAdd(NULL, numberRows_, 0.0, region2, unscale);
1518 multiplyAdd(region2, numberRows_, -1.0, region1 + numberColumns_, 0.0);
1519 CoinZeroN(region1, numberColumns_);
1520 matrix_->transposeTimes(1.0, region2, region1);
1521 for (iColumn = 0; iColumn < numberTotal; iColumn++)
1522 region1[iColumn] = (region1[iColumn] - region1In[iColumn]) * diagonal_[iColumn];
1523 } else {
1524 for (iColumn = 0; iColumn < numberTotal; iColumn++)
1525 region1[iColumn] = region1In[iColumn];
1526 cholesky_->solveKKT(region1, region2, diagonal_, diagonalScaleFactor_);
1527 }
1528 if (saveRegion2) {
1529 //refine?
1530 CoinWorkDouble scaleX = 1.0;
1531 if (gentleRefine)
1532 scaleX = 0.8;
1533 multiplyAdd(saveRegion2, numberRows_, 1.0, region2, scaleX);
1534 assert (saveRegion1);
1535 multiplyAdd(saveRegion1, numberTotal, 1.0, region1, scaleX);
1536 }
1537}
1538// findDirectionVector.
1539CoinWorkDouble ClpPredictorCorrector::findDirectionVector(const int phase)
1540{
1541 CoinWorkDouble projectionTolerance = projectionTolerance_;
1542 //temporary
1543 //projectionTolerance=1.0e-15;
1544 CoinWorkDouble errorCheck = 0.9 * maximumRHSError_ / solutionNorm_;
1545 if (errorCheck > primalTolerance()) {
1546 if (errorCheck < projectionTolerance) {
1547 projectionTolerance = errorCheck;
1548 }
1549 } else {
1550 if (primalTolerance() < projectionTolerance) {
1551 projectionTolerance = primalTolerance();
1552 }
1553 }
1554 CoinWorkDouble * newError = new CoinWorkDouble [numberRows_];
1555 int numberTotal = numberRows_ + numberColumns_;
1556 //if flagged then entries zero so can do
1557 // For KKT separate out
1558 CoinWorkDouble * region1Save = NULL; //for refinement
1559 int iColumn;
1560 if (cholesky_->type() < 20) {
1561 int iColumn;
1562 for (iColumn = 0; iColumn < numberTotal; iColumn++)
1563 deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
1564 multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, deltaY_, 0.0);
1565 matrix_->times(1.0, deltaX_, deltaY_);
1566 } else {
1567 // regions in will be workArray and newError
1568 // regions out deltaX_ and deltaY_
1569 multiplyAdd(solution_ + numberColumns_, numberRows_, 1.0, newError, 0.0);
1570 matrix_->times(-1.0, solution_, newError);
1571 // This is inefficient but just for now get values which will be in deltay
1572 int iColumn;
1573 for (iColumn = 0; iColumn < numberTotal; iColumn++)
1574 deltaX_[iColumn] = workArray_[iColumn] - solution_[iColumn];
1575 multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, deltaY_, 0.0);
1576 matrix_->times(1.0, deltaX_, deltaY_);
1577 }
1578 bool goodSolve = false;
1579 CoinWorkDouble * regionSave = NULL; //for refinement
1580 int numberTries = 0;
1581 CoinWorkDouble relativeError = COIN_DBL_MAX;
1582 CoinWorkDouble tryError = 1.0e31;
1583 CoinWorkDouble saveMaximum = 0.0;
1584 double firstError = 0.0;
1585 double lastError2 = 0.0;
1586 while (!goodSolve && numberTries < 30) {
1587 CoinWorkDouble lastError = relativeError;
1588 goodSolve = true;
1589 CoinWorkDouble maximumRHS;
1590 maximumRHS = CoinMax(maximumAbsElement(deltaY_, numberRows_), 1.0e-12);
1591 if (!numberTries)
1592 saveMaximum = maximumRHS;
1593 if (cholesky_->type() < 20) {
1594 // no kkt
1595 CoinWorkDouble scale = 1.0;
1596 CoinWorkDouble unscale = 1.0;
1597 if (maximumRHS > 1.0e-30) {
1598 if (maximumRHS <= 0.5) {
1599 CoinWorkDouble factor = 2.0;
1600 while (maximumRHS <= 0.5) {
1601 maximumRHS *= factor;
1602 scale *= factor;
1603 } /* endwhile */
1604 } else if (maximumRHS >= 2.0 && maximumRHS <= COIN_DBL_MAX) {
1605 CoinWorkDouble factor = 0.5;
1606 while (maximumRHS >= 2.0) {
1607 maximumRHS *= factor;
1608 scale *= factor;
1609 } /* endwhile */
1610 }
1611 unscale = diagonalScaleFactor_ / scale;
1612 } else {
1613 //effectively zero
1614 scale = 0.0;
1615 unscale = 0.0;
1616 }
1617 //printf("--putting scales to 1.0\n");
1618 //scale=1.0;
1619 //unscale=1.0;
1620 multiplyAdd(NULL, numberRows_, 0.0, deltaY_, scale);
1621 cholesky_->solve(deltaY_);
1622 multiplyAdd(NULL, numberRows_, 0.0, deltaY_, unscale);
1623#if 0
1624 {
1625 printf("deltay\n");
1626 for (int i = 0; i < numberRows_; i++)
1627 printf("%d %.18g\n", i, deltaY_[i]);
1628 }
1629 exit(66);
1630#endif
1631 if (numberTries) {
1632 //refine?
1633 CoinWorkDouble scaleX = 1.0;
1634 if (lastError > 1.0e-5)
1635 scaleX = 0.8;
1636 multiplyAdd(regionSave, numberRows_, 1.0, deltaY_, scaleX);
1637 }
1638 //CoinZeroN(newError,numberRows_);
1639 multiplyAdd(deltaY_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
1640 CoinZeroN(deltaX_, numberColumns_);
1641 matrix_->transposeTimes(1.0, deltaY_, deltaX_);
1642 //if flagged then entries zero so can do
1643 for (iColumn = 0; iColumn < numberTotal; iColumn++)
1644 deltaX_[iColumn] = deltaX_[iColumn] * diagonal_[iColumn]
1645 - workArray_[iColumn];
1646 } else {
1647 // KKT
1648 solveSystem(deltaX_, deltaY_,
1649 workArray_, newError, region1Save, regionSave, lastError > 1.0e-5);
1650 }
1651 multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, newError, 0.0);
1652 matrix_->times(1.0, deltaX_, newError);
1653 numberTries++;
1654
1655 //now add in old Ax - doing extra checking
1656 CoinWorkDouble maximumRHSError = 0.0;
1657 CoinWorkDouble maximumRHSChange = 0.0;
1658 int iRow;
1659 char * dropped = cholesky_->rowsDropped();
1660 for (iRow = 0; iRow < numberRows_; iRow++) {
1661 if (!dropped[iRow]) {
1662 CoinWorkDouble newValue = newError[iRow];
1663 CoinWorkDouble oldValue = errorRegion_[iRow];
1664 //severity of errors depend on signs
1665 //**later */
1666 if (CoinAbs(newValue) > maximumRHSChange) {
1667 maximumRHSChange = CoinAbs(newValue);
1668 }
1669 CoinWorkDouble result = newValue + oldValue;
1670 if (CoinAbs(result) > maximumRHSError) {
1671 maximumRHSError = CoinAbs(result);
1672 }
1673 newError[iRow] = result;
1674 } else {
1675 CoinWorkDouble newValue = newError[iRow];
1676 CoinWorkDouble oldValue = errorRegion_[iRow];
1677 if (CoinAbs(newValue) > maximumRHSChange) {
1678 maximumRHSChange = CoinAbs(newValue);
1679 }
1680 CoinWorkDouble result = newValue + oldValue;
1681 newError[iRow] = result;
1682 //newError[iRow]=0.0;
1683 //assert(deltaY_[iRow]==0.0);
1684 deltaY_[iRow] = 0.0;
1685 }
1686 }
1687 relativeError = maximumRHSError / solutionNorm_;
1688 relativeError = maximumRHSError / saveMaximum;
1689 if (relativeError > tryError)
1690 relativeError = tryError;
1691 if (numberTries == 1)
1692 firstError = relativeError;
1693 if (relativeError < lastError) {
1694 lastError2 = relativeError;
1695 maximumRHSChange_ = maximumRHSChange;
1696 if (relativeError > projectionTolerance && numberTries <= 3) {
1697 //try and refine
1698 goodSolve = false;
1699 }
1700 //*** extra test here
1701 if (!goodSolve) {
1702 if (!regionSave) {
1703 regionSave = new CoinWorkDouble [numberRows_];
1704 if (cholesky_->type() >= 20)
1705 region1Save = new CoinWorkDouble [numberTotal];
1706 }
1707 CoinMemcpyN(deltaY_, numberRows_, regionSave);
1708 if (cholesky_->type() < 20) {
1709 // not KKT
1710 multiplyAdd(newError, numberRows_, -1.0, deltaY_, 0.0);
1711 } else {
1712 // KKT
1713 CoinMemcpyN(deltaX_, numberTotal, region1Save);
1714 // and back to input region
1715 CoinMemcpyN(deltaY_, numberRows_, newError);
1716 }
1717 }
1718 } else {
1719 //std::cout <<" worse residual = "<<relativeError;
1720 //bring back previous
1721 relativeError = lastError;
1722 if (regionSave) {
1723 CoinMemcpyN(regionSave, numberRows_, deltaY_);
1724 if (cholesky_->type() < 20) {
1725 // not KKT
1726 multiplyAdd(deltaY_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
1727 CoinZeroN(deltaX_, numberColumns_);
1728 matrix_->transposeTimes(1.0, deltaY_, deltaX_);
1729 //if flagged then entries zero so can do
1730 for (iColumn = 0; iColumn < numberTotal; iColumn++)
1731 deltaX_[iColumn] = deltaX_[iColumn] * diagonal_[iColumn]
1732 - workArray_[iColumn];
1733 } else {
1734 // KKT
1735 CoinMemcpyN(region1Save, numberTotal, deltaX_);
1736 }
1737 } else {
1738 // disaster
1739 CoinFillN(deltaX_, numberTotal, static_cast<CoinWorkDouble>(1.0));
1740 CoinFillN(deltaY_, numberRows_, static_cast<CoinWorkDouble>(1.0));
1741 COIN_DETAIL_PRINT(printf("bad cholesky\n"));
1742 }
1743 }
1744 } /* endwhile */
1745 if (firstError > 1.0e-8 || numberTries > 1) {
1746 handler_->message(CLP_BARRIER_ACCURACY, messages_)
1747 << phase << numberTries << static_cast<double>(firstError)
1748 << static_cast<double>(lastError2)
1749 << CoinMessageEol;
1750 }
1751 delete [] regionSave;
1752 delete [] region1Save;
1753 delete [] newError;
1754 // now rest
1755 CoinWorkDouble extra = eExtra;
1756 //multiplyAdd(deltaY_,numberRows_,1.0,deltaW_+numberColumns_,0.0);
1757 //CoinZeroN(deltaW_,numberColumns_);
1758 //matrix_->transposeTimes(-1.0,deltaY_,deltaW_);
1759
1760 for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
1761 deltaSU_[iColumn] = 0.0;
1762 deltaSL_[iColumn] = 0.0;
1763 deltaZ_[iColumn] = 0.0;
1764 CoinWorkDouble dd = deltaW_[iColumn];
1765 deltaW_[iColumn] = 0.0;
1766 if (!flagged(iColumn)) {
1767 CoinWorkDouble deltaX = deltaX_[iColumn];
1768 if (lowerBound(iColumn)) {
1769 CoinWorkDouble zValue = rhsZ_[iColumn];
1770 CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
1771 CoinWorkDouble slack = lowerSlack_[iColumn] + extra;
1772 deltaSL_[iColumn] = -rhsL_[iColumn] + deltaX;
1773 deltaZ_[iColumn] = (gHat - zVec_[iColumn] * deltaX) / slack;
1774 }
1775 if (upperBound(iColumn)) {
1776 CoinWorkDouble wValue = rhsW_[iColumn];
1777 CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
1778 CoinWorkDouble slack = upperSlack_[iColumn] + extra;
1779 deltaSU_[iColumn] = rhsU_[iColumn] - deltaX;
1780 deltaW_[iColumn] = (hHat + wVec_[iColumn] * deltaX) / slack;
1781 }
1782 if (0) {
1783 // different way of calculating
1784 CoinWorkDouble gamma2 = gamma_ * gamma_;
1785 CoinWorkDouble dZ = 0.0;
1786 CoinWorkDouble dW = 0.0;
1787 CoinWorkDouble zValue = rhsZ_[iColumn];
1788 CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
1789 CoinWorkDouble slackL = lowerSlack_[iColumn] + extra;
1790 CoinWorkDouble wValue = rhsW_[iColumn];
1791 CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
1792 CoinWorkDouble slackU = upperSlack_[iColumn] + extra;
1793 CoinWorkDouble q = rhsC_[iColumn] + gamma2 * deltaX + dd;
1794 if (primalR_)
1795 q += deltaX * primalR_[iColumn];
1796 dW = (gHat + hHat - slackL * q + (wValue - zValue) * deltaX) / (slackL + slackU);
1797 dZ = dW + q;
1798 //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
1799 //deltaW_[iColumn],dZ,dW);
1800 if (lowerBound(iColumn)) {
1801 if (upperBound(iColumn)) {
1802 //printf("B %d old %g %g new %g %g\n",iColumn,deltaZ_[iColumn],
1803 //deltaW_[iColumn],dZ,dW);
1804 deltaW_[iColumn] = dW;
1805 deltaZ_[iColumn] = dZ;
1806 } else {
1807 // just lower
1808 //printf("L %d old %g new %g\n",iColumn,deltaZ_[iColumn],
1809 //dZ);
1810 }
1811 } else {
1812 assert (upperBound(iColumn));
1813 //printf("U %d old %g new %g\n",iColumn,deltaW_[iColumn],
1814 //dW);
1815 }
1816 }
1817 }
1818 }
1819#if 0
1820 CoinWorkDouble * check = new CoinWorkDouble[numberTotal];
1821 // Check out rhsC_
1822 multiplyAdd(deltaY_, numberRows_, -1.0, check + numberColumns_, 0.0);
1823 CoinZeroN(check, numberColumns_);
1824 matrix_->transposeTimes(1.0, deltaY_, check);
1825 quadraticDjs(check, deltaX_, -1.0);
1826 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1827 check[iColumn] += deltaZ_[iColumn] - deltaW_[iColumn];
1828 if (CoinAbs(check[iColumn] - rhsC_[iColumn]) > 1.0e-3)
1829 printf("rhsC %d %g %g\n", iColumn, check[iColumn], rhsC_[iColumn]);
1830 }
1831 // Check out rhsZ_
1832 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1833 check[iColumn] += lowerSlack_[iColumn] * deltaZ_[iColumn] +
1834 zVec_[iColumn] * deltaSL_[iColumn];
1835 if (CoinAbs(check[iColumn] - rhsZ_[iColumn]) > 1.0e-3)
1836 printf("rhsZ %d %g %g\n", iColumn, check[iColumn], rhsZ_[iColumn]);
1837 }
1838 // Check out rhsW_
1839 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1840 check[iColumn] += upperSlack_[iColumn] * deltaW_[iColumn] +
1841 wVec_[iColumn] * deltaSU_[iColumn];
1842 if (CoinAbs(check[iColumn] - rhsW_[iColumn]) > 1.0e-3)
1843 printf("rhsW %d %g %g\n", iColumn, check[iColumn], rhsW_[iColumn]);
1844 }
1845 delete [] check;
1846#endif
1847 return relativeError;
1848}
1849// createSolution. Creates solution from scratch
1850int ClpPredictorCorrector::createSolution()
1851{
1852 int numberTotal = numberRows_ + numberColumns_;
1853 int iColumn;
1854 CoinWorkDouble tolerance = primalTolerance();
1855 // See if quadratic objective
1856#ifndef NO_RTTI
1857 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
1858#else
1859 ClpQuadraticObjective * quadraticObj = NULL;
1860 if (objective_->type() == 2)
1861 quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
1862#endif
1863 if (!quadraticObj) {
1864 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1865 if (upper_[iColumn] - lower_[iColumn] > tolerance)
1866 clearFixed(iColumn);
1867 else
1868 setFixed(iColumn);
1869 }
1870 } else {
1871 // try leaving fixed
1872 for (iColumn = 0; iColumn < numberTotal; iColumn++)
1873 clearFixed(iColumn);
1874 }
1875
1876 CoinWorkDouble maximumObjective = 0.0;
1877 CoinWorkDouble objectiveNorm2 = 0.0;
1878 getNorms(cost_, numberTotal, maximumObjective, objectiveNorm2);
1879 if (!maximumObjective) {
1880 maximumObjective = 1.0; // objective all zero
1881 }
1882 objectiveNorm2 = CoinSqrt(objectiveNorm2) / static_cast<CoinWorkDouble> (numberTotal);
1883 objectiveNorm_ = maximumObjective;
1884 scaleFactor_ = 1.0;
1885 if (maximumObjective > 0.0) {
1886 if (maximumObjective < 1.0) {
1887 scaleFactor_ = maximumObjective;
1888 } else if (maximumObjective > 1.0e4) {
1889 scaleFactor_ = maximumObjective / 1.0e4;
1890 }
1891 }
1892 if (scaleFactor_ != 1.0) {
1893 objectiveNorm2 *= scaleFactor_;
1894 multiplyAdd(NULL, numberTotal, 0.0, cost_, 1.0 / scaleFactor_);
1895 objectiveNorm_ = maximumObjective / scaleFactor_;
1896 }
1897 // See if quadratic objective
1898 if (quadraticObj) {
1899 // If scaled then really scale matrix
1900 CoinWorkDouble scaleFactor =
1901 scaleFactor_ * optimizationDirection_ * objectiveScale_ *
1902 rhsScale_;
1903 if ((scalingFlag_ > 0 && rowScale_) || scaleFactor != 1.0) {
1904 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1905 const int * columnQuadratic = quadratic->getIndices();
1906 const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1907 const int * columnQuadraticLength = quadratic->getVectorLengths();
1908 double * quadraticElement = quadratic->getMutableElements();
1909 int numberColumns = quadratic->getNumCols();
1910 CoinWorkDouble scale = 1.0 / scaleFactor;
1911 if (scalingFlag_ > 0 && rowScale_) {
1912 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1913 CoinWorkDouble scaleI = columnScale_[iColumn] * scale;
1914 for (CoinBigIndex j = columnQuadraticStart[iColumn];
1915 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1916 int jColumn = columnQuadratic[j];
1917 CoinWorkDouble scaleJ = columnScale_[jColumn];
1918 quadraticElement[j] *= scaleI * scaleJ;
1919 objectiveNorm_ = CoinMax(objectiveNorm_, CoinAbs(quadraticElement[j]));
1920 }
1921 }
1922 } else {
1923 // not scaled
1924 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1925 for (CoinBigIndex j = columnQuadraticStart[iColumn];
1926 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1927 quadraticElement[j] *= scale;
1928 objectiveNorm_ = CoinMax(objectiveNorm_, CoinAbs(quadraticElement[j]));
1929 }
1930 }
1931 }
1932 }
1933 }
1934 baseObjectiveNorm_ = objectiveNorm_;
1935 //accumulate fixed in dj region (as spare)
1936 //accumulate primal solution in primal region
1937 //DZ in lowerDual
1938 //DW in upperDual
1939 CoinWorkDouble infiniteCheck = 1.0e40;
1940 //CoinWorkDouble fakeCheck=1.0e10;
1941 //use deltaX region for work region
1942 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
1943 CoinWorkDouble primalValue = solution_[iColumn];
1944 clearFlagged(iColumn);
1945 clearFixedOrFree(iColumn);
1946 clearLowerBound(iColumn);
1947 clearUpperBound(iColumn);
1948 clearFakeLower(iColumn);
1949 clearFakeUpper(iColumn);
1950 if (!fixed(iColumn)) {
1951 dj_[iColumn] = 0.0;
1952 diagonal_[iColumn] = 1.0;
1953 deltaX_[iColumn] = 1.0;
1954 CoinWorkDouble lowerValue = lower_[iColumn];
1955 CoinWorkDouble upperValue = upper_[iColumn];
1956 if (lowerValue > -infiniteCheck) {
1957 if (upperValue < infiniteCheck) {
1958 //upper and lower bounds
1959 setLowerBound(iColumn);
1960 setUpperBound(iColumn);
1961 if (lowerValue >= 0.0) {
1962 solution_[iColumn] = lowerValue;
1963 } else if (upperValue <= 0.0) {
1964 solution_[iColumn] = upperValue;
1965 } else {
1966 solution_[iColumn] = 0.0;
1967 }
1968 } else {
1969 //just lower bound
1970 setLowerBound(iColumn);
1971 if (lowerValue >= 0.0) {
1972 solution_[iColumn] = lowerValue;
1973 } else {
1974 solution_[iColumn] = 0.0;
1975 }
1976 }
1977 } else {
1978 if (upperValue < infiniteCheck) {
1979 //just upper bound
1980 setUpperBound(iColumn);
1981 if (upperValue <= 0.0) {
1982 solution_[iColumn] = upperValue;
1983 } else {
1984 solution_[iColumn] = 0.0;
1985 }
1986 } else {
1987 //free
1988 setFixedOrFree(iColumn);
1989 solution_[iColumn] = 0.0;
1990 //std::cout<<" free "<<i<<std::endl;
1991 }
1992 }
1993 } else {
1994 setFlagged(iColumn);
1995 setFixedOrFree(iColumn);
1996 setLowerBound(iColumn);
1997 setUpperBound(iColumn);
1998 dj_[iColumn] = primalValue;
1999 solution_[iColumn] = lower_[iColumn];
2000 diagonal_[iColumn] = 0.0;
2001 deltaX_[iColumn] = 0.0;
2002 }
2003 }
2004 // modify fixed RHS
2005 multiplyAdd(dj_ + numberColumns_, numberRows_, -1.0, rhsFixRegion_, 0.0);
2006 // create plausible RHS?
2007 matrix_->times(-1.0, dj_, rhsFixRegion_);
2008 multiplyAdd(solution_ + numberColumns_, numberRows_, 1.0, errorRegion_, 0.0);
2009 matrix_->times(-1.0, solution_, errorRegion_);
2010 rhsNorm_ = maximumAbsElement(errorRegion_, numberRows_);
2011 if (rhsNorm_ < 1.0) {
2012 rhsNorm_ = 1.0;
2013 }
2014 int * rowsDropped = new int [numberRows_];
2015 int returnCode = cholesky_->factorize(diagonal_, rowsDropped);
2016 if (returnCode == -1) {
2017 COIN_DETAIL_PRINT(printf("Out of memory\n"));
2018 problemStatus_ = 4;
2019 return -1;
2020 }
2021 if (cholesky_->status()) {
2022 std::cout << "singular on initial cholesky?" << std::endl;
2023 cholesky_->resetRowsDropped();
2024 //cholesky_->factorize(rowDropped_);
2025 //if (cholesky_->status()) {
2026 //std::cout << "bad cholesky??? (after retry)" <<std::endl;
2027 //abort();
2028 //}
2029 }
2030 delete [] rowsDropped;
2031 if (cholesky_->type() < 20) {
2032 // not KKT
2033 cholesky_->solve(errorRegion_);
2034 //create information for solution
2035 multiplyAdd(errorRegion_, numberRows_, -1.0, deltaX_ + numberColumns_, 0.0);
2036 CoinZeroN(deltaX_, numberColumns_);
2037 matrix_->transposeTimes(1.0, errorRegion_, deltaX_);
2038 } else {
2039 // KKT
2040 // reverse sign on solution
2041 multiplyAdd(NULL, numberRows_ + numberColumns_, 0.0, solution_, -1.0);
2042 solveSystem(deltaX_, errorRegion_, solution_, NULL, NULL, NULL, false);
2043 }
2044 CoinWorkDouble initialValue = 1.0e2;
2045 if (rhsNorm_ * 1.0e-2 > initialValue) {
2046 initialValue = rhsNorm_ * 1.0e-2;
2047 }
2048 //initialValue = CoinMax(1.0,rhsNorm_);
2049 CoinWorkDouble smallestBoundDifference = COIN_DBL_MAX;
2050 CoinWorkDouble * fakeSolution = deltaX_;
2051 for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
2052 if (!flagged(iColumn)) {
2053 if (lower_[iColumn] - fakeSolution[iColumn] > initialValue) {
2054 initialValue = lower_[iColumn] - fakeSolution[iColumn];
2055 }
2056 if (fakeSolution[iColumn] - upper_[iColumn] > initialValue) {
2057 initialValue = fakeSolution[iColumn] - upper_[iColumn];
2058 }
2059 if (upper_[iColumn] - lower_[iColumn] < smallestBoundDifference) {
2060 smallestBoundDifference = upper_[iColumn] - lower_[iColumn];
2061 }
2062 }
2063 }
2064 solutionNorm_ = 1.0e-12;
2065 handler_->message(CLP_BARRIER_SAFE, messages_)
2066 << static_cast<double>(initialValue) << static_cast<double>(objectiveNorm_)
2067 << CoinMessageEol;
2068 CoinWorkDouble extra = 1.0e-10;
2069 CoinWorkDouble largeGap = 1.0e15;
2070 //CoinWorkDouble safeObjectiveValue=2.0*objectiveNorm_;
2071 CoinWorkDouble safeObjectiveValue = objectiveNorm_ + 1.0;
2072 CoinWorkDouble safeFree = 1.0e-1 * initialValue;
2073 //printf("normal safe dual value of %g, primal value of %g\n",
2074 // safeObjectiveValue,initialValue);
2075 //safeObjectiveValue=CoinMax(2.0,1.0e-1*safeObjectiveValue);
2076 //initialValue=CoinMax(100.0,1.0e-1*initialValue);
2077 //printf("temp safe dual value of %g, primal value of %g\n",
2078 // safeObjectiveValue,initialValue);
2079 CoinWorkDouble zwLarge = 1.0e2 * initialValue;
2080 //zwLarge=1.0e40;
2081 if (cholesky_->choleskyCondition() < 0.0 && cholesky_->type() < 20) {
2082 // looks bad - play safe
2083 initialValue *= 10.0;
2084 safeObjectiveValue *= 10.0;
2085 safeFree *= 10.0;
2086 }
2087 CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
2088 // First do primal side
2089 for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
2090 if (!flagged(iColumn)) {
2091 CoinWorkDouble lowerValue = lower_[iColumn];
2092 CoinWorkDouble upperValue = upper_[iColumn];
2093 CoinWorkDouble newValue;
2094 CoinWorkDouble setPrimal = initialValue;
2095 if (quadraticObj) {
2096 // perturb primal solution a bit
2097 //fakeSolution[iColumn] *= 0.002*CoinDrand48()+0.999;
2098 }
2099 if (lowerBound(iColumn)) {
2100 if (upperBound(iColumn)) {
2101 //upper and lower bounds
2102 if (upperValue - lowerValue > 2.0 * setPrimal) {
2103 CoinWorkDouble fakeValue = fakeSolution[iColumn];
2104 if (fakeValue < lowerValue + setPrimal) {
2105 fakeValue = lowerValue + setPrimal;
2106 }
2107 if (fakeValue > upperValue - setPrimal) {
2108 fakeValue = upperValue - setPrimal;
2109 }
2110 newValue = fakeValue;
2111 } else {
2112 newValue = 0.5 * (upperValue + lowerValue);
2113 }
2114 } else {
2115 //just lower bound
2116 CoinWorkDouble fakeValue = fakeSolution[iColumn];
2117 if (fakeValue < lowerValue + setPrimal) {
2118 fakeValue = lowerValue + setPrimal;
2119 }
2120 newValue = fakeValue;
2121 }
2122 } else {
2123 if (upperBound(iColumn)) {
2124 //just upper bound
2125 CoinWorkDouble fakeValue = fakeSolution[iColumn];
2126 if (fakeValue > upperValue - setPrimal) {
2127 fakeValue = upperValue - setPrimal;
2128 }
2129 newValue = fakeValue;
2130 } else {
2131 //free
2132 newValue = fakeSolution[iColumn];
2133 if (newValue >= 0.0) {
2134 if (newValue < safeFree) {
2135 newValue = safeFree;
2136 }
2137 } else {
2138 if (newValue > -safeFree) {
2139 newValue = -safeFree;
2140 }
2141 }
2142 }
2143 }
2144 solution_[iColumn] = newValue;
2145 } else {
2146 // fixed
2147 lowerSlack_[iColumn] = 0.0;
2148 upperSlack_[iColumn] = 0.0;
2149 solution_[iColumn] = lower_[iColumn];
2150 zVec_[iColumn] = 0.0;
2151 wVec_[iColumn] = 0.0;
2152 diagonal_[iColumn] = 0.0;
2153 }
2154 }
2155 solutionNorm_ = maximumAbsElement(solution_, numberTotal);
2156 // Set bounds and do dj including quadratic
2157 largeGap = CoinMax(1.0e7, 1.02 * solutionNorm_);
2158 CoinPackedMatrix * quadratic = NULL;
2159 const int * columnQuadratic = NULL;
2160 const CoinBigIndex * columnQuadraticStart = NULL;
2161 const int * columnQuadraticLength = NULL;
2162 const double * quadraticElement = NULL;
2163 if (quadraticObj) {
2164 quadratic = quadraticObj->quadraticObjective();
2165 columnQuadratic = quadratic->getIndices();
2166 columnQuadraticStart = quadratic->getVectorStarts();
2167 columnQuadraticLength = quadratic->getVectorLengths();
2168 quadraticElement = quadratic->getElements();
2169 }
2170 CoinWorkDouble quadraticNorm = 0.0;
2171 for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
2172 if (!flagged(iColumn)) {
2173 CoinWorkDouble primalValue = solution_[iColumn];
2174 CoinWorkDouble lowerValue = lower_[iColumn];
2175 CoinWorkDouble upperValue = upper_[iColumn];
2176 // Do dj
2177 CoinWorkDouble reducedCost = cost_[iColumn];
2178 if (lowerBound(iColumn)) {
2179 reducedCost += linearPerturbation_;
2180 }
2181 if (upperBound(iColumn)) {
2182 reducedCost -= linearPerturbation_;
2183 }
2184 if (quadraticObj && iColumn < numberColumns_) {
2185 for (CoinBigIndex j = columnQuadraticStart[iColumn];
2186 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
2187 int jColumn = columnQuadratic[j];
2188 CoinWorkDouble valueJ = solution_[jColumn];
2189 CoinWorkDouble elementValue = quadraticElement[j];
2190 reducedCost += valueJ * elementValue;
2191 }
2192 quadraticNorm = CoinMax(quadraticNorm, CoinAbs(reducedCost));
2193 }
2194 dj_[iColumn] = reducedCost;
2195 if (primalValue > lowerValue + largeGap && primalValue < upperValue - largeGap) {
2196 clearFixedOrFree(iColumn);
2197 setLowerBound(iColumn);
2198 setUpperBound(iColumn);
2199 lowerValue = CoinMax(lowerValue, primalValue - largeGap);
2200 upperValue = CoinMin(upperValue, primalValue + largeGap);
2201 lower_[iColumn] = lowerValue;
2202 upper_[iColumn] = upperValue;
2203 }
2204 }
2205 }
2206 safeObjectiveValue = CoinMax(safeObjectiveValue, quadraticNorm);
2207 for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
2208 if (!flagged(iColumn)) {
2209 CoinWorkDouble primalValue = solution_[iColumn];
2210 CoinWorkDouble lowerValue = lower_[iColumn];
2211 CoinWorkDouble upperValue = upper_[iColumn];
2212 CoinWorkDouble reducedCost = dj_[iColumn];
2213 CoinWorkDouble low = 0.0;
2214 CoinWorkDouble high = 0.0;
2215 if (lowerBound(iColumn)) {
2216 if (upperBound(iColumn)) {
2217 //upper and lower bounds
2218 if (upperValue - lowerValue > 2.0 * initialValue) {
2219 low = primalValue - lowerValue;
2220 high = upperValue - primalValue;
2221 } else {
2222 low = initialValue;
2223 high = initialValue;
2224 }
2225 CoinWorkDouble s = low + extra;
2226 CoinWorkDouble ratioZ;
2227 if (s < zwLarge) {
2228 ratioZ = 1.0;
2229 } else {
2230 ratioZ = CoinSqrt(zwLarge / s);
2231 }
2232 CoinWorkDouble t = high + extra;
2233 CoinWorkDouble ratioT;
2234 if (t < zwLarge) {
2235 ratioT = 1.0;
2236 } else {
2237 ratioT = CoinSqrt(zwLarge / t);
2238 }
2239 //modify s and t
2240 if (s > largeGap) {
2241 s = largeGap;
2242 }
2243 if (t > largeGap) {
2244 t = largeGap;
2245 }
2246 //modify if long long way away from bound
2247 if (reducedCost >= 0.0) {
2248 zVec_[iColumn] = reducedCost + safeObjectiveValue * ratioZ;
2249 zVec_[iColumn] = CoinMax(reducedCost, safeObjectiveValue * ratioZ);
2250 wVec_[iColumn] = safeObjectiveValue * ratioT;
2251 } else {
2252 zVec_[iColumn] = safeObjectiveValue * ratioZ;
2253 wVec_[iColumn] = -reducedCost + safeObjectiveValue * ratioT;
2254 wVec_[iColumn] = CoinMax(-reducedCost , safeObjectiveValue * ratioT);
2255 }
2256 CoinWorkDouble gammaTerm = gamma2;
2257 if (primalR_)
2258 gammaTerm += primalR_[iColumn];
2259 diagonal_[iColumn] = (t * s) /
2260 (s * wVec_[iColumn] + t * zVec_[iColumn] + gammaTerm * t * s);
2261 } else {
2262 //just lower bound
2263 low = primalValue - lowerValue;
2264 high = 0.0;
2265 CoinWorkDouble s = low + extra;
2266 CoinWorkDouble ratioZ;
2267 if (s < zwLarge) {
2268 ratioZ = 1.0;
2269 } else {
2270 ratioZ = CoinSqrt(zwLarge / s);
2271 }
2272 //modify s
2273 if (s > largeGap) {
2274 s = largeGap;
2275 }
2276 if (reducedCost >= 0.0) {
2277 zVec_[iColumn] = reducedCost + safeObjectiveValue * ratioZ;
2278 zVec_[iColumn] = CoinMax(reducedCost , safeObjectiveValue * ratioZ);
2279 wVec_[iColumn] = 0.0;
2280 } else {
2281 zVec_[iColumn] = safeObjectiveValue * ratioZ;
2282 wVec_[iColumn] = 0.0;
2283 }
2284 CoinWorkDouble gammaTerm = gamma2;
2285 if (primalR_)
2286 gammaTerm += primalR_[iColumn];
2287 diagonal_[iColumn] = s / (zVec_[iColumn] + s * gammaTerm);
2288 }
2289 } else {
2290 if (upperBound(iColumn)) {
2291 //just upper bound
2292 low = 0.0;
2293 high = upperValue - primalValue;
2294 CoinWorkDouble t = high + extra;
2295 CoinWorkDouble ratioT;
2296 if (t < zwLarge) {
2297 ratioT = 1.0;
2298 } else {
2299 ratioT = CoinSqrt(zwLarge / t);
2300 }
2301 //modify t
2302 if (t > largeGap) {
2303 t = largeGap;
2304 }
2305 if (reducedCost >= 0.0) {
2306 zVec_[iColumn] = 0.0;
2307 wVec_[iColumn] = safeObjectiveValue * ratioT;
2308 } else {
2309 zVec_[iColumn] = 0.0;
2310 wVec_[iColumn] = -reducedCost + safeObjectiveValue * ratioT;
2311 wVec_[iColumn] = CoinMax(-reducedCost , safeObjectiveValue * ratioT);
2312 }
2313 CoinWorkDouble gammaTerm = gamma2;
2314 if (primalR_)
2315 gammaTerm += primalR_[iColumn];
2316 diagonal_[iColumn] = t / (wVec_[iColumn] + t * gammaTerm);
2317 }
2318 }
2319 lowerSlack_[iColumn] = low;
2320 upperSlack_[iColumn] = high;
2321 }
2322 }
2323#if 0
2324 if (solution_[0] > 0.0) {
2325 for (int i = 0; i < numberTotal; i++)
2326 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(solution_[i]),
2327 diagonal_[i], CoinAbs(dj_[i]),
2328 lowerSlack_[i], zVec_[i],
2329 upperSlack_[i], wVec_[i]);
2330 } else {
2331 for (int i = 0; i < numberTotal; i++)
2332 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(solution_[i]),
2333 diagonal_[i], CoinAbs(dj_[i]),
2334 upperSlack_[i], wVec_[i],
2335 lowerSlack_[i], zVec_[i] );
2336 }
2337 exit(66);
2338#endif
2339 return 0;
2340}
2341// complementarityGap. Computes gap
2342//phase 0=as is , 1 = after predictor , 2 after corrector
2343CoinWorkDouble ClpPredictorCorrector::complementarityGap(int & numberComplementarityPairs,
2344 int & numberComplementarityItems,
2345 const int phase)
2346{
2347 CoinWorkDouble gap = 0.0;
2348 //seems to be same coding for phase = 1 or 2
2349 numberComplementarityPairs = 0;
2350 numberComplementarityItems = 0;
2351 int numberTotal = numberRows_ + numberColumns_;
2352 CoinWorkDouble toleranceGap = 0.0;
2353 CoinWorkDouble largestGap = 0.0;
2354 CoinWorkDouble smallestGap = COIN_DBL_MAX;
2355 //seems to be same coding for phase = 1 or 2
2356 int numberNegativeGaps = 0;
2357 CoinWorkDouble sumNegativeGap = 0.0;
2358 CoinWorkDouble largeGap = 1.0e2 * solutionNorm_;
2359 if (largeGap < 1.0e10) {
2360 largeGap = 1.0e10;
2361 }
2362 largeGap = 1.0e30;
2363 CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance];
2364 CoinWorkDouble primalTolerance = dblParam_[ClpPrimalTolerance];
2365 dualTolerance = dualTolerance / scaleFactor_;
2366 for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
2367 if (!fixedOrFree(iColumn)) {
2368 numberComplementarityPairs++;
2369 //can collapse as if no lower bound both zVec and deltaZ 0.0
2370 CoinWorkDouble newZ = 0.0;
2371 CoinWorkDouble newW = 0.0;
2372 if (lowerBound(iColumn)) {
2373 numberComplementarityItems++;
2374 CoinWorkDouble dualValue;
2375 CoinWorkDouble primalValue;
2376 if (!phase) {
2377 dualValue = zVec_[iColumn];
2378 primalValue = lowerSlack_[iColumn];
2379 } else {
2380 CoinWorkDouble change;
2381 change = solution_[iColumn] + deltaX_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn];
2382 dualValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
2383 newZ = dualValue;
2384 primalValue = lowerSlack_[iColumn] + actualPrimalStep_ * change;
2385 }
2386 //reduce primalValue
2387 if (primalValue > largeGap) {
2388 primalValue = largeGap;
2389 }
2390 CoinWorkDouble gapProduct = dualValue * primalValue;
2391 if (gapProduct < 0.0) {
2392 //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
2393 //primalValue<<endl;
2394 numberNegativeGaps++;
2395 sumNegativeGap -= gapProduct;
2396 gapProduct = 0.0;
2397 }
2398 gap += gapProduct;
2399 //printf("l %d prim %g dual %g totalGap %g\n",
2400 // iColumn,primalValue,dualValue,gap);
2401 if (gapProduct > largestGap) {
2402 largestGap = gapProduct;
2403 }
2404 smallestGap = CoinMin(smallestGap, gapProduct);
2405 if (dualValue > dualTolerance && primalValue > primalTolerance) {
2406 toleranceGap += dualValue * primalValue;
2407 }
2408 }
2409 if (upperBound(iColumn)) {
2410 numberComplementarityItems++;
2411 CoinWorkDouble dualValue;
2412 CoinWorkDouble primalValue;
2413 if (!phase) {
2414 dualValue = wVec_[iColumn];
2415 primalValue = upperSlack_[iColumn];
2416 } else {
2417 CoinWorkDouble change;
2418 change = upper_[iColumn] - solution_[iColumn] - deltaX_[iColumn] - upperSlack_[iColumn];
2419 dualValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
2420 newW = dualValue;
2421 primalValue = upperSlack_[iColumn] + actualPrimalStep_ * change;
2422 }
2423 //reduce primalValue
2424 if (primalValue > largeGap) {
2425 primalValue = largeGap;
2426 }
2427 CoinWorkDouble gapProduct = dualValue * primalValue;
2428 if (gapProduct < 0.0) {
2429 //cout<<"negative gap component "<<iColumn<<" "<<dualValue<<" "<<
2430 //primalValue<<endl;
2431 numberNegativeGaps++;
2432 sumNegativeGap -= gapProduct;
2433 gapProduct = 0.0;
2434 }
2435 gap += gapProduct;
2436 //printf("u %d prim %g dual %g totalGap %g\n",
2437 // iColumn,primalValue,dualValue,gap);
2438 if (gapProduct > largestGap) {
2439 largestGap = gapProduct;
2440 }
2441 if (dualValue > dualTolerance && primalValue > primalTolerance) {
2442 toleranceGap += dualValue * primalValue;
2443 }
2444 }
2445 }
2446 }
2447 //if (numberIterations_>4)
2448 //exit(9);
2449 if (!phase && numberNegativeGaps) {
2450 handler_->message(CLP_BARRIER_NEGATIVE_GAPS, messages_)
2451 << numberNegativeGaps << static_cast<double>(sumNegativeGap)
2452 << CoinMessageEol;
2453 }
2454
2455 //in case all free!
2456 if (!numberComplementarityPairs) {
2457 numberComplementarityPairs = 1;
2458 }
2459#ifdef SOME_DEBUG
2460 printf("with d,p steps %g,%g gap %g - smallest %g, largest %g, pairs %d\n",
2461 actualDualStep_, actualPrimalStep_,
2462 gap, smallestGap, largestGap, numberComplementarityPairs);
2463#endif
2464 return gap;
2465}
2466// setupForSolve.
2467//phase 0=affine , 1 = corrector , 2 = primal-dual
2468void ClpPredictorCorrector::setupForSolve(const int phase)
2469{
2470 CoinWorkDouble extra = eExtra;
2471 int numberTotal = numberRows_ + numberColumns_;
2472 int iColumn;
2473#ifdef SOME_DEBUG
2474 printf("phase %d in setupForSolve, mu %.18g\n", phase, mu_);
2475#endif
2476 CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
2477 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
2478 switch (phase) {
2479 case 0:
2480 CoinMemcpyN(errorRegion_, numberRows_, rhsB_);
2481 if (delta_ || dualR_) {
2482 // add in regularization
2483 CoinWorkDouble delta2 = delta_ * delta_;
2484 for (int iRow = 0; iRow < numberRows_; iRow++) {
2485 rhsB_[iRow] -= delta2 * dualArray[iRow];
2486 if (dualR_)
2487 rhsB_[iRow] -= dualR_[iRow] * dualArray[iRow];
2488 }
2489 }
2490 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
2491 rhsC_[iColumn] = 0.0;
2492 rhsU_[iColumn] = 0.0;
2493 rhsL_[iColumn] = 0.0;
2494 rhsZ_[iColumn] = 0.0;
2495 rhsW_[iColumn] = 0.0;
2496 if (!flagged(iColumn)) {
2497 rhsC_[iColumn] = dj_[iColumn] - zVec_[iColumn] + wVec_[iColumn];
2498 rhsC_[iColumn] += gamma2 * solution_[iColumn];
2499 if (primalR_)
2500 rhsC_[iColumn] += primalR_[iColumn] * solution_[iColumn];
2501 if (lowerBound(iColumn)) {
2502 rhsZ_[iColumn] = -zVec_[iColumn] * (lowerSlack_[iColumn] + extra);
2503 rhsL_[iColumn] = CoinMax(0.0, (lower_[iColumn] + lowerSlack_[iColumn]) - solution_[iColumn]);
2504 }
2505 if (upperBound(iColumn)) {
2506 rhsW_[iColumn] = -wVec_[iColumn] * (upperSlack_[iColumn] + extra);
2507 rhsU_[iColumn] = CoinMin(0.0, (upper_[iColumn] - upperSlack_[iColumn]) - solution_[iColumn]);
2508 }
2509 }
2510 }
2511#if 0
2512 for (int i = 0; i < 3; i++) {
2513 if (!CoinAbs(rhsZ_[i]))
2514 rhsZ_[i] = 0.0;
2515 if (!CoinAbs(rhsW_[i]))
2516 rhsW_[i] = 0.0;
2517 if (!CoinAbs(rhsU_[i]))
2518 rhsU_[i] = 0.0;
2519 if (!CoinAbs(rhsL_[i]))
2520 rhsL_[i] = 0.0;
2521 }
2522 if (solution_[0] > 0.0) {
2523 for (int i = 0; i < 3; i++)
2524 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, solution_[i],
2525 diagonal_[i], dj_[i],
2526 lowerSlack_[i], zVec_[i],
2527 upperSlack_[i], wVec_[i]);
2528 for (int i = 0; i < 3; i++)
2529 printf("%d %.18g %.18g %.18g %.18g %.18g\n", i, rhsC_[i],
2530 rhsZ_[i], rhsL_[i],
2531 rhsW_[i], rhsU_[i]);
2532 } else {
2533 for (int i = 0; i < 3; i++)
2534 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, solution_[i],
2535 diagonal_[i], dj_[i],
2536 lowerSlack_[i], zVec_[i],
2537 upperSlack_[i], wVec_[i]);
2538 for (int i = 0; i < 3; i++)
2539 printf("%d %.18g %.18g %.18g %.18g %.18g\n", i, rhsC_[i],
2540 rhsZ_[i], rhsL_[i],
2541 rhsW_[i], rhsU_[i]);
2542 }
2543#endif
2544 break;
2545 case 1:
2546 // could be stored in delta2?
2547 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
2548 rhsZ_[iColumn] = 0.0;
2549 rhsW_[iColumn] = 0.0;
2550 if (!flagged(iColumn)) {
2551 if (lowerBound(iColumn)) {
2552 rhsZ_[iColumn] = mu_ - zVec_[iColumn] * (lowerSlack_[iColumn] + extra)
2553 - deltaZ_[iColumn] * deltaX_[iColumn];
2554 // To bring in line with OSL
2555 rhsZ_[iColumn] += deltaZ_[iColumn] * rhsL_[iColumn];
2556 }
2557 if (upperBound(iColumn)) {
2558 rhsW_[iColumn] = mu_ - wVec_[iColumn] * (upperSlack_[iColumn] + extra)
2559 + deltaW_[iColumn] * deltaX_[iColumn];
2560 // To bring in line with OSL
2561 rhsW_[iColumn] -= deltaW_[iColumn] * rhsU_[iColumn];
2562 }
2563 }
2564 }
2565#if 0
2566 for (int i = 0; i < numberTotal; i++) {
2567 if (!CoinAbs(rhsZ_[i]))
2568 rhsZ_[i] = 0.0;
2569 if (!CoinAbs(rhsW_[i]))
2570 rhsW_[i] = 0.0;
2571 if (!CoinAbs(rhsU_[i]))
2572 rhsU_[i] = 0.0;
2573 if (!CoinAbs(rhsL_[i]))
2574 rhsL_[i] = 0.0;
2575 }
2576 if (solution_[0] > 0.0) {
2577 for (int i = 0; i < numberTotal; i++)
2578 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(solution_[i]),
2579 diagonal_[i], CoinAbs(dj_[i]),
2580 lowerSlack_[i], zVec_[i],
2581 upperSlack_[i], wVec_[i]);
2582 for (int i = 0; i < numberTotal; i++)
2583 printf("%d %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(rhsC_[i]),
2584 rhsZ_[i], rhsL_[i],
2585 rhsW_[i], rhsU_[i]);
2586 } else {
2587 for (int i = 0; i < numberTotal; i++)
2588 printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(solution_[i]),
2589 diagonal_[i], CoinAbs(dj_[i]),
2590 upperSlack_[i], wVec_[i],
2591 lowerSlack_[i], zVec_[i] );
2592 for (int i = 0; i < numberTotal; i++)
2593 printf("%d %.18g %.18g %.18g %.18g %.18g\n", i, CoinAbs(rhsC_[i]),
2594 rhsW_[i], rhsU_[i],
2595 rhsZ_[i], rhsL_[i]);
2596 }
2597 exit(66);
2598#endif
2599 break;
2600 case 2:
2601 CoinMemcpyN(errorRegion_, numberRows_, rhsB_);
2602 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
2603 rhsZ_[iColumn] = 0.0;
2604 rhsW_[iColumn] = 0.0;
2605 if (!flagged(iColumn)) {
2606 if (lowerBound(iColumn)) {
2607 rhsZ_[iColumn] = mu_ - zVec_[iColumn] * (lowerSlack_[iColumn] + extra);
2608 }
2609 if (upperBound(iColumn)) {
2610 rhsW_[iColumn] = mu_ - wVec_[iColumn] * (upperSlack_[iColumn] + extra);
2611 }
2612 }
2613 }
2614 break;
2615 case 3: {
2616 CoinWorkDouble minBeta = 0.1 * mu_;
2617 CoinWorkDouble maxBeta = 10.0 * mu_;
2618 CoinWorkDouble dualStep = CoinMin(1.0, actualDualStep_ + 0.1);
2619 CoinWorkDouble primalStep = CoinMin(1.0, actualPrimalStep_ + 0.1);
2620#ifdef SOME_DEBUG
2621 printf("good complementarity range %g to %g\n", minBeta, maxBeta);
2622#endif
2623 //minBeta=0.0;
2624 //maxBeta=COIN_DBL_MAX;
2625 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
2626 if (!flagged(iColumn)) {
2627 if (lowerBound(iColumn)) {
2628 CoinWorkDouble change = -rhsL_[iColumn] + deltaX_[iColumn];
2629 CoinWorkDouble dualValue = zVec_[iColumn] + dualStep * deltaZ_[iColumn];
2630 CoinWorkDouble primalValue = lowerSlack_[iColumn] + primalStep * change;
2631 CoinWorkDouble gapProduct = dualValue * primalValue;
2632 if (gapProduct > 0.0 && dualValue < 0.0)
2633 gapProduct = - gapProduct;
2634#ifdef FULL_DEBUG
2635 delta2Z_[iColumn] = gapProduct;
2636 if (delta2Z_[iColumn] < minBeta || delta2Z_[iColumn] > maxBeta)
2637 printf("lower %d primal %g, dual %g, gap %g\n",
2638 iColumn, primalValue, dualValue, gapProduct);
2639#endif
2640 CoinWorkDouble value = 0.0;
2641 if (gapProduct < minBeta) {
2642 value = 2.0 * (minBeta - gapProduct);
2643 value = (mu_ - gapProduct);
2644 value = (minBeta - gapProduct);
2645 assert (value > 0.0);
2646 } else if (gapProduct > maxBeta) {
2647 value = CoinMax(maxBeta - gapProduct, -maxBeta);
2648 assert (value < 0.0);
2649 }
2650 rhsZ_[iColumn] += value;
2651 }
2652 if (upperBound(iColumn)) {
2653 CoinWorkDouble change = rhsU_[iColumn] - deltaX_[iColumn];
2654 CoinWorkDouble dualValue = wVec_[iColumn] + dualStep * deltaW_[iColumn];
2655 CoinWorkDouble primalValue = upperSlack_[iColumn] + primalStep * change;
2656 CoinWorkDouble gapProduct = dualValue * primalValue;
2657 if (gapProduct > 0.0 && dualValue < 0.0)
2658 gapProduct = - gapProduct;
2659#ifdef FULL_DEBUG
2660 delta2W_[iColumn] = gapProduct;
2661 if (delta2W_[iColumn] < minBeta || delta2W_[iColumn] > maxBeta)
2662 printf("upper %d primal %g, dual %g, gap %g\n",
2663 iColumn, primalValue, dualValue, gapProduct);
2664#endif
2665 CoinWorkDouble value = 0.0;
2666 if (gapProduct < minBeta) {
2667 value = (minBeta - gapProduct);
2668 assert (value > 0.0);
2669 } else if (gapProduct > maxBeta) {
2670 value = CoinMax(maxBeta - gapProduct, -maxBeta);
2671 assert (value < 0.0);
2672 }
2673 rhsW_[iColumn] += value;
2674 }
2675 }
2676 }
2677 }
2678 break;
2679 } /* endswitch */
2680 if (cholesky_->type() < 20) {
2681 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
2682 CoinWorkDouble value = rhsC_[iColumn];
2683 CoinWorkDouble zValue = rhsZ_[iColumn];
2684 CoinWorkDouble wValue = rhsW_[iColumn];
2685#if 0
2686#if 1
2687 if (phase == 0) {
2688 // more accurate
2689 value = dj[iColumn];
2690 zValue = 0.0;
2691 wValue = 0.0;
2692 } else if (phase == 2) {
2693 // more accurate
2694 value = dj[iColumn];
2695 zValue = mu_;
2696 wValue = mu_;
2697 }
2698#endif
2699 assert (rhsL_[iColumn] >= 0.0);
2700 assert (rhsU_[iColumn] <= 0.0);
2701 if (lowerBound(iColumn)) {
2702 value += (-zVec_[iColumn] * rhsL_[iColumn] - zValue) /
2703 (lowerSlack_[iColumn] + extra);
2704 }
2705 if (upperBound(iColumn)) {
2706 value += (wValue - wVec_[iColumn] * rhsU_[iColumn]) /
2707 (upperSlack_[iColumn] + extra);
2708 }
2709#else
2710 if (lowerBound(iColumn)) {
2711 CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
2712 value -= gHat / (lowerSlack_[iColumn] + extra);
2713 }
2714 if (upperBound(iColumn)) {
2715 CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
2716 value += hHat / (upperSlack_[iColumn] + extra);
2717 }
2718#endif
2719 workArray_[iColumn] = diagonal_[iColumn] * value;
2720 }
2721#if 0
2722 if (solution_[0] > 0.0) {
2723 for (int i = 0; i < numberTotal; i++)
2724 printf("%d %.18g\n", i, workArray_[i]);
2725 } else {
2726 for (int i = 0; i < numberTotal; i++)
2727 printf("%d %.18g\n", i, workArray_[i]);
2728 }
2729 exit(66);
2730#endif
2731 } else {
2732 // KKT
2733 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
2734 CoinWorkDouble value = rhsC_[iColumn];
2735 CoinWorkDouble zValue = rhsZ_[iColumn];
2736 CoinWorkDouble wValue = rhsW_[iColumn];
2737 if (lowerBound(iColumn)) {
2738 CoinWorkDouble gHat = zValue + zVec_[iColumn] * rhsL_[iColumn];
2739 value -= gHat / (lowerSlack_[iColumn] + extra);
2740 }
2741 if (upperBound(iColumn)) {
2742 CoinWorkDouble hHat = wValue - wVec_[iColumn] * rhsU_[iColumn];
2743 value += hHat / (upperSlack_[iColumn] + extra);
2744 }
2745 workArray_[iColumn] = value;
2746 }
2747 }
2748}
2749//method: sees if looks plausible change in complementarity
2750bool ClpPredictorCorrector::checkGoodMove(const bool doCorrector,
2751 CoinWorkDouble & bestNextGap,
2752 bool allowIncreasingGap)
2753{
2754 const CoinWorkDouble beta3 = 0.99997;
2755 bool goodMove = false;
2756 int nextNumber;
2757 int nextNumberItems;
2758 int numberTotal = numberRows_ + numberColumns_;
2759 CoinWorkDouble returnGap = bestNextGap;
2760 CoinWorkDouble nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
2761#ifndef NO_RTTI
2762 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
2763#else
2764 ClpQuadraticObjective * quadraticObj = NULL;
2765 if (objective_->type() == 2)
2766 quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
2767#endif
2768 if (nextGap > bestNextGap && nextGap > 0.9 * complementarityGap_ && doCorrector
2769 && !quadraticObj && !allowIncreasingGap) {
2770#ifdef SOME_DEBUG
2771 printf("checkGood phase 1 next gap %.18g, phase 0 %.18g, old gap %.18g\n",
2772 nextGap, bestNextGap, complementarityGap_);
2773#endif
2774 return false;
2775 } else {
2776 returnGap = nextGap;
2777 }
2778 CoinWorkDouble step;
2779 if (actualDualStep_ > actualPrimalStep_) {
2780 step = actualDualStep_;
2781 } else {
2782 step = actualPrimalStep_;
2783 }
2784 CoinWorkDouble testValue = 1.0 - step * (1.0 - beta3);
2785 //testValue=0.0;
2786 testValue *= complementarityGap_;
2787 if (nextGap < testValue) {
2788 //std::cout <<"predicted duality gap "<<nextGap<<std::endl;
2789 goodMove = true;
2790 } else if(doCorrector) {
2791 //if (actualDualStep_<actualPrimalStep_) {
2792 //step=actualDualStep_;
2793 //} else {
2794 //step=actualPrimalStep_;
2795 //}
2796 CoinWorkDouble gap = bestNextGap;
2797 goodMove = checkGoodMove2(step, gap, allowIncreasingGap);
2798 if (goodMove)
2799 returnGap = gap;
2800 } else {
2801 goodMove = true;
2802 }
2803 if (goodMove)
2804 goodMove = checkGoodMove2(step, bestNextGap, allowIncreasingGap);
2805 // Say good if small
2806 //if (quadraticObj) {
2807 if (CoinMax(actualDualStep_, actualPrimalStep_) < 1.0e-6)
2808 goodMove = true;
2809 if (!goodMove) {
2810 //try smaller of two
2811 if (actualDualStep_ < actualPrimalStep_) {
2812 step = actualDualStep_;
2813 } else {
2814 step = actualPrimalStep_;
2815 }
2816 if (step > 1.0) {
2817 step = 1.0;
2818 }
2819 actualPrimalStep_ = step;
2820 //if (quadraticObj)
2821 //actualPrimalStep_ *=0.5;
2822 actualDualStep_ = step;
2823 goodMove = checkGoodMove2(step, bestNextGap, allowIncreasingGap);
2824 int pass = 0;
2825 while (!goodMove) {
2826 pass++;
2827 CoinWorkDouble gap = bestNextGap;
2828 goodMove = checkGoodMove2(step, gap, allowIncreasingGap);
2829 if (goodMove || pass > 3) {
2830 returnGap = gap;
2831 break;
2832 }
2833 if (step < 1.0e-4) {
2834 break;
2835 }
2836 step *= 0.5;
2837 actualPrimalStep_ = step;
2838 //if (quadraticObj)
2839 //actualPrimalStep_ *=0.5;
2840 actualDualStep_ = step;
2841 } /* endwhile */
2842 if (doCorrector) {
2843 //say bad move if both small
2844 if (numberIterations_ & 1) {
2845 if (actualPrimalStep_ < 1.0e-2 && actualDualStep_ < 1.0e-2) {
2846 goodMove = false;
2847 }
2848 } else {
2849 if (actualPrimalStep_ < 1.0e-5 && actualDualStep_ < 1.0e-5) {
2850 goodMove = false;
2851 }
2852 if (actualPrimalStep_ * actualDualStep_ < 1.0e-20) {
2853 goodMove = false;
2854 }
2855 }
2856 }
2857 }
2858 if (goodMove) {
2859 //compute delta in objectives
2860 CoinWorkDouble deltaObjectivePrimal = 0.0;
2861 CoinWorkDouble deltaObjectiveDual =
2862 innerProduct(deltaY_, numberRows_,
2863 rhsFixRegion_);
2864 CoinWorkDouble error = 0.0;
2865 CoinWorkDouble * workArray = workArray_;
2866 CoinZeroN(workArray, numberColumns_);
2867 CoinMemcpyN(deltaY_, numberRows_, workArray + numberColumns_);
2868 matrix_->transposeTimes(-1.0, deltaY_, workArray);
2869 //CoinWorkDouble sumPerturbCost=0.0;
2870 for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
2871 if (!flagged(iColumn)) {
2872 if (lowerBound(iColumn)) {
2873 //sumPerturbCost+=deltaX_[iColumn];
2874 deltaObjectiveDual += deltaZ_[iColumn] * lower_[iColumn];
2875 }
2876 if (upperBound(iColumn)) {
2877 //sumPerturbCost-=deltaX_[iColumn];
2878 deltaObjectiveDual -= deltaW_[iColumn] * upper_[iColumn];
2879 }
2880 CoinWorkDouble change = CoinAbs(workArray_[iColumn] - deltaZ_[iColumn] + deltaW_[iColumn]);
2881 error = CoinMax(change, error);
2882 }
2883 deltaObjectivePrimal += cost_[iColumn] * deltaX_[iColumn];
2884 }
2885 //deltaObjectivePrimal+=sumPerturbCost*linearPerturbation_;
2886 CoinWorkDouble testValue;
2887 if (error > 0.0) {
2888 testValue = 1.0e1 * CoinMax(maximumDualError_, 1.0e-12) / error;
2889 } else {
2890 testValue = 1.0e1;
2891 }
2892 // If quadratic then primal step may compensate
2893 if (testValue < actualDualStep_ && !quadraticObj) {
2894 handler_->message(CLP_BARRIER_REDUCING, messages_)
2895 << "dual" << static_cast<double>(actualDualStep_)
2896 << static_cast<double>(testValue)
2897 << CoinMessageEol;
2898 actualDualStep_ = testValue;
2899 }
2900 }
2901 if (maximumRHSError_ < 1.0e1 * solutionNorm_ * primalTolerance()
2902 && maximumRHSChange_ > 1.0e-16 * solutionNorm_) {
2903 //check change in AX not too much
2904 //??? could be dropped row going infeasible
2905 CoinWorkDouble ratio = 1.0e1 * CoinMax(maximumRHSError_, 1.0e-12) / maximumRHSChange_;
2906 if (ratio < actualPrimalStep_) {
2907 handler_->message(CLP_BARRIER_REDUCING, messages_)
2908 << "primal" << static_cast<double>(actualPrimalStep_)
2909 << static_cast<double>(ratio)
2910 << CoinMessageEol;
2911 if (ratio > 1.0e-6) {
2912 actualPrimalStep_ = ratio;
2913 } else {
2914 actualPrimalStep_ = ratio;
2915 //std::cout <<"sign we should be stopping"<<std::endl;
2916 }
2917 }
2918 }
2919 if (goodMove)
2920 bestNextGap = returnGap;
2921 return goodMove;
2922}
2923//: checks for one step size
2924bool ClpPredictorCorrector::checkGoodMove2(CoinWorkDouble move,
2925 CoinWorkDouble & bestNextGap,
2926 bool allowIncreasingGap)
2927{
2928 CoinWorkDouble complementarityMultiplier = 1.0 / numberComplementarityPairs_;
2929 const CoinWorkDouble gamma = 1.0e-8;
2930 const CoinWorkDouble gammap = 1.0e-8;
2931 CoinWorkDouble gammad = 1.0e-8;
2932 int nextNumber;
2933 int nextNumberItems;
2934 CoinWorkDouble nextGap = complementarityGap(nextNumber, nextNumberItems, 2);
2935 if (nextGap > bestNextGap && !allowIncreasingGap)
2936 return false;
2937 CoinWorkDouble lowerBoundGap = gamma * nextGap * complementarityMultiplier;
2938 bool goodMove = true;
2939 int iColumn;
2940 for ( iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
2941 if (!flagged(iColumn)) {
2942 if (lowerBound(iColumn)) {
2943 CoinWorkDouble part1 = lowerSlack_[iColumn] + actualPrimalStep_ * deltaSL_[iColumn];
2944 CoinWorkDouble part2 = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
2945 if (part1 * part2 < lowerBoundGap) {
2946 goodMove = false;
2947 break;
2948 }
2949 }
2950 if (upperBound(iColumn)) {
2951 CoinWorkDouble part1 = upperSlack_[iColumn] + actualPrimalStep_ * deltaSU_[iColumn];
2952 CoinWorkDouble part2 = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
2953 if (part1 * part2 < lowerBoundGap) {
2954 goodMove = false;
2955 break;
2956 }
2957 }
2958 }
2959 }
2960 CoinWorkDouble * nextDj = NULL;
2961 CoinWorkDouble maximumDualError = maximumDualError_;
2962#ifndef NO_RTTI
2963 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
2964#else
2965 ClpQuadraticObjective * quadraticObj = NULL;
2966 if (objective_->type() == 2)
2967 quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
2968#endif
2969 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
2970 if (quadraticObj) {
2971 // change gammad
2972 gammad = 1.0e-4;
2973 CoinWorkDouble gamma2 = gamma_ * gamma_;
2974 nextDj = new CoinWorkDouble [numberColumns_];
2975 CoinWorkDouble * nextSolution = new CoinWorkDouble [numberColumns_];
2976 // put next primal into nextSolution
2977 for ( iColumn = 0; iColumn < numberColumns_; iColumn++) {
2978 if (!flagged(iColumn)) {
2979 nextSolution[iColumn] = solution_[iColumn] +
2980 actualPrimalStep_ * deltaX_[iColumn];
2981 } else {
2982 nextSolution[iColumn] = solution_[iColumn];
2983 }
2984 }
2985 // do reduced costs
2986 CoinMemcpyN(cost_, numberColumns_, nextDj);
2987 matrix_->transposeTimes(-1.0, dualArray, nextDj);
2988 matrix_->transposeTimes(-actualDualStep_, deltaY_, nextDj);
2989 quadraticDjs(nextDj, nextSolution, 1.0);
2990 delete [] nextSolution;
2991 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2992 const int * columnQuadraticLength = quadratic->getVectorLengths();
2993 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2994 if (!fixedOrFree(iColumn)) {
2995 CoinWorkDouble newZ = 0.0;
2996 CoinWorkDouble newW = 0.0;
2997 if (lowerBound(iColumn)) {
2998 newZ = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
2999 }
3000 if (upperBound(iColumn)) {
3001 newW = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
3002 }
3003 if (columnQuadraticLength[iColumn]) {
3004 CoinWorkDouble gammaTerm = gamma2;
3005 if (primalR_)
3006 gammaTerm += primalR_[iColumn];
3007 //CoinWorkDouble dualInfeasibility=
3008 //dj_[iColumn]-zVec_[iColumn]+wVec_[iColumn]
3009 //+gammaTerm*solution_[iColumn];
3010 CoinWorkDouble newInfeasibility =
3011 nextDj[iColumn] - newZ + newW
3012 + gammaTerm * (solution_[iColumn] + actualPrimalStep_ * deltaX_[iColumn]);
3013 maximumDualError = CoinMax(maximumDualError, newInfeasibility);
3014 //if (CoinAbs(newInfeasibility)>CoinMax(2000.0*maximumDualError_,1.0e-2)) {
3015 //if (dualInfeasibility*newInfeasibility<0.0) {
3016 // printf("%d current %g next %g\n",iColumn,dualInfeasibility,
3017 // newInfeasibility);
3018 // goodMove=false;
3019 //}
3020 //}
3021 }
3022 }
3023 }
3024 delete [] nextDj;
3025 }
3026// Satisfy g_p(alpha)?
3027 if (rhsNorm_ > solutionNorm_) {
3028 solutionNorm_ = rhsNorm_;
3029 }
3030 CoinWorkDouble errorCheck = maximumRHSError_ / solutionNorm_;
3031 if (errorCheck < maximumBoundInfeasibility_) {
3032 errorCheck = maximumBoundInfeasibility_;
3033 }
3034 // scale back move
3035 move = CoinMin(move, 0.95);
3036 //scale
3037 if ((1.0 - move)*errorCheck > primalTolerance()) {
3038 if (nextGap < gammap*(1.0 - move)*errorCheck) {
3039 goodMove = false;
3040 }
3041 }
3042 // Satisfy g_d(alpha)?
3043 errorCheck = maximumDualError / objectiveNorm_;
3044 if ((1.0 - move)*errorCheck > dualTolerance()) {
3045 if (nextGap < gammad*(1.0 - move)*errorCheck) {
3046 goodMove = false;
3047 }
3048 }
3049 if (goodMove)
3050 bestNextGap = nextGap;
3051 return goodMove;
3052}
3053// updateSolution. Updates solution at end of iteration
3054//returns number fixed
3055int ClpPredictorCorrector::updateSolution(CoinWorkDouble /*nextGap*/)
3056{
3057 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
3058 int numberTotal = numberRows_ + numberColumns_;
3059 //update pi
3060 multiplyAdd(deltaY_, numberRows_, actualDualStep_, dualArray, 1.0);
3061 CoinZeroN(errorRegion_, numberRows_);
3062 CoinZeroN(rhsFixRegion_, numberRows_);
3063 CoinWorkDouble maximumRhsInfeasibility = 0.0;
3064 CoinWorkDouble maximumBoundInfeasibility = 0.0;
3065 CoinWorkDouble maximumDualError = 1.0e-12;
3066 CoinWorkDouble primalObjectiveValue = 0.0;
3067 CoinWorkDouble dualObjectiveValue = 0.0;
3068 CoinWorkDouble solutionNorm = 1.0e-12;
3069 int numberKilled = 0;
3070 CoinWorkDouble freeMultiplier = 1.0e6;
3071 CoinWorkDouble trueNorm = diagonalNorm_ / diagonalScaleFactor_;
3072 if (freeMultiplier < trueNorm) {
3073 freeMultiplier = trueNorm;
3074 }
3075 if (freeMultiplier > 1.0e12) {
3076 freeMultiplier = 1.0e12;
3077 }
3078 freeMultiplier = 0.5 / freeMultiplier;
3079 CoinWorkDouble condition = CoinAbs(cholesky_->choleskyCondition());
3080 bool caution;
3081 if ((condition < 1.0e10 && trueNorm < 1.0e12) || numberIterations_ < 20) {
3082 caution = false;
3083 } else {
3084 caution = true;
3085 }
3086 CoinWorkDouble extra = eExtra;
3087 const CoinWorkDouble largeFactor = 1.0e2;
3088 CoinWorkDouble largeGap = largeFactor * solutionNorm_;
3089 if (largeGap < largeFactor) {
3090 largeGap = largeFactor;
3091 }
3092 CoinWorkDouble dualFake = 0.0;
3093 CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance];
3094 dualTolerance = dualTolerance / scaleFactor_;
3095 if (dualTolerance < 1.0e-12) {
3096 dualTolerance = 1.0e-12;
3097 }
3098 CoinWorkDouble offsetObjective = 0.0;
3099 CoinWorkDouble killTolerance = primalTolerance();
3100 CoinWorkDouble qDiagonal;
3101 if (mu_ < 1.0) {
3102 qDiagonal = 1.0e-8;
3103 } else {
3104 qDiagonal = 1.0e-8 * mu_;
3105 }
3106 //CoinWorkDouble nextMu = nextGap/(static_cast<CoinWorkDouble>(2*numberComplementarityPairs_));
3107 //printf("using gap of %g\n",nextMu);
3108 //qDiagonal *= 1.0e2;
3109 //largest allowable ratio of lowerSlack/zVec (etc)
3110 CoinWorkDouble largestRatio;
3111 CoinWorkDouble epsilonBase;
3112 CoinWorkDouble diagonalLimit;
3113 if (!caution) {
3114 epsilonBase = eBase;
3115 largestRatio = eRatio;
3116 diagonalLimit = eDiagonal;
3117 } else {
3118 epsilonBase = eBaseCaution;
3119 largestRatio = eRatioCaution;
3120 diagonalLimit = eDiagonalCaution;
3121 }
3122 CoinWorkDouble smallGap = 1.0e2;
3123 CoinWorkDouble maximumDJInfeasibility = 0.0;
3124 int numberIncreased = 0;
3125 int numberDecreased = 0;
3126 CoinWorkDouble largestDiagonal = 0.0;
3127 CoinWorkDouble smallestDiagonal = 1.0e50;
3128 CoinWorkDouble largeGap2 = CoinMax(1.0e7, 1.0e2 * solutionNorm_);
3129 //largeGap2 = 1.0e9;
3130 // When to start looking at killing (factor0
3131 CoinWorkDouble killFactor;
3132#ifndef NO_RTTI
3133 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
3134#else
3135 ClpQuadraticObjective * quadraticObj = NULL;
3136 if (objective_->type() == 2)
3137 quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
3138#endif
3139#ifndef CLP_CAUTION
3140#define KILL_ITERATION 50
3141#else
3142#if CLP_CAUTION < 1
3143#define KILL_ITERATION 50
3144#else
3145#define KILL_ITERATION 100
3146#endif
3147#endif
3148 if (!quadraticObj || 1) {
3149 if (numberIterations_ < KILL_ITERATION) {
3150 killFactor = 1.0;
3151 } else if (numberIterations_ < 2 * KILL_ITERATION) {
3152 killFactor = 5.0;
3153 stepLength_ = CoinMax(stepLength_, 0.9995);
3154 } else if (numberIterations_ < 4 * KILL_ITERATION) {
3155 killFactor = 20.0;
3156 stepLength_ = CoinMax(stepLength_, 0.99995);
3157 } else {
3158 killFactor = 1.0e2;
3159 stepLength_ = CoinMax(stepLength_, 0.999995);
3160 }
3161 } else {
3162 killFactor = 1.0;
3163 }
3164 // put next primal into deltaSL_
3165 int iColumn;
3166 int iRow;
3167 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
3168 CoinWorkDouble thisWeight = deltaX_[iColumn];
3169 CoinWorkDouble newPrimal = solution_[iColumn] + 1.0 * actualPrimalStep_ * thisWeight;
3170 deltaSL_[iColumn] = newPrimal;
3171 }
3172#if 0
3173 // nice idea but doesn't work
3174 multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, errorRegion_, 0.0);
3175 matrix_->times(1.0, solution_, errorRegion_);
3176 multiplyAdd(deltaSL_ + numberColumns_, numberRows_, -1.0, rhsFixRegion_, 0.0);
3177 matrix_->times(1.0, deltaSL_, rhsFixRegion_);
3178 CoinWorkDouble newNorm = maximumAbsElement(deltaSL_, numberTotal);
3179 CoinWorkDouble tol = newNorm * primalTolerance();
3180 bool goneInf = false;
3181 for (iRow = 0; iRow < numberRows_; iRow++) {
3182 CoinWorkDouble value = errorRegion_[iRow];
3183 CoinWorkDouble valueNew = rhsFixRegion_[iRow];
3184 if (CoinAbs(value) < tol && CoinAbs(valueNew) > tol) {
3185 printf("row %d old %g new %g\n", iRow, value, valueNew);
3186 goneInf = true;
3187 }
3188 }
3189 if (goneInf) {
3190 actualPrimalStep_ *= 0.5;
3191 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
3192 CoinWorkDouble thisWeight = deltaX_[iColumn];
3193 CoinWorkDouble newPrimal = solution_[iColumn] + 1.0 * actualPrimalStep_ * thisWeight;
3194 deltaSL_[iColumn] = newPrimal;
3195 }
3196 }
3197 CoinZeroN(errorRegion_, numberRows_);
3198 CoinZeroN(rhsFixRegion_, numberRows_);
3199#endif
3200 // do reduced costs
3201 CoinMemcpyN(dualArray, numberRows_, dj_ + numberColumns_);
3202 CoinMemcpyN(cost_, numberColumns_, dj_);
3203 CoinWorkDouble quadraticOffset = quadraticDjs(dj_, deltaSL_, 1.0);
3204 // Save modified costs for fixed variables
3205 CoinMemcpyN(dj_, numberColumns_, deltaSU_);
3206 matrix_->transposeTimes(-1.0, dualArray, dj_);
3207 CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
3208 CoinWorkDouble gammaOffset = 0.0;
3209#if 0
3210 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3211 const int * columnLength = matrix_->getVectorLengths();
3212 const int * row = matrix_->getIndices();
3213 const double * element = matrix_->getElements();
3214#endif
3215 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
3216 if (!flagged(iColumn)) {
3217 CoinWorkDouble reducedCost = dj_[iColumn];
3218 bool thisKilled = false;
3219 CoinWorkDouble zValue = zVec_[iColumn] + actualDualStep_ * deltaZ_[iColumn];
3220 CoinWorkDouble wValue = wVec_[iColumn] + actualDualStep_ * deltaW_[iColumn];
3221 zVec_[iColumn] = zValue;
3222 wVec_[iColumn] = wValue;
3223 CoinWorkDouble thisWeight = deltaX_[iColumn];
3224 CoinWorkDouble oldPrimal = solution_[iColumn];
3225 CoinWorkDouble newPrimal = solution_[iColumn] + actualPrimalStep_ * thisWeight;
3226 CoinWorkDouble dualObjectiveThis = 0.0;
3227 CoinWorkDouble sUpper = extra;
3228 CoinWorkDouble sLower = extra;
3229 CoinWorkDouble kill;
3230 if (CoinAbs(newPrimal) > 1.0e4) {
3231 kill = killTolerance * 1.0e-4 * newPrimal;
3232 } else {
3233 kill = killTolerance;
3234 }
3235 kill *= 1.0e-3; //be conservative
3236 CoinWorkDouble smallerSlack = COIN_DBL_MAX;
3237 bool fakeOldBounds = false;
3238 bool fakeNewBounds = false;
3239 CoinWorkDouble trueLower;
3240 CoinWorkDouble trueUpper;
3241 if (iColumn < numberColumns_) {
3242 trueLower = columnLower_[iColumn];
3243 trueUpper = columnUpper_[iColumn];
3244 } else {
3245 trueLower = rowLower_[iColumn-numberColumns_];
3246 trueUpper = rowUpper_[iColumn-numberColumns_];
3247 }
3248 if (oldPrimal > trueLower + largeGap2 &&
3249 oldPrimal < trueUpper - largeGap2)
3250 fakeOldBounds = true;
3251 if (newPrimal > trueLower + largeGap2 &&
3252 newPrimal < trueUpper - largeGap2)
3253 fakeNewBounds = true;
3254 if (fakeOldBounds) {
3255 if (fakeNewBounds) {
3256 lower_[iColumn] = newPrimal - largeGap2;
3257 lowerSlack_[iColumn] = largeGap2;
3258 upper_[iColumn] = newPrimal + largeGap2;
3259 upperSlack_[iColumn] = largeGap2;
3260 } else {
3261 lower_[iColumn] = trueLower;
3262 setLowerBound(iColumn);
3263 lowerSlack_[iColumn] = CoinMax(newPrimal - trueLower, 1.0);
3264 upper_[iColumn] = trueUpper;
3265 setUpperBound(iColumn);
3266 upperSlack_[iColumn] = CoinMax(trueUpper - newPrimal, 1.0);
3267 }
3268 } else if (fakeNewBounds) {
3269 lower_[iColumn] = newPrimal - largeGap2;
3270 lowerSlack_[iColumn] = largeGap2;
3271 upper_[iColumn] = newPrimal + largeGap2;
3272 upperSlack_[iColumn] = largeGap2;
3273 // so we can just have one test
3274 fakeOldBounds = true;
3275 }
3276 CoinWorkDouble lowerBoundInfeasibility = 0.0;
3277 CoinWorkDouble upperBoundInfeasibility = 0.0;
3278 //double saveNewPrimal = newPrimal;
3279 if (lowerBound(iColumn)) {
3280 CoinWorkDouble oldSlack = lowerSlack_[iColumn];
3281 CoinWorkDouble newSlack;
3282 newSlack =
3283 lowerSlack_[iColumn] + actualPrimalStep_ * (oldPrimal - oldSlack
3284 + thisWeight - lower_[iColumn]);
3285 if (fakeOldBounds)
3286 newSlack = lowerSlack_[iColumn];
3287 CoinWorkDouble epsilon = CoinAbs(newSlack) * epsilonBase;
3288 epsilon = CoinMin(epsilon, 1.0e-5);
3289 //epsilon=1.0e-14;
3290 //make sure reasonable
3291 if (zValue < epsilon) {
3292 zValue = epsilon;
3293 }
3294 CoinWorkDouble feasibleSlack = newPrimal - lower_[iColumn];
3295 if (feasibleSlack > 0.0 && newSlack > 0.0) {
3296 CoinWorkDouble smallGap2 = smallGap;
3297 if (CoinAbs(0.1 * newPrimal) > smallGap) {
3298 smallGap2 = 0.1 * CoinAbs(newPrimal);
3299 }
3300 CoinWorkDouble larger;
3301 if (newSlack > feasibleSlack) {
3302 larger = newSlack;
3303 } else {
3304 larger = feasibleSlack;
3305 }
3306 if (CoinAbs(feasibleSlack - newSlack) < 1.0e-6 * larger) {
3307 newSlack = feasibleSlack;
3308 }
3309 }
3310 if (zVec_[iColumn] > dualTolerance) {
3311 dualObjectiveThis += lower_[iColumn] * zVec_[iColumn];
3312 }
3313 lowerSlack_[iColumn] = newSlack;
3314 if (newSlack < smallerSlack) {
3315 smallerSlack = newSlack;
3316 }
3317 lowerBoundInfeasibility = CoinAbs(newPrimal - lowerSlack_[iColumn] - lower_[iColumn]);
3318 if (lowerSlack_[iColumn] <= kill * killFactor && CoinAbs(newPrimal - lower_[iColumn]) <= kill * killFactor) {
3319 CoinWorkDouble step = CoinMin(actualPrimalStep_ * 1.1, 1.0);
3320 CoinWorkDouble newPrimal2 = solution_[iColumn] + step * thisWeight;
3321 if (newPrimal2 < newPrimal && dj_[iColumn] > 1.0e-5 && numberIterations_ > 50 - 40) {
3322 newPrimal = lower_[iColumn];
3323 lowerSlack_[iColumn] = 0.0;
3324 //printf("fixing %d to lower\n",iColumn);
3325 }
3326 }
3327 if (lowerSlack_[iColumn] <= kill && CoinAbs(newPrimal - lower_[iColumn]) <= kill) {
3328 //may be better to leave at value?
3329 newPrimal = lower_[iColumn];
3330 lowerSlack_[iColumn] = 0.0;
3331 thisKilled = true;
3332 //cout<<j<<" l "<<reducedCost<<" "<<zVec_[iColumn]<<endl;
3333 } else {
3334 sLower += lowerSlack_[iColumn];
3335 }
3336 }
3337 if (upperBound(iColumn)) {
3338 CoinWorkDouble oldSlack = upperSlack_[iColumn];
3339 CoinWorkDouble newSlack;
3340 newSlack =
3341 upperSlack_[iColumn] + actualPrimalStep_ * (-oldPrimal - oldSlack
3342 - thisWeight + upper_[iColumn]);
3343 if (fakeOldBounds)
3344 newSlack = upperSlack_[iColumn];
3345 CoinWorkDouble epsilon = CoinAbs(newSlack) * epsilonBase;
3346 epsilon = CoinMin(epsilon, 1.0e-5);
3347 //make sure reasonable
3348 //epsilon=1.0e-14;
3349 if (wValue < epsilon) {
3350 wValue = epsilon;
3351 }
3352 CoinWorkDouble feasibleSlack = upper_[iColumn] - newPrimal;
3353 if (feasibleSlack > 0.0 && newSlack > 0.0) {
3354 CoinWorkDouble smallGap2 = smallGap;
3355 if (CoinAbs(0.1 * newPrimal) > smallGap) {
3356 smallGap2 = 0.1 * CoinAbs(newPrimal);
3357 }
3358 CoinWorkDouble larger;
3359 if (newSlack > feasibleSlack) {
3360 larger = newSlack;
3361 } else {
3362 larger = feasibleSlack;
3363 }
3364 if (CoinAbs(feasibleSlack - newSlack) < 1.0e-6 * larger) {
3365 newSlack = feasibleSlack;
3366 }
3367 }
3368 if (wVec_[iColumn] > dualTolerance) {
3369 dualObjectiveThis -= upper_[iColumn] * wVec_[iColumn];
3370 }
3371 upperSlack_[iColumn] = newSlack;
3372 if (newSlack < smallerSlack) {
3373 smallerSlack = newSlack;
3374 }
3375 upperBoundInfeasibility = CoinAbs(newPrimal + upperSlack_[iColumn] - upper_[iColumn]);
3376 if (upperSlack_[iColumn] <= kill * killFactor && CoinAbs(newPrimal - upper_[iColumn]) <= kill * killFactor) {
3377 CoinWorkDouble step = CoinMin(actualPrimalStep_ * 1.1, 1.0);
3378 CoinWorkDouble newPrimal2 = solution_[iColumn] + step * thisWeight;
3379 if (newPrimal2 > newPrimal && dj_[iColumn] < -1.0e-5 && numberIterations_ > 50 - 40) {
3380 newPrimal = upper_[iColumn];
3381 upperSlack_[iColumn] = 0.0;
3382 //printf("fixing %d to upper\n",iColumn);
3383 }
3384 }
3385 if (upperSlack_[iColumn] <= kill && CoinAbs(newPrimal - upper_[iColumn]) <= kill) {
3386 //may be better to leave at value?
3387 newPrimal = upper_[iColumn];
3388 upperSlack_[iColumn] = 0.0;
3389 thisKilled = true;
3390 } else {
3391 sUpper += upperSlack_[iColumn];
3392 }
3393 }
3394#if 0
3395 if (newPrimal != saveNewPrimal && iColumn < numberColumns_) {
3396 // adjust slacks
3397 double movement = newPrimal - saveNewPrimal;
3398 for (CoinBigIndex j = columnStart[iColumn];
3399 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3400 int iRow = row[j];
3401 double slackMovement = element[j] * movement;
3402 solution_[iRow+numberColumns_] += slackMovement; // sign?
3403 }
3404 }
3405#endif
3406 solution_[iColumn] = newPrimal;
3407 if (CoinAbs(newPrimal) > solutionNorm) {
3408 solutionNorm = CoinAbs(newPrimal);
3409 }
3410 if (!thisKilled) {
3411 CoinWorkDouble gammaTerm = gamma2;
3412 if (primalR_) {
3413 gammaTerm += primalR_[iColumn];
3414 quadraticOffset += newPrimal * newPrimal * primalR_[iColumn];
3415 }
3416 CoinWorkDouble dualInfeasibility =
3417 reducedCost - zVec_[iColumn] + wVec_[iColumn] + gammaTerm * newPrimal;
3418 if (CoinAbs(dualInfeasibility) > dualTolerance) {
3419#if 0
3420 if (dualInfeasibility > 0.0) {
3421 // To improve we could reduce t and/or increase z
3422 if (lowerBound(iColumn)) {
3423 CoinWorkDouble complementarity = zVec_[iColumn] * lowerSlack_[iColumn];
3424 if (complementarity < nextMu) {
3425 CoinWorkDouble change =
3426 CoinMin(dualInfeasibility,
3427 (nextMu - complementarity) / lowerSlack_[iColumn]);
3428 dualInfeasibility -= change;
3429 COIN_DETAIL_PRINT(printf("%d lb locomp %g - dual inf from %g to %g\n",
3430 iColumn, complementarity, dualInfeasibility + change,
3431 dualInfeasibility));
3432 zVec_[iColumn] += change;
3433 zValue = CoinMax(zVec_[iColumn], 1.0e-12);
3434 }
3435 }
3436 if (upperBound(iColumn)) {
3437 CoinWorkDouble complementarity = wVec_[iColumn] * upperSlack_[iColumn];
3438 if (complementarity > nextMu) {
3439 CoinWorkDouble change =
3440 CoinMin(dualInfeasibility,
3441 (complementarity - nextMu) / upperSlack_[iColumn]);
3442 dualInfeasibility -= change;
3443 COIN_DETAIL_PRINT(printf("%d ub hicomp %g - dual inf from %g to %g\n",
3444 iColumn, complementarity, dualInfeasibility + change,
3445 dualInfeasibility));
3446 wVec_[iColumn] -= change;
3447 wValue = CoinMax(wVec_[iColumn], 1.0e-12);
3448 }
3449 }
3450 } else {
3451 // To improve we could reduce z and/or increase t
3452 if (lowerBound(iColumn)) {
3453 CoinWorkDouble complementarity = zVec_[iColumn] * lowerSlack_[iColumn];
3454 if (complementarity > nextMu) {
3455 CoinWorkDouble change =
3456 CoinMax(dualInfeasibility,
3457 (nextMu - complementarity) / lowerSlack_[iColumn]);
3458 dualInfeasibility -= change;
3459 COIN_DETAIL_PRINT(printf("%d lb hicomp %g - dual inf from %g to %g\n",
3460 iColumn, complementarity, dualInfeasibility + change,
3461 dualInfeasibility));
3462 zVec_[iColumn] += change;
3463 zValue = CoinMax(zVec_[iColumn], 1.0e-12);
3464 }
3465 }
3466 if (upperBound(iColumn)) {
3467 CoinWorkDouble complementarity = wVec_[iColumn] * upperSlack_[iColumn];
3468 if (complementarity < nextMu) {
3469 CoinWorkDouble change =
3470 CoinMax(dualInfeasibility,
3471 (complementarity - nextMu) / upperSlack_[iColumn]);
3472 dualInfeasibility -= change;
3473 COIN_DETAIL_PRINT(printf("%d ub locomp %g - dual inf from %g to %g\n",
3474 iColumn, complementarity, dualInfeasibility + change,
3475 dualInfeasibility));
3476 wVec_[iColumn] -= change;
3477 wValue = CoinMax(wVec_[iColumn], 1.0e-12);
3478 }
3479 }
3480 }
3481#endif
3482 dualFake += newPrimal * dualInfeasibility;
3483 }
3484 if (lowerBoundInfeasibility > maximumBoundInfeasibility) {
3485 maximumBoundInfeasibility = lowerBoundInfeasibility;
3486 }
3487 if (upperBoundInfeasibility > maximumBoundInfeasibility) {
3488 maximumBoundInfeasibility = upperBoundInfeasibility;
3489 }
3490 dualInfeasibility = CoinAbs(dualInfeasibility);
3491 if (dualInfeasibility > maximumDualError) {
3492 //printf("bad dual %d %g\n",iColumn,
3493 // reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal);
3494 maximumDualError = dualInfeasibility;
3495 }
3496 dualObjectiveValue += dualObjectiveThis;
3497 gammaOffset += newPrimal * newPrimal;
3498 if (sLower > largeGap) {
3499 sLower = largeGap;
3500 }
3501 if (sUpper > largeGap) {
3502 sUpper = largeGap;
3503 }
3504#if 1
3505 CoinWorkDouble divisor = sLower * wValue + sUpper * zValue + gammaTerm * sLower * sUpper;
3506 CoinWorkDouble diagonalValue = (sUpper * sLower) / divisor;
3507#else
3508 CoinWorkDouble divisor = sLower * wValue + sUpper * zValue + gammaTerm * sLower * sUpper;
3509 CoinWorkDouble diagonalValue2 = (sUpper * sLower) / divisor;
3510 CoinWorkDouble diagonalValue;
3511 if (!lowerBound(iColumn)) {
3512 diagonalValue = wValue / sUpper + gammaTerm;
3513 } else if (!upperBound(iColumn)) {
3514 diagonalValue = zValue / sLower + gammaTerm;
3515 } else {
3516 diagonalValue = zValue / sLower + wValue / sUpper + gammaTerm;
3517 }
3518 diagonalValue = 1.0 / diagonalValue;
3519#endif
3520 diagonal_[iColumn] = diagonalValue;
3521 //FUDGE
3522 if (diagonalValue > diagonalLimit) {
3523#ifdef COIN_DEVELOP
3524 std::cout << "large diagonal " << diagonalValue << std::endl;
3525#endif
3526 diagonal_[iColumn] = diagonalLimit;
3527 }
3528#ifdef COIN_DEVELOP
3529 if (diagonalValue < 1.0e-10) {
3530 //std::cout<<"small diagonal "<<diagonalValue<<std::endl;
3531 }
3532#endif
3533 if (diagonalValue > largestDiagonal) {
3534 largestDiagonal = diagonalValue;
3535 }
3536 if (diagonalValue < smallestDiagonal) {
3537 smallestDiagonal = diagonalValue;
3538 }
3539 deltaX_[iColumn] = 0.0;
3540 } else {
3541 numberKilled++;
3542 if (solution_[iColumn] != lower_[iColumn] &&
3543 solution_[iColumn] != upper_[iColumn]) {
3544 COIN_DETAIL_PRINT(printf("%d %g %g %g\n", iColumn, static_cast<double>(lower_[iColumn]),
3545 static_cast<double>(solution_[iColumn]), static_cast<double>(upper_[iColumn])));
3546 }
3547 diagonal_[iColumn] = 0.0;
3548 zVec_[iColumn] = 0.0;
3549 wVec_[iColumn] = 0.0;
3550 setFlagged(iColumn);
3551 setFixedOrFree(iColumn);
3552 deltaX_[iColumn] = newPrimal;
3553 offsetObjective += newPrimal * deltaSU_[iColumn];
3554 }
3555 } else {
3556 deltaX_[iColumn] = solution_[iColumn];
3557 diagonal_[iColumn] = 0.0;
3558 offsetObjective += solution_[iColumn] * deltaSU_[iColumn];
3559 if (upper_[iColumn] - lower_[iColumn] > 1.0e-5) {
3560 if (solution_[iColumn] < lower_[iColumn] + 1.0e-8 && dj_[iColumn] < -1.0e-8) {
3561 if (-dj_[iColumn] > maximumDJInfeasibility) {
3562 maximumDJInfeasibility = -dj_[iColumn];
3563 }
3564 }
3565 if (solution_[iColumn] > upper_[iColumn] - 1.0e-8 && dj_[iColumn] > 1.0e-8) {
3566 if (dj_[iColumn] > maximumDJInfeasibility) {
3567 maximumDJInfeasibility = dj_[iColumn];
3568 }
3569 }
3570 }
3571 }
3572 primalObjectiveValue += solution_[iColumn] * cost_[iColumn];
3573 }
3574 handler_->message(CLP_BARRIER_DIAGONAL, messages_)
3575 << static_cast<double>(largestDiagonal) << static_cast<double>(smallestDiagonal)
3576 << CoinMessageEol;
3577#if 0
3578 // If diagonal wild - kill some
3579 if (largestDiagonal > 1.0e17 * smallestDiagonal) {
3580 CoinWorkDouble killValue = largestDiagonal * 1.0e-17;
3581 for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
3582 if (CoinAbs(diagonal_[iColumn]) < killValue)
3583 diagonal_[iolumn] = 0.0;
3584 }
3585 }
3586#endif
3587 // update rhs region
3588 multiplyAdd(deltaX_ + numberColumns_, numberRows_, -1.0, rhsFixRegion_, 1.0);
3589 matrix_->times(1.0, deltaX_, rhsFixRegion_);
3590 primalObjectiveValue += 0.5 * gamma2 * gammaOffset + 0.5 * quadraticOffset;
3591 if (quadraticOffset) {
3592 // printf("gamma offset %g %g, quadoffset %g\n",gammaOffset,gamma2*gammaOffset,quadraticOffset);
3593 }
3594
3595 dualObjectiveValue += offsetObjective + dualFake;
3596 dualObjectiveValue -= 0.5 * gamma2 * gammaOffset + 0.5 * quadraticOffset;
3597 if (numberIncreased || numberDecreased) {
3598 handler_->message(CLP_BARRIER_SLACKS, messages_)
3599 << numberIncreased << numberDecreased
3600 << CoinMessageEol;
3601 }
3602 if (maximumDJInfeasibility) {
3603 handler_->message(CLP_BARRIER_DUALINF, messages_)
3604 << static_cast<double>(maximumDJInfeasibility)
3605 << CoinMessageEol;
3606 }
3607 // Need to rethink (but it is only for printing)
3608 sumPrimalInfeasibilities_ = maximumRhsInfeasibility;
3609 sumDualInfeasibilities_ = maximumDualError;
3610 maximumBoundInfeasibility_ = maximumBoundInfeasibility;
3611 //compute error and fixed RHS
3612 multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, errorRegion_, 0.0);
3613 matrix_->times(1.0, solution_, errorRegion_);
3614 maximumDualError_ = maximumDualError;
3615 maximumBoundInfeasibility_ = maximumBoundInfeasibility;
3616 solutionNorm_ = solutionNorm;
3617 //finish off objective computation
3618 primalObjective_ = primalObjectiveValue * scaleFactor_;
3619 CoinWorkDouble dualValue2 = innerProduct(dualArray, numberRows_,
3620 rhsFixRegion_);
3621 dualObjectiveValue -= dualValue2;
3622 dualObjective_ = dualObjectiveValue * scaleFactor_;
3623 if (numberKilled) {
3624 handler_->message(CLP_BARRIER_KILLED, messages_)
3625 << numberKilled
3626 << CoinMessageEol;
3627 }
3628 CoinWorkDouble maximumRHSError1 = 0.0;
3629 CoinWorkDouble maximumRHSError2 = 0.0;
3630 CoinWorkDouble primalOffset = 0.0;
3631 char * dropped = cholesky_->rowsDropped();
3632 for (iRow = 0; iRow < numberRows_; iRow++) {
3633 CoinWorkDouble value = errorRegion_[iRow];
3634 if (!dropped[iRow]) {
3635 if (CoinAbs(value) > maximumRHSError1) {
3636 maximumRHSError1 = CoinAbs(value);
3637 }
3638 } else {
3639 if (CoinAbs(value) > maximumRHSError2) {
3640 maximumRHSError2 = CoinAbs(value);
3641 }
3642 primalOffset += value * dualArray[iRow];
3643 }
3644 }
3645 primalObjective_ -= primalOffset * scaleFactor_;
3646 if (maximumRHSError1 > maximumRHSError2) {
3647 maximumRHSError_ = maximumRHSError1;
3648 } else {
3649 maximumRHSError_ = maximumRHSError1; //note change
3650 if (maximumRHSError2 > primalTolerance()) {
3651 handler_->message(CLP_BARRIER_ABS_DROPPED, messages_)
3652 << static_cast<double>(maximumRHSError2)
3653 << CoinMessageEol;
3654 }
3655 }
3656 objectiveNorm_ = maximumAbsElement(dualArray, numberRows_);
3657 if (objectiveNorm_ < 1.0e-12) {
3658 objectiveNorm_ = 1.0e-12;
3659 }
3660 if (objectiveNorm_ < baseObjectiveNorm_) {
3661 //std::cout<<" base "<<baseObjectiveNorm_<<" "<<objectiveNorm_<<std::endl;
3662 if (objectiveNorm_ < baseObjectiveNorm_ * 1.0e-4) {
3663 objectiveNorm_ = baseObjectiveNorm_ * 1.0e-4;
3664 }
3665 }
3666 bool primalFeasible = true;
3667 if (maximumRHSError_ > primalTolerance() ||
3668 maximumDualError_ > dualTolerance / scaleFactor_) {
3669 handler_->message(CLP_BARRIER_ABS_ERROR, messages_)
3670 << static_cast<double>(maximumRHSError_) << static_cast<double>(maximumDualError_)
3671 << CoinMessageEol;
3672 }
3673 if (rhsNorm_ > solutionNorm_) {
3674 solutionNorm_ = rhsNorm_;
3675 }
3676 CoinWorkDouble scaledRHSError = maximumRHSError_ / (solutionNorm_ + 10.0);
3677 bool dualFeasible = true;
3678#if KEEP_GOING_IF_FIXED > 5
3679 if (maximumBoundInfeasibility_ > primalTolerance() ||
3680 scaledRHSError > primalTolerance())
3681 primalFeasible = false;
3682#else
3683 if (maximumBoundInfeasibility_ > primalTolerance() ||
3684 scaledRHSError > CoinMax(CoinMin(100.0 * primalTolerance(), 1.0e-5),
3685 primalTolerance()))
3686 primalFeasible = false;
3687#endif
3688 // relax dual test if obj big and gap smallish
3689 CoinWorkDouble gap = CoinAbs(primalObjective_ - dualObjective_);
3690 CoinWorkDouble sizeObj = CoinMin(CoinAbs(primalObjective_), CoinAbs(dualObjective_)) + 1.0e-50;
3691 //printf("gap %g sizeObj %g ratio %g comp %g\n",
3692 // gap,sizeObj,gap/sizeObj,complementarityGap_);
3693 if (numberIterations_ > 100 && gap / sizeObj < 1.0e-9 && complementarityGap_ < 1.0e-7 * sizeObj)
3694 dualTolerance *= 1.0e2;
3695 if (maximumDualError_ > objectiveNorm_ * dualTolerance)
3696 dualFeasible = false;
3697 if (!primalFeasible || !dualFeasible) {
3698 handler_->message(CLP_BARRIER_FEASIBLE, messages_)
3699 << static_cast<double>(maximumBoundInfeasibility_) << static_cast<double>(scaledRHSError)
3700 << static_cast<double>(maximumDualError_ / objectiveNorm_)
3701 << CoinMessageEol;
3702 }
3703 if (!gonePrimalFeasible_) {
3704 gonePrimalFeasible_ = primalFeasible;
3705 } else if (!primalFeasible) {
3706 gonePrimalFeasible_ = primalFeasible;
3707 if (!numberKilled) {
3708 handler_->message(CLP_BARRIER_GONE_INFEASIBLE, messages_)
3709 << CoinMessageEol;
3710 }
3711 }
3712 if (!goneDualFeasible_) {
3713 goneDualFeasible_ = dualFeasible;
3714 } else if (!dualFeasible) {
3715 handler_->message(CLP_BARRIER_GONE_INFEASIBLE, messages_)
3716 << CoinMessageEol;
3717 goneDualFeasible_ = dualFeasible;
3718 }
3719 //objectiveValue();
3720 if (solutionNorm_ > 1.0e40) {
3721 std::cout << "primal off to infinity" << std::endl;
3722 abort();
3723 }
3724 if (objectiveNorm_ > 1.0e40) {
3725 std::cout << "dual off to infinity" << std::endl;
3726 abort();
3727 }
3728 handler_->message(CLP_BARRIER_STEP, messages_)
3729 << static_cast<double>(actualPrimalStep_)
3730 << static_cast<double>(actualDualStep_)
3731 << static_cast<double>(mu_)
3732 << CoinMessageEol;
3733 numberIterations_++;
3734 return numberKilled;
3735}
3736// Save info on products of affine deltaSU*deltaW and deltaSL*deltaZ
3737CoinWorkDouble
3738ClpPredictorCorrector::affineProduct()
3739{
3740 CoinWorkDouble product = 0.0;
3741 //IF zVec starts as 0 then deltaZ always zero
3742 //(remember if free then zVec not 0)
3743 //I think free can be done with careful use of boundSlacks to zero
3744 //out all we want
3745 for (int iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
3746 CoinWorkDouble w3 = deltaZ_[iColumn] * deltaX_[iColumn];
3747 CoinWorkDouble w4 = -deltaW_[iColumn] * deltaX_[iColumn];
3748 if (lowerBound(iColumn)) {
3749 w3 += deltaZ_[iColumn] * (solution_[iColumn] - lowerSlack_[iColumn] - lower_[iColumn]);
3750 product += w3;
3751 }
3752 if (upperBound(iColumn)) {
3753 w4 += deltaW_[iColumn] * (-solution_[iColumn] - upperSlack_[iColumn] + upper_[iColumn]);
3754 product += w4;
3755 }
3756 }
3757 return product;
3758}
3759//See exactly what would happen given current deltas
3760void
3761ClpPredictorCorrector::debugMove(int /*phase*/,
3762 CoinWorkDouble primalStep, CoinWorkDouble dualStep)
3763{
3764#ifndef SOME_DEBUG
3765 return;
3766#endif
3767 CoinWorkDouble * dualArray = reinterpret_cast<CoinWorkDouble *>(dual_);
3768 int numberTotal = numberRows_ + numberColumns_;
3769 CoinWorkDouble * dualNew = ClpCopyOfArray(dualArray, numberRows_);
3770 CoinWorkDouble * errorRegionNew = new CoinWorkDouble [numberRows_];
3771 CoinWorkDouble * rhsFixRegionNew = new CoinWorkDouble [numberRows_];
3772 CoinWorkDouble * primalNew = ClpCopyOfArray(solution_, numberTotal);
3773 CoinWorkDouble * djNew = new CoinWorkDouble[numberTotal];
3774 //update pi
3775 multiplyAdd(deltaY_, numberRows_, dualStep, dualNew, 1.0);
3776 // do reduced costs
3777 CoinMemcpyN(dualNew, numberRows_, djNew + numberColumns_);
3778 CoinMemcpyN(cost_, numberColumns_, djNew);
3779 matrix_->transposeTimes(-1.0, dualNew, djNew);
3780 // update x
3781 int iColumn;
3782 for (iColumn = 0; iColumn < numberTotal; iColumn++) {
3783 if (!flagged(iColumn))
3784 primalNew[iColumn] += primalStep * deltaX_[iColumn];
3785 }
3786 CoinWorkDouble quadraticOffset = quadraticDjs(djNew, primalNew, 1.0);
3787 CoinZeroN(errorRegionNew, numberRows_);
3788 CoinZeroN(rhsFixRegionNew, numberRows_);
3789 CoinWorkDouble maximumBoundInfeasibility = 0.0;
3790 CoinWorkDouble maximumDualError = 1.0e-12;
3791 CoinWorkDouble primalObjectiveValue = 0.0;
3792 CoinWorkDouble dualObjectiveValue = 0.0;
3793 CoinWorkDouble solutionNorm = 1.0e-12;
3794 const CoinWorkDouble largeFactor = 1.0e2;
3795 CoinWorkDouble largeGap = largeFactor * solutionNorm_;
3796 if (largeGap < largeFactor) {
3797 largeGap = largeFactor;
3798 }
3799 CoinWorkDouble dualFake = 0.0;
3800 CoinWorkDouble dualTolerance = dblParam_[ClpDualTolerance];
3801 dualTolerance = dualTolerance / scaleFactor_;
3802 if (dualTolerance < 1.0e-12) {
3803 dualTolerance = 1.0e-12;
3804 }
3805 CoinWorkDouble newGap = 0.0;
3806 CoinWorkDouble offsetObjective = 0.0;
3807 CoinWorkDouble gamma2 = gamma_ * gamma_; // gamma*gamma will be added to diagonal
3808 CoinWorkDouble gammaOffset = 0.0;
3809 CoinWorkDouble maximumDjInfeasibility = 0.0;
3810 for ( iColumn = 0; iColumn < numberTotal; iColumn++) {
3811 if (!flagged(iColumn)) {
3812 CoinWorkDouble reducedCost = djNew[iColumn];
3813 CoinWorkDouble zValue = zVec_[iColumn] + dualStep * deltaZ_[iColumn];
3814 CoinWorkDouble wValue = wVec_[iColumn] + dualStep * deltaW_[iColumn];
3815 CoinWorkDouble thisWeight = deltaX_[iColumn];
3816 CoinWorkDouble oldPrimal = solution_[iColumn];
3817 CoinWorkDouble newPrimal = primalNew[iColumn];
3818 CoinWorkDouble lowerBoundInfeasibility = 0.0;
3819 CoinWorkDouble upperBoundInfeasibility = 0.0;
3820 if (lowerBound(iColumn)) {
3821 CoinWorkDouble oldSlack = lowerSlack_[iColumn];
3822 CoinWorkDouble newSlack =
3823 lowerSlack_[iColumn] + primalStep * (oldPrimal - oldSlack
3824 + thisWeight - lower_[iColumn]);
3825 if (zValue > dualTolerance) {
3826 dualObjectiveValue += lower_[iColumn] * zVec_[iColumn];
3827 }
3828 lowerBoundInfeasibility = CoinAbs(newPrimal - newSlack - lower_[iColumn]);
3829 newGap += newSlack * zValue;
3830 }
3831 if (upperBound(iColumn)) {
3832 CoinWorkDouble oldSlack = upperSlack_[iColumn];
3833 CoinWorkDouble newSlack =
3834 upperSlack_[iColumn] + primalStep * (-oldPrimal - oldSlack
3835 - thisWeight + upper_[iColumn]);
3836 if (wValue > dualTolerance) {
3837 dualObjectiveValue -= upper_[iColumn] * wVec_[iColumn];
3838 }
3839 upperBoundInfeasibility = CoinAbs(newPrimal + newSlack - upper_[iColumn]);
3840 newGap += newSlack * wValue;
3841 }
3842 if (CoinAbs(newPrimal) > solutionNorm) {
3843 solutionNorm = CoinAbs(newPrimal);
3844 }
3845 CoinWorkDouble gammaTerm = gamma2;
3846 if (primalR_) {
3847 gammaTerm += primalR_[iColumn];
3848 quadraticOffset += newPrimal * newPrimal * primalR_[iColumn];
3849 }
3850 CoinWorkDouble dualInfeasibility =
3851 reducedCost - zValue + wValue + gammaTerm * newPrimal;
3852 if (CoinAbs(dualInfeasibility) > dualTolerance) {
3853 dualFake += newPrimal * dualInfeasibility;
3854 }
3855 if (lowerBoundInfeasibility > maximumBoundInfeasibility) {
3856 maximumBoundInfeasibility = lowerBoundInfeasibility;
3857 }
3858 if (upperBoundInfeasibility > maximumBoundInfeasibility) {
3859 maximumBoundInfeasibility = upperBoundInfeasibility;
3860 }
3861 dualInfeasibility = CoinAbs(dualInfeasibility);
3862 if (dualInfeasibility > maximumDualError) {
3863 //printf("bad dual %d %g\n",iColumn,
3864 // reducedCost-zVec_[iColumn]+wVec_[iColumn]+gammaTerm*newPrimal);
3865 maximumDualError = dualInfeasibility;
3866 }
3867 gammaOffset += newPrimal * newPrimal;
3868 djNew[iColumn] = 0.0;
3869 } else {
3870 offsetObjective += primalNew[iColumn] * cost_[iColumn];
3871 if (upper_[iColumn] - lower_[iColumn] > 1.0e-5) {
3872 if (primalNew[iColumn] < lower_[iColumn] + 1.0e-8 && djNew[iColumn] < -1.0e-8) {
3873 if (-djNew[iColumn] > maximumDjInfeasibility) {
3874 maximumDjInfeasibility = -djNew[iColumn];
3875 }
3876 }
3877 if (primalNew[iColumn] > upper_[iColumn] - 1.0e-8 && djNew[iColumn] > 1.0e-8) {
3878 if (djNew[iColumn] > maximumDjInfeasibility) {
3879 maximumDjInfeasibility = djNew[iColumn];
3880 }
3881 }
3882 }
3883 djNew[iColumn] = primalNew[iColumn];
3884 }
3885 primalObjectiveValue += solution_[iColumn] * cost_[iColumn];
3886 }
3887 // update rhs region
3888 multiplyAdd(djNew + numberColumns_, numberRows_, -1.0, rhsFixRegionNew, 1.0);
3889 matrix_->times(1.0, djNew, rhsFixRegionNew);
3890 primalObjectiveValue += 0.5 * gamma2 * gammaOffset + 0.5 * quadraticOffset;
3891 dualObjectiveValue += offsetObjective + dualFake;
3892 dualObjectiveValue -= 0.5 * gamma2 * gammaOffset + 0.5 * quadraticOffset;
3893 // Need to rethink (but it is only for printing)
3894 //compute error and fixed RHS
3895 multiplyAdd(primalNew + numberColumns_, numberRows_, -1.0, errorRegionNew, 0.0);
3896 matrix_->times(1.0, primalNew, errorRegionNew);
3897 //finish off objective computation
3898 CoinWorkDouble primalObjectiveNew = primalObjectiveValue * scaleFactor_;
3899 CoinWorkDouble dualValue2 = innerProduct(dualNew, numberRows_,
3900 rhsFixRegionNew);
3901 dualObjectiveValue -= dualValue2;
3902 //CoinWorkDouble dualObjectiveNew=dualObjectiveValue*scaleFactor_;
3903 CoinWorkDouble maximumRHSError1 = 0.0;
3904 CoinWorkDouble maximumRHSError2 = 0.0;
3905 CoinWorkDouble primalOffset = 0.0;
3906 char * dropped = cholesky_->rowsDropped();
3907 int iRow;
3908 for (iRow = 0; iRow < numberRows_; iRow++) {
3909 CoinWorkDouble value = errorRegionNew[iRow];
3910 if (!dropped[iRow]) {
3911 if (CoinAbs(value) > maximumRHSError1) {
3912 maximumRHSError1 = CoinAbs(value);
3913 }
3914 } else {
3915 if (CoinAbs(value) > maximumRHSError2) {
3916 maximumRHSError2 = CoinAbs(value);
3917 }
3918 primalOffset += value * dualNew[iRow];
3919 }
3920 }
3921 primalObjectiveNew -= primalOffset * scaleFactor_;
3922 CoinWorkDouble maximumRHSError;
3923 if (maximumRHSError1 > maximumRHSError2) {
3924 maximumRHSError = maximumRHSError1;
3925 } else {
3926 maximumRHSError = maximumRHSError1; //note change
3927 if (maximumRHSError2 > primalTolerance()) {
3928 handler_->message(CLP_BARRIER_ABS_DROPPED, messages_)
3929 << static_cast<double>(maximumRHSError2)
3930 << CoinMessageEol;
3931 }
3932 }
3933 /*printf("PH %d %g, %g new comp %g, b %g, p %g, d %g\n",phase,
3934 primalStep,dualStep,newGap,maximumBoundInfeasibility,
3935 maximumRHSError,maximumDualError);
3936 if (handler_->logLevel()>1)
3937 printf(" objs %g %g\n",
3938 primalObjectiveNew,dualObjectiveNew);
3939 if (maximumDjInfeasibility) {
3940 printf(" max dj error on fixed %g\n",
3941 maximumDjInfeasibility);
3942 } */
3943 delete [] dualNew;
3944 delete [] errorRegionNew;
3945 delete [] rhsFixRegionNew;
3946 delete [] primalNew;
3947 delete [] djNew;
3948}
3949