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

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