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

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