source: stable/1.15/Clp/src/AbcSimplexDual.hpp @ 1912

Last change on this file since 1912 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: 11.6 KB
Line 
1/* $Id: AbcSimplexDual.hpp 1910 2013-01-27 02:00:13Z 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 AbcSimplexDual_H
12#define AbcSimplexDual_H
13
14#include "AbcSimplex.hpp"
15#if 0
16#undef ABC_PARALLEL
17#define ABC_PARALLEL 2
18#undef cilk_for
19#undef cilk_spawn
20#undef cilk_sync
21#include <cilk/cilk.h>
22#endif
23typedef struct {
24  double theta;
25  double totalThru;
26  double useThru;
27  double bestEverPivot;
28  double increaseInObjective;
29  double tentativeTheta;
30  double lastPivotValue;
31  double thisPivotValue;
32  double thruThis;
33  double increaseInThis;
34  int lastSequence;
35  int sequence;
36  int block;
37  int numberSwapped;
38  int numberRemaining;
39  int numberLastSwapped;
40  bool modifyCosts;
41} dualColumnResult;
42/** This solves LPs using the dual simplex method
43   
44    It inherits from AbcSimplex.  It has no data of its own and
45    is never created - only cast from a AbcSimplex object at algorithm time.
46   
47*/
48
49class AbcSimplexDual : public AbcSimplex {
50 
51public:
52 
53  /**@name Description of algorithm */
54  //@{
55  /** Dual algorithm
56     
57      Method
58     
59      It tries to be a single phase approach with a weight of 1.0 being
60      given to getting optimal and a weight of updatedDualBound_ being
61      given to getting dual feasible.  In this version I have used the
62      idea that this weight can be thought of as a fake bound.  If the
63      distance between the lower and upper bounds on a variable is less
64      than the feasibility weight then we are always better off flipping
65      to other bound to make dual feasible.  If the distance is greater
66      then we make up a fake bound updatedDualBound_ away from one bound.
67      If we end up optimal or primal infeasible, we check to see if
68      bounds okay.  If so we have finished, if not we increase updatedDualBound_
69      and continue (after checking if unbounded). I am undecided about
70      free variables - there is coding but I am not sure about it.  At
71      present I put them in basis anyway.
72     
73      The code is designed to take advantage of sparsity so arrays are
74      seldom zeroed out from scratch or gone over in their entirety.
75      The only exception is a full scan to find outgoing variable for
76      Dantzig row choice.  For steepest edge we keep an updated list
77      of infeasibilities (actually squares).
78      On easy problems we don't need full scan - just
79      pick first reasonable.
80     
81      One problem is how to tackle degeneracy and accuracy.  At present
82      I am using the modification of costs which I put in OSL and some
83      of what I think is the dual analog of Gill et al.
84      I am still not sure of the exact details.
85     
86      The flow of dual is three while loops as follows:
87     
88      while (not finished) {
89     
90      while (not clean solution) {
91     
92      Factorize and/or clean up solution by flipping variables so
93      dual feasible.  If looks finished check fake dual bounds.
94      Repeat until status is iterating (-1) or finished (0,1,2)
95     
96      }
97     
98      while (status==-1) {
99     
100      Iterate until no pivot in or out or time to re-factorize.
101     
102      Flow is:
103     
104      choose pivot row (outgoing variable).  if none then
105      we are primal feasible so looks as if done but we need to
106      break and check bounds etc.
107     
108      Get pivot row in tableau
109     
110      Choose incoming column.  If we don't find one then we look
111      primal infeasible so break and check bounds etc.  (Also the
112      pivot tolerance is larger after any iterations so that may be
113      reason)
114     
115      If we do find incoming column, we may have to adjust costs to
116      keep going forwards (anti-degeneracy).  Check pivot will be stable
117      and if unstable throw away iteration and break to re-factorize.
118      If minor error re-factorize after iteration.
119     
120      Update everything (this may involve flipping variables to stay
121      dual feasible.
122     
123      }
124     
125      }
126     
127      TODO's (or maybe not)
128     
129      At present we never check we are going forwards.  I overdid that in
130      OSL so will try and make a last resort.
131     
132      Needs partial scan pivot out option.
133     
134      May need other anti-degeneracy measures, especially if we try and use
135      loose tolerances as a way to solve in fewer iterations.
136     
137      I like idea of dynamic scaling.  This gives opportunity to decouple
138      different implications of scaling for accuracy, iteration count and
139      feasibility tolerance.
140     
141      for use of exotic parameter startFinishoptions see Abcsimplex.hpp
142  */
143 
144  int dual();
145  /** For strong branching.  On input lower and upper are new bounds
146      while on output they are change in objective function values
147      (>1.0e50 infeasible).
148      Return code is 0 if nothing interesting, -1 if infeasible both
149      ways and +1 if infeasible one way (check values to see which one(s))
150      Solutions are filled in as well - even down, odd up - also
151      status and number of iterations
152  */
153  int strongBranching(int numberVariables, const int * variables,
154                      double * newLower, double * newUpper,
155                      double ** outputSolution,
156                      int * outputStatus, int * outputIterations,
157                      bool stopOnFirstInfeasible = true,
158                      bool alwaysFinish = false,
159                      int startFinishOptions = 0);
160  /// This does first part of StrongBranching
161  AbcSimplexFactorization * setupForStrongBranching(char * arrays, int numberRows,
162                                                    int numberColumns, bool solveLp = false);
163  /// This cleans up after strong branching
164  void cleanupAfterStrongBranching(AbcSimplexFactorization * factorization);
165  //@}
166 
167  /**@name Functions used in dual */
168  //@{
169  /** This has the flow between re-factorizations
170      Broken out for clarity and will be used by strong branching
171     
172      Reasons to come out:
173      -1 iterations etc
174      -2 inaccuracy
175      -3 slight inaccuracy (and done iterations)
176      +0 looks optimal (might be unbounded - but we will investigate)
177      +1 looks infeasible
178      +3 max iterations
179     
180      If givenPi not NULL then in values pass (copy from ClpSimplexDual)
181  */
182  int whileIteratingSerial();
183#if ABC_PARALLEL==1
184  int whileIteratingThread();
185#endif
186#if ABC_PARALLEL==2
187  int whileIteratingCilk();
188#endif
189  void whileIterating2();
190  int whileIteratingParallel(int numberIterations);
191  int whileIterating3();
192  void updatePrimalSolution();
193  int noPivotRow();
194  int noPivotColumn();
195  void dualPivotColumn();
196  /// Create dual pricing vector
197  void createDualPricingVectorSerial();
198  int getTableauColumnFlipAndStartReplaceSerial();
199  void getTableauColumnPart1Serial();
200#if ABC_PARALLEL==1
201  /// Create dual pricing vector
202  void createDualPricingVectorThread();
203  int getTableauColumnFlipAndStartReplaceThread();
204  void getTableauColumnPart1Thread();
205#endif
206#if ABC_PARALLEL==2
207  /// Create dual pricing vector
208  void createDualPricingVectorCilk();
209  int getTableauColumnFlipAndStartReplaceCilk();
210  void getTableauColumnPart1Cilk();
211#endif
212  void getTableauColumnPart2();
213  int checkReplace();
214  void replaceColumnPart3();
215  void checkReplacePart1();
216  void checkReplacePart1a();
217  void checkReplacePart1b();
218  /// The duals are updated
219  void updateDualsInDual();
220  /** The duals are updated by the given arrays.
221      This is in values pass - so no changes to primal is made
222  */
223  //void updateDualsInValuesPass(CoinIndexedVector * array,
224  //                           double theta);
225  /** While dualColumn gets flips
226      this does actual flipping.
227      returns number flipped
228  */
229  int flipBounds();
230  /** Undo a flip
231  */
232  void flipBack(int number);
233  /**
234     Array has tableau row (row section)
235     Puts candidates for rows in list
236     Returns guess at upper theta (infinite if no pivot) and may set sequenceIn_ if free
237     Can do all (if tableauRow created)
238  */
239  void dualColumn1(bool doAll=false);
240  /**
241     Array has tableau row (row section)
242     Just does slack part
243     Returns guess at upper theta (infinite if no pivot) and may set sequenceIn_ if free
244  */
245  double dualColumn1A();
246  /// Do all given tableau row
247  double dualColumn1B();
248  /**
249     Chooses incoming
250     Puts flipped ones in list
251     If necessary will modify costs
252  */
253  void dualColumn2();
254  void dualColumn2Most(dualColumnResult & result);
255  void dualColumn2First(dualColumnResult & result);
256  /**
257     Chooses part of incoming
258     Puts flipped ones in list
259     If necessary will modify costs
260  */
261  void dualColumn2(dualColumnResult & result);
262  /**
263     This sees what is best thing to do in branch and bound cleanup
264     If sequenceIn_ < 0 then can't do anything
265  */
266  void checkPossibleCleanup(CoinIndexedVector * array);
267  /**
268     Chooses dual pivot row
269     Would be faster with separate region to scan
270     and will have this (with square of infeasibility) when steepest
271     For easy problems we can just choose one of the first rows we look at
272  */
273  void dualPivotRow();
274  /** Checks if any fake bounds active - if so returns number and modifies
275      updatedDualBound_ and everything.
276      Free variables will be left as free
277      Returns number of bounds changed if >=0
278      Returns -1 if not initialize and no effect
279      fills cost of change vector
280  */
281  int changeBounds(int initialize, double & changeCost);
282  /** As changeBounds but just changes new bounds for a single variable.
283      Returns true if change */
284  bool changeBound( int iSequence);
285  /// Restores bound to original bound
286  void originalBound(int iSequence);
287  /** Checks if tentative optimal actually means unbounded in dual
288      Returns -3 if not, 2 if is unbounded */
289  int checkUnbounded(CoinIndexedVector & ray,  double changeCost);
290  /**  Refactorizes if necessary
291       Checks if finished.  Updates status.
292       lastCleaned refers to iteration at which some objective/feasibility
293       cleaning too place.
294       
295       type - 0 initial so set up save arrays etc
296       - 1 normal -if good update save
297       - 2 restoring from saved
298  */
299  void statusOfProblemInDual(int type);
300  /** Fast iterations.  Misses out a lot of initialization.
301      Normally stops on maximum iterations, first re-factorization
302      or tentative optimum.  If looks interesting then continues as
303      normal.  Returns 0 if finished properly, 1 otherwise.
304  */
305  /// Gets tableau column - does flips and checks what to do next
306  /// Knows tableau column in 1, flips in 2 and gets an array for flips (as serial here)
307  int whatNext();
308  /// see if cutoff reached
309  bool checkCutoff(bool computeObjective);
310  /// Does something about fake tolerances
311  int bounceTolerances(int type);
312  /// Perturbs problem
313  void perturb(double factor);
314  /// Perturbs problem B
315  void perturbB(double factor,int type);
316  /// Make non free variables dual feasible by moving to a bound
317  int makeNonFreeVariablesDualFeasible(bool changeCosts=false);
318  int fastDual(bool alwaysFinish = false);
319  /** Checks number of variables at fake bounds.  This is used by fastDual
320      so can exit gracefully before end */
321  int numberAtFakeBound();
322 
323  /** Pivot in a variable and choose an outgoing one.  Assumes dual
324      feasible - will not go through a reduced cost.  Returns step length in theta
325      Return codes as before but -1 means no acceptable pivot
326  */
327  int pivotResultPart1();
328  /** Get next free , -1 if none */
329  int nextSuperBasic();
330  /// Startup part of dual
331  void startupSolve();
332  /// Ending part of dual
333  void finishSolve();
334  void gutsOfDual();
335  //int dual2(int ifValuesPass,int startFinishOptions=0);
336  int resetFakeBounds(int type);
337 
338  //@}
339};
340#if ABC_PARALLEL==1
341void * abc_parallelManager(void * simplex);
342#endif
343#endif
Note: See TracBrowser for help on using the repository browser.