source: trunk/include/ClpSimplex.hpp @ 383

Last change on this file since 383 was 383, checked in by forrest, 16 years ago

add barrier

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 40.4 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5   Authors
6 
7   John Forrest
8
9 */
10#ifndef ClpSimplex_H
11#define ClpSimplex_H
12
13#include <iostream>
14#include <cfloat>
15#include "ClpModel.hpp"
16#include "ClpMatrixBase.hpp"
17#include "ClpSolve.hpp"
18class ClpDualRowPivot;
19class ClpPrimalColumnPivot;
20class ClpFactorization;
21class CoinIndexedVector;
22class ClpNonLinearCost;
23class ClpSimplexProgress;
24
25/** This solves LPs using the simplex method
26
27    It inherits from ClpModel and all its arrays are created at
28    algorithm time. Originally I tried to work with model arrays
29    but for simplicity of coding I changed to single arrays with
30    structural variables then row variables.  Some coding is still
31    based on old style and needs cleaning up.
32
33    For a description of algorithms:
34
35    for dual see ClpSimplexDual.hpp and at top of ClpSimplexDual.cpp
36    for primal see ClpSimplexPrimal.hpp and at top of ClpSimplexPrimal.cpp
37
38    There is an algorithm data member.  + for primal variations
39    and - for dual variations
40
41    This file also includes (at end) a very simple class ClpSimplexProgress
42    which is where anti-looping stuff should migrate to
43
44*/
45
46class ClpSimplex : public ClpModel {
47   friend void ClpSimplexUnitTest(const std::string & mpsDir,
48                                  const std::string & netlibDir);
49
50public:
51  /** enums for status of various sorts.
52      First 4 match CoinWarmStartBasis,
53      isFixed means fixed at lower bound and out of basis
54  */
55  enum Status {
56    isFree = 0x00,
57    basic = 0x01,
58    atUpperBound = 0x02,
59    atLowerBound = 0x03,
60    superBasic = 0x04,
61    isFixed = 0x05
62  };
63  // For Dual
64  enum FakeBound {
65    noFake = 0x00,
66    bothFake = 0x01,
67    upperFake = 0x02,
68    lowerFake = 0x03
69  };
70
71  /**@name Constructors and destructor and copy */
72  //@{
73  /// Default constructor
74    ClpSimplex (  );
75
76  /// Copy constructor.
77  ClpSimplex(const ClpSimplex &);
78  /// Copy constructor from model.
79  ClpSimplex(const ClpModel &);
80  /** Subproblem constructor.  A subset of whole model is created from the
81      row and column lists given.  The new order is given by list order and
82      duplicates are allowed.  Name and integer information can be dropped
83  */
84  ClpSimplex (const ClpModel * wholeModel,
85              int numberRows, const int * whichRows,
86              int numberColumns, const int * whichColumns,
87              bool dropNames=true, bool dropIntegers=true);
88  /** This constructor modifies original ClpSimplex and stores
89      original stuff in created ClpSimplex.  It is only to be used in
90      conjunction with originalModel */
91  ClpSimplex (ClpSimplex * wholeModel,
92              int numberColumns, const int * whichColumns);
93  /** This copies back stuff from miniModel and then deletes miniModel.
94      Only to be used with mini constructor */
95  void originalModel(ClpSimplex * miniModel);
96  /// Assignment operator. This copies the data
97    ClpSimplex & operator=(const ClpSimplex & rhs);
98  /// Destructor
99   ~ClpSimplex (  );
100  // Ones below are just ClpModel with some changes
101  /** Loads a problem (the constraints on the
102        rows are given by lower and upper bounds). If a pointer is 0 then the
103        following values are the default:
104        <ul>
105          <li> <code>colub</code>: all columns have upper bound infinity
106          <li> <code>collb</code>: all columns have lower bound 0
107          <li> <code>rowub</code>: all rows have upper bound infinity
108          <li> <code>rowlb</code>: all rows have lower bound -infinity
109          <li> <code>obj</code>: all variables have 0 objective coefficient
110        </ul>
111    */
112  void loadProblem (  const ClpMatrixBase& matrix,
113                     const double* collb, const double* colub,   
114                     const double* obj,
115                     const double* rowlb, const double* rowub,
116                      const double * rowObjective=NULL);
117  void loadProblem (  const CoinPackedMatrix& matrix,
118                     const double* collb, const double* colub,   
119                     const double* obj,
120                     const double* rowlb, const double* rowub,
121                      const double * rowObjective=NULL);
122
123  /** Just like the other loadProblem() method except that the matrix is
124        given in a standard column major ordered format (without gaps). */
125  void loadProblem (  const int numcols, const int numrows,
126                     const CoinBigIndex* start, const int* index,
127                     const double* value,
128                     const double* collb, const double* colub,   
129                     const double* obj,
130                      const double* rowlb, const double* rowub,
131                      const double * rowObjective=NULL);
132  /// This one is for after presolve to save memory
133  void loadProblem (  const int numcols, const int numrows,
134                     const CoinBigIndex* start, const int* index,
135                      const double* value,const int * length,
136                     const double* collb, const double* colub,   
137                     const double* obj,
138                      const double* rowlb, const double* rowub,
139                      const double * rowObjective=NULL);
140  /// Read an mps file from the given filename
141  int readMps(const char *filename,
142              bool keepNames=false,
143              bool ignoreErrors = false);
144  /** Borrow model.  This is so we dont have to copy large amounts
145      of data around.  It assumes a derived class wants to overwrite
146      an empty model with a real one - while it does an algorithm.
147      This is same as ClpModel one, but sets scaling on etc. */
148  void borrowModel(ClpModel & otherModel);
149  void borrowModel(ClpSimplex & otherModel);
150   /// Pass in Event handler (cloned and deleted at end)
151   void passInEventHandler(const ClpEventHandler * eventHandler);
152  //@}
153
154  /**@name Functions most useful to user */
155  //@{
156  /** General solve algorithm which can do presolve.
157      See  ClpSolve.hpp for options
158   */
159  int initialSolve(ClpSolve & options);
160  /// Default initial solve
161  int initialSolve();
162  /// Dual initial solve
163  int initialDualSolve();
164  /// Primal initial solve
165  int initialPrimalSolve();
166  /** Dual algorithm - see ClpSimplexDual.hpp for method.
167      ifValuesPass==2 just does values pass and then stops */
168  int dual(int ifValuesPass=0);
169  /** Primal algorithm - see ClpSimplexPrimal.hpp for method.
170      ifValuesPass==2 just does values pass and then stops */
171  int primal(int ifValuesPass=0);
172  /** Solves quadratic problem using SLP - may be used as crash
173      for other algorithms when number of iterations small.
174      Also exits if all problematical variables are changing
175      less than deltaTolerance
176  */
177  int quadraticSLP(int numberPasses,double deltaTolerance);
178  /// Solves quadratic using Dantzig's algorithm - primal
179  int quadraticPrimal(int phase=2);
180  /** Solves using barrier (assumes you have good cholesky factor code).
181      Does crossover to simplex if asked*/
182  int barrier(bool crossover=true);
183  /** Dual ranging.
184      This computes increase/decrease in cost for each given variable and corresponding
185      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
186      and numberColumns.. for artificials/slacks.
187      For non-basic variables the sequence number will be that of the non-basic variables.
188      The increase/decrease value is always >= 0.0
189
190      Up to user to provide correct length arrays.
191
192      Returns non-zero if infeasible unbounded etc
193  */
194  int dualRanging(int numberCheck,const int * which,
195                  double * costIncrease, int * sequenceIncrease,
196                  double * costDecrease, int * sequenceDecrease);
197  /// Passes in factorization
198  void setFactorization( ClpFactorization & factorization);
199  /** Tightens primal bounds to make dual faster.  Unless
200      fixed, bounds are slightly looser than they could be.
201      This is to make dual go faster and is probably not needed
202      with a presolve.  Returns non-zero if problem infeasible.
203
204      Fudge for branch and bound - put bounds on columns of factor *
205      largest value (at continuous) - should improve stability
206      in branch and bound on infeasible branches (0.0 is off)
207  */
208  int tightenPrimalBounds(double factor=0.0);
209  /** Crash - at present just aimed at dual, returns
210      -2 if dual preferred and crash basis created
211      -1 if dual preferred and all slack basis preferred
212       0 if basis going in was not all slack
213       1 if primal preferred and all slack basis preferred
214       2 if primal preferred and crash basis created.
215       
216       if gap between bounds <="gap" variables can be flipped
217
218       If "pivot" is
219       0 No pivoting (so will just be choice of algorithm)
220       1 Simple pivoting e.g. gub
221       2 Mini iterations
222  */
223  int crash(double gap,int pivot);
224  /// Sets row pivot choice algorithm in dual
225  void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
226  /// Sets column pivot choice algorithm in primal
227  void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
228  /** For strong branching.  On input lower and upper are new bounds
229      while on output they are change in objective function values
230      (>1.0e50 infeasible).
231      Return code is 0 if nothing interesting, -1 if infeasible both
232      ways and +1 if infeasible one way (check values to see which one(s))
233      Solutions are filled in as well - even down, odd up - also
234      status and number of iterations
235  */
236  int strongBranching(int numberVariables,const int * variables,
237                      double * newLower, double * newUpper,
238                      double ** outputSolution,
239                      int * outputStatus, int * outputIterations,
240                      bool stopOnFirstInfeasible=true,
241                      bool alwaysFinish=false);
242  //@}
243
244  /**@name Needed for functionality of OsiSimplexInterface */
245  //@{
246  /** Pivot in a variable and out a variable.  Returns 0 if okay,
247      1 if inaccuracy forced re-factorization, -1 if would be singular.
248      Also updates primal/dual infeasibilities.
249      Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
250  */
251  int pivot();
252
253  /** Pivot in a variable and choose an outgoing one.  Assumes primal
254      feasible - will not go through a bound.  Returns step length in theta
255      Returns ray in ray_ (or NULL if no pivot)
256      Return codes as before but -1 means no acceptable pivot
257  */
258  int primalPivotResult();
259 
260  /** Pivot out a variable and choose an incoing one.  Assumes dual
261      feasible - will not go through a reduced cost. 
262      Returns step length in theta
263      Returns ray in ray_ (or NULL if no pivot)
264      Return codes as before but -1 means no acceptable pivot
265  */
266  int dualPivotResult();
267
268  /** Common bits of coding for dual and primal.  Return s0 if okay,
269      1 if bad matrix, 2 if very bad factorization
270  */
271  int startup(int ifValuesPass);
272  void finish();
273 
274  /** Factorizes and returns true if optimal.  Used by user */
275  bool statusOfProblem();
276  //@}
277
278  /**@name most useful gets and sets */
279  //@{
280  /// If problem is primal feasible
281  inline bool primalFeasible() const
282         { return (numberPrimalInfeasibilities_==0);};
283  /// If problem is dual feasible
284  inline bool dualFeasible() const
285         { return (numberDualInfeasibilities_==0);};
286  /// factorization
287  inline ClpFactorization * factorization() const 
288          { return factorization_;};
289  /// Sparsity on or off
290  bool sparseFactorization() const;
291  void setSparseFactorization(bool value);
292  /// Factorization frequency
293  int factorizationFrequency() const;
294  void setFactorizationFrequency(int value);
295  /// Dual bound
296  inline double dualBound() const
297          { return dualBound_;};
298  void setDualBound(double value);
299  /// Infeasibility cost
300  inline double infeasibilityCost() const
301          { return infeasibilityCost_;};
302  void setInfeasibilityCost(double value);
303  /** Amount of print out:
304      0 - none
305      1 - just final
306      2 - just factorizations
307      3 - as 2 plus a bit more
308      4 - verbose
309      above that 8,16,32 etc just for selective debug
310  */
311  /** Perturbation:
312      50  - switch on perturbation
313      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
314      101 - we are perturbed
315      102 - don't try perturbing again
316      default is 100
317      others are for playing
318  */
319  inline int perturbation() const
320    { return perturbation_;};
321  void setPerturbation(int value);
322  /// Current (or last) algorithm
323  inline int algorithm() const 
324  {return algorithm_; } ;
325  /// Set algorithm
326  inline void setAlgorithm(int value)
327  {algorithm_=value; } ;
328  /// Sum of dual infeasibilities
329  inline double sumDualInfeasibilities() const 
330          { return sumDualInfeasibilities_;} ;
331  inline void setSumDualInfeasibilities(double value)
332          { sumDualInfeasibilities_=value;} ;
333  /// Sum of relaxed dual infeasibilities
334  inline double sumOfRelaxedDualInfeasibilities() const 
335          { return sumOfRelaxedDualInfeasibilities_;} ;
336  inline void setSumOfRelaxedDualInfeasibilities(double value)
337          { sumOfRelaxedDualInfeasibilities_=value;} ;
338  /// Number of dual infeasibilities
339  inline int numberDualInfeasibilities() const 
340          { return numberDualInfeasibilities_;} ;
341  inline void setNumberDualInfeasibilities(int value)
342          { numberDualInfeasibilities_=value;} ;
343  /// Sum of primal infeasibilities
344  inline double sumPrimalInfeasibilities() const 
345          { return sumPrimalInfeasibilities_;} ;
346  inline void setSumPrimalInfeasibilities(double value)
347          { sumPrimalInfeasibilities_=value;} ;
348  /// Sum of relaxed primal infeasibilities
349  inline double sumOfRelaxedPrimalInfeasibilities() const 
350          { return sumOfRelaxedPrimalInfeasibilities_;} ;
351  inline void setSumOfRelaxedPrimalInfeasibilities(double value)
352          { sumOfRelaxedPrimalInfeasibilities_=value;} ;
353  /// Number of primal infeasibilities
354  inline int numberPrimalInfeasibilities() const 
355          { return numberPrimalInfeasibilities_;} ;
356  inline void setNumberPrimalInfeasibilities(int value)
357          { numberPrimalInfeasibilities_=value;} ;
358  /** Save model to file, returns 0 if success.  This is designed for
359      use outside algorithms so does not save iterating arrays etc.
360  It does not save any messaging information.
361  Does not save scaling values.
362  It does not know about all types of virtual functions.
363  */
364  int saveModel(const char * fileName);
365  /** Restore model from file, returns 0 if success,
366      deletes current model */
367  int restoreModel(const char * fileName);
368 
369  /** Just check solution (for external use) - sets sum of
370      infeasibilities etc */
371  void checkSolution();
372  /// Useful row length arrays (0,1,2,3,4,5)
373  inline CoinIndexedVector * rowArray(int index) const
374  { return rowArray_[index];};
375  /// Useful column length arrays (0,1,2,3,4,5)
376  inline CoinIndexedVector * columnArray(int index) const
377  { return columnArray_[index];};
378  //@}
379
380  /******************** End of most useful part **************/
381  /**@name Functions less likely to be useful to casual user */
382  //@{
383  /** Given an existing factorization computes and checks
384      primal and dual solutions.  Uses input arrays for variables at
385      bounds.  Returns feasibility states */
386  int getSolution (  const double * rowActivities,
387                     const double * columnActivities);
388  /** Given an existing factorization computes and checks
389      primal and dual solutions.  Uses current problem arrays for
390      bounds.  Returns feasibility states */
391  int getSolution ();
392  /** Constructs a non linear cost from list of non-linearities (columns only)
393      First lower of each column is taken as real lower
394      Last lower is taken as real upper and cost ignored
395
396      Returns nonzero if bad data e.g. lowers not monotonic
397  */
398  int createPiecewiseLinearCosts(const int * starts,
399                   const double * lower, const double * gradient);
400  /** Return model - updates any scalars */
401  void returnModel(ClpSimplex & otherModel);
402  /** Factorizes using current basis. 
403      solveType - 1 iterating, 0 initial, -1 external
404      If 10 added then in primal values pass
405      Return codes are as from ClpFactorization unless initial factorization
406      when total number of singularities is returned
407  */
408  int internalFactorize(int solveType);
409  /// Save data
410  ClpDataSave saveData() ;
411  /// Restore data
412    void restoreData(ClpDataSave saved);
413  /// Clean up status
414  void cleanStatus();
415  /// Factorizes using current basis. For external use
416  int factorize();
417  /** Computes duals from scratch. If givenDjs then
418      allows for nonzero basic djs */
419  void computeDuals(double * givenDjs);
420  /// Computes primals from scratch
421  void computePrimals (  const double * rowActivities,
422                     const double * columnActivities);
423  /** Adds multiple of a column into an array */
424  void add(double * array,
425                   int column, double multiplier) const;
426  /**
427     Unpacks one column of the matrix into indexed array
428     Uses sequenceIn_
429     Also applies scaling if needed
430  */
431  void unpack(CoinIndexedVector * rowArray) const ;
432  /**
433     Unpacks one column of the matrix into indexed array
434     Slack if sequence>= numberColumns
435     Also applies scaling if needed
436  */
437  void unpack(CoinIndexedVector * rowArray,int sequence) const;
438  /**
439     Unpacks one column of the matrix into indexed array
440     ** as packed vector
441     Uses sequenceIn_
442     Also applies scaling if needed
443  */
444  void unpackPacked(CoinIndexedVector * rowArray) ;
445  /**
446     Unpacks one column of the matrix into indexed array
447     ** as packed vector
448     Slack if sequence>= numberColumns
449     Also applies scaling if needed
450  */
451  void unpackPacked(CoinIndexedVector * rowArray,int sequence);
452 
453  /**
454      This does basis housekeeping and does values for in/out variables.
455      Can also decide to re-factorize
456  */
457  int housekeeping(double objectiveChange);
458  /** This sets largest infeasibility and most infeasible and sum
459      and number of infeasibilities (Primal) */
460  void checkPrimalSolution(const double * rowActivities=NULL,
461                           const double * columnActivies=NULL);
462  /** This sets largest infeasibility and most infeasible and sum
463      and number of infeasibilities (Dual) */
464  void checkDualSolution();
465  /** For advanced use.  When doing iterative solves things can get
466      nasty so on values pass if incoming solution has largest
467      infeasibility < incomingInfeasibility throw out variables
468      from basis until largest infeasibility < allowedInfeasibility
469      or incoming largest infeasibility.
470      If allowedInfeasibility>= incomingInfeasibility this is
471      always possible altough you may end up with an all slack basis.
472
473      Defaults are 1.0,10.0
474  */
475  void setValuesPassAction(float incomingInfeasibility,
476                           float allowedInfeasibility);
477  //@}
478  /**@name most useful gets and sets */
479  //@{
480  /// Worst column primal infeasibility
481  inline double columnPrimalInfeasibility() const 
482          { return columnPrimalInfeasibility_;} ;
483  /// Sequence of worst (-1 if feasible)
484  inline int columnPrimalSequence() const 
485          { return columnPrimalSequence_;} ;
486  /// Worst row primal infeasibility
487  inline double rowPrimalInfeasibility() const 
488          { return rowPrimalInfeasibility_;} ;
489  /// Sequence of worst (-1 if feasible)
490  inline int rowPrimalSequence() const 
491          { return rowPrimalSequence_;} ;
492  /** Worst column dual infeasibility (note - these may not be as meaningful
493      if the problem is primal infeasible */
494  inline double columnDualInfeasibility() const 
495          { return columnDualInfeasibility_;} ;
496  /// Sequence of worst (-1 if feasible)
497  inline int columnDualSequence() const 
498          { return columnDualSequence_;} ;
499  /// Worst row dual infeasibility
500  inline double rowDualInfeasibility() const 
501          { return rowDualInfeasibility_;} ;
502  /// Sequence of worst (-1 if feasible)
503  inline int rowDualSequence() const 
504          { return rowDualSequence_;} ;
505  /// Primal tolerance needed to make dual feasible (<largeTolerance)
506  inline double primalToleranceToGetOptimal() const 
507          { return primalToleranceToGetOptimal_;} ;
508  /// Remaining largest dual infeasibility
509  inline double remainingDualInfeasibility() const 
510          { return remainingDualInfeasibility_;} ;
511  /// Large bound value (for complementarity etc)
512  inline double largeValue() const 
513          { return largeValue_;} ;
514  void setLargeValue( double value) ;
515  /// Largest error on Ax-b
516  inline double largestPrimalError() const
517          { return largestPrimalError_;} ;
518  /// Largest error on basic duals
519  inline double largestDualError() const
520          { return largestDualError_;} ;
521  /// Largest difference between input primal solution and computed
522  inline double largestSolutionError() const
523          { return largestSolutionError_;} ;
524  /// Basic variables pivoting on which rows
525  inline int * pivotVariable() const
526          { return pivotVariable_;};
527  /// Scaling of objective
528  inline double objectiveScale() const 
529          { return objectiveScale_;} ;
530  inline void setObjectiveScale(double value)
531          { objectiveScale_ = value;} ;
532  /// Scaling of rhs and bounds
533  inline double rhsScale() const 
534          { return rhsScale_;} ;
535  inline void setRhsScale(double value)
536          { rhsScale_ = value;} ;
537  /// If automatic scaling on
538  inline bool automaticScaling() const
539  { return automaticScale_!=0;};
540  inline void setAutomaticScaling(bool onOff)
541  { automaticScale_ = onOff ? 1: 0;}; 
542  /// Current dual tolerance
543  inline double currentDualTolerance() const 
544          { return dualTolerance_;} ;
545  inline void setCurrentDualTolerance(double value)
546          { dualTolerance_ = value;} ;
547  /// Current primal tolerance
548  inline double currentPrimalTolerance() const 
549          { return primalTolerance_;} ;
550  inline void setCurrentPrimalTolerance(double value)
551          { primalTolerance_ = value;} ;
552  /// How many iterative refinements to do
553  inline int numberRefinements() const 
554          { return numberRefinements_;} ;
555  void setNumberRefinements( int value) ;
556  /// Alpha (pivot element) for use by classes e.g. steepestedge
557  inline double alpha() const { return alpha_;};
558  inline void setAlpha(double value) { alpha_ = value;};
559  /// Reduced cost of last incoming for use by classes e.g. steepestedge
560  inline double dualIn() const { return dualIn_;};
561  /// Pivot Row for use by classes e.g. steepestedge
562  inline int pivotRow() const{ return pivotRow_;};
563  inline void setPivotRow(int value) { pivotRow_=value;};
564  /// value of incoming variable (in Dual)
565  double valueIncomingDual() const;
566  //@}
567
568  protected:
569  /**@name protected methods */
570  //@{
571  /** May change basis and then returns number changed.
572      Computation of solutions may be overriden by given pi and solution
573  */
574  int gutsOfSolution ( double * givenDuals,
575                       const double * givenPrimals,
576                       bool valuesPass=false);
577  /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
578  void gutsOfDelete(int type);
579  /// Does most of copying
580  void gutsOfCopy(const ClpSimplex & rhs);
581  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
582      1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
583      8 bit does solution scaling in
584      16 bit does rowArray and columnArray indexed vectors
585      and makes row copy if wanted, also sets columnStart_ etc
586      Also creates scaling arrays if needed.  It does scaling if needed.
587      16 also moves solutions etc in to work arrays
588      On 16 returns false if problem "bad" i.e. matrix or bounds bad
589  */
590  bool createRim(int what,bool makeRowCopy=false);
591  /** releases above arrays and does solution scaling out.  May also
592      get rid of factorization data -
593      0 get rid of nothing, 1 get rid of arrays, 2 also factorization
594  */
595  void deleteRim(int getRidOfFactorizationData=2);
596  /// Sanity check on input rim data (after scaling) - returns true if okay
597  bool sanityCheck();
598  //@}
599  public:
600  /**@name public methods */
601  //@{
602  /** Return row or column sections - not as much needed as it
603      once was.  These just map into single arrays */
604  inline double * solutionRegion(int section)
605  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
606  inline double * djRegion(int section)
607  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
608  inline double * lowerRegion(int section)
609  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
610  inline double * upperRegion(int section)
611  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
612  inline double * costRegion(int section)
613  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
614  /// Return region as single array
615  inline double * solutionRegion()
616  { return solution_;};
617  inline double * djRegion()
618  { return dj_;};
619  inline double * lowerRegion()
620  { return lower_;};
621  inline double * upperRegion()
622  { return upper_;};
623  inline double * costRegion()
624  { return cost_;};
625  inline Status getStatus(int sequence) const
626  {return static_cast<Status> (status_[sequence]&7);};
627  inline void setStatus(int sequence, Status status)
628  {
629    unsigned char & st_byte = status_[sequence];
630    st_byte &= ~7;
631    st_byte |= status;
632  };
633  /** Normally the first factorization does sparse coding because
634      the factorization could be singular.  This allows initial dense
635      factorization when it is known to be safe
636  */
637  void setInitialDenseFactorization(bool onOff);
638  bool  initialDenseFactorization() const;
639  /** Return sequence In or Out */
640  inline int sequenceIn() const
641  {return sequenceIn_;};
642  inline int sequenceOut() const
643  {return sequenceOut_;};
644  /** Set sequenceIn or Out */
645  inline void  setSequenceIn(int sequence)
646  { sequenceIn_=sequence;};
647  inline void  setSequenceOut(int sequence)
648  { sequenceOut_=sequence;};
649  /** Return direction In or Out */
650  inline int directionIn() const
651  {return directionIn_;};
652  inline int directionOut() const
653  {return directionOut_;};
654  /** Set directionIn or Out */
655  inline void  setDirectionIn(int direction)
656  { directionIn_=direction;};
657  inline void  setDirectionOut(int direction)
658  { directionOut_=direction;};
659  /// Value of Out variable
660  inline double valueOut() const
661  { return valueOut_;};
662  /// Returns 1 if sequence indicates column
663  inline int isColumn(int sequence) const
664  { return sequence<numberColumns_ ? 1 : 0;};
665  /// Returns sequence number within section
666  inline int sequenceWithin(int sequence) const
667  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
668  /// Return row or column values
669  inline double solution(int sequence)
670  { return solution_[sequence];};
671  /// Return address of row or column values
672  inline double & solutionAddress(int sequence)
673  { return solution_[sequence];};
674  inline double reducedCost(int sequence)
675   { return dj_[sequence];};
676  inline double & reducedCostAddress(int sequence)
677   { return dj_[sequence];};
678  inline double lower(int sequence)
679  { return lower_[sequence];};
680  /// Return address of row or column lower bound
681  inline double & lowerAddress(int sequence)
682  { return lower_[sequence];};
683  inline double upper(int sequence)
684  { return upper_[sequence];};
685  /// Return address of row or column upper bound
686  inline double & upperAddress(int sequence)
687  { return upper_[sequence];};
688  inline double cost(int sequence)
689  { return cost_[sequence];};
690  /// Return address of row or column cost
691  inline double & costAddress(int sequence)
692  { return cost_[sequence];};
693  /// Return original lower bound
694  inline double originalLower(int iSequence) const
695  { if (iSequence<numberColumns_) return columnLower_[iSequence]; else
696    return rowLower_[iSequence-numberColumns_];};
697  /// Return original lower bound
698  inline double originalUpper(int iSequence) const
699  { if (iSequence<numberColumns_) return columnUpper_[iSequence]; else
700    return rowUpper_[iSequence-numberColumns_];};
701  /// Theta (pivot change)
702  inline double theta() const
703  { return theta_;};
704  /// Return pointer to details of costs
705  inline ClpNonLinearCost * nonLinearCost() const
706  { return nonLinearCost_;};
707  //@}
708  /**@name status methods */
709  //@{
710  inline void setFakeBound(int sequence, FakeBound fakeBound)
711  {
712    unsigned char & st_byte = status_[sequence];
713    st_byte &= ~24;
714    st_byte |= fakeBound<<3;
715  };
716  inline FakeBound getFakeBound(int sequence) const
717  {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
718  inline void setRowStatus(int sequence, Status status)
719  {
720    unsigned char & st_byte = status_[sequence+numberColumns_];
721    st_byte &= ~7;
722    st_byte |= status;
723  };
724  inline Status getRowStatus(int sequence) const
725  {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
726  inline void setColumnStatus(int sequence, Status status)
727  {
728    unsigned char & st_byte = status_[sequence];
729    st_byte &= ~7;
730    st_byte |= status;
731  };
732  inline Status getColumnStatus(int sequence) const
733  {return static_cast<Status> (status_[sequence]&7);};
734  inline void setPivoted( int sequence)
735  { status_[sequence] |= 32;};
736  inline void clearPivoted( int sequence)
737  { status_[sequence] &= ~32; };
738  inline bool pivoted(int sequence) const
739  {return (((status_[sequence]>>5)&1)!=0);};
740  /// To flag a variable (not inline to allow for column generation)
741  void setFlagged( int sequence);
742  inline void clearFlagged( int sequence)
743  {
744    status_[sequence] &= ~64;
745  };
746  inline bool flagged(int sequence) const
747  {return ((status_[sequence]&64)!=0);};
748  /// To say row active in primal pivot row choice
749  inline void setActive( int iRow)
750  {
751    status_[iRow] |= 128;
752  };
753  inline void clearActive( int iRow)
754  {
755    status_[iRow] &= ~128;
756  };
757  inline bool active(int iRow) const
758  {return ((status_[iRow]&128)!=0);};
759  /** Set up status array (can be used by OsiClp).
760      Also can be used to set up all slack basis */
761  void createStatus() ;
762  inline void allSlackBasis()
763  { createStatus();};
764   
765  /// So we know when to be cautious
766  inline int lastBadIteration() const
767  {return lastBadIteration_;};
768  /// Progress flag - at present 0 bit says artificials out
769  inline int progressFlag() const
770  {return progressFlag_;};
771  /// Force re-factorization early
772  inline void forceFactorization(int value)
773  { forceFactorization_ = value;};
774  /// Raw objective value (so always minimize in primal)
775  inline double rawObjectiveValue() const
776  { return objectiveValue_;};
777  /** Number of extra rows.  These are ones which will be dynamically created
778      each iteration.  This is for GUB but may have other uses.
779  */
780  inline int numberExtraRows() const
781  { return numberExtraRows_;};
782  /** Maximum number of basic variables - can be more than number of rows if GUB
783  */
784  inline int maximumBasic() const
785  { return maximumBasic_;};
786  /** For advanced options
787      1 - Don't keep changing infeasibility weight
788      2 - Keep nonLinearCost round solves
789      4 - Force outgoing variables to exact bound (primal)
790      8 - Safe to use dense initial factorization
791      16 -Just use basic variables for operation if column generation
792      32 -Clean up with primal before strong branching
793  */
794  inline unsigned int specialOptions() const
795  { return specialOptions_;};
796  inline void setSpecialOptions(unsigned int value)
797  { specialOptions_=value;};
798  //@}
799
800////////////////// data //////////////////
801protected:
802
803  /**@name data.  Many arrays have a row part and a column part.
804   There is a single array with both - columns then rows and
805   then normally two arrays pointing to rows and columns.  The
806   single array is the owner of memory
807  */
808  //@{
809  /// Worst column primal infeasibility
810  double columnPrimalInfeasibility_;
811  /// Worst row primal infeasibility
812  double rowPrimalInfeasibility_;
813  /// Sequence of worst (-1 if feasible)
814  int columnPrimalSequence_;
815  /// Sequence of worst (-1 if feasible)
816  int rowPrimalSequence_;
817  /// Worst column dual infeasibility
818  double columnDualInfeasibility_;
819  /// Worst row dual infeasibility
820  double rowDualInfeasibility_;
821  /// Sequence of worst (-1 if feasible)
822  int columnDualSequence_;
823  /// Sequence of worst (-1 if feasible)
824  int rowDualSequence_;
825  /// Primal tolerance needed to make dual feasible (<largeTolerance)
826  double primalToleranceToGetOptimal_;
827  /// Remaining largest dual infeasibility
828  double remainingDualInfeasibility_;
829  /// Large bound value (for complementarity etc)
830  double largeValue_;
831  /// Scaling of objective
832  double objectiveScale_;
833  /// Scaling of rhs and bounds
834  double rhsScale_;
835  /// Largest error on Ax-b
836  double largestPrimalError_;
837  /// Largest error on basic duals
838  double largestDualError_;
839  /// Largest difference between input primal solution and computed
840  double largestSolutionError_;
841  /// Dual bound
842  double dualBound_;
843  /// Alpha (pivot element)
844  double alpha_;
845  /// Theta (pivot change)
846  double theta_;
847  /// Lower Bound on In variable
848  double lowerIn_;
849  /// Value of In variable
850  double valueIn_;
851  /// Upper Bound on In variable
852  double upperIn_;
853  /// Reduced cost of In variable
854  double dualIn_;
855  /// Lower Bound on Out variable
856  double lowerOut_;
857  /// Value of Out variable
858  double valueOut_;
859  /// Upper Bound on Out variable
860  double upperOut_;
861  /// Infeasibility (dual) or ? (primal) of Out variable
862  double dualOut_;
863  /// Current dual tolerance for algorithm
864  double dualTolerance_;
865  /// Current primal tolerance for algorithm
866  double primalTolerance_;
867  /// Sum of dual infeasibilities
868  double sumDualInfeasibilities_;
869  /// Sum of primal infeasibilities
870  double sumPrimalInfeasibilities_;
871  /// Weight assigned to being infeasible in primal
872  double infeasibilityCost_;
873  /// Sum of Dual infeasibilities using tolerance based on error in duals
874  double sumOfRelaxedDualInfeasibilities_;
875  /// Sum of Primal infeasibilities using tolerance based on error in primals
876  double sumOfRelaxedPrimalInfeasibilities_;
877  /// Working copy of lower bounds (Owner of arrays below)
878  double * lower_;
879  /// Row lower bounds - working copy
880  double * rowLowerWork_;
881  /// Column lower bounds - working copy
882  double * columnLowerWork_;
883  /// Working copy of upper bounds (Owner of arrays below)
884  double * upper_;
885  /// Row upper bounds - working copy
886  double * rowUpperWork_;
887  /// Column upper bounds - working copy
888  double * columnUpperWork_;
889  /// Working copy of objective (Owner of arrays below)
890  double * cost_;
891  /// Row objective - working copy
892  double * rowObjectiveWork_;
893  /// Column objective - working copy
894  double * objectiveWork_;
895  /// Useful row length arrays
896  CoinIndexedVector * rowArray_[6];
897  /// Useful column length arrays
898  CoinIndexedVector * columnArray_[6];
899  /// Sequence of In variable
900  int sequenceIn_;
901  /// Direction of In, 1 going up, -1 going down, 0 not a clude
902  int directionIn_;
903  /// Sequence of Out variable
904  int sequenceOut_;
905  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
906  int directionOut_;
907  /// Pivot Row
908  int pivotRow_;
909  /// Last good iteration (immediately after a re-factorization)
910  int lastGoodIteration_;
911  /// Working copy of reduced costs (Owner of arrays below)
912  double * dj_;
913  /// Reduced costs of slacks not same as duals (or - duals)
914  double * rowReducedCost_;
915  /// Possible scaled reduced costs
916  double * reducedCostWork_;
917  /// Working copy of primal solution (Owner of arrays below)
918  double * solution_;
919  /// Row activities - working copy
920  double * rowActivityWork_;
921  /// Column activities - working copy
922  double * columnActivityWork_;
923  /// Number of dual infeasibilities
924  int numberDualInfeasibilities_;
925  /// Number of dual infeasibilities (without free)
926  int numberDualInfeasibilitiesWithoutFree_;
927  /// Number of primal infeasibilities
928  int numberPrimalInfeasibilities_;
929  /// How many iterative refinements to do
930  int numberRefinements_;
931  /// dual row pivot choice
932  ClpDualRowPivot * dualRowPivot_;
933  /// primal column pivot choice
934  ClpPrimalColumnPivot * primalColumnPivot_;
935  /// Basic variables pivoting on which rows
936  int * pivotVariable_;
937  /// factorization
938  ClpFactorization * factorization_;
939  /// Saved version of solution
940  double * savedSolution_;
941  /// Number of times code has tentatively thought optimal
942  int numberTimesOptimal_;
943  /// If change has been made (first attempt at stopping looping)
944  int changeMade_;
945  /// Algorithm >0 == Primal, <0 == Dual
946  int algorithm_;
947  /** Now for some reliability aids
948      This forces re-factorization early */
949  int forceFactorization_;
950  /** Perturbation:
951      -50 to +50 - perturb by this power of ten (-6 sounds good)
952      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
953      101 - we are perturbed
954      102 - don't try perturbing again
955      default is 100
956  */
957  int perturbation_;
958  /// Saved status regions
959  unsigned char * saveStatus_;
960  /** Very wasteful way of dealing with infeasibilities in primal.
961      However it will allow non-linearities and use of dual
962      analysis.  If it doesn't work it can easily be replaced.
963  */
964  ClpNonLinearCost * nonLinearCost_;
965  /** For advanced options
966      1 - Don't keep changing infeasibility weight
967      2 - Keep nonLinearCost round solves
968      4 - Force outgoing variables to exact bound (primal)
969      8 - Safe to use dense initial factorization
970      16 -Just use basic variables for operation
971  */
972  unsigned int specialOptions_;
973  /// So we know when to be cautious
974  int lastBadIteration_;
975  /// Can be used for count of fake bounds (dual) or fake costs (primal)
976  int numberFake_;
977  /// Progress flag - at present 0 bit says artificials out, 1 free in
978  int progressFlag_;
979  /// First free/super-basic variable (-1 if none)
980  int firstFree_;
981  /** Number of extra rows.  These are ones which will be dynamically created
982      each iteration.  This is for GUB but may have other uses.
983  */
984  int numberExtraRows_;
985  /** Maximum number of basic variables - can be more than number of rows if GUB
986  */
987  int maximumBasic_;
988  /** For advanced use.  When doing iterative solves things can get
989      nasty so on values pass if incoming solution has largest
990      infeasibility < incomingInfeasibility throw out variables
991      from basis until largest infeasibility < allowedInfeasibility.
992      if allowedInfeasibility>= incomingInfeasibility this is
993      always possible altough you may end up with an all slack basis.
994
995      Defaults are 1.0,10.0
996  */
997  float incomingInfeasibility_;
998  float allowedInfeasibility_;
999  /// Automatic scaling of objective and rhs and bounds
1000  int automaticScale_;
1001  /// For dealing with all issues of cycling etc
1002  ClpSimplexProgress * progress_;
1003  //@}
1004};
1005//#############################################################################
1006/** A function that tests the methods in the ClpSimplex class. The
1007    only reason for it not to be a member method is that this way it doesn't
1008    have to be compiled into the library. And that's a gain, because the
1009    library should be compiled with optimization on, but this method should be
1010    compiled with debugging.
1011
1012    It also does some testing of ClpFactorization class
1013 */
1014void
1015ClpSimplexUnitTest(const std::string & mpsDir,
1016                   const std::string & netlibDir);
1017
1018
1019/// For saving extra information to see if looping.
1020class ClpSimplexProgress {
1021
1022public:
1023
1024
1025  /**@name Constructors and destructor and copy */
1026  //@{
1027  /// Default constructor
1028    ClpSimplexProgress (  );
1029
1030  /// Constructor from model
1031    ClpSimplexProgress ( ClpSimplex * model );
1032
1033  /// Copy constructor.
1034  ClpSimplexProgress(const ClpSimplexProgress &);
1035
1036  /// Assignment operator. This copies the data
1037    ClpSimplexProgress & operator=(const ClpSimplexProgress & rhs);
1038  /// Destructor
1039   ~ClpSimplexProgress (  );
1040  //@}
1041
1042  /**@name Check progress */
1043  //@{
1044  /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
1045      >=0 if give up and use as problem status
1046  */
1047    int looping (  );
1048  /// Start check at beginning of whileIterating
1049  void startCheck();
1050  /// Returns cycle length in whileIterating
1051  int cycle(int in, int out,int wayIn,int wayOut); 
1052
1053  /// Returns previous objective (if -1) - current if (0)
1054  double lastObjective(int back=1) const;
1055  /// Modify objective e.g. if dual infeasible in dual
1056  void modifyObjective(double value);
1057  /// Returns previous iteration number (if -1) - current if (0)
1058  int lastIterationNumber(int back=1) const;
1059  /// Odd state
1060  inline void newOddState()
1061  { oddState_= - oddState_-1;};
1062  inline void endOddState()
1063  { oddState_=abs(oddState_);};
1064  inline void clearOddState() 
1065  { oddState_=0;};
1066  inline int oddState() const
1067  { return oddState_;};
1068  /// number of bad times
1069  inline int badTimes() const
1070  { return numberBadTimes_;};
1071  inline void clearBadTimes()
1072  { numberBadTimes_=0;};
1073
1074  //@}
1075  /**@name Data  */
1076#define CLP_PROGRESS 5
1077  //@{
1078  /// Objective values
1079  double objective_[CLP_PROGRESS];
1080  /// Sum of infeasibilities for algorithm
1081  double infeasibility_[CLP_PROGRESS];
1082#define CLP_CYCLE 12
1083  /// For cycle checking
1084  //double obj_[CLP_CYCLE];
1085  int in_[CLP_CYCLE];
1086  int out_[CLP_CYCLE];
1087  char way_[CLP_CYCLE];
1088  /// Pointer back to model so we can get information
1089  ClpSimplex * model_;
1090  /// Number of infeasibilities
1091  int numberInfeasibilities_[CLP_PROGRESS];
1092  /// Iteration number at which occurred
1093  int iterationNumber_[CLP_PROGRESS];
1094  /// Number of times checked (so won't stop too early)
1095  int numberTimes_;
1096  /// Number of times it looked like loop
1097  int numberBadTimes_;
1098  /// If things are in an odd state
1099  int oddState_;
1100  //@}
1101};
1102#endif
Note: See TracBrowser for help on using the repository browser.