1/* $Id: ClpDualRowSteepest.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#include "ClpSimplex.hpp"
8#include "ClpDualRowSteepest.hpp"
9#include "CoinIndexedVector.hpp"
10#include "ClpFactorization.hpp"
11#include "CoinHelperFunctions.hpp"
12#include <cstdio>
13//#############################################################################
14// Constructors / Destructor / Assignment
15//#############################################################################
16//#define CLP_DEBUG 4
17//-------------------------------------------------------------------
18// Default Constructor
19//-------------------------------------------------------------------
20ClpDualRowSteepest::ClpDualRowSteepest (int mode)
21 : ClpDualRowPivot(),
22 state_(-1),
23 mode_(mode),
24 persistence_(normal),
25 weights_(NULL),
26 infeasible_(NULL),
27 alternateWeights_(NULL),
28 savedWeights_(NULL),
29 dubiousWeights_(NULL)
30{
31 type_ = 2 + 64 * mode;
32}
33
34//-------------------------------------------------------------------
35// Copy constructor
36//-------------------------------------------------------------------
37ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs)
38 : ClpDualRowPivot(rhs)
39{
40 state_ = rhs.state_;
41 mode_ = rhs.mode_;
42 persistence_ = rhs.persistence_;
43 model_ = rhs.model_;
44 if ((model_ && model_->whatsChanged() & 1) != 0) {
45 int number = model_->numberRows();
46 if (rhs.savedWeights_)
47 number = CoinMin(number, rhs.savedWeights_->capacity());
48 if (rhs.infeasible_) {
49 infeasible_ = new CoinIndexedVector(rhs.infeasible_);
50 } else {
51 infeasible_ = NULL;
52 }
53 if (rhs.weights_) {
54 weights_ = new double[number];
55 ClpDisjointCopyN(rhs.weights_, number, weights_);
56 } else {
57 weights_ = NULL;
58 }
59 if (rhs.alternateWeights_) {
60 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
61 } else {
62 alternateWeights_ = NULL;
63 }
64 if (rhs.savedWeights_) {
65 savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
66 } else {
67 savedWeights_ = NULL;
68 }
69 if (rhs.dubiousWeights_) {
70 assert(model_);
71 int number = model_->numberRows();
72 dubiousWeights_ = new int[number];
73 ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
74 } else {
75 dubiousWeights_ = NULL;
76 }
77 } else {
78 infeasible_ = NULL;
79 weights_ = NULL;
80 alternateWeights_ = NULL;
81 savedWeights_ = NULL;
82 dubiousWeights_ = NULL;
83 }
84}
85
86//-------------------------------------------------------------------
87// Destructor
88//-------------------------------------------------------------------
89ClpDualRowSteepest::~ClpDualRowSteepest ()
90{
91 delete [] weights_;
92 delete [] dubiousWeights_;
93 delete infeasible_;
94 delete alternateWeights_;
95 delete savedWeights_;
96}
97
98//----------------------------------------------------------------
99// Assignment operator
100//-------------------------------------------------------------------
101ClpDualRowSteepest &
102ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
103{
104 if (this != &rhs) {
105 ClpDualRowPivot::operator=(rhs);
106 state_ = rhs.state_;
107 mode_ = rhs.mode_;
108 persistence_ = rhs.persistence_;
109 model_ = rhs.model_;
110 delete [] weights_;
111 delete [] dubiousWeights_;
112 delete infeasible_;
113 delete alternateWeights_;
114 delete savedWeights_;
115 assert(model_);
116 int number = model_->numberRows();
117 if (rhs.savedWeights_)
118 number = CoinMin(number, rhs.savedWeights_->capacity());
119 if (rhs.infeasible_ != NULL) {
120 infeasible_ = new CoinIndexedVector(rhs.infeasible_);
121 } else {
122 infeasible_ = NULL;
123 }
124 if (rhs.weights_ != NULL) {
125 weights_ = new double[number];
126 ClpDisjointCopyN(rhs.weights_, number, weights_);
127 } else {
128 weights_ = NULL;
129 }
130 if (rhs.alternateWeights_ != NULL) {
131 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
132 } else {
133 alternateWeights_ = NULL;
134 }
135 if (rhs.savedWeights_ != NULL) {
136 savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
137 } else {
138 savedWeights_ = NULL;
139 }
140 if (rhs.dubiousWeights_) {
141 assert(model_);
142 int number = model_->numberRows();
143 dubiousWeights_ = new int[number];
144 ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
145 } else {
146 dubiousWeights_ = NULL;
147 }
148 }
149 return *this;
150}
151// Fill most values
152void
153ClpDualRowSteepest::fill(const ClpDualRowSteepest& rhs)
154{
155 state_ = rhs.state_;
156 mode_ = rhs.mode_;
157 persistence_ = rhs.persistence_;
158 assert (model_->numberRows() == rhs.model_->numberRows());
159 model_ = rhs.model_;
160 assert(model_);
161 int number = model_->numberRows();
162 if (rhs.savedWeights_)
163 number = CoinMin(number, rhs.savedWeights_->capacity());
164 if (rhs.infeasible_ != NULL) {
165 if (!infeasible_)
166 infeasible_ = new CoinIndexedVector(rhs.infeasible_);
167 else
168 *infeasible_ = *rhs.infeasible_;
169 } else {
170 delete infeasible_;
171 infeasible_ = NULL;
172 }
173 if (rhs.weights_ != NULL) {
174 if (!weights_)
175 weights_ = new double[number];
176 ClpDisjointCopyN(rhs.weights_, number, weights_);
177 } else {
178 delete [] weights_;
179 weights_ = NULL;
180 }
181 if (rhs.alternateWeights_ != NULL) {
182 if (!alternateWeights_)
183 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
184 else
185 *alternateWeights_ = *rhs.alternateWeights_;
186 } else {
187 delete alternateWeights_;
188 alternateWeights_ = NULL;
189 }
190 if (rhs.savedWeights_ != NULL) {
191 if (!savedWeights_)
192 savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
193 else
194 *savedWeights_ = *rhs.savedWeights_;
195 } else {
196 delete savedWeights_;
197 savedWeights_ = NULL;
198 }
199 if (rhs.dubiousWeights_) {
200 assert(model_);
201 int number = model_->numberRows();
202 if (!dubiousWeights_)
203 dubiousWeights_ = new int[number];
204 ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
205 } else {
206 delete [] dubiousWeights_;
207 dubiousWeights_ = NULL;
208 }
209}
210// Returns pivot row, -1 if none
211int
212ClpDualRowSteepest::pivotRow()
213{
214 assert(model_);
215 int i, iRow;
216 double * infeas = infeasible_->denseVector();
217 double largest = 0.0;
218 int * index = infeasible_->getIndices();
219 int number = infeasible_->getNumElements();
220 const int * pivotVariable = model_->pivotVariable();
221 int chosenRow = -1;
222 int lastPivotRow = model_->pivotRow();
223 assert (lastPivotRow < model_->numberRows());
224 double tolerance = model_->currentPrimalTolerance();
225 // we can't really trust infeasibilities if there is primal error
226 // this coding has to mimic coding in checkPrimalSolution
227 double error = CoinMin(1.0e-2, model_->largestPrimalError());
228 // allow tolerance at least slightly bigger than standard
229 tolerance = tolerance + error;
230 // But cap
231 tolerance = CoinMin(1000.0, tolerance);
232 tolerance *= tolerance; // as we are using squares
233 bool toleranceChanged = false;
234 double * solution = model_->solutionRegion();
235 double * lower = model_->lowerRegion();
236 double * upper = model_->upperRegion();
237 // do last pivot row one here
238 //#define CLP_DUAL_FIXED_COLUMN_MULTIPLIER 10.0
239 if (lastPivotRow >= 0 && lastPivotRow < model_->numberRows()) {
240#ifdef CLP_DUAL_COLUMN_MULTIPLIER
241 int numberColumns = model_->numberColumns();
242#endif
243 int iPivot = pivotVariable[lastPivotRow];
244 double value = solution[iPivot];
245 double lower = model_->lower(iPivot);
246 double upper = model_->upper(iPivot);
247 if (value > upper + tolerance) {
248 value -= upper;
249 value *= value;
250#ifdef CLP_DUAL_COLUMN_MULTIPLIER
251 if (iPivot < numberColumns)
252 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
253#endif
254 // store square in list
255 if (infeas[lastPivotRow])
256 infeas[lastPivotRow] = value; // already there
257 else
258 infeasible_->quickAdd(lastPivotRow, value);
259 } else if (value < lower - tolerance) {
260 value -= lower;
261 value *= value;
262#ifdef CLP_DUAL_COLUMN_MULTIPLIER
263 if (iPivot < numberColumns)
264 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
265#endif
266 // store square in list
267 if (infeas[lastPivotRow])
268 infeas[lastPivotRow] = value; // already there
269 else
270 infeasible_->add(lastPivotRow, value);
271 } else {
272 // feasible - was it infeasible - if so set tiny
273 if (infeas[lastPivotRow])
274 infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
275 }
276 number = infeasible_->getNumElements();
277 }
278 if(model_->numberIterations() < model_->lastBadIteration() + 200) {
279 // we can't really trust infeasibilities if there is dual error
280 if (model_->largestDualError() > model_->largestPrimalError()) {
281 tolerance *= CoinMin(model_->largestDualError() / model_->largestPrimalError(), 1000.0);
282 toleranceChanged = true;
283 }
284 }
285 int numberWanted;
286 if (mode_ < 2 ) {
287 numberWanted = number + 1;
288 } else if (mode_ == 2) {
289 numberWanted = CoinMax(2000, number / 8);
290 } else {
291 int numberElements = model_->factorization()->numberElements();
292 double ratio = static_cast<double> (numberElements) /
293 static_cast<double> (model_->numberRows());
294 numberWanted = CoinMax(2000, number / 8);
295 if (ratio < 1.0) {
296 numberWanted = CoinMax(2000, number / 20);
297 } else if (ratio > 10.0) {
298 ratio = number * (ratio / 80.0);
299 if (ratio > number)
300 numberWanted = number + 1;
301 else
302 numberWanted = CoinMax(2000, static_cast<int> (ratio));
303 }
304 }
305 if (model_->largestPrimalError() > 1.0e-3)
306 numberWanted = number + 1; // be safe
307 int iPass;
308 // Setup two passes
309 int start[4];
310 start[1] = number;
311 start[2] = 0;
312 double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
313 start[0] = static_cast<int> (dstart);
314 start[3] = start[0];
315 //double largestWeight=0.0;
316 //double smallestWeight=1.0e100;
317 for (iPass = 0; iPass < 2; iPass++) {
318 int end = start[2*iPass+1];
319 for (i = start[2*iPass]; i < end; i++) {
320 iRow = index[i];
321 double value = infeas[iRow];
322 if (value > tolerance) {
323 //#define OUT_EQ
324#ifdef OUT_EQ
325 {
326 int iSequence = pivotVariable[iRow];
327 if (upper[iSequence] == lower[iSequence])
328 value *= 2.0;
329 }
330#endif
331 double weight = CoinMin(weights_[iRow], 1.0e50);
332 //largestWeight = CoinMax(largestWeight,weight);
333 //smallestWeight = CoinMin(smallestWeight,weight);
334 //double dubious = dubiousWeights_[iRow];
335 //weight *= dubious;
336 //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
337 if (value > largest * weight) {
338 // make last pivot row last resort choice
339 if (iRow == lastPivotRow) {
340 if (value * 1.0e-10 < largest * weight)
341 continue;
342 else
343 value *= 1.0e-10;
344 }
345 int iSequence = pivotVariable[iRow];
346 if (!model_->flagged(iSequence)) {
347 //#define CLP_DEBUG 3
348#ifdef CLP_DEBUG
349 double value2 = 0.0;
350 if (solution[iSequence] > upper[iSequence] + tolerance)
351 value2 = solution[iSequence] - upper[iSequence];
352 else if (solution[iSequence] < lower[iSequence] - tolerance)
353 value2 = solution[iSequence] - lower[iSequence];
354 assert(fabs(value2 * value2 - infeas[iRow]) < 1.0e-8 * CoinMin(value2 * value2, infeas[iRow]));
355#endif
356 if (solution[iSequence] > upper[iSequence] + tolerance ||
357 solution[iSequence] < lower[iSequence] - tolerance) {
358 chosenRow = iRow;
359 largest = value / weight;
360 //largestWeight = dubious;
361 }
362 } else {
363 // just to make sure we don't exit before got something
364 numberWanted++;
365 }
366 }
367 numberWanted--;
368 if (!numberWanted)
369 break;
370 }
371 }
372 if (!numberWanted)
373 break;
374 }
375 //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
376 if (chosenRow < 0 && toleranceChanged) {
377 // won't line up with checkPrimalSolution - do again
378 double saveError = model_->largestDualError();
379 model_->setLargestDualError(0.0);
380 // can't loop
381 chosenRow = pivotRow();
382 model_->setLargestDualError(saveError);
383 }
384 if (chosenRow < 0 && lastPivotRow < 0) {
385 int nLeft = 0;
386 for (int i = 0; i < number; i++) {
387 int iRow = index[i];
388 if (fabs(infeas[iRow]) > 1.0e-50) {
389 index[nLeft++] = iRow;
390 } else {
391 infeas[iRow] = 0.0;
392 }
393 }
394 infeasible_->setNumElements(nLeft);
395 model_->setNumberPrimalInfeasibilities(nLeft);
396 }
397 return chosenRow;
398}
399#if 0
400static double ft_count = 0.0;
401static double up_count = 0.0;
402static double ft_count_in = 0.0;
403static double up_count_in = 0.0;
404static int xx_count = 0;
405#endif
406/* Updates weights and returns pivot alpha.
407 Also does FT update */
408double
409ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
410 CoinIndexedVector * spare,
411 CoinIndexedVector * spare2,
412 CoinIndexedVector * updatedColumn)
413{
414 //#define CLP_DEBUG 3
415#if CLP_DEBUG>2
416 // Very expensive debug
417 {
418 int numberRows = model_->numberRows();
419 CoinIndexedVector * temp = new CoinIndexedVector();
420 temp->reserve(numberRows +
421 model_->factorization()->maximumPivots());
422 double * array = alternateWeights_->denseVector();
423 int * which = alternateWeights_->getIndices();
424 alternateWeights_->clear();
425 int i;
426 for (i = 0; i < numberRows; i++) {
427 double value = 0.0;
428 array[i] = 1.0;
429 which[0] = i;
430 alternateWeights_->setNumElements(1);
431 model_->factorization()->updateColumnTranspose(temp,
432 alternateWeights_);
433 int number = alternateWeights_->getNumElements();
434 int j;
435 for (j = 0; j < number; j++) {
436 int iRow = which[j];
437 value += array[iRow] * array[iRow];
438 array[iRow] = 0.0;
439 }
440 alternateWeights_->setNumElements(0);
441 double w = CoinMax(weights_[i], value) * .1;
442 if (fabs(weights_[i] - value) > w) {
443 printf("%d old %g, true %g\n", i, weights_[i], value);
444 weights_[i] = value; // to reduce printout
445 }
446 //else
447 //printf("%d matches %g\n",i,value);
448 }
449 delete temp;
450 }
451#endif
452 assert (input->packedMode());
453 if (!updatedColumn->packedMode()) {
454 // I think this means empty
455#ifdef COIN_DEVELOP
456 printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
457#endif
458 return 0.0;
459 }
460 double alpha = 0.0;
461 if (!model_->factorization()->networkBasis()) {
462 // clear other region
463 alternateWeights_->clear();
464 double norm = 0.0;
465 int i;
466 double * work = input->denseVector();
467 int numberNonZero = input->getNumElements();
468 int * which = input->getIndices();
469 double * work2 = spare->denseVector();
470 int * which2 = spare->getIndices();
471 // ftran
472 //permute and move indices into index array
473 //also compute norm
474 //int *regionIndex = alternateWeights_->getIndices ( );
475 const int *permute = model_->factorization()->permute();
476 //double * region = alternateWeights_->denseVector();
477 if (permute) {
478 for ( i = 0; i < numberNonZero; i ++ ) {
479 int iRow = which[i];
480 double value = work[i];
481 norm += value * value;
482 iRow = permute[iRow];
483 work2[iRow] = value;
484 which2[i] = iRow;
485 }
486 } else {
487 for ( i = 0; i < numberNonZero; i ++ ) {
488 int iRow = which[i];
489 double value = work[i];
490 norm += value * value;
491 //iRow = permute[iRow];
492 work2[iRow] = value;
493 which2[i] = iRow;
494 }
495 }
496 spare->setNumElements ( numberNonZero );
497 // Do FT update
498#if 0
499 ft_count_in += updatedColumn->getNumElements();
500 up_count_in += spare->getNumElements();
501#endif
502 if (permute || true) {
503#if CLP_DEBUG>2
504 printf("REGION before %d els\n", spare->getNumElements());
505 spare->print();
506#endif
507 model_->factorization()->updateTwoColumnsFT(spare2, updatedColumn,
508 spare, permute != NULL);
509#if CLP_DEBUG>2
510 printf("REGION after %d els\n", spare->getNumElements());
511 spare->print();
512#endif
513 } else {
514 // Leave as old way
515 model_->factorization()->updateColumnFT(spare2, updatedColumn);
516 model_->factorization()->updateColumn(spare2, spare, false);
517 }
518#undef CLP_DEBUG
519#if 0
520 ft_count += updatedColumn->getNumElements();
521 up_count += spare->getNumElements();
522 xx_count++;
523 if ((xx_count % 1000) == 0)
524 printf("zz %d ft %g up %g (in %g %g)\n", xx_count, ft_count, up_count,
525 ft_count_in, up_count_in);
526#endif
527 numberNonZero = spare->getNumElements();
528 // alternateWeights_ should still be empty
529 int pivotRow = model_->pivotRow();
530#ifdef CLP_DEBUG
531 if ( model_->logLevel ( ) > 4 &&
532 fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
533 printf("on row %d, true weight %g, old %g\n",
534 pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
535#endif
536 // could re-initialize here (could be expensive)
537 norm /= model_->alpha() * model_->alpha();
538 assert(model_->alpha());
539 assert(norm);
540 // pivot element
541 alpha = 0.0;
542 double multiplier = 2.0 / model_->alpha();
543 // look at updated column
544 work = updatedColumn->denseVector();
545 numberNonZero = updatedColumn->getNumElements();
546 which = updatedColumn->getIndices();
547
548 int nSave = 0;
549 double * work3 = alternateWeights_->denseVector();
550 int * which3 = alternateWeights_->getIndices();
551 const int * pivotColumn = model_->factorization()->pivotColumn();
552 for (i = 0; i < numberNonZero; i++) {
553 int iRow = which[i];
554 double theta = work[i];
555 if (iRow == pivotRow)
556 alpha = theta;
557 double devex = weights_[iRow];
558 work3[nSave] = devex; // save old
559 which3[nSave++] = iRow;
560 // transform to match spare
561 int jRow = permute ? pivotColumn[iRow] : iRow;
562 double value = work2[jRow];
563 devex += theta * (theta * norm + value * multiplier);
564 if (devex < DEVEX_TRY_NORM)
565 devex = DEVEX_TRY_NORM;
566 weights_[iRow] = devex;
567 }
568 alternateWeights_->setPackedMode(true);
569 alternateWeights_->setNumElements(nSave);
570 if (norm < DEVEX_TRY_NORM)
571 norm = DEVEX_TRY_NORM;
572 // Try this to make less likely will happen again and stop cycling
573 //norm *= 1.02;
574 weights_[pivotRow] = norm;
575 spare->clear();
576#ifdef CLP_DEBUG
577 spare->checkClear();
578#endif
579 } else {
580 // Do FT update
581 model_->factorization()->updateColumnFT(spare, updatedColumn);
582 // clear other region
583 alternateWeights_->clear();
584 double norm = 0.0;
585 int i;
586 double * work = input->denseVector();
587 int number = input->getNumElements();
588 int * which = input->getIndices();
589 double * work2 = spare->denseVector();
590 int * which2 = spare->getIndices();
591 for (i = 0; i < number; i++) {
592 int iRow = which[i];
593 double value = work[i];
594 norm += value * value;
595 work2[iRow] = value;
596 which2[i] = iRow;
597 }
598 spare->setNumElements(number);
599 // ftran
600#ifndef NDEBUG
601 alternateWeights_->checkClear();
602#endif
603 model_->factorization()->updateColumn(alternateWeights_, spare);
604 // alternateWeights_ should still be empty
605#ifndef NDEBUG
606 alternateWeights_->checkClear();
607#endif
608 int pivotRow = model_->pivotRow();
609#ifdef CLP_DEBUG
610 if ( model_->logLevel ( ) > 4 &&
611 fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
612 printf("on row %d, true weight %g, old %g\n",
613 pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
614#endif
615 // could re-initialize here (could be expensive)
616 norm /= model_->alpha() * model_->alpha();
617
618 assert(norm);
619 //if (norm < DEVEX_TRY_NORM)
620 //norm = DEVEX_TRY_NORM;
621 // pivot element
622 alpha = 0.0;
623 double multiplier = 2.0 / model_->alpha();
624 // look at updated column
625 work = updatedColumn->denseVector();
626 number = updatedColumn->getNumElements();
627 which = updatedColumn->getIndices();
628
629 int nSave = 0;
630 double * work3 = alternateWeights_->denseVector();
631 int * which3 = alternateWeights_->getIndices();
632 for (i = 0; i < number; i++) {
633 int iRow = which[i];
634 double theta = work[i];
635 if (iRow == pivotRow)
636 alpha = theta;
637 double devex = weights_[iRow];
638 work3[nSave] = devex; // save old
639 which3[nSave++] = iRow;
640 double value = work2[iRow];
641 devex += theta * (theta * norm + value * multiplier);
642 if (devex < DEVEX_TRY_NORM)
643 devex = DEVEX_TRY_NORM;
644 weights_[iRow] = devex;
645 }
646 if (!alpha) {
647 // error - but carry on
648 alpha = 1.0e-50;
649 }
650 alternateWeights_->setPackedMode(true);
651 alternateWeights_->setNumElements(nSave);
652 if (norm < DEVEX_TRY_NORM)
653 norm = DEVEX_TRY_NORM;
654 weights_[pivotRow] = norm;
655 spare->clear();
656 }
657#ifdef CLP_DEBUG
658 spare->checkClear();
659#endif
660 return alpha;
661}
662
663/* Updates primal solution (and maybe list of candidates)
664 Uses input vector which it deletes
665 Computes change in objective function
666*/
667void
668ClpDualRowSteepest::updatePrimalSolution(
669 CoinIndexedVector * primalUpdate,
670 double primalRatio,
671 double & objectiveChange)
672{
673 double * COIN_RESTRICT work = primalUpdate->denseVector();
674 int number = primalUpdate->getNumElements();
675 int * COIN_RESTRICT which = primalUpdate->getIndices();
676 int i;
677 double changeObj = 0.0;
678 double tolerance = model_->currentPrimalTolerance();
679 const int * COIN_RESTRICT pivotVariable = model_->pivotVariable();
680 double * COIN_RESTRICT infeas = infeasible_->denseVector();
681 double * COIN_RESTRICT solution = model_->solutionRegion();
682 const double * COIN_RESTRICT costModel = model_->costRegion();
683 const double * COIN_RESTRICT lowerModel = model_->lowerRegion();
684 const double * COIN_RESTRICT upperModel = model_->upperRegion();
685#ifdef CLP_DUAL_COLUMN_MULTIPLIER
686 int numberColumns = model_->numberColumns();
687#endif
688 if (primalUpdate->packedMode()) {
689 for (i = 0; i < number; i++) {
690 int iRow = which[i];
691 int iPivot = pivotVariable[iRow];
692 double value = solution[iPivot];
693 double cost = costModel[iPivot];
694 double change = primalRatio * work[i];
695 work[i] = 0.0;
696 value -= change;
697 changeObj -= change * cost;
698 double lower = lowerModel[iPivot];
699 double upper = upperModel[iPivot];
700 solution[iPivot] = value;
701 if (value < lower - tolerance) {
702 value -= lower;
703 value *= value;
704#ifdef CLP_DUAL_COLUMN_MULTIPLIER
705 if (iPivot < numberColumns)
706 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
707#endif
708#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
709 if (lower == upper)
710 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
711#endif
712 // store square in list
713 if (infeas[iRow])
714 infeas[iRow] = value; // already there
715 else
716 infeasible_->quickAdd(iRow, value);
717 } else if (value > upper + tolerance) {
718 value -= upper;
719 value *= value;
720#ifdef CLP_DUAL_COLUMN_MULTIPLIER
721 if (iPivot < numberColumns)
722 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
723#endif
724#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
725 if (lower == upper)
726 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
727#endif
728 // store square in list
729 if (infeas[iRow])
730 infeas[iRow] = value; // already there
731 else
732 infeasible_->quickAdd(iRow, value);
733 } else {
734 // feasible - was it infeasible - if so set tiny
735 if (infeas[iRow])
736 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
737 }
738 }
739 } else {
740 for (i = 0; i < number; i++) {
741 int iRow = which[i];
742 int iPivot = pivotVariable[iRow];
743 double value = solution[iPivot];
744 double cost = costModel[iPivot];
745 double change = primalRatio * work[iRow];
746 value -= change;
747 changeObj -= change * cost;
748 double lower = lowerModel[iPivot];
749 double upper = upperModel[iPivot];
750 solution[iPivot] = value;
751 if (value < lower - tolerance) {
752 value -= lower;
753 value *= value;
754#ifdef CLP_DUAL_COLUMN_MULTIPLIER
755 if (iPivot < numberColumns)
756 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
757#endif
758#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
759 if (lower == upper)
760 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
761#endif
762 // store square in list
763 if (infeas[iRow])
764 infeas[iRow] = value; // already there
765 else
766 infeasible_->quickAdd(iRow, value);
767 } else if (value > upper + tolerance) {
768 value -= upper;
769 value *= value;
770#ifdef CLP_DUAL_COLUMN_MULTIPLIER
771 if (iPivot < numberColumns)
772 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
773#endif
774#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
775 if (lower == upper)
776 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
777#endif
778 // store square in list
779 if (infeas[iRow])
780 infeas[iRow] = value; // already there
781 else
782 infeasible_->quickAdd(iRow, value);
783 } else {
784 // feasible - was it infeasible - if so set tiny
785 if (infeas[iRow])
786 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
787 }
788 work[iRow] = 0.0;
789 }
790 }
791 // Do pivot row
792 {
793 int iRow = model_->pivotRow();
794 // feasible - was it infeasible - if so set tiny
795 //assert (infeas[iRow]);
796 if (infeas[iRow])
797 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
798 }
799 primalUpdate->setNumElements(0);
800 objectiveChange += changeObj;
801}
802/* Saves any weights round factorization as pivot rows may change
803 1) before factorization
804 2) after factorization
805 3) just redo infeasibilities
806 4) restore weights
807*/
808void
809ClpDualRowSteepest::saveWeights(ClpSimplex * model, int mode)
810{
811 // alternateWeights_ is defined as indexed but is treated oddly
812 model_ = model;
813 int numberRows = model_->numberRows();
814 int numberColumns = model_->numberColumns();
815 const int * pivotVariable = model_->pivotVariable();
816 int i;
817 if (mode == 1) {
818 if(weights_) {
819 // Check if size has changed
820 if (infeasible_->capacity() == numberRows) {
821 alternateWeights_->clear();
822 // change from row numbers to sequence numbers
823 int * which = alternateWeights_->getIndices();
824 for (i = 0; i < numberRows; i++) {
825 int iPivot = pivotVariable[i];
826 which[i] = iPivot;
827 }
828 state_ = 1;
829 } else {
830 // size has changed - clear everything
831 delete [] weights_;
832 weights_ = NULL;
833 delete [] dubiousWeights_;
834 dubiousWeights_ = NULL;
835 delete infeasible_;
836 infeasible_ = NULL;
837 delete alternateWeights_;
838 alternateWeights_ = NULL;
839 delete savedWeights_;
840 savedWeights_ = NULL;
841 state_ = -1;
842 }
843 }
844 } else if (mode == 2 || mode == 4 || mode >= 5) {
845 // restore
846 if (!weights_ || state_ == -1 || mode == 5) {
847 // initialize weights
848 delete [] weights_;
849 delete alternateWeights_;
850 weights_ = new double[numberRows];
851 alternateWeights_ = new CoinIndexedVector();
852 // enough space so can use it for factorization
853 alternateWeights_->reserve(numberRows +
854 model_->factorization()->maximumPivots());
855 if (mode_ != 1 || mode == 5) {
856 // initialize to 1.0 (can we do better?)
857 for (i = 0; i < numberRows; i++) {
858 weights_[i] = 1.0;
859 }
860 } else {
861 CoinIndexedVector * temp = new CoinIndexedVector();
862 temp->reserve(numberRows +
863 model_->factorization()->maximumPivots());
864 double * array = alternateWeights_->denseVector();
865 int * which = alternateWeights_->getIndices();
866 for (i = 0; i < numberRows; i++) {
867 double value = 0.0;
868 array[0] = 1.0;
869 which[0] = i;
870 alternateWeights_->setNumElements(1);
871 alternateWeights_->setPackedMode(true);
872 model_->factorization()->updateColumnTranspose(temp,
873 alternateWeights_);
874 int number = alternateWeights_->getNumElements();
875 int j;
876 for (j = 0; j < number; j++) {
877 value += array[j] * array[j];
878 array[j] = 0.0;
879 }
880 alternateWeights_->setNumElements(0);
881 weights_[i] = value;
882 }
883 delete temp;
884 }
885 // create saved weights (not really indexedvector)
886 savedWeights_ = new CoinIndexedVector();
887 savedWeights_->reserve(numberRows);
888
889 double * array = savedWeights_->denseVector();
890 int * which = savedWeights_->getIndices();
891 for (i = 0; i < numberRows; i++) {
892 array[i] = weights_[i];
893 which[i] = pivotVariable[i];
894 }
895 } else if (mode != 6) {
896 int * which = alternateWeights_->getIndices();
897 CoinIndexedVector * rowArray3 = model_->rowArray(3);
898 rowArray3->clear();
899 int * back = rowArray3->getIndices();
900 // In case something went wrong
901 for (i = 0; i < numberRows + numberColumns; i++)
902 back[i] = -1;
903 if (mode != 4) {
904 // save
905 CoinMemcpyN(which, numberRows, savedWeights_->getIndices());
906 CoinMemcpyN(weights_, numberRows, savedWeights_->denseVector());
907 } else {
908 // restore
909 //memcpy(which,savedWeights_->getIndices(),
910 // numberRows*sizeof(int));
911 //memcpy(weights_,savedWeights_->denseVector(),
912 // numberRows*sizeof(double));
913 which = savedWeights_->getIndices();
914 }
915 // restore (a bit slow - but only every re-factorization)
916 double * array = savedWeights_->denseVector();
917 for (i = 0; i < numberRows; i++) {
918 int iSeq = which[i];
919 back[iSeq] = i;
920 }
921 for (i = 0; i < numberRows; i++) {
922 int iPivot = pivotVariable[i];
923 iPivot = back[iPivot];
924 if (iPivot >= 0) {
925 weights_[i] = array[iPivot];
926 if (weights_[i] < DEVEX_TRY_NORM)
927 weights_[i] = DEVEX_TRY_NORM; // may need to check more
928 } else {
929 // odd
930 weights_[i] = 1.0;
931 }
932 }
933 } else {
934 // mode 6 - scale back weights as primal errors
935 double primalError = model_->largestPrimalError();
936 double allowed ;
937 if (primalError > 1.0e3)
938 allowed = 10.0;
939 else if (primalError > 1.0e2)
940 allowed = 50.0;
941 else if (primalError > 1.0e1)
942 allowed = 100.0;
943 else
944 allowed = 1000.0;
945 double allowedInv = 1.0 / allowed;
946 for (i = 0; i < numberRows; i++) {
947 double value = weights_[i];
948 if (value < allowedInv)
949 value = allowedInv;
950 else if (value > allowed)
951 value = allowed;
952 weights_[i] = allowed;
953 }
954 }
955 state_ = 0;
956 // set up infeasibilities
957 if (!infeasible_) {
958 infeasible_ = new CoinIndexedVector();
959 infeasible_->reserve(numberRows);
960 }
961 }
962 if (mode >= 2) {
963 // Get dubious weights
964 //if (!dubiousWeights_)
965 //dubiousWeights_=new int[numberRows];
966 //model_->factorization()->getWeights(dubiousWeights_);
967 infeasible_->clear();
968 int iRow;
969 const int * pivotVariable = model_->pivotVariable();
970 double tolerance = model_->currentPrimalTolerance();
971 for (iRow = 0; iRow < numberRows; iRow++) {
972 int iPivot = pivotVariable[iRow];
973 double value = model_->solution(iPivot);
974 double lower = model_->lower(iPivot);
975 double upper = model_->upper(iPivot);
976 if (value < lower - tolerance) {
977 value -= lower;
978 value *= value;
979#ifdef CLP_DUAL_COLUMN_MULTIPLIER
980 if (iPivot < numberColumns)
981 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
982#endif
983#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
984 if (lower == upper)
985 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
986#endif
987 // store square in list
988 infeasible_->quickAdd(iRow, value);
989 } else if (value > upper + tolerance) {
990 value -= upper;
991 value *= value;
992#ifdef CLP_DUAL_COLUMN_MULTIPLIER
993 if (iPivot < numberColumns)
994 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
995#endif
996#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
997 if (lower == upper)
998 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
999#endif
1000 // store square in list
1001 infeasible_->quickAdd(iRow, value);
1002 }
1003 }
1004 }
1005}
1006// Gets rid of last update
1007void
1008ClpDualRowSteepest::unrollWeights()
1009{
1010 double * saved = alternateWeights_->denseVector();
1011 int number = alternateWeights_->getNumElements();
1012 int * which = alternateWeights_->getIndices();
1013 int i;
1014 if (alternateWeights_->packedMode()) {
1015 for (i = 0; i < number; i++) {
1016 int iRow = which[i];
1017 weights_[iRow] = saved[i];
1018 saved[i] = 0.0;
1019 }
1020 } else {
1021 for (i = 0; i < number; i++) {
1022 int iRow = which[i];
1023 weights_[iRow] = saved[iRow];
1024 saved[iRow] = 0.0;
1025 }
1026 }
1027 alternateWeights_->setNumElements(0);
1028}
1029//-------------------------------------------------------------------
1030// Clone
1031//-------------------------------------------------------------------
1032ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
1033{
1034 if (CopyData) {
1035 return new ClpDualRowSteepest(*this);
1036 } else {
1037 return new ClpDualRowSteepest();
1038 }
1039}
1040// Gets rid of all arrays
1041void
1042ClpDualRowSteepest::clearArrays()
1043{
1044 if (persistence_ == normal) {
1045 delete [] weights_;
1046 weights_ = NULL;
1047 delete [] dubiousWeights_;
1048 dubiousWeights_ = NULL;
1049 delete infeasible_;
1050 infeasible_ = NULL;
1051 delete alternateWeights_;
1052 alternateWeights_ = NULL;
1053 delete savedWeights_;
1054 savedWeights_ = NULL;
1055 }
1056 state_ = -1;
1057}
1058// Returns true if would not find any row
1059bool
1060ClpDualRowSteepest::looksOptimal() const
1061{
1062 int iRow;
1063 const int * pivotVariable = model_->pivotVariable();
1064 double tolerance = model_->currentPrimalTolerance();
1065 // we can't really trust infeasibilities if there is primal error
1066 // this coding has to mimic coding in checkPrimalSolution
1067 double error = CoinMin(1.0e-2, model_->largestPrimalError());
1068 // allow tolerance at least slightly bigger than standard
1069 tolerance = tolerance + error;
1070 // But cap
1071 tolerance = CoinMin(1000.0, tolerance);
1072 int numberRows = model_->numberRows();
1073 int numberInfeasible = 0;
1074 for (iRow = 0; iRow < numberRows; iRow++) {
1075 int iPivot = pivotVariable[iRow];
1076 double value = model_->solution(iPivot);
1077 double lower = model_->lower(iPivot);
1078 double upper = model_->upper(iPivot);
1079 if (value < lower - tolerance) {
1080 numberInfeasible++;
1081 } else if (value > upper + tolerance) {
1082 numberInfeasible++;
1083 }
1084 }
1085 return (numberInfeasible == 0);
1086}
1087// Called when maximum pivots changes
1088void
1089ClpDualRowSteepest::maximumPivotsChanged()
1090{
1091 if (alternateWeights_ &&
1092 alternateWeights_->capacity() != model_->numberRows() +
1093 model_->factorization()->maximumPivots()) {
1094 delete alternateWeights_;
1095 alternateWeights_ = new CoinIndexedVector();
1096 // enough space so can use it for factorization
1097 alternateWeights_->reserve(model_->numberRows() +
1098 model_->factorization()->maximumPivots());
1099 }
1100}
1101