1/* $Id: ClpSimplexPrimal.hpp 1665 2011-01-04 17:55:54Z lou $ */
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 Authors
7
8 John Forrest
9
10 */
11#ifndef ClpSimplexPrimal_H
12#define ClpSimplexPrimal_H
13
14#include "ClpSimplex.hpp"
15
16/** This solves LPs using the primal simplex method
17
18 It inherits from ClpSimplex. It has no data of its own and
19 is never created - only cast from a ClpSimplex object at algorithm time.
20
21*/
22
23class ClpSimplexPrimal : public ClpSimplex {
24
25public:
26
27 /**@name Description of algorithm */
28 //@{
29 /** Primal algorithm
30
31 Method
32
33 It tries to be a single phase approach with a weight of 1.0 being
34 given to getting optimal and a weight of infeasibilityCost_ being
35 given to getting primal feasible. In this version I have tried to
36 be clever in a stupid way. The idea of fake bounds in dual
37 seems to work so the primal analogue would be that of getting
38 bounds on reduced costs (by a presolve approach) and using
39 these for being above or below feasible region. I decided to waste
40 memory and keep these explicitly. This allows for non-linear
41 costs! I have not tested non-linear costs but will be glad
42 to do something if a reasonable example is provided.
43
44 The code is designed to take advantage of sparsity so arrays are
45 seldom zeroed out from scratch or gone over in their entirety.
46 The only exception is a full scan to find incoming variable for
47 Dantzig row choice. For steepest edge we keep an updated list
48 of dual infeasibilities (actually squares).
49 On easy problems we don't need full scan - just
50 pick first reasonable. This method has not been coded.
51
52 One problem is how to tackle degeneracy and accuracy. At present
53 I am using the modification of costs which I put in OSL and which was
54 extended by Gill et al. I am still not sure whether we will also
55 need explicit perturbation.
56
57 The flow of primal is three while loops as follows:
58
59 while (not finished) {
60
61 while (not clean solution) {
62
63 Factorize and/or clean up solution by changing bounds so
64 primal feasible. If looks finished check fake primal bounds.
65 Repeat until status is iterating (-1) or finished (0,1,2)
66
67 }
68
69 while (status==-1) {
70
71 Iterate until no pivot in or out or time to re-factorize.
72
73 Flow is:
74
75 choose pivot column (incoming variable). if none then
76 we are primal feasible so looks as if done but we need to
77 break and check bounds etc.
78
79 Get pivot column in tableau
80
81 Choose outgoing row. If we don't find one then we look
82 primal unbounded so break and check bounds etc. (Also the
83 pivot tolerance is larger after any iterations so that may be
84 reason)
85
86 If we do find outgoing row, we may have to adjust costs to
87 keep going forwards (anti-degeneracy). Check pivot will be stable
88 and if unstable throw away iteration and break to re-factorize.
89 If minor error re-factorize after iteration.
90
91 Update everything (this may involve changing bounds on
92 variables to stay primal feasible.
93
94 }
95
96 }
97
98 TODO's (or maybe not)
99
100 At present we never check we are going forwards. I overdid that in
101 OSL so will try and make a last resort.
102
103 Needs partial scan pivot in option.
104
105 May need other anti-degeneracy measures, especially if we try and use
106 loose tolerances as a way to solve in fewer iterations.
107
108 I like idea of dynamic scaling. This gives opportunity to decouple
109 different implications of scaling for accuracy, iteration count and
110 feasibility tolerance.
111
112 for use of exotic parameter startFinishoptions see Clpsimplex.hpp
113 */
114
115 int primal(int ifValuesPass = 0, int startFinishOptions = 0);
116 //@}
117
118 /**@name For advanced users */
119 //@{
120 /// Do not change infeasibility cost and always say optimal
121 void alwaysOptimal(bool onOff);
122 bool alwaysOptimal() const;
123 /** Normally outgoing variables can go out to slightly negative
124 values (but within tolerance) - this is to help stability and
125 and degeneracy. This can be switched off
126 */
127 void exactOutgoing(bool onOff);
128 bool exactOutgoing() const;
129 //@}
130
131 /**@name Functions used in primal */
132 //@{
133 /** This has the flow between re-factorizations
134
135 Returns a code to say where decision to exit was made
136 Problem status set to:
137
138 -2 re-factorize
139 -4 Looks optimal/infeasible
140 -5 Looks unbounded
141 +3 max iterations
142
143 valuesOption has original value of valuesPass
144 */
145 int whileIterating(int valuesOption);
146
147 /** Do last half of an iteration. This is split out so people can
148 force incoming variable. If solveType_ is 2 then this may
149 re-factorize while normally it would exit to re-factorize.
150 Return codes
151 Reasons to come out (normal mode/user mode):
152 -1 normal
153 -2 factorize now - good iteration/ NA
154 -3 slight inaccuracy - refactorize - iteration done/ same but factor done
155 -4 inaccuracy - refactorize - no iteration/ NA
156 -5 something flagged - go round again/ pivot not possible
157 +2 looks unbounded
158 +3 max iterations (iteration done)
159
160 With solveType_ ==2 this should
161 Pivot in a variable and choose an outgoing one. Assumes primal
162 feasible - will not go through a bound. Returns step length in theta
163 Returns ray in ray_
164 */
165 int pivotResult(int ifValuesPass = 0);
166
167
168 /** The primals are updated by the given array.
169 Returns number of infeasibilities.
170 After rowArray will have cost changes for use next iteration
171 */
172 int updatePrimalsInPrimal(CoinIndexedVector * rowArray,
173 double theta,
174 double & objectiveChange,
175 int valuesPass);
176 /**
177 Row array has pivot column
178 This chooses pivot row.
179 Rhs array is used for distance to next bound (for speed)
180 For speed, we may need to go to a bucket approach when many
181 variables go through bounds
182 If valuesPass non-zero then compute dj for direction
183 */
184 void primalRow(CoinIndexedVector * rowArray,
185 CoinIndexedVector * rhsArray,
186 CoinIndexedVector * spareArray,
187 int valuesPass);
188 /**
189 Chooses primal pivot column
190 updateArray has cost updates (also use pivotRow_ from last iteration)
191 Would be faster with separate region to scan
192 and will have this (with square of infeasibility) when steepest
193 For easy problems we can just choose one of the first columns we look at
194 */
195 void primalColumn(CoinIndexedVector * updateArray,
196 CoinIndexedVector * spareRow1,
197 CoinIndexedVector * spareRow2,
198 CoinIndexedVector * spareColumn1,
199 CoinIndexedVector * spareColumn2);
200
201 /** Checks if tentative optimal actually means unbounded in primal
202 Returns -3 if not, 2 if is unbounded */
203 int checkUnbounded(CoinIndexedVector * ray, CoinIndexedVector * spare,
204 double changeCost);
205 /** Refactorizes if necessary
206 Checks if finished. Updates status.
207 lastCleaned refers to iteration at which some objective/feasibility
208 cleaning too place.
209
210 type - 0 initial so set up save arrays etc
211 - 1 normal -if good update save
212 - 2 restoring from saved
213 saveModel is normally NULL but may not be if doing Sprint
214 */
215 void statusOfProblemInPrimal(int & lastCleaned, int type,
216 ClpSimplexProgress * progress,
217 bool doFactorization,
218 int ifValuesPass,
219 ClpSimplex * saveModel = nullptr);
220 /// Perturbs problem (method depends on perturbation())
221 void perturb(int type);
222 /// Take off effect of perturbation and say whether to try dual
223 bool unPerturb();
224 /// Unflag all variables and return number unflagged
225 int unflag();
226 /** Get next superbasic -1 if none,
227 Normal type is 1
228 If type is 3 then initializes sorted list
229 if 2 uses list.
230 */
231 int nextSuperBasic(int superBasicType, CoinIndexedVector * columnArray);
232
233 /// Create primal ray
234 void primalRay(CoinIndexedVector * rowArray);
235 /// Clears all bits and clears rowArray[1] etc
236 void clearAll();
237
238 /// Sort of lexicographic resolve
239 int lexSolve();
240
241 //@}
242};
243#endif
244