source: trunk/Clp/src/AbcSimplexPrimal.hpp @ 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: 10.0 KB
Line 
1/* $Id: AbcSimplexPrimal.hpp 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   Authors
7
8   John Forrest
9
10 */
11#ifndef AbcSimplexPrimal_H
12#define AbcSimplexPrimal_H
13
14#include "AbcSimplex.hpp"
15
16/** This solves LPs using the primal simplex method
17
18    It inherits from AbcSimplex.  It has no data of its own and
19    is never created - only cast from a AbcSimplex object at algorithm time.
20
21*/
22
23class AbcSimplexPrimal : public AbcSimplex {
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     int pivotResult4(int ifValuesPass = 0);
167
168
169     /** The primals are updated by the given array.
170         Returns number of infeasibilities.
171         After rowArray will have cost changes for use next iteration
172     */
173     int updatePrimalsInPrimal(CoinIndexedVector * rowArray,
174                               double theta,
175                               double & objectiveChange,
176                               int valuesPass);
177     /** The primals are updated by the given array.
178         costs are changed
179     */
180     void updatePrimalsInPrimal(CoinIndexedVector & rowArray,
181                                double theta,bool valuesPass);
182     /** After rowArray will have cost changes for use next iteration
183     */
184     void createUpdateDuals(CoinIndexedVector & rowArray,
185                           const double * originalCost,
186                            const double extraCost[4],
187                           double & objectiveChange,
188                           int valuesPass);
189     /** Update minor candidate vector - new reduced cost returned
190         later try and get change in reduced cost (then may not need sequence in)*/
191     double updateMinorCandidate(const CoinIndexedVector & updateBy,
192                                 CoinIndexedVector & candidate,
193                                 int sequenceIn);
194     /// Update partial Ftran by R update
195     void updatePartialUpdate(CoinIndexedVector & partialUpdate);
196  /// Do FT update as separate function for minor iterations (nonzero return code on problems)
197  int doFTUpdate(CoinIndexedVector * vector[4]);
198     /**
199         Row array has pivot column
200         This chooses pivot row.
201         Rhs array is used for distance to next bound (for speed)
202         For speed, we may need to go to a bucket approach when many
203         variables go through bounds
204         If valuesPass non-zero then compute dj for direction
205     */
206     void primalRow(CoinIndexedVector * rowArray,
207                    CoinIndexedVector * rhsArray,
208                    CoinIndexedVector * spareArray,
209                    int valuesPass);
210  typedef struct {
211    double theta_;
212    double alpha_;
213    double saveDualIn_;
214    double dualIn_;
215    double lowerIn_;
216    double upperIn_;
217    double valueIn_;
218    int sequenceIn_;
219    int directionIn_;
220    double dualOut_;
221    double lowerOut_;
222    double upperOut_;
223    double valueOut_;
224    int sequenceOut_;
225    int directionOut_;
226    int pivotRow_;
227    int valuesPass_;
228  } pivotStruct;
229     void primalRow(CoinIndexedVector * rowArray,
230                    CoinIndexedVector * rhsArray,
231                    CoinIndexedVector * spareArray,
232                    pivotStruct & stuff);
233     /**
234         Chooses primal pivot column
235         updateArray has cost updates (also use pivotRow_ from last iteration)
236         Would be faster with separate region to scan
237         and will have this (with square of infeasibility) when steepest
238         For easy problems we can just choose one of the first columns we look at
239     */
240     void primalColumn(CoinPartitionedVector * updateArray,
241                       CoinPartitionedVector * spareRow2,
242                       CoinPartitionedVector * spareColumn1);
243
244     /** Checks if tentative optimal actually means unbounded in primal
245         Returns -3 if not, 2 if is unbounded */
246     int checkUnbounded(CoinIndexedVector * ray, CoinIndexedVector * spare,
247                        double changeCost);
248     /**  Refactorizes if necessary
249          Checks if finished.  Updates status.
250          lastCleaned refers to iteration at which some objective/feasibility
251          cleaning too place.
252
253          type - 0 initial so set up save arrays etc
254               - 1 normal -if good update save
255           - 2 restoring from saved
256     */
257  void statusOfProblemInPrimal(int type);
258     /// Perturbs problem (method depends on perturbation())
259     void perturb(int type);
260     /// Take off effect of perturbation and say whether to try dual
261     bool unPerturb();
262     /// Unflag all variables and return number unflagged
263     int unflag();
264     /** Get next superbasic -1 if none,
265         Normal type is 1
266         If type is 3 then initializes sorted list
267         if 2 uses list.
268     */
269     int nextSuperBasic(int superBasicType, CoinIndexedVector * columnArray);
270
271     /// Create primal ray
272     void primalRay(CoinIndexedVector * rowArray);
273     /// Clears all bits and clears rowArray[1] etc
274     void clearAll();
275
276     /// Sort of lexicographic resolve
277     int lexSolve();
278
279     //@}
280};
281#endif
282
Note: See TracBrowser for help on using the repository browser.