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

Last change on this file since 2470 was 2385, checked in by unxusr, 11 months ago

formatting

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