source: trunk/Clp/src/ClpDualRowSteepest.cpp @ 2030

Last change on this file since 2030 was 1732, checked in by forrest, 8 years ago

various fixes plus slightly weighted pricing plus lagomory switches

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 42.2 KB
Line 
1/* $Id: ClpDualRowSteepest.cpp 1732 2011-05-31 08:09:41Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#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
Note: See TracBrowser for help on using the repository browser.