1/* $Id: ClpDualRowDantzig.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#include "ClpSimplex.hpp"
8#include "ClpDualRowDantzig.hpp"
9#include "ClpFactorization.hpp"
10#include "CoinIndexedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12#ifndef CLP_DUAL_COLUMN_MULTIPLIER
13#define CLP_DUAL_COLUMN_MULTIPLIER 1.01
14#endif
15
16//#############################################################################
17// Constructors / Destructor / Assignment
18//#############################################################################
19
20//-------------------------------------------------------------------
21// Default Constructor
22//-------------------------------------------------------------------
23ClpDualRowDantzig::ClpDualRowDantzig ()
24 : ClpDualRowPivot()
25{
26 type_ = 1;
27}
28
29//-------------------------------------------------------------------
30// Copy constructor
31//-------------------------------------------------------------------
32ClpDualRowDantzig::ClpDualRowDantzig (const ClpDualRowDantzig & source)
33 : ClpDualRowPivot(source)
34{
35
36}
37
38//-------------------------------------------------------------------
39// Destructor
40//-------------------------------------------------------------------
41ClpDualRowDantzig::~ClpDualRowDantzig ()
42{
43
44}
45
46//----------------------------------------------------------------
47// Assignment operator
48//-------------------------------------------------------------------
49ClpDualRowDantzig &
50ClpDualRowDantzig::operator=(const ClpDualRowDantzig& rhs)
51{
52 if (this != &rhs) {
53 ClpDualRowPivot::operator=(rhs);
54 }
55 return *this;
56}
57
58// Returns pivot row, -1 if none
59int
60ClpDualRowDantzig::pivotRow()
61{
62 assert(model_);
63 int iRow;
64 const int * pivotVariable = model_->pivotVariable();
65 double tolerance = model_->currentPrimalTolerance();
66 // we can't really trust infeasibilities if there is primal error
67 if (model_->largestPrimalError() > 1.0e-8)
68 tolerance *= model_->largestPrimalError() / 1.0e-8;
69 double largest = 0.0;
70 int chosenRow = -1;
71 int numberRows = model_->numberRows();
72#ifdef CLP_DUAL_COLUMN_MULTIPLIER
73 int numberColumns = model_->numberColumns();
74#endif
75 for (iRow = 0; iRow < numberRows; iRow++) {
76 int iSequence = pivotVariable[iRow];
77 double value = model_->solution(iSequence);
78 double lower = model_->lower(iSequence);
79 double upper = model_->upper(iSequence);
80 double infeas = CoinMax(value - upper , lower - value);
81 if (infeas > tolerance) {
82#ifdef CLP_DUAL_COLUMN_MULTIPLIER
83 if (iSequence < numberColumns)
84 infeas *= CLP_DUAL_COLUMN_MULTIPLIER;
85#endif
86 if (infeas > largest) {
87 if (!model_->flagged(iSequence)) {
88 chosenRow = iRow;
89 largest = infeas;
90 }
91 }
92 }
93 }
94 return chosenRow;
95}
96// FT update and returns pivot alpha
97double
98ClpDualRowDantzig::updateWeights(CoinIndexedVector * /*input*/,
99 CoinIndexedVector * spare,
100 CoinIndexedVector * /*spare2*/,
101 CoinIndexedVector * updatedColumn)
102{
103 // Do FT update
104 model_->factorization()->updateColumnFT(spare, updatedColumn);
105 // pivot element
106 double alpha = 0.0;
107 // look at updated column
108 double * work = updatedColumn->denseVector();
109 int number = updatedColumn->getNumElements();
110 int * which = updatedColumn->getIndices();
111 int i;
112 int pivotRow = model_->pivotRow();
113
114 if (updatedColumn->packedMode()) {
115 for (i = 0; i < number; i++) {
116 int iRow = which[i];
117 if (iRow == pivotRow) {
118 alpha = work[i];
119 break;
120 }
121 }
122 } else {
123 alpha = work[pivotRow];
124 }
125 return alpha;
126}
127
128/* Updates primal solution (and maybe list of candidates)
129 Uses input vector which it deletes
130 Computes change in objective function
131*/
132void
133ClpDualRowDantzig::updatePrimalSolution(CoinIndexedVector * primalUpdate,
134 double primalRatio,
135 double & objectiveChange)
136{
137 double * work = primalUpdate->denseVector();
138 int number = primalUpdate->getNumElements();
139 int * which = primalUpdate->getIndices();
140 int i;
141 double changeObj = 0.0;
142 const int * pivotVariable = model_->pivotVariable();
143 if (primalUpdate->packedMode()) {
144 for (i = 0; i < number; i++) {
145 int iRow = which[i];
146 int iPivot = pivotVariable[iRow];
147 double & value = model_->solutionAddress(iPivot);
148 double cost = model_->cost(iPivot);
149 double change = primalRatio * work[i];
150 value -= change;
151 changeObj -= change * cost;
152 work[i] = 0.0;
153 }
154 } else {
155 for (i = 0; i < number; i++) {
156 int iRow = which[i];
157 int iPivot = pivotVariable[iRow];
158 double & value = model_->solutionAddress(iPivot);
159 double cost = model_->cost(iPivot);
160 double change = primalRatio * work[iRow];
161 value -= change;
162 changeObj -= change * cost;
163 work[iRow] = 0.0;
164 }
165 }
166 primalUpdate->setNumElements(0);
167 objectiveChange += changeObj;
168}
169//-------------------------------------------------------------------
170// Clone
171//-------------------------------------------------------------------
172ClpDualRowPivot * ClpDualRowDantzig::clone(bool CopyData) const
173{
174 if (CopyData) {
175 return new ClpDualRowDantzig(*this);
176 } else {
177 return new ClpDualRowDantzig();
178 }
179}
180