source: trunk/Clp/src/AbcPrimalColumnDantzig.cpp

Last change on this file was 2385, checked in by unxusr, 10 months ago

formatting

  • Property svn:keywords set to Id
File size: 8.3 KB
Line 
1/* $Id: AbcPrimalColumnDantzig.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
8#include <cstdio>
9
10#include "CoinIndexedVector.hpp"
11
12#include "AbcSimplex.hpp"
13#include "AbcPrimalColumnDantzig.hpp"
14#include "AbcSimplexFactorization.hpp"
15#include "AbcMatrix.hpp"
16
17//#############################################################################
18// Constructors / Destructor / Assignment
19//#############################################################################
20
21//-------------------------------------------------------------------
22// Default Constructor
23//-------------------------------------------------------------------
24AbcPrimalColumnDantzig::AbcPrimalColumnDantzig()
25  : AbcPrimalColumnPivot()
26{
27  type_ = 1;
28}
29
30//-------------------------------------------------------------------
31// Copy constructor
32//-------------------------------------------------------------------
33AbcPrimalColumnDantzig::AbcPrimalColumnDantzig(const AbcPrimalColumnDantzig &source)
34  : AbcPrimalColumnPivot(source)
35{
36}
37
38//-------------------------------------------------------------------
39// Destructor
40//-------------------------------------------------------------------
41AbcPrimalColumnDantzig::~AbcPrimalColumnDantzig()
42{
43}
44
45//----------------------------------------------------------------
46// Assignment operator
47//-------------------------------------------------------------------
48AbcPrimalColumnDantzig &
49AbcPrimalColumnDantzig::operator=(const AbcPrimalColumnDantzig &rhs)
50{
51  if (this != &rhs) {
52    AbcPrimalColumnPivot::operator=(rhs);
53  }
54  return *this;
55}
56// Returns pivot column, -1 if none
57int AbcPrimalColumnDantzig::pivotColumn(CoinPartitionedVector *updates,
58  CoinPartitionedVector * /*spareRow2*/,
59  CoinPartitionedVector *spareColumn1)
60{
61  assert(model_);
62  int iSection, j;
63  int number;
64  double multiplier;
65  int *index;
66  double *updateBy;
67  double *reducedCost;
68
69  bool anyUpdates;
70
71  if (updates->getNumElements()) {
72    anyUpdates = true;
73  } else {
74    // sub flip - nothing to do
75    anyUpdates = false;
76  }
77  if (anyUpdates) {
78    model_->factorization()->updateColumnTranspose(*updates);
79    int iVector = model_->getAvailableArray();
80    int bestSequence = model_->abcMatrix()->pivotColumnDantzig(*updates, *model_->usefulArray(iVector));
81    model_->setAvailableArray(iVector);
82    int pivotRow = model_->pivotRow();
83    if (pivotRow >= 0) {
84      // make sure infeasibility on incoming is 0.0
85      int sequenceIn = model_->sequenceIn();
86      double *reducedCost = model_->djRegion();
87      reducedCost[sequenceIn] = 0.0;
88    }
89#if 1
90    if (model_->logLevel() == 127) {
91      double *reducedCost = model_->djRegion();
92      int numberRows = model_->numberRows();
93      printf("Best sequence %d\n", bestSequence);
94      for (int i = 0; i < numberRows; i++)
95        printf("row %d dj %g\n", i, reducedCost[i]);
96      for (int i = 0; i < model_->numberColumns(); i++)
97        printf("column %d dj %g\n", i, reducedCost[i + numberRows]);
98    }
99#endif
100    looksOptimal_ = bestSequence < 0;
101    return bestSequence;
102    // put row of tableau in rowArray and columnArray
103    model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
104    for (iSection = 0; iSection < 2; iSection++) {
105
106      reducedCost = model_->djRegion(iSection);
107
108      if (!iSection) {
109        number = updates->getNumElements();
110        index = updates->getIndices();
111        updateBy = updates->denseVector();
112        multiplier = -1.0;
113      } else {
114        number = spareColumn1->getNumElements();
115        index = spareColumn1->getIndices();
116        updateBy = spareColumn1->denseVector();
117        multiplier = 1.0;
118      }
119
120      for (j = 0; j < number; j++) {
121        int iSequence = index[j];
122        double value = reducedCost[iSequence];
123        value -= multiplier * updateBy[iSequence];
124        updateBy[iSequence] = 0.0;
125        reducedCost[iSequence] = value;
126      }
127    }
128    updates->setNumElements(0);
129    spareColumn1->setNumElements(0);
130  }
131
132  // update of duals finished - now do pricing
133#if 0
134     double largest = model_->currentPrimalTolerance();
135     // we can't really trust infeasibilities if there is primal error
136     if (model_->largestDualError() > 1.0e-8)
137          largest *= model_->largestDualError() / 1.0e-8;
138#endif
139
140  double bestDj = model_->dualTolerance();
141  int bestSequence = -1;
142
143  double bestFreeDj = model_->dualTolerance();
144  int bestFreeSequence = -1;
145
146  number = model_->numberTotal();
147  int iSequence;
148  reducedCost = model_->djRegion();
149
150#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
151  for (iSequence = 0; iSequence < number; iSequence++) {
152    // check flagged variable
153    if (!model_->flagged(iSequence)) {
154      double value = reducedCost[iSequence];
155      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
156
157      switch (status) {
158
159      case AbcSimplex::basic:
160      case AbcSimplex::isFixed:
161        break;
162      case AbcSimplex::isFree:
163      case AbcSimplex::superBasic:
164        if (fabs(value) > bestFreeDj) {
165          bestFreeDj = fabs(value);
166          bestFreeSequence = iSequence;
167        }
168        break;
169      case AbcSimplex::atUpperBound:
170        if (value > bestDj) {
171          bestDj = value;
172          bestSequence = iSequence;
173        }
174        break;
175      case AbcSimplex::atLowerBound:
176        if (value < -bestDj) {
177          bestDj = -value;
178          bestSequence = iSequence;
179        }
180      }
181    }
182  }
183#else
184  int numberColumns = model_->numberColumns();
185  int maximumRows = model_->maximumAbcNumberRows();
186  for (iSequence = maximumRows; iSequence < numberColumns + maximumRows; iSequence++) {
187    // check flagged variable
188    if (!model_->flagged(iSequence)) {
189      double value = reducedCost[iSequence];
190      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
191
192      switch (status) {
193
194      case AbcSimplex::basic:
195      case AbcSimplex::isFixed:
196        break;
197      case AbcSimplex::isFree:
198      case AbcSimplex::superBasic:
199        if (fabs(value) > bestFreeDj) {
200          bestFreeDj = fabs(value);
201          bestFreeSequence = iSequence;
202        }
203        break;
204      case AbcSimplex::atUpperBound:
205        if (value > bestDj) {
206          bestDj = value;
207          bestSequence = iSequence;
208        }
209        break;
210      case AbcSimplex::atLowerBound:
211        if (value < -bestDj) {
212          bestDj = -value;
213          bestSequence = iSequence;
214        }
215      }
216    }
217  }
218  // Rows
219  number = model_->numberRows();
220  for (iSequence = 0; iSequence < number; iSequence++) {
221    // check flagged variable
222    if (!model_->flagged(iSequence)) {
223      double value = reducedCost[iSequence] * CLP_PRIMAL_SLACK_MULTIPLIER;
224      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
225
226      switch (status) {
227
228      case AbcSimplex::basic:
229      case AbcSimplex::isFixed:
230        break;
231      case AbcSimplex::isFree:
232      case AbcSimplex::superBasic:
233        if (fabs(value) > bestFreeDj) {
234          bestFreeDj = fabs(value);
235          bestFreeSequence = iSequence;
236        }
237        break;
238      case AbcSimplex::atUpperBound:
239        if (value > bestDj) {
240          bestDj = value;
241          bestSequence = iSequence;
242        }
243        break;
244      case AbcSimplex::atLowerBound:
245        if (value < -bestDj) {
246          bestDj = -value;
247          bestSequence = iSequence;
248        }
249      }
250    }
251  }
252#endif
253  // bias towards free
254  if (bestFreeSequence >= 0 && bestFreeDj > 0.1 * bestDj)
255    bestSequence = bestFreeSequence;
256#if 1
257  if (model_->logLevel() == 127) {
258    double *reducedCost = model_->djRegion();
259    int numberRows = model_->numberRows();
260    printf("Best sequence %d\n", bestSequence);
261    for (int i = 0; i < numberRows; i++)
262      printf("row %d dj %g\n", i, reducedCost[i]);
263    for (int i = 0; i < model_->numberColumns(); i++)
264      printf("column %d dj %g\n", i, reducedCost[i + numberRows]);
265  }
266#endif
267  looksOptimal_ = bestSequence < 0;
268  return bestSequence;
269}
270
271//-------------------------------------------------------------------
272// Clone
273//-------------------------------------------------------------------
274AbcPrimalColumnPivot *AbcPrimalColumnDantzig::clone(bool CopyData) const
275{
276  if (CopyData) {
277    return new AbcPrimalColumnDantzig(*this);
278  } else {
279    return new AbcPrimalColumnDantzig();
280  }
281}
282
283/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
284*/
Note: See TracBrowser for help on using the repository browser.