1/* $Id: ClpPrimalColumnDantzig.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others. 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 "ClpSimplex.hpp"
13#include "ClpPrimalColumnDantzig.hpp"
14#include "ClpFactorization.hpp"
15#include "ClpPackedMatrix.hpp"
16
17//#############################################################################
18// Constructors / Destructor / Assignment
19//#############################################################################
20
21//-------------------------------------------------------------------
22// Default Constructor
23//-------------------------------------------------------------------
24ClpPrimalColumnDantzig::ClpPrimalColumnDantzig ()
25 : ClpPrimalColumnPivot()
26{
27 type_ = 1;
28}
29
30//-------------------------------------------------------------------
31// Copy constructor
32//-------------------------------------------------------------------
33ClpPrimalColumnDantzig::ClpPrimalColumnDantzig (const ClpPrimalColumnDantzig & source)
34 : ClpPrimalColumnPivot(source)
35{
36
37}
38
39//-------------------------------------------------------------------
40// Destructor
41//-------------------------------------------------------------------
42ClpPrimalColumnDantzig::~ClpPrimalColumnDantzig ()
43{
44
45}
46
47//----------------------------------------------------------------
48// Assignment operator
49//-------------------------------------------------------------------
50ClpPrimalColumnDantzig &
51ClpPrimalColumnDantzig::operator=(const ClpPrimalColumnDantzig& rhs)
52{
53 if (this != &rhs) {
54 ClpPrimalColumnPivot::operator=(rhs);
55 }
56 return *this;
57}
58
59// Returns pivot column, -1 if none
60int
61ClpPrimalColumnDantzig::pivotColumn(CoinIndexedVector * updates,
62 CoinIndexedVector * /*spareRow1*/,
63 CoinIndexedVector * spareRow2,
64 CoinIndexedVector * spareColumn1,
65 CoinIndexedVector * spareColumn2)
66{
67 assert(model_);
68 int iSection, j;
69 int number;
70 int * index;
71 double * updateBy;
72 double * reducedCost;
73
74 bool anyUpdates;
75
76 if (updates->getNumElements()) {
77 anyUpdates = true;
78 } else {
79 // sub flip - nothing to do
80 anyUpdates = false;
81 }
82 if (anyUpdates) {
83 model_->factorization()->updateColumnTranspose(spareRow2, updates);
84 // put row of tableau in rowArray and columnArray
85 model_->clpMatrix()->transposeTimes(model_, -1.0,
86 updates, spareColumn2, spareColumn1);
87 for (iSection = 0; iSection < 2; iSection++) {
88
89 reducedCost = model_->djRegion(iSection);
90
91 if (!iSection) {
92 number = updates->getNumElements();
93 index = updates->getIndices();
94 updateBy = updates->denseVector();
95 } else {
96 number = spareColumn1->getNumElements();
97 index = spareColumn1->getIndices();
98 updateBy = spareColumn1->denseVector();
99 }
100
101 for (j = 0; j < number; j++) {
102 int iSequence = index[j];
103 double value = reducedCost[iSequence];
104 value -= updateBy[j];
105 updateBy[j] = 0.0;
106 reducedCost[iSequence] = value;
107 }
108
109 }
110 updates->setNumElements(0);
111 spareColumn1->setNumElements(0);
112 }
113
114
115 // update of duals finished - now do pricing
116
117 double largest = model_->currentPrimalTolerance();
118 // we can't really trust infeasibilities if there is primal error
119 if (model_->largestDualError() > 1.0e-8)
120 largest *= model_->largestDualError() / 1.0e-8;
121
122
123
124 double bestDj = model_->dualTolerance();
125 int bestSequence = -1;
126
127 double bestFreeDj = model_->dualTolerance();
128 int bestFreeSequence = -1;
129
130 number = model_->numberRows() + model_->numberColumns();
131 int iSequence;
132 reducedCost = model_->djRegion();
133
134#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
135 for (iSequence = 0; iSequence < number; iSequence++) {
136 // check flagged variable
137 if (!model_->flagged(iSequence)) {
138 double value = reducedCost[iSequence];
139 ClpSimplex::Status status = model_->getStatus(iSequence);
140
141 switch(status) {
142
143 case ClpSimplex::basic:
144 case ClpSimplex::isFixed:
145 break;
146 case ClpSimplex::isFree:
147 case ClpSimplex::superBasic:
148 if (fabs(value) > bestFreeDj) {
149 bestFreeDj = fabs(value);
150 bestFreeSequence = iSequence;
151 }
152 break;
153 case ClpSimplex::atUpperBound:
154 if (value > bestDj) {
155 bestDj = value;
156 bestSequence = iSequence;
157 }
158 break;
159 case ClpSimplex::atLowerBound:
160 if (value < -bestDj) {
161 bestDj = -value;
162 bestSequence = iSequence;
163 }
164 }
165 }
166 }
167#else
168 // Columns
169 int numberColumns = model_->numberColumns();
170 for (iSequence = 0; iSequence < numberColumns; iSequence++) {
171 // check flagged variable
172 if (!model_->flagged(iSequence)) {
173 double value = reducedCost[iSequence];
174 ClpSimplex::Status status = model_->getStatus(iSequence);
175
176 switch(status) {
177
178 case ClpSimplex::basic:
179 case ClpSimplex::isFixed:
180 break;
181 case ClpSimplex::isFree:
182 case ClpSimplex::superBasic:
183 if (fabs(value) > bestFreeDj) {
184 bestFreeDj = fabs(value);
185 bestFreeSequence = iSequence;
186 }
187 break;
188 case ClpSimplex::atUpperBound:
189 if (value > bestDj) {
190 bestDj = value;
191 bestSequence = iSequence;
192 }
193 break;
194 case ClpSimplex::atLowerBound:
195 if (value < -bestDj) {
196 bestDj = -value;
197 bestSequence = iSequence;
198 }
199 }
200 }
201 }
202 // Rows
203 for ( ; iSequence < number; iSequence++) {
204 // check flagged variable
205 if (!model_->flagged(iSequence)) {
206 double value = reducedCost[iSequence] * CLP_PRIMAL_SLACK_MULTIPLIER;
207 ClpSimplex::Status status = model_->getStatus(iSequence);
208
209 switch(status) {
210
211 case ClpSimplex::basic:
212 case ClpSimplex::isFixed:
213 break;
214 case ClpSimplex::isFree:
215 case ClpSimplex::superBasic:
216 if (fabs(value) > bestFreeDj) {
217 bestFreeDj = fabs(value);
218 bestFreeSequence = iSequence;
219 }
220 break;
221 case ClpSimplex::atUpperBound:
222 if (value > bestDj) {
223 bestDj = value;
224 bestSequence = iSequence;
225 }
226 break;
227 case ClpSimplex::atLowerBound:
228 if (value < -bestDj) {
229 bestDj = -value;
230 bestSequence = iSequence;
231 }
232 }
233 }
234 }
235#endif
236 // bias towards free
237 if (bestFreeSequence >= 0 && bestFreeDj > 0.1 * bestDj)
238 bestSequence = bestFreeSequence;
239 return bestSequence;
240}
241
242//-------------------------------------------------------------------
243// Clone
244//-------------------------------------------------------------------
245ClpPrimalColumnPivot * ClpPrimalColumnDantzig::clone(bool CopyData) const
246{
247 if (CopyData) {
248 return new ClpPrimalColumnDantzig(*this);
249 } else {
250 return new ClpPrimalColumnDantzig();
251 }
252}
253