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

Last change on this file since 1878 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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