source: trunk/Clp/src/AbcSimplexPrimal.hpp @ 2470

Last change on this file since 2470 was 2385, checked in by unxusr, 10 months ago

formatting

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