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

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

formatting
  • Property svn:keywords set to Id
File size: 13.0 KB
Line 
1/* $Id: AbcDualRowDantzig.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include "AbcSimplex.hpp"
8#include "AbcDualRowDantzig.hpp"
9#include "AbcSimplexFactorization.hpp"
10#include "CoinIndexedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12#ifndef CLP_DUAL_COLUMN_MULTIPLIER
13#define CLP_DUAL_COLUMN_MULTIPLIER 1.01
14#endif
15
16//#############################################################################
17// Constructors / Destructor / Assignment
18//#############################################################################
19
20//-------------------------------------------------------------------
21// Default Constructor
22//-------------------------------------------------------------------
23AbcDualRowDantzig::AbcDualRowDantzig()
24  : AbcDualRowPivot()
25  , infeasible_(NULL)
26{
27  type_ = 1;
28}
29
30//-------------------------------------------------------------------
31// Copy constructor
32//-------------------------------------------------------------------
33AbcDualRowDantzig::AbcDualRowDantzig(const AbcDualRowDantzig &rhs)
34  : AbcDualRowPivot(rhs)
35  , infeasible_(NULL)
36{
37  model_ = rhs.model_;
38  if ((model_ && model_->whatsChanged() & 1) != 0) {
39    if (rhs.infeasible_) {
40      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
41    } else {
42      infeasible_ = NULL;
43    }
44  }
45}
46
47//-------------------------------------------------------------------
48// Destructor
49//-------------------------------------------------------------------
50AbcDualRowDantzig::~AbcDualRowDantzig()
51{
52  delete infeasible_;
53}
54
55//----------------------------------------------------------------
56// Assignment operator
57//-------------------------------------------------------------------
58AbcDualRowDantzig &
59AbcDualRowDantzig::operator=(const AbcDualRowDantzig &rhs)
60{
61  if (this != &rhs) {
62    AbcDualRowPivot::operator=(rhs);
63    delete infeasible_;
64    if (rhs.infeasible_ != NULL) {
65      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
66    } else {
67      infeasible_ = NULL;
68    }
69  }
70  return *this;
71}
72
73/*
74   1) before factorization
75   2) after factorization
76   3) just redo infeasibilities
77   4) restore weights
78*/
79void AbcDualRowDantzig::saveWeights(AbcSimplex *model, int mode)
80{
81  model_ = model;
82  int numberRows = model_->numberRows();
83  if (mode == 1) {
84    // Check if size has changed
85    if (infeasible_ && infeasible_->capacity() != numberRows) {
86      // size has changed - clear everything
87      delete infeasible_;
88      infeasible_ = NULL;
89    }
90  } else if (mode != 3 && !infeasible_) {
91    infeasible_ = new CoinIndexedVector();
92    infeasible_->reserve(numberRows);
93  }
94  if (mode >= 2) {
95    recomputeInfeasibilities();
96  }
97}
98// Recompute infeasibilities
99void AbcDualRowDantzig::recomputeInfeasibilities()
100{
101  int numberRows = model_->numberRows();
102  infeasible_->clear();
103  double tolerance = model_->currentPrimalTolerance();
104  const double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
105  const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
106  const double *COIN_RESTRICT upperBasic = model_->upperBasic();
107  for (int iRow = 0; iRow < numberRows; iRow++) {
108    double value = solutionBasic[iRow];
109    double lower = lowerBasic[iRow];
110    double upper = upperBasic[iRow];
111    if (value < lower - tolerance) {
112      value -= lower;
113#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
114      if (lower == upper)
115        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
116#endif
117      // store in list
118      infeasible_->quickAdd(iRow, fabs(value));
119    } else if (value > upper + tolerance) {
120      value -= upper;
121#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
122      if (lower == upper)
123        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
124#endif
125      // store in list
126      infeasible_->quickAdd(iRow, fabs(value));
127    }
128  }
129}
130#if ABC_PARALLEL == 2
131static void choose(CoinIndexedVector *infeasible,
132  int &chosenRowSave, double &largestSave, int first, int last,
133  double tolerance)
134{
135  if (last - first > 256) {
136    int mid = (last + first) >> 1;
137    int chosenRow2 = chosenRowSave;
138    double largest2 = largestSave;
139    cilk_spawn choose(infeasible, chosenRow2, largest2, first, mid,
140      tolerance);
141    choose(infeasible, chosenRowSave, largestSave, mid, last,
142      tolerance);
143    cilk_sync;
144    if (largest2 > largestSave) {
145      largestSave = largest2;
146      chosenRowSave = chosenRow2;
147    }
148  } else {
149    const int *index = infeasible->getIndices();
150    const double *infeas = infeasible->denseVector();
151    double largest = largestSave;
152    int chosenRow = chosenRowSave;
153    for (int i = first; i < last; i++) {
154      int iRow = index[i];
155      double value = infeas[iRow];
156      if (value > largest) {
157        largest = value;
158        chosenRow = iRow;
159      }
160    }
161    chosenRowSave = chosenRow;
162    largestSave = largest;
163  }
164}
165#endif
166// Returns pivot row, -1 if none
167int AbcDualRowDantzig::pivotRow()
168{
169  assert(model_);
170  double *COIN_RESTRICT infeas = infeasible_->denseVector();
171  int *COIN_RESTRICT index = infeasible_->getIndices();
172  int number = infeasible_->getNumElements();
173  double tolerance = model_->currentPrimalTolerance();
174  // we can't really trust infeasibilities if there is primal error
175  if (model_->largestPrimalError() > 1.0e-8)
176    tolerance *= model_->largestPrimalError() / 1.0e-8;
177  int numberRows = model_->numberRows();
178  const double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
179  const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
180  const double *COIN_RESTRICT upperBasic = model_->upperBasic();
181  const int *pivotVariable = model_->pivotVariable();
182  // code so has list of infeasibilities (like steepest)
183  int numberWanted = CoinMax(2000, numberRows >> 4);
184  numberWanted = CoinMax(numberWanted, number >> 2);
185  if (model_->largestPrimalError() > 1.0e-3)
186    numberWanted = number + 1; // be safe
187  // Setup two passes
188  int start[4];
189  start[1] = number;
190  start[2] = 0;
191  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
192  start[0] = static_cast< int >(dstart);
193  start[3] = start[0];
194  double largest = tolerance;
195  int chosenRow = -1;
196  int saveNumberWanted = numberWanted;
197#ifdef DO_REDUCE
198  bool doReduce = true;
199  int lastChosen = -1;
200  double lastLargest = 0.0;
201#endif
202  for (int iPass = 0; iPass < 2; iPass++) {
203    int endThis = start[2 * iPass + 1];
204    int startThis = start[2 * iPass];
205    while (startThis < endThis) {
206      int end = CoinMin(endThis, startThis + numberWanted);
207#ifdef DO_REDUCE
208      if (doReduce) {
209        choose(infeasible, chosenRow, largest, startThis, end, tolerance);
210        if (chosenRow != lastChosen) {
211          assert(chosenRow >= 0);
212          if (model_->flagged(pivotVariable[chosenRow]) || (solutionBasic[chosenRow] <= upperBasic[chosenRow] + tolerance && solutionBasic[chosenRow] >= lowerBasic[chosenRow] - tolerance)) {
213            doReduce = false;
214            chosenRow = lastChosen;
215            largest = lastLargest;
216          } else {
217            lastChosen = chosenRow;
218            lastLargest = largest;
219          }
220        }
221      }
222      if (!doReduce) {
223#endif
224        for (int i = startThis; i < end; i++) {
225          int iRow = index[i];
226          double value = infeas[iRow];
227          if (value > largest) {
228            if (!model_->flagged(pivotVariable[iRow])) {
229              if (solutionBasic[iRow] > upperBasic[iRow] + tolerance || solutionBasic[iRow] < lowerBasic[iRow] - tolerance) {
230                chosenRow = iRow;
231                largest = value;
232              }
233            }
234          }
235        }
236#ifdef DO_REDUCE
237      }
238#endif
239      numberWanted -= (end - startThis);
240      if (!numberWanted) {
241        if (chosenRow >= 0)
242          break;
243        else
244          numberWanted = (saveNumberWanted + 1) >> 1;
245      }
246      startThis = end;
247    }
248    if (!numberWanted) {
249      if (chosenRow >= 0)
250        break;
251      else
252        numberWanted = (saveNumberWanted + 1) >> 1;
253    }
254  }
255  return chosenRow;
256}
257// FT update and returns pivot alpha
258double
259AbcDualRowDantzig::updateWeights(CoinIndexedVector &input, CoinIndexedVector &updatedColumn)
260{
261  // Do FT update
262  model_->factorization()->updateColumnFT(updatedColumn);
263  // pivot element
264  double alpha = 0.0;
265  // look at updated column
266  double *work = updatedColumn.denseVector();
267  int pivotRow = model_->lastPivotRow();
268  assert(pivotRow == model_->pivotRow());
269
270  assert(!updatedColumn.packedMode());
271  alpha = work[pivotRow];
272  return alpha;
273}
274double
275AbcDualRowDantzig::updateWeights1(CoinIndexedVector &input, CoinIndexedVector &updateColumn)
276{
277  return updateWeights(input, updateColumn);
278}
279#if ABC_PARALLEL == 2
280static void update(int first, int last,
281  const int *COIN_RESTRICT which, double *COIN_RESTRICT work,
282  const double *COIN_RESTRICT lowerBasic, double *COIN_RESTRICT solutionBasic,
283  const double *COIN_RESTRICT upperBasic, double theta, double tolerance)
284{
285  if (last - first > 256) {
286    int mid = (last + first) >> 1;
287    cilk_spawn update(first, mid, which, work, lowerBasic, solutionBasic,
288      upperBasic, theta, tolerance);
289    update(mid, last, which, work, lowerBasic, solutionBasic,
290      upperBasic, theta, tolerance);
291    cilk_sync;
292  } else {
293    for (int i = first; i < last; i++) {
294      int iRow = which[i];
295      double updateValue = work[iRow];
296      double value = solutionBasic[iRow];
297      double change = theta * updateValue;
298      value -= change;
299      double lower = lowerBasic[iRow];
300      double upper = upperBasic[iRow];
301      solutionBasic[iRow] = value;
302      if (value < lower - tolerance) {
303        value -= lower;
304#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
305        if (lower == upper)
306          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
307#endif
308      } else if (value > upper + tolerance) {
309        value -= upper;
310#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
311        if (lower == upper)
312          value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
313#endif
314      } else {
315        // feasible
316        value = 0.0;
317      }
318      // store
319      work[iRow] = fabs(value);
320    }
321  }
322}
323#endif
324/* Updates primal solution (and maybe list of candidates)
325   Uses input vector which it deletes
326   Computes change in objective function
327*/
328void AbcDualRowDantzig::updatePrimalSolution(CoinIndexedVector &primalUpdate,
329  double theta)
330{
331  double *COIN_RESTRICT work = primalUpdate.denseVector();
332  int numberNonZero = primalUpdate.getNumElements();
333  int *which = primalUpdate.getIndices();
334  double tolerance = model_->currentPrimalTolerance();
335  double *COIN_RESTRICT infeas = infeasible_->denseVector();
336  double *COIN_RESTRICT solutionBasic = model_->solutionBasic();
337  const double *COIN_RESTRICT lowerBasic = model_->lowerBasic();
338  const double *COIN_RESTRICT upperBasic = model_->upperBasic();
339  assert(!primalUpdate.packedMode());
340#if 0 //ABC_PARALLEL==2
341  update(0,numberNonZero,which,work,
342         lowerBasic,solutionBasic,upperBasic,
343         theta,tolerance);
344  for (int i = 0; i < numberNonZero; i++) {
345    int iRow = which[i];
346    double infeasibility=work[iRow];
347    work[iRow]=0.0;
348    if (infeasibility) {
349      if (infeas[iRow])
350        infeas[iRow] = infeasibility; // already there
351      else
352        infeasible_->quickAdd(iRow, infeasibility);
353    } else {
354      // feasible - was it infeasible - if so set tiny
355      if (infeas[iRow])
356        infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
357    }
358  }
359#else
360  for (int i = 0; i < numberNonZero; i++) {
361    int iRow = which[i];
362    double updateValue = work[iRow];
363    work[iRow] = 0.0;
364    double value = solutionBasic[iRow];
365    double change = theta * updateValue;
366    value -= change;
367    double lower = lowerBasic[iRow];
368    double upper = upperBasic[iRow];
369    solutionBasic[iRow] = value;
370    if (value < lower - tolerance) {
371      value -= lower;
372#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
373      if (lower == upper)
374        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
375#endif
376      // store in list
377      if (infeas[iRow])
378        infeas[iRow] = fabs(value); // already there
379      else
380        infeasible_->quickAdd(iRow, fabs(value));
381    } else if (value > upper + tolerance) {
382      value -= upper;
383#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
384      if (lower == upper)
385        value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
386#endif
387      // store in list
388      if (infeas[iRow])
389        infeas[iRow] = fabs(value); // already there
390      else
391        infeasible_->quickAdd(iRow, fabs(value));
392    } else {
393      // feasible - was it infeasible - if so set tiny
394      if (infeas[iRow])
395        infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
396    }
397  }
398#endif
399  primalUpdate.setNumElements(0);
400}
401//-------------------------------------------------------------------
402// Clone
403//-------------------------------------------------------------------
404AbcDualRowPivot *AbcDualRowDantzig::clone(bool CopyData) const
405{
406  if (CopyData) {
407    return new AbcDualRowDantzig(*this);
408  } else {
409    return new AbcDualRowDantzig();
410  }
411}
412
413/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
414*/
Note: See TracBrowser for help on using the repository browser.