1/* $Id: ClpPrimalColumnSteepest.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7
8#include "ClpSimplex.hpp"
9#include "ClpPrimalColumnSteepest.hpp"
10#include "CoinIndexedVector.hpp"
11#include "ClpFactorization.hpp"
12#include "ClpNonLinearCost.hpp"
13#include "ClpMessage.hpp"
14#include "CoinHelperFunctions.hpp"
15#include <stdio.h>
16//#define CLP_DEBUG
17//#############################################################################
18// Constructors / Destructor / Assignment
19//#############################################################################
20
21//-------------------------------------------------------------------
22// Default Constructor
23//-------------------------------------------------------------------
24ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (int mode)
25 : ClpPrimalColumnPivot(),
26 devex_(0.0),
27 weights_(NULL),
28 infeasible_(NULL),
29 alternateWeights_(NULL),
30 savedWeights_(NULL),
31 reference_(NULL),
32 state_(-1),
33 mode_(mode),
34 persistence_(normal),
35 numberSwitched_(0),
36 pivotSequence_(-1),
37 savedPivotSequence_(-1),
38 savedSequenceOut_(-1),
39 sizeFactorization_(0)
40{
41 type_ = 2 + 64 * mode;
42}
43//-------------------------------------------------------------------
44// Copy constructor
45//-------------------------------------------------------------------
46ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (const ClpPrimalColumnSteepest & rhs)
47 : ClpPrimalColumnPivot(rhs)
48{
49 state_ = rhs.state_;
50 mode_ = rhs.mode_;
51 persistence_ = rhs.persistence_;
52 numberSwitched_ = rhs.numberSwitched_;
53 model_ = rhs.model_;
54 pivotSequence_ = rhs.pivotSequence_;
55 savedPivotSequence_ = rhs.savedPivotSequence_;
56 savedSequenceOut_ = rhs.savedSequenceOut_;
57 sizeFactorization_ = rhs.sizeFactorization_;
58 devex_ = rhs.devex_;
59 if ((model_ && model_->whatsChanged() & 1) != 0) {
60 if (rhs.infeasible_) {
61 infeasible_ = new CoinIndexedVector(rhs.infeasible_);
62 } else {
63 infeasible_ = NULL;
64 }
65 reference_ = NULL;
66 if (rhs.weights_) {
67 assert(model_);
68 int number = model_->numberRows() + model_->numberColumns();
69 assert (number == rhs.model_->numberRows() + rhs.model_->numberColumns());
70 weights_ = new double[number];
71 CoinMemcpyN(rhs.weights_, number, weights_);
72 savedWeights_ = new double[number];
73 CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
74 if (mode_ != 1) {
75 reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
76 }
77 } else {
78 weights_ = NULL;
79 savedWeights_ = NULL;
80 }
81 if (rhs.alternateWeights_) {
82 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
83 } else {
84 alternateWeights_ = NULL;
85 }
86 } else {
87 infeasible_ = NULL;
88 reference_ = NULL;
89 weights_ = NULL;
90 savedWeights_ = NULL;
91 alternateWeights_ = NULL;
92 }
93}
94
95//-------------------------------------------------------------------
96// Destructor
97//-------------------------------------------------------------------
98ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest ()
99{
100 delete [] weights_;
101 delete infeasible_;
102 delete alternateWeights_;
103 delete [] savedWeights_;
104 delete [] reference_;
105}
106
107//----------------------------------------------------------------
108// Assignment operator
109//-------------------------------------------------------------------
110ClpPrimalColumnSteepest &
111ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest& rhs)
112{
113 if (this != &rhs) {
114 ClpPrimalColumnPivot::operator=(rhs);
115 state_ = rhs.state_;
116 mode_ = rhs.mode_;
117 persistence_ = rhs.persistence_;
118 numberSwitched_ = rhs.numberSwitched_;
119 model_ = rhs.model_;
120 pivotSequence_ = rhs.pivotSequence_;
121 savedPivotSequence_ = rhs.savedPivotSequence_;
122 savedSequenceOut_ = rhs.savedSequenceOut_;
123 sizeFactorization_ = rhs.sizeFactorization_;
124 devex_ = rhs.devex_;
125 delete [] weights_;
126 delete [] reference_;
127 reference_ = NULL;
128 delete infeasible_;
129 delete alternateWeights_;
130 delete [] savedWeights_;
131 savedWeights_ = NULL;
132 if (rhs.infeasible_ != NULL) {
133 infeasible_ = new CoinIndexedVector(rhs.infeasible_);
134 } else {
135 infeasible_ = NULL;
136 }
137 if (rhs.weights_ != NULL) {
138 assert(model_);
139 int number = model_->numberRows() + model_->numberColumns();
140 assert (number == rhs.model_->numberRows() + rhs.model_->numberColumns());
141 weights_ = new double[number];
142 CoinMemcpyN(rhs.weights_, number, weights_);
143 savedWeights_ = new double[number];
144 CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
145 if (mode_ != 1) {
146 reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
147 }
148 } else {
149 weights_ = NULL;
150 }
151 if (rhs.alternateWeights_ != NULL) {
152 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
153 } else {
154 alternateWeights_ = NULL;
155 }
156 }
157 return *this;
158}
159// These have to match ClpPackedMatrix version
160#define TRY_NORM 1.0e-4
161#define ADD_ONE 1.0
162// Returns pivot column, -1 if none
163/* The Packed CoinIndexedVector updates has cost updates - for normal LP
164 that is just +-weight where a feasibility changed. It also has
165 reduced cost from last iteration in pivot row*/
166int
167ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
168 CoinIndexedVector * spareRow1,
169 CoinIndexedVector * spareRow2,
170 CoinIndexedVector * spareColumn1,
171 CoinIndexedVector * spareColumn2)
172{
173 assert(model_);
174 if (model_->nonLinearCost()->lookBothWays() || model_->algorithm() == 2) {
175 // Do old way
176 updates->expand();
177 return pivotColumnOldMethod(updates, spareRow1, spareRow2,
178 spareColumn1, spareColumn2);
179 }
180 int number = 0;
181 int * index;
182 double tolerance = model_->currentDualTolerance();
183 // we can't really trust infeasibilities if there is dual error
184 // this coding has to mimic coding in checkDualSolution
185 double error = CoinMin(1.0e-2, model_->largestDualError());
186 // allow tolerance at least slightly bigger than standard
187 tolerance = tolerance + error;
188 int pivotRow = model_->pivotRow();
189 int anyUpdates;
190 double * infeas = infeasible_->denseVector();
191
192 // Local copy of mode so can decide what to do
193 int switchType;
194 if (mode_ == 4)
195 switchType = 5 - numberSwitched_;
196 else if (mode_ >= 10)
197 switchType = 3;
198 else
199 switchType = mode_;
200 /* switchType -
201 0 - all exact devex
202 1 - all steepest
203 2 - some exact devex
204 3 - auto some exact devex
205 4 - devex
206 5 - dantzig
207 10 - can go to mini-sprint
208 */
209 // Look at gub
210#if 1
211 model_->clpMatrix()->dualExpanded(model_, updates, NULL, 4);
212#else
213 updates->clear();
214 model_->computeDuals(NULL);
215#endif
216 if (updates->getNumElements() > 1) {
217 // would have to have two goes for devex, three for steepest
218 anyUpdates = 2;
219 } else if (updates->getNumElements()) {
220 if (updates->getIndices()[0] == pivotRow && fabs(updates->denseVector()[0]) > 1.0e-6) {
221 // reasonable size
222 anyUpdates = 1;
223 //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
224 //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
225 } else {
226 // too small
227 anyUpdates = 2;
228 }
229 } else if (pivotSequence_ >= 0) {
230 // just after re-factorization
231 anyUpdates = -1;
232 } else {
233 // sub flip - nothing to do
234 anyUpdates = 0;
235 }
236 int sequenceOut = model_->sequenceOut();
237 if (switchType == 5) {
238 // If known matrix then we will do partial pricing
239 if (model_->clpMatrix()->canDoPartialPricing()) {
240 pivotSequence_ = -1;
241 pivotRow = -1;
242 // See if to switch
243 int numberRows = model_->numberRows();
244 int numberWanted = 10;
245 int numberColumns = model_->numberColumns();
246 int numberHiddenRows = model_->clpMatrix()->hiddenRows();
247 double ratio = static_cast<double> (sizeFactorization_ + numberHiddenRows) /
248 static_cast<double> (numberRows + 2 * numberHiddenRows);
249 // Number of dual infeasibilities at last invert
250 int numberDual = model_->numberDualInfeasibilities();
251 int numberLook = CoinMin(numberDual, numberColumns / 10);
252 if (ratio < 1.0) {
253 numberWanted = 100;
254 numberLook /= 20;
255 numberWanted = CoinMax(numberWanted, numberLook);
256 } else if (ratio < 3.0) {
257 numberWanted = 500;
258 numberLook /= 15;
259 numberWanted = CoinMax(numberWanted, numberLook);
260 } else if (ratio < 4.0 || mode_ == 5) {
261 numberWanted = 1000;
262 numberLook /= 10;
263 numberWanted = CoinMax(numberWanted, numberLook);
264 } else if (mode_ != 5) {
265 switchType = 4;
266 // initialize
267 numberSwitched_++;
268 // Make sure will re-do
269 delete [] weights_;
270 weights_ = NULL;
271 model_->computeDuals(NULL);
272 saveWeights(model_, 4);
273 anyUpdates = 0;
274 COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
275 }
276 if (switchType == 5) {
277 numberLook *= 5; // needs tuning for gub
278 if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
279 COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d look %d\n",
280 sizeFactorization_, ratio, numberWanted, numberLook));
281 }
282 // Update duals and row djs
283 // Do partial pricing
284 return partialPricing(updates, spareRow2,
285 numberWanted, numberLook);
286 }
287 }
288 }
289 if (switchType == 5) {
290 if (anyUpdates > 0) {
291 justDjs(updates, spareRow2,
292 spareColumn1, spareColumn2);
293 }
294 } else if (anyUpdates == 1) {
295 if (switchType < 4) {
296 // exact etc when can use dj
297 djsAndSteepest(updates, spareRow2,
298 spareColumn1, spareColumn2);
299 } else {
300 // devex etc when can use dj
301 djsAndDevex(updates, spareRow2,
302 spareColumn1, spareColumn2);
303 }
304 } else if (anyUpdates == -1) {
305 if (switchType < 4) {
306 // exact etc when djs okay
307 justSteepest(updates, spareRow2,
308 spareColumn1, spareColumn2);
309 } else {
310 // devex etc when djs okay
311 justDevex(updates, spareRow2,
312 spareColumn1, spareColumn2);
313 }
314 } else if (anyUpdates == 2) {
315 if (switchType < 4) {
316 // exact etc when have to use pivot
317 djsAndSteepest2(updates, spareRow2,
318 spareColumn1, spareColumn2);
319 } else {
320 // devex etc when have to use pivot
321 djsAndDevex2(updates, spareRow2,
322 spareColumn1, spareColumn2);
323 }
324 }
325#ifdef CLP_DEBUG
326 alternateWeights_->checkClear();
327#endif
328 // make sure outgoing from last iteration okay
329 if (sequenceOut >= 0) {
330 ClpSimplex::Status status = model_->getStatus(sequenceOut);
331 double value = model_->reducedCost(sequenceOut);
332
333 switch(status) {
334
335 case ClpSimplex::basic:
336 case ClpSimplex::isFixed:
337 break;
338 case ClpSimplex::isFree:
339 case ClpSimplex::superBasic:
340 if (fabs(value) > FREE_ACCEPT * tolerance) {
341 // we are going to bias towards free (but only if reasonable)
342 value *= FREE_BIAS;
343 // store square in list
344 if (infeas[sequenceOut])
345 infeas[sequenceOut] = value * value; // already there
346 else
347 infeasible_->quickAdd(sequenceOut, value * value);
348 } else {
349 infeasible_->zero(sequenceOut);
350 }
351 break;
352 case ClpSimplex::atUpperBound:
353 if (value > tolerance) {
354 // store square in list
355 if (infeas[sequenceOut])
356 infeas[sequenceOut] = value * value; // already there
357 else
358 infeasible_->quickAdd(sequenceOut, value * value);
359 } else {
360 infeasible_->zero(sequenceOut);
361 }
362 break;
363 case ClpSimplex::atLowerBound:
364 if (value < -tolerance) {
365 // store square in list
366 if (infeas[sequenceOut])
367 infeas[sequenceOut] = value * value; // already there
368 else
369 infeasible_->quickAdd(sequenceOut, value * value);
370 } else {
371 infeasible_->zero(sequenceOut);
372 }
373 }
374 }
375 // update of duals finished - now do pricing
376 // See what sort of pricing
377 int numberWanted = 10;
378 number = infeasible_->getNumElements();
379 int numberColumns = model_->numberColumns();
380 if (switchType == 5) {
381 pivotSequence_ = -1;
382 pivotRow = -1;
383 // See if to switch
384 int numberRows = model_->numberRows();
385 // ratio is done on number of columns here
386 //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
387 double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
388 //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
389 if (ratio < 1.0) {
390 numberWanted = CoinMax(100, number / 200);
391 } else if (ratio < 2.0 - 1.0) {
392 numberWanted = CoinMax(500, number / 40);
393 } else if (ratio < 4.0 - 3.0 || mode_ == 5) {
394 numberWanted = CoinMax(2000, number / 10);
395 numberWanted = CoinMax(numberWanted, numberColumns / 30);
396 } else if (mode_ != 5) {
397 switchType = 4;
398 // initialize
399 numberSwitched_++;
400 // Make sure will re-do
401 delete [] weights_;
402 weights_ = NULL;
403 saveWeights(model_, 4);
404 COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
405 }
406 //if (model_->numberIterations()%1000==0)
407 //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
408 }
409 int numberRows = model_->numberRows();
410 // ratio is done on number of rows here
411 double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
412 if(switchType == 4) {
413 // Still in devex mode
414 // Go to steepest if lot of iterations?
415 if (ratio < 5.0) {
416 numberWanted = CoinMax(2000, number / 10);
417 numberWanted = CoinMax(numberWanted, numberColumns / 20);
418 } else if (ratio < 7.0) {
419 numberWanted = CoinMax(2000, number / 5);
420 numberWanted = CoinMax(numberWanted, numberColumns / 10);
421 } else {
422 // we can zero out
423 updates->clear();
424 spareColumn1->clear();
425 switchType = 3;
426 // initialize
427 pivotSequence_ = -1;
428 pivotRow = -1;
429 numberSwitched_++;
430 // Make sure will re-do
431 delete [] weights_;
432 weights_ = NULL;
433 saveWeights(model_, 4);
434 COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
435 updates->clear();
436 }
437 if (model_->numberIterations() % 1000 == 0)
438 COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
439 }
440 if (switchType < 4) {
441 if (switchType < 2 ) {
442 numberWanted = number + 1;
443 } else if (switchType == 2) {
444 numberWanted = CoinMax(2000, number / 8);
445 } else {
446 if (ratio < 1.0) {
447 numberWanted = CoinMax(2000, number / 20);
448 } else if (ratio < 5.0) {
449 numberWanted = CoinMax(2000, number / 10);
450 numberWanted = CoinMax(numberWanted, numberColumns / 40);
451 } else if (ratio < 10.0) {
452 numberWanted = CoinMax(2000, number / 8);
453 numberWanted = CoinMax(numberWanted, numberColumns / 20);
454 } else {
455 ratio = number * (ratio / 80.0);
456 if (ratio > number) {
457 numberWanted = number + 1;
458 } else {
459 numberWanted = CoinMax(2000, static_cast<int> (ratio));
460 numberWanted = CoinMax(numberWanted, numberColumns / 10);
461 }
462 }
463 }
464 //if (model_->numberIterations()%1000==0)
465 //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
466 //switchType);
467 }
468
469
470 double bestDj = 1.0e-30;
471 int bestSequence = -1;
472
473 int i, iSequence;
474 index = infeasible_->getIndices();
475 number = infeasible_->getNumElements();
476 // Re-sort infeasible every 100 pivots
477 // Not a good idea
478 if (0 && model_->factorization()->pivots() > 0 &&
479 (model_->factorization()->pivots() % 100) == 0) {
480 int nLook = model_->numberRows() + numberColumns;
481 number = 0;
482 for (i = 0; i < nLook; i++) {
483 if (infeas[i]) {
484 if (fabs(infeas[i]) > COIN_INDEXED_TINY_ELEMENT)
485 index[number++] = i;
486 else
487 infeas[i] = 0.0;
488 }
489 }
490 infeasible_->setNumElements(number);
491 }
492 if(model_->numberIterations() < model_->lastBadIteration() + 200 &&
493 model_->factorization()->pivots() > 10) {
494 // we can't really trust infeasibilities if there is dual error
495 double checkTolerance = 1.0e-8;
496 if (model_->largestDualError() > checkTolerance)
497 tolerance *= model_->largestDualError() / checkTolerance;
498 // But cap
499 tolerance = CoinMin(1000.0, tolerance);
500 }
501#ifdef CLP_DEBUG
502 if (model_->numberDualInfeasibilities() == 1)
503 printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
504 model_->largestDualError(), model_, model_->messageHandler(),
505 number);
506#endif
507 // stop last one coming immediately
508 double saveOutInfeasibility = 0.0;
509 if (sequenceOut >= 0) {
510 saveOutInfeasibility = infeas[sequenceOut];
511 infeas[sequenceOut] = 0.0;
512 }
513 if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
514 tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
515 tolerance *= tolerance; // as we are using squares
516
517 int iPass;
518 // Setup two passes
519 int start[4];
520 start[1] = number;
521 start[2] = 0;
522 double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
523 start[0] = static_cast<int> (dstart);
524 start[3] = start[0];
525 //double largestWeight=0.0;
526 //double smallestWeight=1.0e100;
527 for (iPass = 0; iPass < 2; iPass++) {
528 int end = start[2*iPass+1];
529 if (switchType < 5) {
530 for (i = start[2*iPass]; i < end; i++) {
531 iSequence = index[i];
532 double value = infeas[iSequence];
533 double weight = weights_[iSequence];
534 if (value > tolerance) {
535 //weight=1.0;
536 if (value > bestDj * weight) {
537 // check flagged variable and correct dj
538 if (!model_->flagged(iSequence)) {
539 bestDj = value / weight;
540 bestSequence = iSequence;
541 } else {
542 // just to make sure we don't exit before got something
543 numberWanted++;
544 }
545 }
546 numberWanted--;
547 }
548 if (!numberWanted)
549 break;
550 }
551 } else {
552 // Dantzig
553 for (i = start[2*iPass]; i < end; i++) {
554 iSequence = index[i];
555 double value = infeas[iSequence];
556 if (value > tolerance) {
557 if (value > bestDj) {
558 // check flagged variable and correct dj
559 if (!model_->flagged(iSequence)) {
560 bestDj = value;
561 bestSequence = iSequence;
562 } else {
563 // just to make sure we don't exit before got something
564 numberWanted++;
565 }
566 }
567 numberWanted--;
568 }
569 if (!numberWanted)
570 break;
571 }
572 }
573 if (!numberWanted)
574 break;
575 }
576 model_->clpMatrix()->setSavedBestSequence(bestSequence);
577 if (bestSequence >= 0)
578 model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
579 if (sequenceOut >= 0) {
580 infeas[sequenceOut] = saveOutInfeasibility;
581 }
582 /*if (model_->numberIterations()%100==0)
583 printf("%d best %g\n",bestSequence,bestDj);*/
584
585#ifndef NDEBUG
586 if (bestSequence >= 0) {
587 if (model_->getStatus(bestSequence) == ClpSimplex::atLowerBound)
588 assert(model_->reducedCost(bestSequence) < 0.0);
589 if (model_->getStatus(bestSequence) == ClpSimplex::atUpperBound) {
590 assert(model_->reducedCost(bestSequence) > 0.0);
591 }
592 }
593#endif
594 return bestSequence;
595}
596// Just update djs
597void
598ClpPrimalColumnSteepest::justDjs(CoinIndexedVector * updates,
599 CoinIndexedVector * spareRow2,
600 CoinIndexedVector * spareColumn1,
601 CoinIndexedVector * spareColumn2)
602{
603 int iSection, j;
604 int number = 0;
605 int * index;
606 double * updateBy;
607 double * reducedCost;
608 double tolerance = model_->currentDualTolerance();
609 // we can't really trust infeasibilities if there is dual error
610 // this coding has to mimic coding in checkDualSolution
611 double error = CoinMin(1.0e-2, model_->largestDualError());
612 // allow tolerance at least slightly bigger than standard
613 tolerance = tolerance + error;
614 int pivotRow = model_->pivotRow();
615 double * infeas = infeasible_->denseVector();
616 //updates->scanAndPack();
617 model_->factorization()->updateColumnTranspose(spareRow2, updates);
618
619 // put row of tableau in rowArray and columnArray (packed mode)
620 model_->clpMatrix()->transposeTimes(model_, -1.0,
621 updates, spareColumn2, spareColumn1);
622 // normal
623 for (iSection = 0; iSection < 2; iSection++) {
624
625 reducedCost = model_->djRegion(iSection);
626 int addSequence;
627#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
628 double slack_multiplier;
629#endif
630
631 if (!iSection) {
632 number = updates->getNumElements();
633 index = updates->getIndices();
634 updateBy = updates->denseVector();
635 addSequence = model_->numberColumns();
636#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
637 slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
638#endif
639 } else {
640 number = spareColumn1->getNumElements();
641 index = spareColumn1->getIndices();
642 updateBy = spareColumn1->denseVector();
643 addSequence = 0;
644#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
645 slack_multiplier = 1.0;
646#endif
647 }
648
649 for (j = 0; j < number; j++) {
650 int iSequence = index[j];
651 double value = reducedCost[iSequence];
652 value -= updateBy[j];
653 updateBy[j] = 0.0;
654 reducedCost[iSequence] = value;
655 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
656
657 switch(status) {
658
659 case ClpSimplex::basic:
660 infeasible_->zero(iSequence + addSequence);
661 case ClpSimplex::isFixed:
662 break;
663 case ClpSimplex::isFree:
664 case ClpSimplex::superBasic:
665 if (fabs(value) > FREE_ACCEPT * tolerance) {
666 // we are going to bias towards free (but only if reasonable)
667 value *= FREE_BIAS;
668 // store square in list
669 if (infeas[iSequence+addSequence])
670 infeas[iSequence+addSequence] = value * value; // already there
671 else
672 infeasible_->quickAdd(iSequence + addSequence, value * value);
673 } else {
674 infeasible_->zero(iSequence + addSequence);
675 }
676 break;
677 case ClpSimplex::atUpperBound:
678 iSequence += addSequence;
679 if (value > tolerance) {
680#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
681 value *= value*slack_multiplier;
682#else
683 value *= value;
684#endif
685 // store square in list
686 if (infeas[iSequence])
687 infeas[iSequence] = value; // already there
688 else
689 infeasible_->quickAdd(iSequence, value);
690 } else {
691 infeasible_->zero(iSequence);
692 }
693 break;
694 case ClpSimplex::atLowerBound:
695 iSequence += addSequence;
696 if (value < -tolerance) {
697#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
698 value *= value*slack_multiplier;
699#else
700 value *= value;
701#endif
702 // store square in list
703 if (infeas[iSequence])
704 infeas[iSequence] = value; // already there
705 else
706 infeasible_->quickAdd(iSequence, value);
707 } else {
708 infeasible_->zero(iSequence);
709 }
710 }
711 }
712 }
713 updates->setNumElements(0);
714 spareColumn1->setNumElements(0);
715 if (pivotRow >= 0) {
716 // make sure infeasibility on incoming is 0.0
717 int sequenceIn = model_->sequenceIn();
718 infeasible_->zero(sequenceIn);
719 }
720}
721// Update djs, weights for Devex
722void
723ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates,
724 CoinIndexedVector * spareRow2,
725 CoinIndexedVector * spareColumn1,
726 CoinIndexedVector * spareColumn2)
727{
728 int j;
729 int number = 0;
730 int * index;
731 double * updateBy;
732 double * reducedCost;
733 double tolerance = model_->currentDualTolerance();
734 // we can't really trust infeasibilities if there is dual error
735 // this coding has to mimic coding in checkDualSolution
736 double error = CoinMin(1.0e-2, model_->largestDualError());
737 // allow tolerance at least slightly bigger than standard
738 tolerance = tolerance + error;
739 // for weights update we use pivotSequence
740 // unset in case sub flip
741 assert (pivotSequence_ >= 0);
742 assert (model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
743 pivotSequence_ = -1;
744 double * infeas = infeasible_->denseVector();
745 //updates->scanAndPack();
746 model_->factorization()->updateColumnTranspose(spareRow2, updates);
747 // and we can see if reference
748 //double referenceIn = 0.0;
749 int sequenceIn = model_->sequenceIn();
750 //if (mode_ != 1 && reference(sequenceIn))
751 // referenceIn = 1.0;
752 // save outgoing weight round update
753 double outgoingWeight = 0.0;
754 int sequenceOut = model_->sequenceOut();
755 if (sequenceOut >= 0)
756 outgoingWeight = weights_[sequenceOut];
757
758 double scaleFactor = 1.0 / updates->denseVector()[0]; // as formula is with 1.0
759 // put row of tableau in rowArray and columnArray (packed mode)
760 model_->clpMatrix()->transposeTimes(model_, -1.0,
761 updates, spareColumn2, spareColumn1);
762 // update weights
763 double * weight;
764 int numberColumns = model_->numberColumns();
765 // rows
766 reducedCost = model_->djRegion(0);
767 int addSequence = model_->numberColumns();
768
769 number = updates->getNumElements();
770 index = updates->getIndices();
771 updateBy = updates->denseVector();
772 weight = weights_ + numberColumns;
773 // Devex
774 for (j = 0; j < number; j++) {
775 double thisWeight;
776 double pivot;
777 double value3;
778 int iSequence = index[j];
779 double value = reducedCost[iSequence];
780 double value2 = updateBy[j];
781 updateBy[j] = 0.0;
782 value -= value2;
783 reducedCost[iSequence] = value;
784 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
785
786 switch(status) {
787
788 case ClpSimplex::basic:
789 infeasible_->zero(iSequence + addSequence);
790 case ClpSimplex::isFixed:
791 break;
792 case ClpSimplex::isFree:
793 case ClpSimplex::superBasic:
794 thisWeight = weight[iSequence];
795 // row has -1
796 pivot = value2 * scaleFactor;
797 value3 = pivot * pivot * devex_;
798 if (reference(iSequence + numberColumns))
799 value3 += 1.0;
800 weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
801 if (fabs(value) > FREE_ACCEPT * tolerance) {
802 // we are going to bias towards free (but only if reasonable)
803 value *= FREE_BIAS;
804 // store square in list
805 if (infeas[iSequence+addSequence])
806 infeas[iSequence+addSequence] = value * value; // already there
807 else
808 infeasible_->quickAdd(iSequence + addSequence, value * value);
809 } else {
810 infeasible_->zero(iSequence + addSequence);
811 }
812 break;
813 case ClpSimplex::atUpperBound:
814 thisWeight = weight[iSequence];
815 // row has -1
816 pivot = value2 * scaleFactor;
817 value3 = pivot * pivot * devex_;
818 if (reference(iSequence + numberColumns))
819 value3 += 1.0;
820 weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
821 iSequence += addSequence;
822 if (value > tolerance) {
823 // store square in list
824#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
825 value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
826#else
827 value *= value;
828#endif
829 if (infeas[iSequence])
830 infeas[iSequence] = value; // already there
831 else
832 infeasible_->quickAdd(iSequence , value);
833 } else {
834 infeasible_->zero(iSequence);
835 }
836 break;
837 case ClpSimplex::atLowerBound:
838 thisWeight = weight[iSequence];
839 // row has -1
840 pivot = value2 * scaleFactor;
841 value3 = pivot * pivot * devex_;
842 if (reference(iSequence + numberColumns))
843 value3 += 1.0;
844 weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
845 iSequence += addSequence;
846 if (value < -tolerance) {
847 // store square in list
848#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
849 value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
850#else
851 value *= value;
852#endif
853 if (infeas[iSequence])
854 infeas[iSequence] = value; // already there
855 else
856 infeasible_->quickAdd(iSequence , value);
857 } else {
858 infeasible_->zero(iSequence);
859 }
860 }
861 }
862
863 // columns
864 weight = weights_;
865
866 scaleFactor = -scaleFactor;
867 reducedCost = model_->djRegion(1);
868 number = spareColumn1->getNumElements();
869 index = spareColumn1->getIndices();
870 updateBy = spareColumn1->denseVector();
871
872 // Devex
873
874 for (j = 0; j < number; j++) {
875 double thisWeight;
876 double pivot;
877 double value3;
878 int iSequence = index[j];
879 double value = reducedCost[iSequence];
880 double value2 = updateBy[j];
881 value -= value2;
882 updateBy[j] = 0.0;
883 reducedCost[iSequence] = value;
884 ClpSimplex::Status status = model_->getStatus(iSequence);
885
886 switch(status) {
887
888 case ClpSimplex::basic:
889 infeasible_->zero(iSequence);
890 case ClpSimplex::isFixed:
891 break;
892 case ClpSimplex::isFree:
893 case ClpSimplex::superBasic:
894 thisWeight = weight[iSequence];
895 // row has -1
896 pivot = value2 * scaleFactor;
897 value3 = pivot * pivot * devex_;
898 if (reference(iSequence))
899 value3 += 1.0;
900 weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
901 if (fabs(value) > FREE_ACCEPT * tolerance) {
902 // we are going to bias towards free (but only if reasonable)
903 value *= FREE_BIAS;
904 // store square in list
905 if (infeas[iSequence])
906 infeas[iSequence] = value * value; // already there
907 else
908 infeasible_->quickAdd(iSequence, value * value);
909 } else {
910 infeasible_->zero(iSequence);
911 }
912 break;
913 case ClpSimplex::atUpperBound:
914 thisWeight = weight[iSequence];
915 // row has -1
916 pivot = value2 * scaleFactor;
917 value3 = pivot * pivot * devex_;
918 if (reference(iSequence))
919 value3 += 1.0;
920 weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
921 if (value > tolerance) {
922 // store square in list
923 if (infeas[iSequence])
924 infeas[iSequence] = value * value; // already there
925 else
926 infeasible_->quickAdd(iSequence, value * value);
927 } else {
928 infeasible_->zero(iSequence);
929 }
930 break;
931 case ClpSimplex::atLowerBound:
932 thisWeight = weight[iSequence];
933 // row has -1
934 pivot = value2 * scaleFactor;
935 value3 = pivot * pivot * devex_;
936 if (reference(iSequence))
937 value3 += 1.0;
938 weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
939 if (value < -tolerance) {
940 // store square in list
941 if (infeas[iSequence])
942 infeas[iSequence] = value * value; // already there
943 else
944 infeasible_->quickAdd(iSequence, value * value);
945 } else {
946 infeasible_->zero(iSequence);
947 }
948 }
949 }
950 // restore outgoing weight
951 if (sequenceOut >= 0)
952 weights_[sequenceOut] = outgoingWeight;
953 // make sure infeasibility on incoming is 0.0
954 infeasible_->zero(sequenceIn);
955 spareRow2->setNumElements(0);
956 //#define SOME_DEBUG_1
957#ifdef SOME_DEBUG_1
958 // check for accuracy
959 int iCheck = 892;
960 //printf("weight for iCheck is %g\n",weights_[iCheck]);
961 int numberRows = model_->numberRows();
962 //int numberColumns = model_->numberColumns();
963 for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
964 if (model_->getStatus(iCheck) != ClpSimplex::basic &&
965 !model_->getStatus(iCheck) != ClpSimplex::isFixed)
966 checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
967 }
968#endif
969 updates->setNumElements(0);
970 spareColumn1->setNumElements(0);
971}
972// Update djs, weights for Steepest
973void
974ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector * updates,
975 CoinIndexedVector * spareRow2,
976 CoinIndexedVector * spareColumn1,
977 CoinIndexedVector * spareColumn2)
978{
979 int j;
980 int number = 0;
981 int * index;
982 double * updateBy;
983 double * reducedCost;
984 double tolerance = model_->currentDualTolerance();
985 // we can't really trust infeasibilities if there is dual error
986 // this coding has to mimic coding in checkDualSolution
987 double error = CoinMin(1.0e-2, model_->largestDualError());
988 // allow tolerance at least slightly bigger than standard
989 tolerance = tolerance + error;
990 // for weights update we use pivotSequence
991 // unset in case sub flip
992 assert (pivotSequence_ >= 0);
993 assert (model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
994 double * infeas = infeasible_->denseVector();
995 double scaleFactor = 1.0 / updates->denseVector()[0]; // as formula is with 1.0
996 assert (updates->getIndices()[0] == pivotSequence_);
997 pivotSequence_ = -1;
998 //updates->scanAndPack();
999 model_->factorization()->updateColumnTranspose(spareRow2, updates);
1000 //alternateWeights_->scanAndPack();
1001 model_->factorization()->updateColumnTranspose(spareRow2,
1002 alternateWeights_);
1003 // and we can see if reference
1004 int sequenceIn = model_->sequenceIn();
1005 double referenceIn;
1006 if (mode_ != 1) {
1007 if(reference(sequenceIn))
1008 referenceIn = 1.0;
1009 else
1010 referenceIn = 0.0;
1011 } else {
1012 referenceIn = -1.0;
1013 }
1014 // save outgoing weight round update
1015 double outgoingWeight = 0.0;
1016 int sequenceOut = model_->sequenceOut();
1017 if (sequenceOut >= 0)
1018 outgoingWeight = weights_[sequenceOut];
1019 // update row weights here so we can scale alternateWeights_
1020 // update weights
1021 double * weight;
1022 double * other = alternateWeights_->denseVector();
1023 int numberColumns = model_->numberColumns();
1024 // rows
1025 reducedCost = model_->djRegion(0);
1026 int addSequence = model_->numberColumns();
1027
1028 number = updates->getNumElements();
1029 index = updates->getIndices();
1030 updateBy = updates->denseVector();
1031 weight = weights_ + numberColumns;
1032
1033 for (j = 0; j < number; j++) {
1034 double thisWeight;
1035 double pivot;
1036 double modification;
1037 double pivotSquared;
1038 int iSequence = index[j];
1039 double value2 = updateBy[j];
1040 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
1041 double value;
1042
1043 switch(status) {
1044
1045 case ClpSimplex::basic:
1046 infeasible_->zero(iSequence + addSequence);
1047 reducedCost[iSequence] = 0.0;
1048 case ClpSimplex::isFixed:
1049 break;
1050 case ClpSimplex::isFree:
1051 case ClpSimplex::superBasic:
1052 value = reducedCost[iSequence] - value2;
1053 modification = other[iSequence];
1054 thisWeight = weight[iSequence];
1055 // row has -1
1056 pivot = value2 * scaleFactor;
1057 pivotSquared = pivot * pivot;
1058
1059 thisWeight += pivotSquared * devex_ + pivot * modification;
1060 reducedCost[iSequence] = value;
1061 if (thisWeight < TRY_NORM) {
1062 if (mode_ == 1) {
1063 // steepest
1064 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1065 } else {
1066 // exact
1067 thisWeight = referenceIn * pivotSquared;
1068 if (reference(iSequence + numberColumns))
1069 thisWeight += 1.0;
1070 thisWeight = CoinMax(thisWeight, TRY_NORM);
1071 }
1072 }
1073 weight[iSequence] = thisWeight;
1074 if (fabs(value) > FREE_ACCEPT * tolerance) {
1075 // we are going to bias towards free (but only if reasonable)
1076 value *= FREE_BIAS;
1077 // store square in list
1078 if (infeas[iSequence+addSequence])
1079 infeas[iSequence+addSequence] = value * value; // already there
1080 else
1081 infeasible_->quickAdd(iSequence + addSequence, value * value);
1082 } else {
1083 infeasible_->zero(iSequence + addSequence);
1084 }
1085 break;
1086 case ClpSimplex::atUpperBound:
1087 value = reducedCost[iSequence] - value2;
1088 modification = other[iSequence];
1089 thisWeight = weight[iSequence];
1090 // row has -1
1091 pivot = value2 * scaleFactor;
1092 pivotSquared = pivot * pivot;
1093
1094 thisWeight += pivotSquared * devex_ + pivot * modification;
1095 reducedCost[iSequence] = value;
1096 if (thisWeight < TRY_NORM) {
1097 if (mode_ == 1) {
1098 // steepest
1099 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1100 } else {
1101 // exact
1102 thisWeight = referenceIn * pivotSquared;
1103 if (reference(iSequence + numberColumns))
1104 thisWeight += 1.0;
1105 thisWeight = CoinMax(thisWeight, TRY_NORM);
1106 }
1107 }
1108 weight[iSequence] = thisWeight;
1109 iSequence += addSequence;
1110 if (value > tolerance) {
1111 // store square in list
1112#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1113 value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
1114#else
1115 value *= value;
1116#endif
1117 if (infeas[iSequence])
1118 infeas[iSequence] = value; // already there
1119 else
1120 infeasible_->quickAdd(iSequence , value);
1121 } else {
1122 infeasible_->zero(iSequence);
1123 }
1124 break;
1125 case ClpSimplex::atLowerBound:
1126 value = reducedCost[iSequence] - value2;
1127 modification = other[iSequence];
1128 thisWeight = weight[iSequence];
1129 // row has -1
1130 pivot = value2 * scaleFactor;
1131 pivotSquared = pivot * pivot;
1132
1133 thisWeight += pivotSquared * devex_ + pivot * modification;
1134 reducedCost[iSequence] = value;
1135 if (thisWeight < TRY_NORM) {
1136 if (mode_ == 1) {
1137 // steepest
1138 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1139 } else {
1140 // exact
1141 thisWeight = referenceIn * pivotSquared;
1142 if (reference(iSequence + numberColumns))
1143 thisWeight += 1.0;
1144 thisWeight = CoinMax(thisWeight, TRY_NORM);
1145 }
1146 }
1147 weight[iSequence] = thisWeight;
1148 iSequence += addSequence;
1149 if (value < -tolerance) {
1150 // store square in list
1151#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1152 value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
1153#else
1154 value *= value;
1155#endif
1156 if (infeas[iSequence])
1157 infeas[iSequence] = value; // already there
1158 else
1159 infeasible_->quickAdd(iSequence, value);
1160 } else {
1161 infeasible_->zero(iSequence);
1162 }
1163 }
1164 }
1165 // put row of tableau in rowArray and columnArray (packed)
1166 // get subset which have nonzero tableau elements
1167 transposeTimes2(updates, spareColumn1, alternateWeights_, spareColumn2, spareRow2,
1168 -scaleFactor);
1169 // zero updateBy
1170 CoinZeroN(updateBy, number);
1171 alternateWeights_->clear();
1172 // columns
1173 assert (scaleFactor);
1174 reducedCost = model_->djRegion(1);
1175 number = spareColumn1->getNumElements();
1176 index = spareColumn1->getIndices();
1177 updateBy = spareColumn1->denseVector();
1178
1179 for (j = 0; j < number; j++) {
1180 int iSequence = index[j];
1181 double value = reducedCost[iSequence];
1182 double value2 = updateBy[j];
1183 updateBy[j] = 0.0;
1184 value -= value2;
1185 reducedCost[iSequence] = value;
1186 ClpSimplex::Status status = model_->getStatus(iSequence);
1187
1188 switch(status) {
1189
1190 case ClpSimplex::basic:
1191 case ClpSimplex::isFixed:
1192 break;
1193 case ClpSimplex::isFree:
1194 case ClpSimplex::superBasic:
1195 if (fabs(value) > FREE_ACCEPT * tolerance) {
1196 // we are going to bias towards free (but only if reasonable)
1197 value *= FREE_BIAS;
1198 // store square in list
1199 if (infeas[iSequence])
1200 infeas[iSequence] = value * value; // already there
1201 else
1202 infeasible_->quickAdd(iSequence, value * value);
1203 } else {
1204 infeasible_->zero(iSequence);
1205 }
1206 break;
1207 case ClpSimplex::atUpperBound:
1208 if (value > tolerance) {
1209 // store square in list
1210 if (infeas[iSequence])
1211 infeas[iSequence] = value * value; // already there
1212 else
1213 infeasible_->quickAdd(iSequence, value * value);
1214 } else {
1215 infeasible_->zero(iSequence);
1216 }
1217 break;
1218 case ClpSimplex::atLowerBound:
1219 if (value < -tolerance) {
1220 // store square in list
1221 if (infeas[iSequence])
1222 infeas[iSequence] = value * value; // already there
1223 else
1224 infeasible_->quickAdd(iSequence, value * value);
1225 } else {
1226 infeasible_->zero(iSequence);
1227 }
1228 }
1229 }
1230 // restore outgoing weight
1231 if (sequenceOut >= 0)
1232 weights_[sequenceOut] = outgoingWeight;
1233 // make sure infeasibility on incoming is 0.0
1234 infeasible_->zero(sequenceIn);
1235 spareColumn2->setNumElements(0);
1236 //#define SOME_DEBUG_1
1237#ifdef SOME_DEBUG_1
1238 // check for accuracy
1239 int iCheck = 892;
1240 //printf("weight for iCheck is %g\n",weights_[iCheck]);
1241 int numberRows = model_->numberRows();
1242 //int numberColumns = model_->numberColumns();
1243 for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1244 if (model_->getStatus(iCheck) != ClpSimplex::basic &&
1245 !model_->getStatus(iCheck) != ClpSimplex::isFixed)
1246 checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
1247 }
1248#endif
1249 updates->setNumElements(0);
1250 spareColumn1->setNumElements(0);
1251}
1252// Update djs, weights for Devex
1253void
1254ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates,
1255 CoinIndexedVector * spareRow2,
1256 CoinIndexedVector * spareColumn1,
1257 CoinIndexedVector * spareColumn2)
1258{
1259 int iSection, j;
1260 int number = 0;
1261 int * index;
1262 double * updateBy;
1263 double * reducedCost;
1264 // dj could be very small (or even zero - take care)
1265 double dj = model_->dualIn();
1266 double tolerance = model_->currentDualTolerance();
1267 // we can't really trust infeasibilities if there is dual error
1268 // this coding has to mimic coding in checkDualSolution
1269 double error = CoinMin(1.0e-2, model_->largestDualError());
1270 // allow tolerance at least slightly bigger than standard
1271 tolerance = tolerance + error;
1272 int pivotRow = model_->pivotRow();
1273 double * infeas = infeasible_->denseVector();
1274 //updates->scanAndPack();
1275 model_->factorization()->updateColumnTranspose(spareRow2, updates);
1276
1277 // put row of tableau in rowArray and columnArray
1278 model_->clpMatrix()->transposeTimes(model_, -1.0,
1279 updates, spareColumn2, spareColumn1);
1280 // normal
1281 for (iSection = 0; iSection < 2; iSection++) {
1282
1283 reducedCost = model_->djRegion(iSection);
1284 int addSequence;
1285#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1286 double slack_multiplier;
1287#endif
1288
1289 if (!iSection) {
1290 number = updates->getNumElements();
1291 index = updates->getIndices();
1292 updateBy = updates->denseVector();
1293 addSequence = model_->numberColumns();
1294#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1295 slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
1296#endif
1297 } else {
1298 number = spareColumn1->getNumElements();
1299 index = spareColumn1->getIndices();
1300 updateBy = spareColumn1->denseVector();
1301 addSequence = 0;
1302#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1303 slack_multiplier = 1.0;
1304#endif
1305 }
1306
1307 for (j = 0; j < number; j++) {
1308 int iSequence = index[j];
1309 double value = reducedCost[iSequence];
1310 value -= updateBy[j];
1311 updateBy[j] = 0.0;
1312 reducedCost[iSequence] = value;
1313 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
1314
1315 switch(status) {
1316
1317 case ClpSimplex::basic:
1318 infeasible_->zero(iSequence + addSequence);
1319 case ClpSimplex::isFixed:
1320 break;
1321 case ClpSimplex::isFree:
1322 case ClpSimplex::superBasic:
1323 if (fabs(value) > FREE_ACCEPT * tolerance) {
1324 // we are going to bias towards free (but only if reasonable)
1325 value *= FREE_BIAS;
1326 // store square in list
1327 if (infeas[iSequence+addSequence])
1328 infeas[iSequence+addSequence] = value * value; // already there
1329 else
1330 infeasible_->quickAdd(iSequence + addSequence, value * value);
1331 } else {
1332 infeasible_->zero(iSequence + addSequence);
1333 }
1334 break;
1335 case ClpSimplex::atUpperBound:
1336 iSequence += addSequence;
1337 if (value > tolerance) {
1338#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1339 value *= value*slack_multiplier;
1340#else
1341 value *= value;
1342#endif
1343 // store square in list
1344 if (infeas[iSequence])
1345 infeas[iSequence] = value; // already there
1346 else
1347 infeasible_->quickAdd(iSequence, value);
1348 } else {
1349 infeasible_->zero(iSequence);
1350 }
1351 break;
1352 case ClpSimplex::atLowerBound:
1353 iSequence += addSequence;
1354 if (value < -tolerance) {
1355#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1356 value *= value*slack_multiplier;
1357#else
1358 value *= value;
1359#endif
1360 // store square in list
1361 if (infeas[iSequence])
1362 infeas[iSequence] = value; // already there
1363 else
1364 infeasible_->quickAdd(iSequence, value);
1365 } else {
1366 infeasible_->zero(iSequence);
1367 }
1368 }
1369 }
1370 }
1371 // They are empty
1372 updates->setNumElements(0);
1373 spareColumn1->setNumElements(0);
1374 // make sure infeasibility on incoming is 0.0
1375 int sequenceIn = model_->sequenceIn();
1376 infeasible_->zero(sequenceIn);
1377 // for weights update we use pivotSequence
1378 if (pivotSequence_ >= 0) {
1379 pivotRow = pivotSequence_;
1380 // unset in case sub flip
1381 pivotSequence_ = -1;
1382 // make sure infeasibility on incoming is 0.0
1383 const int * pivotVariable = model_->pivotVariable();
1384 sequenceIn = pivotVariable[pivotRow];
1385 infeasible_->zero(sequenceIn);
1386 // and we can see if reference
1387 //double referenceIn = 0.0;
1388 //if (mode_ != 1 && reference(sequenceIn))
1389 // referenceIn = 1.0;
1390 // save outgoing weight round update
1391 double outgoingWeight = 0.0;
1392 int sequenceOut = model_->sequenceOut();
1393 if (sequenceOut >= 0)
1394 outgoingWeight = weights_[sequenceOut];
1395 // update weights
1396 updates->setNumElements(0);
1397 spareColumn1->setNumElements(0);
1398 // might as well set dj to 1
1399 dj = 1.0;
1400 updates->insert(pivotRow, -dj);
1401 model_->factorization()->updateColumnTranspose(spareRow2, updates);
1402 // put row of tableau in rowArray and columnArray
1403 model_->clpMatrix()->transposeTimes(model_, -1.0,
1404 updates, spareColumn2, spareColumn1);
1405 double * weight;
1406 int numberColumns = model_->numberColumns();
1407 // rows
1408 number = updates->getNumElements();
1409 index = updates->getIndices();
1410 updateBy = updates->denseVector();
1411 weight = weights_ + numberColumns;
1412
1413 assert (devex_ > 0.0);
1414 for (j = 0; j < number; j++) {
1415 int iSequence = index[j];
1416 double thisWeight = weight[iSequence];
1417 // row has -1
1418 double pivot = - updateBy[iSequence];
1419 updateBy[iSequence] = 0.0;
1420 double value = pivot * pivot * devex_;
1421 if (reference(iSequence + numberColumns))
1422 value += 1.0;
1423 weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1424 }
1425
1426 // columns
1427 weight = weights_;
1428
1429 number = spareColumn1->getNumElements();
1430 index = spareColumn1->getIndices();
1431 updateBy = spareColumn1->denseVector();
1432 for (j = 0; j < number; j++) {
1433 int iSequence = index[j];
1434 double thisWeight = weight[iSequence];
1435 // row has -1
1436 double pivot = updateBy[iSequence];
1437 updateBy[iSequence] = 0.0;
1438 double value = pivot * pivot * devex_;
1439 if (reference(iSequence))
1440 value += 1.0;
1441 weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1442 }
1443 // restore outgoing weight
1444 if (sequenceOut >= 0)
1445 weights_[sequenceOut] = outgoingWeight;
1446 spareColumn2->setNumElements(0);
1447 //#define SOME_DEBUG_1
1448#ifdef SOME_DEBUG_1
1449 // check for accuracy
1450 int iCheck = 892;
1451 //printf("weight for iCheck is %g\n",weights_[iCheck]);
1452 int numberRows = model_->numberRows();
1453 //int numberColumns = model_->numberColumns();
1454 for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1455 if (model_->getStatus(iCheck) != ClpSimplex::basic &&
1456 !model_->getStatus(iCheck) != ClpSimplex::isFixed)
1457 checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
1458 }
1459#endif
1460 updates->setNumElements(0);
1461 spareColumn1->setNumElements(0);
1462 }
1463}
1464// Update djs, weights for Steepest
1465void
1466ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector * updates,
1467 CoinIndexedVector * spareRow2,
1468 CoinIndexedVector * spareColumn1,
1469 CoinIndexedVector * spareColumn2)
1470{
1471 int iSection, j;
1472 int number = 0;
1473 int * index;
1474 double * updateBy;
1475 double * reducedCost;
1476 // dj could be very small (or even zero - take care)
1477 double dj = model_->dualIn();
1478 double tolerance = model_->currentDualTolerance();
1479 // we can't really trust infeasibilities if there is dual error
1480 // this coding has to mimic coding in checkDualSolution
1481 double error = CoinMin(1.0e-2, model_->largestDualError());
1482 // allow tolerance at least slightly bigger than standard
1483 tolerance = tolerance + error;
1484 int pivotRow = model_->pivotRow();
1485 double * infeas = infeasible_->denseVector();
1486 //updates->scanAndPack();
1487 model_->factorization()->updateColumnTranspose(spareRow2, updates);
1488
1489 // put row of tableau in rowArray and columnArray
1490 model_->clpMatrix()->transposeTimes(model_, -1.0,
1491 updates, spareColumn2, spareColumn1);
1492 // normal
1493 for (iSection = 0; iSection < 2; iSection++) {
1494
1495 reducedCost = model_->djRegion(iSection);
1496 int addSequence;
1497#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1498 double slack_multiplier;
1499#endif
1500
1501 if (!iSection) {
1502 number = updates->getNumElements();
1503 index = updates->getIndices();
1504 updateBy = updates->denseVector();
1505 addSequence = model_->numberColumns();
1506#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1507 slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
1508#endif
1509 } else {
1510 number = spareColumn1->getNumElements();
1511 index = spareColumn1->getIndices();
1512 updateBy = spareColumn1->denseVector();
1513 addSequence = 0;
1514#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1515 slack_multiplier = 1.0;
1516#endif
1517 }
1518
1519 for (j = 0; j < number; j++) {
1520 int iSequence = index[j];
1521 double value = reducedCost[iSequence];
1522 value -= updateBy[j];
1523 updateBy[j] = 0.0;
1524 reducedCost[iSequence] = value;
1525 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
1526
1527 switch(status) {
1528
1529 case ClpSimplex::basic:
1530 infeasible_->zero(iSequence + addSequence);
1531 case ClpSimplex::isFixed:
1532 break;
1533 case ClpSimplex::isFree:
1534 case ClpSimplex::superBasic:
1535 if (fabs(value) > FREE_ACCEPT * tolerance) {
1536 // we are going to bias towards free (but only if reasonable)
1537 value *= FREE_BIAS;
1538 // store square in list
1539 if (infeas[iSequence+addSequence])
1540 infeas[iSequence+addSequence] = value * value; // already there
1541 else
1542 infeasible_->quickAdd(iSequence + addSequence, value * value);
1543 } else {
1544 infeasible_->zero(iSequence + addSequence);
1545 }
1546 break;
1547 case ClpSimplex::atUpperBound:
1548 iSequence += addSequence;
1549 if (value > tolerance) {
1550#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1551 value *= value*slack_multiplier;
1552#else
1553 value *= value;
1554#endif
1555 // store square in list
1556 if (infeas[iSequence])
1557 infeas[iSequence] = value; // already there
1558 else
1559 infeasible_->quickAdd(iSequence, value);
1560 } else {
1561 infeasible_->zero(iSequence);
1562 }
1563 break;
1564 case ClpSimplex::atLowerBound:
1565 iSequence += addSequence;
1566 if (value < -tolerance) {
1567#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1568 value *= value*slack_multiplier;
1569#else
1570 value *= value;
1571#endif
1572 // store square in list
1573 if (infeas[iSequence])
1574 infeas[iSequence] = value; // already there
1575 else
1576 infeasible_->quickAdd(iSequence, value);
1577 } else {
1578 infeasible_->zero(iSequence);
1579 }
1580 }
1581 }
1582 }
1583 // we can zero out as will have to get pivot row
1584 // ***** move
1585 updates->setNumElements(0);
1586 spareColumn1->setNumElements(0);
1587 if (pivotRow >= 0) {
1588 // make sure infeasibility on incoming is 0.0
1589 int sequenceIn = model_->sequenceIn();
1590 infeasible_->zero(sequenceIn);
1591 }
1592 // for weights update we use pivotSequence
1593 pivotRow = pivotSequence_;
1594 // unset in case sub flip
1595 pivotSequence_ = -1;
1596 if (pivotRow >= 0) {
1597 // make sure infeasibility on incoming is 0.0
1598 const int * pivotVariable = model_->pivotVariable();
1599 int sequenceIn = pivotVariable[pivotRow];
1600 assert (sequenceIn == model_->sequenceIn());
1601 infeasible_->zero(sequenceIn);
1602 // and we can see if reference
1603 double referenceIn;
1604 if (mode_ != 1) {
1605 if(reference(sequenceIn))
1606 referenceIn = 1.0;
1607 else
1608 referenceIn = 0.0;
1609 } else {
1610 referenceIn = -1.0;
1611 }
1612 // save outgoing weight round update
1613 double outgoingWeight = 0.0;
1614 int sequenceOut = model_->sequenceOut();
1615 if (sequenceOut >= 0)
1616 outgoingWeight = weights_[sequenceOut];
1617 // update weights
1618 updates->setNumElements(0);
1619 spareColumn1->setNumElements(0);
1620 // might as well set dj to 1
1621 dj = -1.0;
1622 updates->createPacked(1, &pivotRow, &dj);
1623 model_->factorization()->updateColumnTranspose(spareRow2, updates);
1624 bool needSubset = (mode_ < 4 || numberSwitched_ > 1 || mode_ >= 10);
1625
1626 double * weight;
1627 double * other = alternateWeights_->denseVector();
1628 int numberColumns = model_->numberColumns();
1629 // rows
1630 number = updates->getNumElements();
1631 index = updates->getIndices();
1632 updateBy = updates->denseVector();
1633 weight = weights_ + numberColumns;
1634 if (needSubset) {
1635 // now update weight update array
1636 model_->factorization()->updateColumnTranspose(spareRow2,
1637 alternateWeights_);
1638 // do alternateWeights_ here so can scale
1639 for (j = 0; j < number; j++) {
1640 int iSequence = index[j];
1641 assert (iSequence >= 0 && iSequence < model_->numberRows());
1642 double thisWeight = weight[iSequence];
1643 // row has -1
1644 double pivot = - updateBy[j];
1645 double modification = other[iSequence];
1646 double pivotSquared = pivot * pivot;
1647
1648 thisWeight += pivotSquared * devex_ + pivot * modification;
1649 if (thisWeight < TRY_NORM) {
1650 if (mode_ == 1) {
1651 // steepest
1652 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1653 } else {
1654 // exact
1655 thisWeight = referenceIn * pivotSquared;
1656 if (reference(iSequence + numberColumns))
1657 thisWeight += 1.0;
1658 thisWeight = CoinMax(thisWeight, TRY_NORM);
1659 }
1660 }
1661 weight[iSequence] = thisWeight;
1662 }
1663 transposeTimes2(updates, spareColumn1, alternateWeights_, spareColumn2, spareRow2, 0.0);
1664 } else {
1665 // put row of tableau in rowArray and columnArray
1666 model_->clpMatrix()->transposeTimes(model_, -1.0,
1667 updates, spareColumn2, spareColumn1);
1668 }
1669
1670 if (needSubset) {
1671 CoinZeroN(updateBy, number);
1672 } else if (mode_ == 4) {
1673 // Devex
1674 assert (devex_ > 0.0);
1675 for (j = 0; j < number; j++) {
1676 int iSequence = index[j];
1677 double thisWeight = weight[iSequence];
1678 // row has -1
1679 double pivot = -updateBy[j];
1680 updateBy[j] = 0.0;
1681 double value = pivot * pivot * devex_;
1682 if (reference(iSequence + numberColumns))
1683 value += 1.0;
1684 weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1685 }
1686 }
1687
1688 // columns
1689 weight = weights_;
1690
1691 number = spareColumn1->getNumElements();
1692 index = spareColumn1->getIndices();
1693 updateBy = spareColumn1->denseVector();
1694 if (needSubset) {
1695 // Exact - already done
1696 } else if (mode_ == 4) {
1697 // Devex
1698 for (j = 0; j < number; j++) {
1699 int iSequence = index[j];
1700 double thisWeight = weight[iSequence];
1701 // row has -1
1702 double pivot = updateBy[j];
1703 updateBy[j] = 0.0;
1704 double value = pivot * pivot * devex_;
1705 if (reference(iSequence))
1706 value += 1.0;
1707 weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1708 }
1709 }
1710 // restore outgoing weight
1711 if (sequenceOut >= 0)
1712 weights_[sequenceOut] = outgoingWeight;
1713 alternateWeights_->clear();
1714 spareColumn2->setNumElements(0);
1715 //#define SOME_DEBUG_1
1716#ifdef SOME_DEBUG_1
1717 // check for accuracy
1718 int iCheck = 892;
1719 //printf("weight for iCheck is %g\n",weights_[iCheck]);
1720 int numberRows = model_->numberRows();
1721 //int numberColumns = model_->numberColumns();
1722 for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1723 if (model_->getStatus(iCheck) != ClpSimplex::basic &&
1724 !model_->getStatus(iCheck) != ClpSimplex::isFixed)
1725 checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
1726 }
1727#endif
1728 }
1729 updates->setNumElements(0);
1730 spareColumn1->setNumElements(0);
1731}
1732// Updates two arrays for steepest
1733void
1734ClpPrimalColumnSteepest::transposeTimes2(const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
1735 const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
1736 CoinIndexedVector * spare,
1737 double scaleFactor)
1738{
1739 // see if reference
1740 int sequenceIn = model_->sequenceIn();
1741 double referenceIn;
1742 if (mode_ != 1) {
1743 if(reference(sequenceIn))
1744 referenceIn = 1.0;
1745 else
1746 referenceIn = 0.0;
1747 } else {
1748 referenceIn = -1.0;
1749 }
1750 if (model_->clpMatrix()->canCombine(model_, pi1)) {
1751 // put row of tableau in rowArray and columnArray
1752 model_->clpMatrix()->transposeTimes2(model_, pi1, dj1, pi2, spare, referenceIn, devex_,
1753 reference_,
1754 weights_, scaleFactor);
1755 } else {
1756 // put row of tableau in rowArray and columnArray
1757 model_->clpMatrix()->transposeTimes(model_, -1.0,
1758 pi1, dj2, dj1);
1759 // get subset which have nonzero tableau elements
1760 model_->clpMatrix()->subsetTransposeTimes(model_, pi2, dj1, dj2);
1761 bool killDjs = (scaleFactor == 0.0);
1762 if (!scaleFactor)
1763 scaleFactor = 1.0;
1764 // columns
1765 double * weight = weights_;
1766
1767 int number = dj1->getNumElements();
1768 const int * index = dj1->getIndices();
1769 double * updateBy = dj1->denseVector();
1770 double * updateBy2 = dj2->denseVector();
1771
1772 for (int j = 0; j < number; j++) {
1773 double thisWeight;
1774 double pivot;
1775 double pivotSquared;
1776 int iSequence = index[j];
1777 double value2 = updateBy[j];
1778 if (killDjs)
1779 updateBy[j] = 0.0;
1780 double modification = updateBy2[j];
1781 updateBy2[j] = 0.0;
1782 ClpSimplex::Status status = model_->getStatus(iSequence);
1783
1784 if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) {
1785 thisWeight = weight[iSequence];
1786 pivot = value2 * scaleFactor;
1787 pivotSquared = pivot * pivot;
1788
1789 thisWeight += pivotSquared * devex_ + pivot * modification;
1790 if (thisWeight < TRY_NORM) {
1791 if (referenceIn < 0.0) {
1792 // steepest
1793 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1794 } else {
1795 // exact
1796 thisWeight = referenceIn * pivotSquared;
1797 if (reference(iSequence))
1798 thisWeight += 1.0;
1799 thisWeight = CoinMax(thisWeight, TRY_NORM);
1800 }
1801 }
1802 weight[iSequence] = thisWeight;
1803 }
1804 }
1805 }
1806 dj2->setNumElements(0);
1807}
1808// Update weights for Devex
1809void
1810ClpPrimalColumnSteepest::justDevex(CoinIndexedVector * updates,
1811 CoinIndexedVector * spareRow2,
1812 CoinIndexedVector * spareColumn1,
1813 CoinIndexedVector * spareColumn2)
1814{
1815 int j;
1816 int number = 0;
1817 int * index;
1818 double * updateBy;
1819 // dj could be very small (or even zero - take care)
1820 double dj = model_->dualIn();
1821 double tolerance = model_->currentDualTolerance();
1822 // we can't really trust infeasibilities if there is dual error
1823 // this coding has to mimic coding in checkDualSolution
1824 double error = CoinMin(1.0e-2, model_->largestDualError());
1825 // allow tolerance at least slightly bigger than standard
1826 tolerance = tolerance + error;
1827 int pivotRow = model_->pivotRow();
1828 // for weights update we use pivotSequence
1829 pivotRow = pivotSequence_;
1830 assert (pivotRow >= 0);
1831 // make sure infeasibility on incoming is 0.0
1832 const int * pivotVariable = model_->pivotVariable();
1833 int sequenceIn = pivotVariable[pivotRow];
1834 infeasible_->zero(sequenceIn);
1835 // and we can see if reference
1836 //double referenceIn = 0.0;
1837 //if (mode_ != 1 && reference(sequenceIn))
1838 // referenceIn = 1.0;
1839 // save outgoing weight round update
1840 double outgoingWeight = 0.0;
1841 int sequenceOut = model_->sequenceOut();
1842 if (sequenceOut >= 0)
1843 outgoingWeight = weights_[sequenceOut];
1844 assert (!updates->getNumElements());
1845 assert (!spareColumn1->getNumElements());
1846 // unset in case sub flip
1847 pivotSequence_ = -1;
1848 // might as well set dj to 1
1849 dj = -1.0;
1850 updates->createPacked(1, &pivotRow, &dj);
1851 model_->factorization()->updateColumnTranspose(spareRow2, updates);
1852 // put row of tableau in rowArray and columnArray
1853 model_->clpMatrix()->transposeTimes(model_, -1.0,
1854 updates, spareColumn2, spareColumn1);
1855 double * weight;
1856 int numberColumns = model_->numberColumns();
1857 // rows
1858 number = updates->getNumElements();
1859 index = updates->getIndices();
1860 updateBy = updates->denseVector();
1861 weight = weights_ + numberColumns;
1862
1863 // Devex
1864 assert (devex_ > 0.0);
1865 for (j = 0; j < number; j++) {
1866 int iSequence = index[j];
1867 double thisWeight = weight[iSequence];
1868 // row has -1
1869 double pivot = - updateBy[j];
1870 updateBy[j] = 0.0;
1871 double value = pivot * pivot * devex_;
1872 if (reference(iSequence + numberColumns))
1873 value += 1.0;
1874 weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1875 }
1876
1877 // columns
1878 weight = weights_;
1879
1880 number = spareColumn1->getNumElements();
1881 index = spareColumn1->getIndices();
1882 updateBy = spareColumn1->denseVector();
1883 // Devex
1884 for (j = 0; j < number; j++) {
1885 int iSequence = index[j];
1886 double thisWeight = weight[iSequence];
1887 // row has -1
1888 double pivot = updateBy[j];
1889 updateBy[j] = 0.0;
1890 double value = pivot * pivot * devex_;
1891 if (reference(iSequence))
1892 value += 1.0;
1893 weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1894 }
1895 // restore outgoing weight
1896 if (sequenceOut >= 0)
1897 weights_[sequenceOut] = outgoingWeight;
1898 spareColumn2->setNumElements(0);
1899 //#define SOME_DEBUG_1
1900#ifdef SOME_DEBUG_1
1901 // check for accuracy
1902 int iCheck = 892;
1903 //printf("weight for iCheck is %g\n",weights_[iCheck]);
1904 int numberRows = model_->numberRows();
1905 //int numberColumns = model_->numberColumns();
1906 for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1907 if (model_->getStatus(iCheck) != ClpSimplex::basic &&
1908 !model_->getStatus(iCheck) != ClpSimplex::isFixed)
1909 checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
1910 }
1911#endif
1912 updates->setNumElements(0);
1913 spareColumn1->setNumElements(0);
1914}
1915// Update weights for Steepest
1916void
1917ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector * updates,
1918 CoinIndexedVector * spareRow2,
1919 CoinIndexedVector * spareColumn1,
1920 CoinIndexedVector * spareColumn2)
1921{
1922 int j;
1923 int number = 0;
1924 int * index;
1925 double * updateBy;
1926 // dj could be very small (or even zero - take care)
1927 double dj = model_->dualIn();
1928 double tolerance = model_->currentDualTolerance();
1929 // we can't really trust infeasibilities if there is dual error
1930 // this coding has to mimic coding in checkDualSolution
1931 double error = CoinMin(1.0e-2, model_->largestDualError());
1932 // allow tolerance at least slightly bigger than standard
1933 tolerance = tolerance + error;
1934 int pivotRow = model_->pivotRow();
1935 // for weights update we use pivotSequence
1936 pivotRow = pivotSequence_;
1937 // unset in case sub flip
1938 pivotSequence_ = -1;
1939 assert (pivotRow >= 0);
1940 // make sure infeasibility on incoming is 0.0
1941 const int * pivotVariable = model_->pivotVariable();
1942 int sequenceIn = pivotVariable[pivotRow];
1943 infeasible_->zero(sequenceIn);
1944 // and we can see if reference
1945 double referenceIn = 0.0;
1946 if (mode_ != 1 && reference(sequenceIn))
1947 referenceIn = 1.0;
1948 // save outgoing weight round update
1949 double outgoingWeight = 0.0;
1950 int sequenceOut = model_->sequenceOut();
1951 if (sequenceOut >= 0)
1952 outgoingWeight = weights_[sequenceOut];
1953 assert (!updates->getNumElements());
1954 assert (!spareColumn1->getNumElements());
1955 // update weights
1956 //updates->setNumElements(0);
1957 //spareColumn1->setNumElements(0);
1958 // might as well set dj to 1
1959 dj = -1.0;
1960 updates->createPacked(1, &pivotRow, &dj);
1961 model_->factorization()->updateColumnTranspose(spareRow2, updates);
1962 // put row of tableau in rowArray and columnArray
1963 model_->clpMatrix()->transposeTimes(model_, -1.0,
1964 updates, spareColumn2, spareColumn1);
1965 double * weight;
1966 double * other = alternateWeights_->denseVector();
1967 int numberColumns = model_->numberColumns();
1968 // rows
1969 number = updates->getNumElements();
1970 index = updates->getIndices();
1971 updateBy = updates->denseVector();
1972 weight = weights_ + numberColumns;
1973
1974 // Exact
1975 // now update weight update array
1976 //alternateWeights_->scanAndPack();
1977 model_->factorization()->updateColumnTranspose(spareRow2,
1978 alternateWeights_);
1979 // get subset which have nonzero tableau elements
1980 model_->clpMatrix()->subsetTransposeTimes(model_, alternateWeights_,
1981 spareColumn1,
1982 spareColumn2);
1983 for (j = 0; j < number; j++) {
1984 int iSequence = index[j];
1985 double thisWeight = weight[iSequence];
1986 // row has -1
1987 double pivot = -updateBy[j];
1988 updateBy[j] = 0.0;
1989 double modification = other[iSequence];
1990 double pivotSquared = pivot * pivot;
1991
1992 thisWeight += pivotSquared * devex_ + pivot * modification;
1993 if (thisWeight < TRY_NORM) {
1994 if (mode_ == 1) {
1995 // steepest
1996 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1997 } else {
1998 // exact
1999 thisWeight = referenceIn * pivotSquared;
2000 if (reference(iSequence + numberColumns))
2001 thisWeight += 1.0;
2002 thisWeight = CoinMax(thisWeight, TRY_NORM);
2003 }
2004 }
2005 weight[iSequence] = thisWeight;
2006 }
2007
2008 // columns
2009 weight = weights_;
2010
2011 number = spareColumn1->getNumElements();
2012 index = spareColumn1->getIndices();
2013 updateBy = spareColumn1->denseVector();
2014 // Exact
2015 double * updateBy2 = spareColumn2->denseVector();
2016 for (j = 0; j < number; j++) {
2017 int iSequence = index[j];
2018 double thisWeight = weight[iSequence];
2019 double pivot = updateBy[j];
2020 updateBy[j] = 0.0;
2021 double modification = updateBy2[j];
2022 updateBy2[j] = 0.0;
2023 double pivotSquared = pivot * pivot;
2024
2025 thisWeight += pivotSquared * devex_ + pivot * modification;
2026 if (thisWeight < TRY_NORM) {
2027 if (mode_ == 1) {
2028 // steepest
2029 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
2030 } else {
2031 // exact
2032 thisWeight = referenceIn * pivotSquared;
2033 if (reference(iSequence))
2034 thisWeight += 1.0;
2035 thisWeight = CoinMax(thisWeight, TRY_NORM);
2036 }
2037 }
2038 weight[iSequence] = thisWeight;
2039 }
2040 // restore outgoing weight
2041 if (sequenceOut >= 0)
2042 weights_[sequenceOut] = outgoingWeight;
2043 alternateWeights_->clear();
2044 spareColumn2->setNumElements(0);
2045 //#define SOME_DEBUG_1
2046#ifdef SOME_DEBUG_1
2047 // check for accuracy
2048 int iCheck = 892;
2049 //printf("weight for iCheck is %g\n",weights_[iCheck]);
2050 int numberRows = model_->numberRows();
2051 //int numberColumns = model_->numberColumns();
2052 for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
2053 if (model_->getStatus(iCheck) != ClpSimplex::basic &&
2054 !model_->getStatus(iCheck) != ClpSimplex::isFixed)
2055 checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
2056 }
2057#endif
2058 updates->setNumElements(0);
2059 spareColumn1->setNumElements(0);
2060}
2061// Returns pivot column, -1 if none
2062int
2063ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector * updates,
2064 CoinIndexedVector * ,
2065 CoinIndexedVector * spareRow2,
2066 CoinIndexedVector * spareColumn1,
2067 CoinIndexedVector * spareColumn2)
2068{
2069 assert(model_);
2070 int iSection, j;
2071 int number = 0;
2072 int * index;
2073 double * updateBy;
2074 double * reducedCost;
2075 // dj could be very small (or even zero - take care)
2076 double dj = model_->dualIn();
2077 double tolerance = model_->currentDualTolerance();
2078 // we can't really trust infeasibilities if there is dual error
2079 // this coding has to mimic coding in checkDualSolution
2080 double error = CoinMin(1.0e-2, model_->largestDualError());
2081 // allow tolerance at least slightly bigger than standard
2082 tolerance = tolerance + error;
2083 int pivotRow = model_->pivotRow();
2084 int anyUpdates;
2085 double * infeas = infeasible_->denseVector();
2086
2087 // Local copy of mode so can decide what to do
2088 int switchType;
2089 if (mode_ == 4)
2090 switchType = 5 - numberSwitched_;
2091 else if (mode_ >= 10)
2092 switchType = 3;
2093 else
2094 switchType = mode_;
2095 /* switchType -
2096 0 - all exact devex
2097 1 - all steepest
2098 2 - some exact devex
2099 3 - auto some exact devex
2100 4 - devex
2101 5 - dantzig
2102 */
2103 if (updates->getNumElements()) {
2104 // would have to have two goes for devex, three for steepest
2105 anyUpdates = 2;
2106 // add in pivot contribution
2107 if (pivotRow >= 0)
2108 updates->add(pivotRow, -dj);
2109 } else if (pivotRow >= 0) {
2110 if (fabs(dj) > 1.0e-15) {
2111 // some dj
2112 updates->insert(pivotRow, -dj);
2113 if (fabs(dj) > 1.0e-6) {
2114 // reasonable size
2115 anyUpdates = 1;
2116 } else {
2117 // too small
2118 anyUpdates = 2;
2119 }
2120 } else {
2121 // zero dj
2122 anyUpdates = -1;
2123 }
2124 } else if (pivotSequence_ >= 0) {
2125 // just after re-factorization
2126 anyUpdates = -1;
2127 } else {
2128 // sub flip - nothing to do
2129 anyUpdates = 0;
2130 }
2131
2132 if (anyUpdates > 0) {
2133 model_->factorization()->updateColumnTranspose(spareRow2, updates);
2134
2135 // put row of tableau in rowArray and columnArray
2136 model_->clpMatrix()->transposeTimes(model_, -1.0,
2137 updates, spareColumn2, spareColumn1);
2138 // normal
2139 for (iSection = 0; iSection < 2; iSection++) {
2140
2141 reducedCost = model_->djRegion(iSection);
2142 int addSequence;
2143
2144 if (!iSection) {
2145 number = updates->getNumElements();
2146 index = updates->getIndices();
2147 updateBy = updates->denseVector();
2148 addSequence = model_->numberColumns();
2149 } else {
2150 number = spareColumn1->getNumElements();
2151 index = spareColumn1->getIndices();
2152 updateBy = spareColumn1->denseVector();
2153 addSequence = 0;
2154 }
2155 if (!model_->nonLinearCost()->lookBothWays()) {
2156
2157 for (j = 0; j < number; j++) {
2158 int iSequence = index[j];
2159 double value = reducedCost[iSequence];
2160 value -= updateBy[iSequence];
2161 reducedCost[iSequence] = value;
2162 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
2163
2164 switch(status) {
2165
2166 case ClpSimplex::basic:
2167 infeasible_->zero(iSequence + addSequence);
2168 case ClpSimplex::isFixed:
2169 break;
2170 case ClpSimplex::isFree:
2171 case ClpSimplex::superBasic:
2172 if (fabs(value) > FREE_ACCEPT * tolerance) {
2173 // we are going to bias towards free (but only if reasonable)
2174 value *= FREE_BIAS;
2175 // store square in list
2176 if (infeas[iSequence+addSequence])
2177 infeas[iSequence+addSequence] = value * value; // already there
2178 else
2179 infeasible_->quickAdd(iSequence + addSequence, value * value);
2180 } else {
2181 infeasible_->zero(iSequence + addSequence);
2182 }
2183 break;
2184 case ClpSimplex::atUpperBound:
2185 if (value > tolerance) {
2186 // store square in list
2187 if (infeas[iSequence+addSequence])
2188 infeas[iSequence+addSequence] = value * value; // already there
2189 else
2190 infeasible_->quickAdd(iSequence + addSequence, value * value);
2191 } else {
2192 infeasible_->zero(iSequence + addSequence);
2193 }
2194 break;
2195 case ClpSimplex::atLowerBound:
2196 if (value < -tolerance) {
2197 // store square in list
2198 if (infeas[iSequence+addSequence])
2199 infeas[iSequence+addSequence] = value * value; // already there
2200 else
2201 infeasible_->quickAdd(iSequence + addSequence, value * value);
2202 } else {
2203 infeasible_->zero(iSequence + addSequence);
2204 }
2205 }
2206 }
2207 } else {
2208 ClpNonLinearCost * nonLinear = model_->nonLinearCost();
2209 // We can go up OR down
2210 for (j = 0; j < number; j++) {
2211 int iSequence = index[j];
2212 double value = reducedCost[iSequence];
2213 value -= updateBy[iSequence];
2214 reducedCost[iSequence] = value;
2215 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
2216
2217 switch(status) {
2218
2219 case ClpSimplex::basic:
2220 infeasible_->zero(iSequence + addSequence);
2221 case ClpSimplex::isFixed:
2222 break;
2223 case ClpSimplex::isFree:
2224 case ClpSimplex::superBasic:
2225 if (fabs(value) > FREE_ACCEPT * tolerance) {
2226 // we are going to bias towards free (but only if reasonable)
2227 value *= FREE_BIAS;
2228 // store square in list
2229 if (infeas[iSequence+addSequence])
2230 infeas[iSequence+addSequence] = value * value; // already there
2231 else
2232 infeasible_->quickAdd(iSequence + addSequence, value * value);
2233 } else {
2234 infeasible_->zero(iSequence + addSequence);
2235 }
2236 break;
2237 case ClpSimplex::atUpperBound:
2238 if (value > tolerance) {
2239 // store square in list
2240 if (infeas[iSequence+addSequence])
2241 infeas[iSequence+addSequence] = value * value; // already there
2242 else
2243 infeasible_->quickAdd(iSequence + addSequence, value * value);
2244 } else {
2245 // look other way - change up should be negative
2246 value -= nonLinear->changeUpInCost(iSequence + addSequence);
2247 if (value > -tolerance) {
2248 infeasible_->zero(iSequence + addSequence);
2249 } else {
2250 // store square in list
2251 if (infeas[iSequence+addSequence])
2252 infeas[iSequence+addSequence] = value * value; // already there
2253 else
2254 infeasible_->quickAdd(iSequence + addSequence, value * value);
2255 }
2256 }
2257 break;
2258 case ClpSimplex::atLowerBound:
2259 if (value < -tolerance) {
2260 // store square in list
2261 if (infeas[iSequence+addSequence])
2262 infeas[iSequence+addSequence] = value * value; // already there
2263 else
2264 infeasible_->quickAdd(iSequence + addSequence, value * value);
2265 } else {
2266 // look other way - change down should be positive
2267 value -= nonLinear->changeDownInCost(iSequence + addSequence);
2268 if (value < tolerance) {
2269 infeasible_->zero(iSequence + addSequence);
2270 } else {
2271 // store square in list
2272 if (infeas[iSequence+addSequence])
2273 infeas[iSequence+addSequence] = value * value; // already there
2274 else
2275 infeasible_->quickAdd(iSequence + addSequence, value * value);
2276 }
2277 }
2278 }
2279 }
2280 }
2281 }
2282 if (anyUpdates == 2) {
2283 // we can zero out as will have to get pivot row
2284 updates->clear();
2285 spareColumn1->clear();
2286 }
2287 if (pivotRow >= 0) {
2288 // make sure infeasibility on incoming is 0.0
2289 int sequenceIn = model_->sequenceIn();
2290 infeasible_->zero(sequenceIn);
2291 }
2292 }
2293 // make sure outgoing from last iteration okay
2294 int sequenceOut = model_->sequenceOut();
2295 if (sequenceOut >= 0) {
2296 ClpSimplex::Status status = model_->getStatus(sequenceOut);
2297 double value = model_->reducedCost(sequenceOut);
2298
2299 switch(status) {
2300
2301 case ClpSimplex::basic:
2302 case ClpSimplex::isFixed:
2303 break;
2304 case ClpSimplex::isFree:
2305 case ClpSimplex::superBasic:
2306 if (fabs(value) > FREE_ACCEPT * tolerance) {
2307 // we are going to bias towards free (but only if reasonable)
2308 value *= FREE_BIAS;
2309 // store square in list
2310 if (infeas[sequenceOut])
2311 infeas[sequenceOut] = value * value; // already there
2312 else
2313 infeasible_->quickAdd(sequenceOut, value * value);
2314 } else {
2315 infeasible_->zero(sequenceOut);
2316 }
2317 break;
2318 case ClpSimplex::atUpperBound:
2319 if (value > tolerance) {
2320 // store square in list
2321 if (infeas[sequenceOut])
2322 infeas[sequenceOut] = value * value; // already there
2323 else
2324 infeasible_->quickAdd(sequenceOut, value * value);
2325 } else {
2326 infeasible_->zero(sequenceOut);
2327 }
2328 break;
2329 case ClpSimplex::atLowerBound:
2330 if (value < -tolerance) {
2331 // store square in list
2332 if (infeas[sequenceOut])
2333 infeas[sequenceOut] = value * value; // already there
2334 else
2335 infeasible_->quickAdd(sequenceOut, value * value);
2336 } else {
2337 infeasible_->zero(sequenceOut);
2338 }
2339 }
2340 }
2341
2342 // If in quadratic re-compute all
2343 if (model_->algorithm() == 2) {
2344 for (iSection = 0; iSection < 2; iSection++) {
2345
2346 reducedCost = model_->djRegion(iSection);
2347 int addSequence;
2348 int iSequence;
2349
2350 if (!iSection) {
2351 number = model_->numberRows();
2352 addSequence = model_->numberColumns();
2353 } else {
2354 number = model_->numberColumns();
2355 addSequence = 0;
2356 }
2357
2358 if (!model_->nonLinearCost()->lookBothWays()) {
2359 for (iSequence = 0; iSequence < number; iSequence++) {
2360 double value = reducedCost[iSequence];
2361 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
2362
2363 switch(status) {
2364
2365 case ClpSimplex::basic:
2366 infeasible_->zero(iSequence + addSequence);
2367 case ClpSimplex::isFixed:
2368 break;
2369 case ClpSimplex::isFree:
2370 case ClpSimplex::superBasic:
2371 if (fabs(value) > tolerance) {
2372 // we are going to bias towards free (but only if reasonable)
2373 value *= FREE_BIAS;
2374 // store square in list
2375 if (infeas[iSequence+addSequence])
2376 infeas[iSequence+addSequence] = value * value; // already there
2377 else
2378 infeasible_->quickAdd(iSequence + addSequence, value * value);
2379 } else {
2380 infeasible_->zero(iSequence + addSequence);
2381 }
2382 break;
2383 case ClpSimplex::atUpperBound:
2384 if (value > tolerance) {
2385 // store square in list
2386 if (infeas[iSequence+addSequence])
2387 infeas[iSequence+addSequence] = value * value; // already there
2388 else
2389 infeasible_->quickAdd(iSequence + addSequence, value * value);
2390 } else {
2391 infeasible_->zero(iSequence + addSequence);
2392 }
2393 break;
2394 case ClpSimplex::atLowerBound:
2395 if (value < -tolerance) {
2396 // store square in list
2397 if (infeas[iSequence+addSequence])
2398 infeas[iSequence+addSequence] = value * value; // already there
2399 else
2400 infeasible_->quickAdd(iSequence + addSequence, value * value);
2401 } else {
2402 infeasible_->zero(iSequence + addSequence);
2403 }
2404 }
2405 }
2406 } else {
2407 // we can go both ways
2408 ClpNonLinearCost * nonLinear = model_->nonLinearCost();
2409 for (iSequence = 0; iSequence < number; iSequence++) {
2410 double value = reducedCost[iSequence];
2411 ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
2412
2413 switch(status) {
2414
2415 case ClpSimplex::basic:
2416 infeasible_->zero(iSequence + addSequence);
2417 case ClpSimplex::isFixed:
2418 break;
2419 case ClpSimplex::isFree:
2420 case ClpSimplex::superBasic:
2421 if (fabs(value) > tolerance) {
2422 // we are going to bias towards free (but only if reasonable)
2423 value *= FREE_BIAS;
2424 // store square in list
2425 if (infeas[iSequence+addSequence])
2426 infeas[iSequence+addSequence] = value * value; // already there
2427 else
2428 infeasible_->quickAdd(iSequence + addSequence, value * value);
2429 } else {
2430 infeasible_->zero(iSequence + addSequence);
2431 }
2432 break;
2433 case ClpSimplex::atUpperBound:
2434 if (value > tolerance) {
2435 // store square in list
2436 if (infeas[iSequence+addSequence])
2437 infeas[iSequence+addSequence] = value * value; // already there
2438 else
2439 infeasible_->quickAdd(iSequence + addSequence, value * value);
2440 } else {
2441 // look other way - change up should be negative
2442 value -= nonLinear->changeUpInCost(iSequence + addSequence);
2443 if (value > -tolerance) {
2444 infeasible_->zero(iSequence + addSequence);
2445 } else {
2446 // store square in list
2447 if (infeas[iSequence+addSequence])
2448 infeas[iSequence+addSequence] = value * value; // already there
2449 else
2450 infeasible_->quickAdd(iSequence + addSequence, value * value);
2451 }
2452 }
2453 break;
2454 case ClpSimplex::atLowerBound:
2455 if (value < -tolerance) {
2456 // store square in list
2457 if (infeas[iSequence+addSequence])
2458 infeas[iSequence+addSequence] = value * value; // already there
2459 else
2460 infeasible_->quickAdd(iSequence + addSequence, value * value);
2461 } else {
2462 // look other way - change down should be positive
2463 value -= nonLinear->changeDownInCost(iSequence + addSequence);
2464 if (value < tolerance) {
2465 infeasible_->zero(iSequence + addSequence);
2466 } else {
2467 // store square in list
2468 if (infeas[iSequence+addSequence])
2469 infeas[iSequence+addSequence] = value * value; // already there
2470 else
2471 infeasible_->quickAdd(iSequence + addSequence, value * value);
2472 }
2473 }
2474 }
2475 }
2476 }
2477 }
2478 }
2479 // See what sort of pricing
2480 int numberWanted = 10;
2481 number = infeasible_->getNumElements();
2482 int numberColumns = model_->numberColumns();
2483 if (switchType == 5) {
2484 // we can zero out
2485 updates->clear();
2486 spareColumn1->clear();
2487 pivotSequence_ = -1;
2488 pivotRow = -1;
2489 // See if to switch
2490 int numberRows = model_->numberRows();
2491 // ratio is done on number of columns here
2492 //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
2493 double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
2494 //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
2495 if (ratio < 0.1) {
2496 numberWanted = CoinMax(100, number / 200);
2497 } else if (ratio < 0.3) {
2498 numberWanted = CoinMax(500, number / 40);
2499 } else if (ratio < 0.5 || mode_ == 5) {
2500 numberWanted = CoinMax(2000, number / 10);
2501 numberWanted = CoinMax(numberWanted, numberColumns / 30);
2502 } else if (mode_ != 5) {
2503 switchType = 4;
2504 // initialize
2505 numberSwitched_++;
2506 // Make sure will re-do
2507 delete [] weights_;
2508 weights_ = NULL;
2509 saveWeights(model_, 4);
2510 COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
2511 updates->clear();
2512 }
2513 if (model_->numberIterations() % 1000 == 0)
2514 COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted));
2515 }
2516 if(switchType == 4) {
2517 // Still in devex mode
2518 int numberRows = model_->numberRows();
2519 // ratio is done on number of rows here
2520 double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
2521 // Go to steepest if lot of iterations?
2522 if (ratio < 1.0) {
2523 numberWanted = CoinMax(2000, number / 20);
2524 } else if (ratio < 5.0) {
2525 numberWanted = CoinMax(2000, number / 10);
2526 numberWanted = CoinMax(numberWanted, numberColumns / 20);
2527 } else {
2528 // we can zero out
2529 updates->clear();
2530 spareColumn1->clear();
2531 switchType = 3;
2532 // initialize
2533 pivotSequence_ = -1;
2534 pivotRow = -1;
2535 numberSwitched_++;
2536 // Make sure will re-do
2537 delete [] weights_;
2538 weights_ = NULL;
2539 saveWeights(model_, 4);
2540 COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
2541 updates->clear();
2542 }
2543 if (model_->numberIterations() % 1000 == 0)
2544 COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted));
2545 }
2546 if (switchType < 4) {
2547 if (switchType < 2 ) {
2548 numberWanted = number + 1;
2549 } else if (switchType == 2) {
2550 numberWanted = CoinMax(2000, number / 8);
2551 } else {
2552 double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (model_->numberRows());
2553 if (ratio < 1.0) {
2554 numberWanted = CoinMax(2000, number / 20);
2555 } else if (ratio < 5.0) {
2556 numberWanted = CoinMax(2000, number / 10);
2557 numberWanted = CoinMax(numberWanted, numberColumns / 20);
2558 } else if (ratio < 10.0) {
2559 numberWanted = CoinMax(2000, number / 8);
2560 numberWanted = CoinMax(numberWanted, numberColumns / 20);
2561 } else {
2562 ratio = number * (ratio / 80.0);
2563 if (ratio > number) {
2564 numberWanted = number + 1;
2565 } else {
2566 numberWanted = CoinMax(2000, static_cast<int> (ratio));
2567 numberWanted = CoinMax(numberWanted, numberColumns / 10);
2568 }
2569 }
2570 }
2571 }
2572 // for weights update we use pivotSequence
2573 pivotRow = pivotSequence_;
2574 // unset in case sub flip
2575 pivotSequence_ = -1;
2576 if (pivotRow >= 0) {
2577 // make sure infeasibility on incoming is 0.0
2578 const int * pivotVariable = model_->pivotVariable();
2579 int sequenceIn = pivotVariable[pivotRow];
2580 infeasible_->zero(sequenceIn);
2581 // and we can see if reference
2582 double referenceIn = 0.0;
2583 if (switchType != 1 && reference(sequenceIn))
2584 referenceIn = 1.0;
2585 // save outgoing weight round update
2586 double outgoingWeight = 0.0;
2587 if (sequenceOut >= 0)
2588 outgoingWeight = weights_[sequenceOut];
2589 // update weights
2590 if (anyUpdates != 1) {
2591 updates->setNumElements(0);
2592 spareColumn1->setNumElements(0);
2593 // might as well set dj to 1
2594 dj = 1.0;
2595 updates->insert(pivotRow, -dj);
2596 model_->factorization()->updateColumnTranspose(spareRow2, updates);
2597 // put row of tableau in rowArray and columnArray
2598 model_->clpMatrix()->transposeTimes(model_, -1.0,
2599 updates, spareColumn2, spareColumn1);
2600 }
2601 double * weight;
2602 double * other = alternateWeights_->denseVector();
2603 int numberColumns = model_->numberColumns();
2604 double scaleFactor = -1.0 / dj; // as formula is with 1.0
2605 // rows
2606 number = updates->getNumElements();
2607 index = updates->getIndices();
2608 updateBy = updates->denseVector();
2609 weight = weights_ + numberColumns;
2610
2611 if (switchType < 4) {
2612 // Exact
2613 // now update weight update array
2614 model_->factorization()->updateColumnTranspose(spareRow2,
2615 alternateWeights_);
2616 for (j = 0; j < number; j++) {
2617 int iSequence = index[j];
2618 double thisWeight = weight[iSequence];
2619 // row has -1
2620 double pivot = updateBy[iSequence] * scaleFactor;
2621 updateBy[iSequence] = 0.0;
2622 double modification = other[iSequence];
2623 double pivotSquared = pivot * pivot;
2624
2625 thisWeight += pivotSquared * devex_ + pivot * modification;
2626 if (thisWeight < TRY_NORM) {
2627 if (switchType == 1) {
2628 // steepest
2629 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
2630 } else {
2631 // exact
2632 thisWeight = referenceIn * pivotSquared;
2633 if (reference(iSequence + numberColumns))
2634 thisWeight += 1.0;
2635 thisWeight = CoinMax(thisWeight, TRY_NORM);
2636 }
2637 }
2638 weight[iSequence] = thisWeight;
2639 }
2640 } else if (switchType == 4) {
2641 // Devex
2642 assert (devex_ > 0.0);
2643 for (j = 0; j < number; j++) {
2644 int iSequence = index[j];
2645 double thisWeight = weight[iSequence];
2646 // row has -1
2647 double pivot = updateBy[iSequence] * scaleFactor;
2648 updateBy[iSequence] = 0.0;
2649 double value = pivot * pivot * devex_;
2650 if (reference(iSequence + numberColumns))
2651 value += 1.0;
2652 weight[iSequence] = CoinMax(0.99 * thisWeight, value);
2653 }
2654 }
2655
2656 // columns
2657 weight = weights_;
2658
2659 scaleFactor = -scaleFactor;
2660
2661 number = spareColumn1->getNumElements();
2662 index = spareColumn1->getIndices();
2663 updateBy = spareColumn1->denseVector();
2664 if (switchType < 4) {
2665 // Exact
2666 // get subset which have nonzero tableau elements
2667 model_->clpMatrix()->subsetTransposeTimes(model_, alternateWeights_,
2668 spareColumn1,
2669 spareColumn2);
2670 double * updateBy2 = spareColumn2->denseVector();
2671 for (j = 0; j < number; j++) {
2672 int iSequence = index[j];
2673 double thisWeight = weight[iSequence];
2674 double pivot = updateBy[iSequence] * scaleFactor;
2675 updateBy[iSequence] = 0.0;
2676 double modification = updateBy2[j];
2677 updateBy2[j] = 0.0;
2678 double pivotSquared = pivot * pivot;
2679
2680 thisWeight += pivotSquared * devex_ + pivot * modification;
2681 if (thisWeight < TRY_NORM) {
2682 if (switchType == 1) {
2683 // steepest
2684 thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
2685 } else {
2686 // exact
2687 thisWeight = referenceIn * pivotSquared;
2688 if (reference(iSequence))
2689 thisWeight += 1.0;
2690 thisWeight = CoinMax(thisWeight, TRY_NORM);
2691 }
2692 }
2693 weight[iSequence] = thisWeight;
2694 }
2695 } else if (switchType == 4) {
2696 // Devex
2697 for (j = 0; j < number; j++) {
2698 int iSequence = index[j];
2699 double thisWeight = weight[iSequence];
2700 // row has -1
2701 double pivot = updateBy[iSequence] * scaleFactor;
2702 updateBy[iSequence] = 0.0;
2703 double value = pivot * pivot * devex_;
2704 if (reference(iSequence))
2705 value += 1.0;
2706 weight[iSequence] = CoinMax(0.99 * thisWeight, value);
2707 }
2708 }
2709 // restore outgoing weight
2710 if (sequenceOut >= 0)
2711 weights_[sequenceOut] = outgoingWeight;
2712 alternateWeights_->clear();
2713 spareColumn2->setNumElements(0);
2714 //#define SOME_DEBUG_1
2715#ifdef SOME_DEBUG_1
2716 // check for accuracy
2717 int iCheck = 892;
2718 //printf("weight for iCheck is %g\n",weights_[iCheck]);
2719 int numberRows = model_->numberRows();
2720 //int numberColumns = model_->numberColumns();
2721 for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
2722 if (model_->getStatus(iCheck) != ClpSimplex::basic &&
2723 !model_->getStatus(iCheck) != ClpSimplex::isFixed)
2724 checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
2725 }
2726#endif
2727 updates->setNumElements(0);
2728 spareColumn1->setNumElements(0);
2729 }
2730
2731 // update of duals finished - now do pricing
2732
2733
2734 double bestDj = 1.0e-30;
2735 int bestSequence = -1;
2736
2737 int i, iSequence;
2738 index = infeasible_->getIndices();
2739 number = infeasible_->getNumElements();
2740 if(model_->numberIterations() < model_->lastBadIteration() + 200) {
2741 // we can't really trust infeasibilities if there is dual error
2742 double checkTolerance = 1.0e-8;
2743 if (!model_->factorization()->pivots())
2744 checkTolerance = 1.0e-6;
2745 if (model_->largestDualError() > checkTolerance)
2746 tolerance *= model_->largestDualError() / checkTolerance;
2747 // But cap
2748 tolerance = CoinMin(1000.0, tolerance);
2749 }
2750#ifdef CLP_DEBUG
2751 if (model_->numberDualInfeasibilities() == 1)
2752 printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
2753 model_->largestDualError(), model_, model_->messageHandler(),
2754 number);
2755#endif
2756 // stop last one coming immediately
2757 double saveOutInfeasibility = 0.0;
2758 if (sequenceOut >= 0) {
2759 saveOutInfeasibility = infeas[sequenceOut];
2760 infeas[sequenceOut] = 0.0;
2761 }
2762 tolerance *= tolerance; // as we are using squares
2763
2764 int iPass;
2765 // Setup two passes
2766 int start[4];
2767 start[1] = number;
2768 start[2] = 0;
2769 double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
2770 start[0] = static_cast<int> (dstart);
2771 start[3] = start[0];
2772 //double largestWeight=0.0;
2773 //double smallestWeight=1.0e100;
2774 for (iPass = 0; iPass < 2; iPass++) {
2775 int end = start[2*iPass+1];
2776 if (switchType < 5) {
2777 for (i = start[2*iPass]; i < end; i++) {
2778 iSequence = index[i];
2779 double value = infeas[iSequence];
2780 if (value > tolerance) {
2781 double weight = weights_[iSequence];
2782 //weight=1.0;
2783 if (value > bestDj * weight) {
2784 // check flagged variable and correct dj
2785 if (!model_->flagged(iSequence)) {
2786 bestDj = value / weight;
2787 bestSequence = iSequence;
2788 } else {
2789 // just to make sure we don't exit before got something
2790 numberWanted++;
2791 }
2792 }
2793 }
2794 numberWanted--;
2795 if (!numberWanted)
2796 break;
2797 }
2798 } else {
2799 // Dantzig
2800 for (i = start[2*iPass]; i < end; i++) {
2801 iSequence = index[i];
2802 double value = infeas[iSequence];
2803 if (value > tolerance) {
2804 if (value > bestDj) {
2805 // check flagged variable and correct dj
2806 if (!model_->flagged(iSequence)) {
2807 bestDj = value;
2808 bestSequence = iSequence;
2809 } else {
2810 // just to make sure we don't exit before got something
2811 numberWanted++;
2812 }
2813 }
2814 }
2815 numberWanted--;
2816 if (!numberWanted)
2817 break;
2818 }
2819 }
2820 if (!numberWanted)
2821 break;
2822 }
2823 if (sequenceOut >= 0) {
2824 infeas[sequenceOut] = saveOutInfeasibility;
2825 }
2826 /*if (model_->numberIterations()%100==0)
2827 printf("%d best %g\n",bestSequence,bestDj);*/
2828 reducedCost = model_->djRegion();
2829 model_->clpMatrix()->setSavedBestSequence(bestSequence);
2830 if (bestSequence >= 0)
2831 model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
2832
2833#ifdef CLP_DEBUG
2834 if (bestSequence >= 0) {
2835 if (model_->getStatus(bestSequence) == ClpSimplex::atLowerBound)
2836 assert(model_->reducedCost(bestSequence) < 0.0);
2837 if (model_->getStatus(bestSequence) == ClpSimplex::atUpperBound)
2838 assert(model_->reducedCost(bestSequence) > 0.0);
2839 }
2840#endif
2841 return bestSequence;
2842}
2843// Called when maximum pivots changes
2844void
2845ClpPrimalColumnSteepest::maximumPivotsChanged()
2846{
2847 if (alternateWeights_ &&
2848 alternateWeights_->capacity() != model_->numberRows() +
2849 model_->factorization()->maximumPivots()) {
2850 delete alternateWeights_;
2851 alternateWeights_ = new CoinIndexedVector();
2852 // enough space so can use it for factorization
2853 alternateWeights_->reserve(model_->numberRows() +
2854 model_->factorization()->maximumPivots());
2855 }
2856}
2857/*
2858 1) before factorization
2859 2) after factorization
2860 3) just redo infeasibilities
2861 4) restore weights
2862 5) at end of values pass (so need initialization)
2863*/
2864void
2865ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model, int mode)
2866{
2867 model_ = model;
2868 if (mode_ == 4 || mode_ == 5) {
2869 if (mode == 1 && !weights_)
2870 numberSwitched_ = 0; // Reset
2871 }
2872 // alternateWeights_ is defined as indexed but is treated oddly
2873 // at times
2874 int numberRows = model_->numberRows();
2875 int numberColumns = model_->numberColumns();
2876 const int * pivotVariable = model_->pivotVariable();
2877 bool doInfeasibilities = true;
2878 if (mode == 1) {
2879 if(weights_) {
2880 // Check if size has changed
2881 if (infeasible_->capacity() == numberRows + numberColumns &&
2882 alternateWeights_->capacity() == numberRows +
2883 model_->factorization()->maximumPivots()) {
2884 //alternateWeights_->clear();
2885 if (pivotSequence_ >= 0 && pivotSequence_ < numberRows) {
2886 // save pivot order
2887 CoinMemcpyN(pivotVariable,
2888 numberRows, alternateWeights_->getIndices());
2889 // change from pivot row number to sequence number
2890 pivotSequence_ = pivotVariable[pivotSequence_];
2891 } else {
2892 pivotSequence_ = -1;
2893 }
2894 state_ = 1;
2895 } else {
2896 // size has changed - clear everything
2897 delete [] weights_;
2898 weights_ = NULL;
2899 delete infeasible_;
2900 infeasible_ = NULL;
2901 delete alternateWeights_;
2902 alternateWeights_ = NULL;
2903 delete [] savedWeights_;
2904 savedWeights_ = NULL;
2905 delete [] reference_;
2906 reference_ = NULL;
2907 state_ = -1;
2908 pivotSequence_ = -1;
2909 }
2910 }
2911 } else if (mode == 2 || mode == 4 || mode == 5) {
2912 // restore
2913 if (!weights_ || state_ == -1 || mode == 5) {
2914 // Partial is only allowed with certain types of matrix
2915 if ((mode_ != 4 && mode_ != 5) || numberSwitched_ || !model_->clpMatrix()->canDoPartialPricing()) {
2916 // initialize weights
2917 delete [] weights_;
2918 delete alternateWeights_;
2919 weights_ = new double[numberRows+numberColumns];
2920 alternateWeights_ = new CoinIndexedVector();
2921 // enough space so can use it for factorization
2922 alternateWeights_->reserve(numberRows +
2923 model_->factorization()->maximumPivots());
2924 initializeWeights();
2925 // create saved weights
2926 delete [] savedWeights_;
2927 savedWeights_ = CoinCopyOfArray(weights_, numberRows + numberColumns);
2928 // just do initialization
2929 mode = 3;
2930 } else {
2931 // Partial pricing
2932 // use region as somewhere to save non-fixed slacks
2933 // set up infeasibilities
2934 if (!infeasible_) {
2935 infeasible_ = new CoinIndexedVector();
2936 infeasible_->reserve(numberColumns + numberRows);
2937 }
2938 infeasible_->clear();
2939 int number = model_->numberRows() + model_->numberColumns();
2940 int iSequence;
2941 int numberLook = 0;
2942 int * which = infeasible_->getIndices();
2943 for (iSequence = model_->numberColumns(); iSequence < number; iSequence++) {
2944 ClpSimplex::Status status = model_->getStatus(iSequence);
2945 if (status != ClpSimplex::isFixed)
2946 which[numberLook++] = iSequence;
2947 }
2948 infeasible_->setNumElements(numberLook);
2949 doInfeasibilities = false;
2950 }
2951 savedPivotSequence_ = -2;
2952 savedSequenceOut_ = -2;
2953
2954 } else {
2955 if (mode != 4) {
2956 // save
2957 CoinMemcpyN(weights_, (numberRows + numberColumns), savedWeights_);
2958 savedPivotSequence_ = pivotSequence_;
2959 savedSequenceOut_ = model_->sequenceOut();
2960 } else {
2961 // restore
2962 CoinMemcpyN(savedWeights_, (numberRows + numberColumns), weights_);
2963 // was - but I think should not be
2964 //pivotSequence_= savedPivotSequence_;
2965 //model_->setSequenceOut(savedSequenceOut_);
2966 pivotSequence_ = -1;
2967 model_->setSequenceOut(-1);
2968 // indices are wrong so clear by hand
2969 //alternateWeights_->clear();
2970 CoinZeroN(alternateWeights_->denseVector(),
2971 alternateWeights_->capacity());
2972 alternateWeights_->setNumElements(0);
2973 }
2974 }
2975 state_ = 0;
2976 // set up infeasibilities
2977 if (!infeasible_) {
2978 infeasible_ = new CoinIndexedVector();
2979 infeasible_->reserve(numberColumns + numberRows);
2980 }
2981 }
2982 if (mode >= 2 && mode != 5) {
2983 if (mode != 3) {
2984 if (pivotSequence_ >= 0) {
2985 // restore pivot row
2986 int iRow;
2987 // permute alternateWeights
2988 double * temp = model_->rowArray(3)->denseVector();
2989 double * work = alternateWeights_->denseVector();
2990 int * savePivotOrder = model_->rowArray(3)->getIndices();
2991 int * oldPivotOrder = alternateWeights_->getIndices();
2992 for (iRow = 0; iRow < numberRows; iRow++) {
2993 int iPivot = oldPivotOrder[iRow];
2994 temp[iPivot] = work[iRow];
2995 savePivotOrder[iRow] = iPivot;
2996 }
2997 int number = 0;
2998 int found = -1;
2999 int * which = oldPivotOrder;
3000 // find pivot row and re-create alternateWeights
3001 for (iRow = 0; iRow < numberRows; iRow++) {
3002 int iPivot = pivotVariable[iRow];
3003 if (iPivot == pivotSequence_)
3004 found = iRow;
3005 work[iRow] = temp[iPivot];
3006 if (work[iRow])
3007 which[number++] = iRow;
3008 }
3009 alternateWeights_->setNumElements(number);
3010#ifdef CLP_DEBUG
3011 // Can happen but I should clean up
3012 assert(found >= 0);
3013#endif
3014 pivotSequence_ = found;
3015 for (iRow = 0; iRow < numberRows; iRow++) {
3016 int iPivot = savePivotOrder[iRow];
3017 temp[iPivot] = 0.0;
3018 }
3019 } else {
3020 // Just clean up
3021 if (alternateWeights_)
3022 alternateWeights_->clear();
3023 }
3024 }
3025 // Save size of factorization
3026 if (!model->factorization()->pivots())
3027 sizeFactorization_ = model_->factorization()->numberElements();
3028 if(!doInfeasibilities)
3029 return; // don't disturb infeasibilities
3030 infeasible_->clear();
3031 double tolerance = model_->currentDualTolerance();
3032 int number = model_->numberRows() + model_->numberColumns();
3033 int iSequence;
3034
3035 double * reducedCost = model_->djRegion();
3036
3037 if (!model_->nonLinearCost()->lookBothWays()) {
3038#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
3039 for (iSequence = 0; iSequence < number; iSequence++) {
3040 double value = reducedCost[iSequence];
3041 ClpSimplex::Status status = model_->getStatus(iSequence);
3042
3043 switch(status) {
3044
3045 case ClpSimplex::basic:
3046 case ClpSimplex::isFixed:
3047 break;
3048 case ClpSimplex::isFree:
3049 case ClpSimplex::superBasic:
3050 if (fabs(value) > FREE_ACCEPT * tolerance) {
3051 // we are going to bias towards free (but only if reasonable)
3052 value *= FREE_BIAS;
3053 // store square in list
3054 infeasible_->quickAdd(iSequence, value * value);
3055 }
3056 break;
3057 case ClpSimplex::atUpperBound:
3058 if (value > tolerance) {
3059 infeasible_->quickAdd(iSequence, value * value);
3060 }
3061 break;
3062 case ClpSimplex::atLowerBound:
3063 if (value < -tolerance) {
3064 infeasible_->quickAdd(iSequence, value * value);
3065 }
3066 }
3067 }
3068#else
3069 // Columns
3070 int numberColumns = model_->numberColumns();
3071 for (iSequence = 0; iSequence < numberColumns; iSequence++) {
3072 double value = reducedCost[iSequence];
3073 ClpSimplex::Status status = model_->getStatus(iSequence);
3074
3075 switch(status) {
3076
3077 case ClpSimplex::basic:
3078 case ClpSimplex::isFixed:
3079 break;
3080 case ClpSimplex::isFree:
3081 case ClpSimplex::superBasic:
3082 if (fabs(value) > FREE_ACCEPT * tolerance) {
3083 // we are going to bias towards free (but only if reasonable)
3084 value *= FREE_BIAS;
3085 // store square in list
3086 infeasible_->quickAdd(iSequence, value * value);
3087 }
3088 break;
3089 case ClpSimplex::atUpperBound:
3090 if (value > tolerance) {
3091 infeasible_->quickAdd(iSequence, value * value);
3092 }
3093 break;
3094 case ClpSimplex::atLowerBound:
3095 if (value < -tolerance) {
3096 infeasible_->quickAdd(iSequence, value * value);
3097 }
3098 }
3099 }
3100 // Rows
3101 for ( ; iSequence < number; iSequence++) {
3102 double value = reducedCost[iSequence];
3103 ClpSimplex::Status status = model_->getStatus(iSequence);
3104
3105 switch(status) {
3106
3107 case ClpSimplex::basic:
3108 case ClpSimplex::isFixed:
3109 break;
3110 case ClpSimplex::isFree:
3111 case ClpSimplex::superBasic:
3112 if (fabs(value) > FREE_ACCEPT * tolerance) {
3113 // we are going to bias towards free (but only if reasonable)
3114 value *= FREE_BIAS;
3115 // store square in list
3116 infeasible_->quickAdd(iSequence, value * value);
3117 }
3118 break;
3119 case ClpSimplex::atUpperBound:
3120 if (value > tolerance) {
3121 infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
3122 }
3123 break;
3124 case ClpSimplex::atLowerBound:
3125 if (value < -tolerance) {
3126 infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
3127 }
3128 }
3129 }
3130#endif
3131 } else {
3132 ClpNonLinearCost * nonLinear = model_->nonLinearCost();
3133 // can go both ways
3134 for (iSequence = 0; iSequence < number; iSequence++) {
3135 double value = reducedCost[iSequence];
3136 ClpSimplex::Status status = model_->getStatus(iSequence);
3137
3138 switch(status) {
3139
3140 case ClpSimplex::basic:
3141 case ClpSimplex::isFixed:
3142 break;
3143 case ClpSimplex::isFree:
3144 case ClpSimplex::superBasic:
3145 if (fabs(value) > FREE_ACCEPT * tolerance) {
3146 // we are going to bias towards free (but only if reasonable)
3147 value *= FREE_BIAS;
3148 // store square in list
3149 infeasible_->quickAdd(iSequence, value * value);
3150 }
3151 break;
3152 case ClpSimplex::atUpperBound:
3153 if (value > tolerance) {
3154 infeasible_->quickAdd(iSequence, value * value);
3155 } else {
3156 // look other way - change up should be negative
3157 value -= nonLinear->changeUpInCost(iSequence);
3158 if (value < -tolerance) {
3159 // store square in list
3160 infeasible_->quickAdd(iSequence, value * value);
3161 }
3162 }
3163 break;
3164 case ClpSimplex::atLowerBound:
3165 if (value < -tolerance) {
3166 infeasible_->quickAdd(iSequence, value * value);
3167 } else {
3168 // look other way - change down should be positive
3169 value -= nonLinear->changeDownInCost(iSequence);
3170 if (value > tolerance) {
3171 // store square in list
3172 infeasible_->quickAdd(iSequence, value * value);
3173 }
3174 }
3175 }
3176 }
3177 }
3178 }
3179}
3180// Gets rid of last update
3181void
3182ClpPrimalColumnSteepest::unrollWeights()
3183{
3184 if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
3185 return;
3186 double * saved = alternateWeights_->denseVector();
3187 int number = alternateWeights_->getNumElements();
3188 int * which = alternateWeights_->getIndices();
3189 int i;
3190 for (i = 0; i < number; i++) {
3191 int iRow = which[i];
3192 weights_[iRow] = saved[iRow];
3193 saved[iRow] = 0.0;
3194 }
3195 alternateWeights_->setNumElements(0);
3196}
3197
3198//-------------------------------------------------------------------
3199// Clone
3200//-------------------------------------------------------------------
3201ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
3202{
3203 if (CopyData) {
3204 return new ClpPrimalColumnSteepest(*this);
3205 } else {
3206 return new ClpPrimalColumnSteepest();
3207 }
3208}
3209void
3210ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
3211{
3212 // Local copy of mode so can decide what to do
3213 int switchType = mode_;
3214 if (mode_ == 4 && numberSwitched_)
3215 switchType = 3;
3216 else if (mode_ == 4 || mode_ == 5)
3217 return;
3218 int number = input->getNumElements();
3219 int * which = input->getIndices();
3220 double * work = input->denseVector();
3221 int newNumber = 0;
3222 int * newWhich = alternateWeights_->getIndices();
3223 double * newWork = alternateWeights_->denseVector();
3224 int i;
3225 int sequenceIn = model_->sequenceIn();
3226 int sequenceOut = model_->sequenceOut();
3227 const int * pivotVariable = model_->pivotVariable();
3228
3229 int pivotRow = model_->pivotRow();
3230 pivotSequence_ = pivotRow;
3231
3232 devex_ = 0.0;
3233 // Can't create alternateWeights_ as packed as needed unpacked
3234 if (!input->packedMode()) {
3235 if (pivotRow >= 0) {
3236 if (switchType == 1) {
3237 for (i = 0; i < number; i++) {
3238 int iRow = which[i];
3239 devex_ += work[iRow] * work[iRow];
3240 newWork[iRow] = -2.0 * work[iRow];
3241 }
3242 newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3243 devex_ += ADD_ONE;
3244 weights_[sequenceOut] = 1.0 + ADD_ONE;
3245 CoinMemcpyN(which, number, newWhich);
3246 alternateWeights_->setNumElements(number);
3247 } else {
3248 if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
3249 for (i = 0; i < number; i++) {
3250 int iRow = which[i];
3251 int iPivot = pivotVariable[iRow];
3252 if (reference(iPivot)) {
3253 devex_ += work[iRow] * work[iRow];
3254 newWork[iRow] = -2.0 * work[iRow];
3255 newWhich[newNumber++] = iRow;
3256 }
3257 }
3258 if (!newWork[pivotRow] && devex_ > 0.0)
3259 newWhich[newNumber++] = pivotRow; // add if not already in
3260 newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3261 } else {
3262 for (i = 0; i < number; i++) {
3263 int iRow = which[i];
3264 int iPivot = pivotVariable[iRow];
3265 if (reference(iPivot))
3266 devex_ += work[iRow] * work[iRow];
3267 }
3268 }
3269 if (reference(sequenceIn)) {
3270 devex_ += 1.0;
3271 } else {
3272 }
3273 if (reference(sequenceOut)) {
3274 weights_[sequenceOut] = 1.0 + 1.0;
3275 } else {
3276 weights_[sequenceOut] = 1.0;
3277 }
3278 alternateWeights_->setNumElements(newNumber);
3279 }
3280 } else {
3281 if (switchType == 1) {
3282 for (i = 0; i < number; i++) {
3283 int iRow = which[i];
3284 devex_ += work[iRow] * work[iRow];
3285 }
3286 devex_ += ADD_ONE;
3287 } else {
3288 for (i = 0; i < number; i++) {
3289 int iRow = which[i];
3290 int iPivot = pivotVariable[iRow];
3291 if (reference(iPivot)) {
3292 devex_ += work[iRow] * work[iRow];
3293 }
3294 }
3295 if (reference(sequenceIn))
3296 devex_ += 1.0;
3297 }
3298 }
3299 } else {
3300 // packed input
3301 if (pivotRow >= 0) {
3302 if (switchType == 1) {
3303 for (i = 0; i < number; i++) {
3304 int iRow = which[i];
3305 devex_ += work[i] * work[i];
3306 newWork[iRow] = -2.0 * work[i];
3307 }
3308 newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3309 devex_ += ADD_ONE;
3310 weights_[sequenceOut] = 1.0 + ADD_ONE;
3311 CoinMemcpyN(which, number, newWhich);
3312 alternateWeights_->setNumElements(number);
3313 } else {
3314 if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
3315 for (i = 0; i < number; i++) {
3316 int iRow = which[i];
3317 int iPivot = pivotVariable[iRow];
3318 if (reference(iPivot)) {
3319 devex_ += work[i] * work[i];
3320 newWork[iRow] = -2.0 * work[i];
3321 newWhich[newNumber++] = iRow;
3322 }
3323 }
3324 if (!newWork[pivotRow] && devex_ > 0.0)
3325 newWhich[newNumber++] = pivotRow; // add if not already in
3326 newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3327 } else {
3328 for (i = 0; i < number; i++) {
3329 int iRow = which[i];
3330 int iPivot = pivotVariable[iRow];
3331 if (reference(iPivot))
3332 devex_ += work[i] * work[i];
3333 }
3334 }
3335 if (reference(sequenceIn)) {
3336 devex_ += 1.0;
3337 } else {
3338 }
3339 if (reference(sequenceOut)) {
3340 weights_[sequenceOut] = 1.0 + 1.0;
3341 } else {
3342 weights_[sequenceOut] = 1.0;
3343 }
3344 alternateWeights_->setNumElements(newNumber);
3345 }
3346 } else {
3347 if (switchType == 1) {
3348 for (i = 0; i < number; i++) {
3349 devex_ += work[i] * work[i];
3350 }
3351 devex_ += ADD_ONE;
3352 } else {
3353 for (i = 0; i < number; i++) {
3354 int iRow = which[i];
3355 int iPivot = pivotVariable[iRow];
3356 if (reference(iPivot)) {
3357 devex_ += work[i] * work[i];
3358 }
3359 }
3360 if (reference(sequenceIn))
3361 devex_ += 1.0;
3362 }
3363 }
3364 }
3365 double oldDevex = weights_[sequenceIn];
3366#ifdef CLP_DEBUG
3367 if ((model_->messageHandler()->logLevel() & 32))
3368 printf("old weight %g, new %g\n", oldDevex, devex_);
3369#endif
3370 double check = CoinMax(devex_, oldDevex) + 0.1;
3371 weights_[sequenceIn] = devex_;
3372 double testValue = 0.1;
3373 if (mode_ == 4 && numberSwitched_ == 1)
3374 testValue = 0.5;
3375 if ( fabs ( devex_ - oldDevex ) > testValue * check ) {
3376#ifdef CLP_DEBUG
3377 if ((model_->messageHandler()->logLevel() & 48) == 16)
3378 printf("old weight %g, new %g\n", oldDevex, devex_);
3379#endif
3380 //printf("old weight %g, new %g\n",oldDevex,devex_);
3381 testValue = 0.99;
3382 if (mode_ == 1)
3383 testValue = 1.01e1; // make unlikely to do if steepest
3384 else if (mode_ == 4 && numberSwitched_ == 1)
3385 testValue = 0.9;
3386 double difference = fabs(devex_ - oldDevex);
3387 if ( difference > testValue * check ) {
3388 // need to redo
3389 model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
3390 *model_->messagesPointer())
3391 << oldDevex << devex_
3392 << CoinMessageEol;
3393 initializeWeights();
3394 }
3395 }
3396 if (pivotRow >= 0) {
3397 // set outgoing weight here
3398 weights_[model_->sequenceOut()] = devex_ / (model_->alpha() * model_->alpha());
3399 }
3400}
3401// Checks accuracy - just for debug
3402void
3403ClpPrimalColumnSteepest::checkAccuracy(int sequence,
3404 double relativeTolerance,
3405 CoinIndexedVector * rowArray1,
3406 CoinIndexedVector * rowArray2)
3407{
3408 if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
3409 return;
3410 model_->unpack(rowArray1, sequence);
3411 model_->factorization()->updateColumn(rowArray2, rowArray1);
3412 int number = rowArray1->getNumElements();
3413 int * which = rowArray1->getIndices();
3414 double * work = rowArray1->denseVector();
3415 const int * pivotVariable = model_->pivotVariable();
3416
3417 double devex = 0.0;
3418 int i;
3419
3420 if (mode_ == 1) {
3421 for (i = 0; i < number; i++) {
3422 int iRow = which[i];
3423 devex += work[iRow] * work[iRow];
3424 work[iRow] = 0.0;
3425 }
3426 devex += ADD_ONE;
3427 } else {
3428 for (i = 0; i < number; i++) {
3429 int iRow = which[i];
3430 int iPivot = pivotVariable[iRow];
3431 if (reference(iPivot)) {
3432 devex += work[iRow] * work[iRow];
3433 }
3434 work[iRow] = 0.0;
3435 }
3436 if (reference(sequence))
3437 devex += 1.0;
3438 }
3439
3440 double oldDevex = weights_[sequence];
3441 double check = CoinMax(devex, oldDevex);
3442 if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
3443 COIN_DETAIL_PRINT(printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex));
3444 // update so won't print again
3445 weights_[sequence] = devex;
3446 }
3447 rowArray1->setNumElements(0);
3448}
3449
3450// Initialize weights
3451void
3452ClpPrimalColumnSteepest::initializeWeights()
3453{
3454 int numberRows = model_->numberRows();
3455 int numberColumns = model_->numberColumns();
3456 int number = numberRows + numberColumns;
3457 int iSequence;
3458 if (mode_ != 1) {
3459 // initialize to 1.0
3460 // and set reference framework
3461 if (!reference_) {
3462 int nWords = (number + 31) >> 5;
3463 reference_ = new unsigned int[nWords];
3464 CoinZeroN(reference_, nWords);
3465 }
3466
3467 for (iSequence = 0; iSequence < number; iSequence++) {
3468 weights_[iSequence] = 1.0;
3469 if (model_->getStatus(iSequence) == ClpSimplex::basic) {
3470 setReference(iSequence, false);
3471 } else {
3472 setReference(iSequence, true);
3473 }
3474 }
3475 } else {
3476 CoinIndexedVector * temp = new CoinIndexedVector();
3477 temp->reserve(numberRows +
3478 model_->factorization()->maximumPivots());
3479 double * array = alternateWeights_->denseVector();
3480 int * which = alternateWeights_->getIndices();
3481
3482 for (iSequence = 0; iSequence < number; iSequence++) {
3483 weights_[iSequence] = 1.0 + ADD_ONE;
3484 if (model_->getStatus(iSequence) != ClpSimplex::basic &&
3485 model_->getStatus(iSequence) != ClpSimplex::isFixed) {
3486 model_->unpack(alternateWeights_, iSequence);
3487 double value = ADD_ONE;
3488 model_->factorization()->updateColumn(temp, alternateWeights_);
3489 int number = alternateWeights_->getNumElements();
3490 int j;
3491 for (j = 0; j < number; j++) {
3492 int iRow = which[j];
3493 value += array[iRow] * array[iRow];
3494 array[iRow] = 0.0;
3495 }
3496 alternateWeights_->setNumElements(0);
3497 weights_[iSequence] = value;
3498 }
3499 }
3500 delete temp;
3501 }
3502}
3503// Gets rid of all arrays
3504void
3505ClpPrimalColumnSteepest::clearArrays()
3506{
3507 if (persistence_ == normal) {
3508 delete [] weights_;
3509 weights_ = NULL;
3510 delete infeasible_;
3511 infeasible_ = NULL;
3512 delete alternateWeights_;
3513 alternateWeights_ = NULL;
3514 delete [] savedWeights_;
3515 savedWeights_ = NULL;
3516 delete [] reference_;
3517 reference_ = NULL;
3518 }
3519 pivotSequence_ = -1;
3520 state_ = -1;
3521 savedPivotSequence_ = -1;
3522 savedSequenceOut_ = -1;
3523 devex_ = 0.0;
3524}
3525// Returns true if would not find any column
3526bool
3527ClpPrimalColumnSteepest::looksOptimal() const
3528{
3529 if (looksOptimal_)
3530 return true; // user overrode
3531 //**** THIS MUST MATCH the action coding above
3532 double tolerance = model_->currentDualTolerance();
3533 // we can't really trust infeasibilities if there is dual error
3534 // this coding has to mimic coding in checkDualSolution
3535 double error = CoinMin(1.0e-2, model_->largestDualError());
3536 // allow tolerance at least slightly bigger than standard
3537 tolerance = tolerance + error;
3538 if(model_->numberIterations() < model_->lastBadIteration() + 200) {
3539 // we can't really trust infeasibilities if there is dual error
3540 double checkTolerance = 1.0e-8;
3541 if (!model_->factorization()->pivots())
3542 checkTolerance = 1.0e-6;
3543 if (model_->largestDualError() > checkTolerance)
3544 tolerance *= model_->largestDualError() / checkTolerance;
3545 // But cap
3546 tolerance = CoinMin(1000.0, tolerance);
3547 }
3548 int number = model_->numberRows() + model_->numberColumns();
3549 int iSequence;
3550
3551 double * reducedCost = model_->djRegion();
3552 int numberInfeasible = 0;
3553 if (!model_->nonLinearCost()->lookBothWays()) {
3554 for (iSequence = 0; iSequence < number; iSequence++) {
3555 double value = reducedCost[iSequence];
3556 ClpSimplex::Status status = model_->getStatus(iSequence);
3557
3558 switch(status) {
3559
3560 case ClpSimplex::basic:
3561 case ClpSimplex::isFixed:
3562 break;
3563 case ClpSimplex::isFree:
3564 case ClpSimplex::superBasic:
3565 if (fabs(value) > FREE_ACCEPT * tolerance)
3566 numberInfeasible++;
3567 break;
3568 case ClpSimplex::atUpperBound:
3569 if (value > tolerance)
3570 numberInfeasible++;
3571 break;
3572 case ClpSimplex::atLowerBound:
3573 if (value < -tolerance)
3574 numberInfeasible++;
3575 }
3576 }
3577 } else {
3578 ClpNonLinearCost * nonLinear = model_->nonLinearCost();
3579 // can go both ways
3580 for (iSequence = 0; iSequence < number; iSequence++) {
3581 double value = reducedCost[iSequence];
3582 ClpSimplex::Status status = model_->getStatus(iSequence);
3583
3584 switch(status) {
3585
3586 case ClpSimplex::basic:
3587 case ClpSimplex::isFixed:
3588 break;
3589 case ClpSimplex::isFree:
3590 case ClpSimplex::superBasic:
3591 if (fabs(value) > FREE_ACCEPT * tolerance)
3592 numberInfeasible++;
3593 break;
3594 case ClpSimplex::atUpperBound:
3595 if (value > tolerance) {
3596 numberInfeasible++;
3597 } else {
3598 // look other way - change up should be negative
3599 value -= nonLinear->changeUpInCost(iSequence);
3600 if (value < -tolerance)
3601 numberInfeasible++;
3602 }
3603 break;
3604 case ClpSimplex::atLowerBound:
3605 if (value < -tolerance) {
3606 numberInfeasible++;
3607 } else {
3608 // look other way - change down should be positive
3609 value -= nonLinear->changeDownInCost(iSequence);
3610 if (value > tolerance)
3611 numberInfeasible++;
3612 }
3613 }
3614 }
3615 }
3616 return numberInfeasible == 0;
3617}
3618/* Returns number of extra columns for sprint algorithm - 0 means off.
3619 Also number of iterations before recompute
3620*/
3621int
3622ClpPrimalColumnSteepest::numberSprintColumns(int & numberIterations) const
3623{
3624 numberIterations = 0;
3625 int numberAdd = 0;
3626 if (!numberSwitched_ && mode_ >= 10) {
3627 numberIterations = CoinMin(2000, model_->numberRows() / 5);
3628 numberIterations = CoinMax(numberIterations, model_->factorizationFrequency());
3629 numberIterations = CoinMax(numberIterations, 500);
3630 if (mode_ == 10) {
3631 numberAdd = CoinMax(300, model_->numberColumns() / 10);
3632 numberAdd = CoinMax(numberAdd, model_->numberRows() / 5);
3633 // fake all
3634 //numberAdd=1000000;
3635 numberAdd = CoinMin(numberAdd, model_->numberColumns());
3636 } else {
3637 abort();
3638 }
3639 }
3640 return numberAdd;
3641}
3642// Switch off sprint idea
3643void
3644ClpPrimalColumnSteepest::switchOffSprint()
3645{
3646 numberSwitched_ = 10;
3647}
3648// Update djs doing partial pricing (dantzig)
3649int
3650ClpPrimalColumnSteepest::partialPricing(CoinIndexedVector * updates,
3651 CoinIndexedVector * spareRow2,
3652 int numberWanted,
3653 int numberLook)
3654{
3655 int number = 0;
3656 int * index;
3657 double * updateBy;
3658 double * reducedCost;
3659 double saveTolerance = model_->currentDualTolerance();
3660 double tolerance = model_->currentDualTolerance();
3661 // we can't really trust infeasibilities if there is dual error
3662 // this coding has to mimic coding in checkDualSolution
3663 double error = CoinMin(1.0e-2, model_->largestDualError());
3664 // allow tolerance at least slightly bigger than standard
3665 tolerance = tolerance + error;
3666 if(model_->numberIterations() < model_->lastBadIteration() + 200) {
3667 // we can't really trust infeasibilities if there is dual error
3668 double checkTolerance = 1.0e-8;
3669 if (!model_->factorization()->pivots())
3670 checkTolerance = 1.0e-6;
3671 if (model_->largestDualError() > checkTolerance)
3672 tolerance *= model_->largestDualError() / checkTolerance;
3673 // But cap
3674 tolerance = CoinMin(1000.0, tolerance);
3675 }
3676 if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
3677 tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
3678 // So partial pricing can use
3679 model_->setCurrentDualTolerance(tolerance);
3680 model_->factorization()->updateColumnTranspose(spareRow2, updates);
3681 int numberColumns = model_->numberColumns();
3682
3683 // Rows
3684 reducedCost = model_->djRegion(0);
3685
3686 number = updates->getNumElements();
3687 index = updates->getIndices();
3688 updateBy = updates->denseVector();
3689 int j;
3690 double * duals = model_->dualRowSolution();
3691 for (j = 0; j < number; j++) {
3692 int iSequence = index[j];
3693 double value = duals[iSequence];
3694 value -= updateBy[j];
3695 updateBy[j] = 0.0;
3696 duals[iSequence] = value;
3697 }
3698 //#define CLP_DEBUG
3699#ifdef CLP_DEBUG
3700 // check duals
3701 {
3702 int numberRows = model_->numberRows();
3703 //work space
3704 CoinIndexedVector arrayVector;
3705 arrayVector.reserve(numberRows + 1000);
3706 CoinIndexedVector workSpace;
3707 workSpace.reserve(numberRows + 1000);
3708
3709
3710 int iRow;
3711 double * array = arrayVector.denseVector();
3712 int * index = arrayVector.getIndices();
3713 int number = 0;
3714 int * pivotVariable = model_->pivotVariable();
3715 double * cost = model_->costRegion();
3716 for (iRow = 0; iRow < numberRows; iRow++) {
3717 int iPivot = pivotVariable[iRow];
3718 double value = cost[iPivot];
3719 if (value) {
3720 array[iRow] = value;
3721 index[number++] = iRow;
3722 }
3723 }
3724 arrayVector.setNumElements(number);
3725 // Extended duals before "updateTranspose"
3726 model_->clpMatrix()->dualExpanded(model_, &arrayVector, NULL, 0);
3727
3728 // Btran basic costs
3729 model_->factorization()->updateColumnTranspose(&workSpace, &arrayVector);
3730
3731 // now look at dual solution
3732 for (iRow = 0; iRow < numberRows; iRow++) {
3733 // slack
3734 double value = array[iRow];
3735 if (fabs(duals[iRow] - value) > 1.0e-3)
3736 printf("bad row %d old dual %g new %g\n", iRow, duals[iRow], value);
3737 //duals[iRow]=value;
3738 }
3739 }
3740#endif
3741#undef CLP_DEBUG
3742 double bestDj = tolerance;
3743 int bestSequence = -1;
3744
3745 const double * cost = model_->costRegion(1);
3746
3747 model_->clpMatrix()->setOriginalWanted(numberWanted);
3748 model_->clpMatrix()->setCurrentWanted(numberWanted);
3749 int iPassR = 0, iPassC = 0;
3750 // Setup two passes
3751 // This biases towards picking row variables
3752 // This probably should be fixed
3753 int startR[4];
3754 const int * which = infeasible_->getIndices();
3755 int nSlacks = infeasible_->getNumElements();
3756 startR[1] = nSlacks;
3757 startR[2] = 0;
3758 double randomR = model_->randomNumberGenerator()->randomDouble();
3759 double dstart = static_cast<double> (nSlacks) * randomR;
3760 startR[0] = static_cast<int> (dstart);
3761 startR[3] = startR[0];
3762 double startC[4];
3763 startC[1] = 1.0;
3764 startC[2] = 0;
3765 double randomC = model_->randomNumberGenerator()->randomDouble();
3766 startC[0] = randomC;
3767 startC[3] = randomC;
3768 reducedCost = model_->djRegion(1);
3769 int sequenceOut = model_->sequenceOut();
3770 double * duals2 = duals - numberColumns;
3771 int chunk = CoinMin(1024, (numberColumns + nSlacks) / 32);
3772#ifdef COIN_DETAIL
3773 if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
3774 printf("%d wanted, nSlacks %d, chunk %d\n", numberWanted, nSlacks, chunk);
3775 int i;
3776 for (i = 0; i < 4; i++)
3777 printf("start R %d C %g ", startR[i], startC[i]);
3778 printf("\n");
3779 }
3780#endif
3781 chunk = CoinMax(chunk, 256);
3782 bool finishedR = false, finishedC = false;
3783 bool doingR = randomR > randomC;
3784 //doingR=false;
3785 int saveNumberWanted = numberWanted;
3786 while (!finishedR || !finishedC) {
3787 if (finishedR)
3788 doingR = false;
3789 if (doingR) {
3790 int saveSequence = bestSequence;
3791 int start = startR[iPassR];
3792 int end = CoinMin(startR[iPassR+1], start + chunk / 2);
3793 int jSequence;
3794 for (jSequence = start; jSequence < end; jSequence++) {
3795 int iSequence = which[jSequence];
3796 if (iSequence != sequenceOut) {
3797 double value;
3798 ClpSimplex::Status status = model_->getStatus(iSequence);
3799
3800 switch(status) {
3801
3802 case ClpSimplex::basic:
3803 case ClpSimplex::isFixed:
3804 break;
3805 case ClpSimplex::isFree:
3806 case ClpSimplex::superBasic:
3807 value = fabs(cost[iSequence] + duals2[iSequence]);
3808 if (value > FREE_ACCEPT * tolerance) {
3809 numberWanted--;
3810 // we are going to bias towards free (but only if reasonable)
3811 value *= FREE_BIAS;
3812 if (value > bestDj) {
3813 // check flagged variable and correct dj
3814 if (!model_->flagged(iSequence)) {
3815 bestDj = value;
3816 bestSequence = iSequence;
3817 } else {
3818 // just to make sure we don't exit before got something
3819 numberWanted++;
3820 }
3821 }
3822 }
3823 break;
3824 case ClpSimplex::atUpperBound:
3825 value = cost[iSequence] + duals2[iSequence];
3826 if (value > tolerance) {
3827 numberWanted--;
3828 if (value > bestDj) {
3829 // check flagged variable and correct dj
3830 if (!model_->flagged(iSequence)) {
3831 bestDj = value;
3832 bestSequence = iSequence;
3833 } else {
3834 // just to make sure we don't exit before got something
3835 numberWanted++;
3836 }
3837 }
3838 }
3839 break;
3840 case ClpSimplex::atLowerBound:
3841 value = -(cost[iSequence] + duals2[iSequence]);
3842 if (value > tolerance) {
3843 numberWanted--;
3844 if (value > bestDj) {
3845 // check flagged variable and correct dj
3846 if (!model_->flagged(iSequence)) {
3847 bestDj = value;
3848 bestSequence = iSequence;
3849 } else {
3850 // just to make sure we don't exit before got something
3851 numberWanted++;
3852 }
3853 }
3854 }
3855 break;
3856 }
3857 }
3858 if (!numberWanted)
3859 break;
3860 }
3861 numberLook -= (end - start);
3862 if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
3863 numberWanted = 0; // give up
3864 if (saveSequence != bestSequence) {
3865 // dj
3866 reducedCost[bestSequence] = cost[bestSequence] + duals[bestSequence-numberColumns];
3867 bestDj = fabs(reducedCost[bestSequence]);
3868 model_->clpMatrix()->setSavedBestSequence(bestSequence);
3869 model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
3870 }
3871 model_->clpMatrix()->setCurrentWanted(numberWanted);
3872 if (!numberWanted)
3873 break;
3874 doingR = false;
3875 // update start
3876 startR[iPassR] = jSequence;
3877 if (jSequence >= startR[iPassR+1]) {
3878 if (iPassR)
3879 finishedR = true;
3880 else
3881 iPassR = 2;
3882 }
3883 }
3884 if (finishedC)
3885 doingR = true;
3886 if (!doingR) {
3887 int saveSequence = bestSequence;
3888 // Columns
3889 double start = startC[iPassC];
3890 // If we put this idea back then each function needs to update endFraction **
3891#if 0
3892 double dchunk = (static_cast<double> chunk) / (static_cast<double> numberColumns);
3893 double end = CoinMin(startC[iPassC+1], start + dchunk);
3894#else
3895 double end = startC[iPassC+1]; // force end
3896#endif
3897 model_->clpMatrix()->partialPricing(model_, start, end, bestSequence, numberWanted);
3898 numberWanted = model_->clpMatrix()->currentWanted();
3899 numberLook -= static_cast<int> ((end - start) * numberColumns);
3900 if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
3901 numberWanted = 0; // give up
3902 if (saveSequence != bestSequence) {
3903 // dj
3904 bestDj = fabs(model_->clpMatrix()->reducedCost(model_, bestSequence));
3905 }
3906 if (!numberWanted)
3907 break;
3908 doingR = true;
3909 // update start
3910 startC[iPassC] = end;
3911 if (end >= startC[iPassC+1] - 1.0e-8) {
3912 if (iPassC)
3913 finishedC = true;
3914 else
3915 iPassC = 2;
3916 }
3917 }
3918 }
3919 updates->setNumElements(0);
3920
3921 // Restore tolerance
3922 model_->setCurrentDualTolerance(saveTolerance);
3923 // Now create variable if column generation
3924 model_->clpMatrix()->createVariable(model_, bestSequence);
3925#ifndef NDEBUG
3926 if (bestSequence >= 0) {
3927 if (model_->getStatus(bestSequence) == ClpSimplex::atLowerBound)
3928 assert(model_->reducedCost(bestSequence) < 0.0);
3929 if (model_->getStatus(bestSequence) == ClpSimplex::atUpperBound)
3930 assert(model_->reducedCost(bestSequence) > 0.0);
3931 }
3932#endif
3933 return bestSequence;
3934}
3935