source: branches/pre/include/ClpSimplex.hpp @ 222

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

Start of mini sprint

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.2 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
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  //@}
150
151  /**@name Functions most useful to user */
152  //@{
153  /** General solve algorithm which can do presolve.
154      See  ClpSolve.hpp for options
155   */
156  int initialSolve(ClpSolve & options);
157  /// Default initial solve
158  int initialSolve();
159  /// Dual initial solve
160  int initialDualSolve();
161  /// Primal initial solve
162  int initialPrimalSolve();
163  /** Dual algorithm - see ClpSimplexDual.hpp for method */
164  int dual(int ifValuesPass=0);
165  /** Primal algorithm - see ClpSimplexPrimal.hpp for method */
166  int primal(int ifValuesPass=0);
167  /** Solves quadratic problem using SLP - may be used as crash
168      for other algorithms when number of iterations small.
169      Also exits if all problematical variables are changing
170      less than deltaTolerance
171  */
172  int quadraticSLP(int numberPasses,double deltaTolerance);
173  /// Solves quadratic using Dantzig's algorithm - primal
174  int quadraticPrimal(int phase=2);
175  /// Passes in factorization
176  void setFactorization( ClpFactorization & factorization);
177  /// Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
178  void scaling(int mode=1);
179  /// Gets scalingFlag
180  inline int scalingFlag() const {return scalingFlag_;};
181  /** Tightens primal bounds to make dual faster.  Unless
182      fixed, bounds are slightly looser than they could be.
183      This is to make dual go faster and is probably not needed
184      with a presolve.  Returns non-zero if problem infeasible.
185
186      Fudge for branch and bound - put bounds on columns of factor *
187      largest value (at continuous) - should improve stability
188      in branch and bound on infeasible branches (0.0 is off)
189  */
190  int tightenPrimalBounds(double factor=0.0);
191  /** Crash - at present just aimed at dual, returns
192      -2 if dual preferred and crash basis created
193      -1 if dual preferred and all slack basis preferred
194       0 if basis going in was not all slack
195       1 if primal preferred and all slack basis preferred
196       2 if primal preferred and crash basis created.
197       
198       if gap between bounds <="gap" variables can be flipped
199
200       If "pivot" is
201       0 No pivoting (so will just be choice of algorithm)
202       1 Simple pivoting e.g. gub
203       2 Mini iterations
204  */
205  int crash(double gap,int pivot);
206  /// Sets row pivot choice algorithm in dual
207  void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
208  /// Sets column pivot choice algorithm in primal
209  void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
210  /** For strong branching.  On input lower and upper are new bounds
211      while on output they are change in objective function values
212      (>1.0e50 infeasible).
213      Return code is 0 if nothing interesting, -1 if infeasible both
214      ways and +1 if infeasible one way (check values to see which one(s))
215      Solutions are filled in as well - even down, odd up - also
216      status and number of iterations
217  */
218  int strongBranching(int numberVariables,const int * variables,
219                      double * newLower, double * newUpper,
220                      double ** outputSolution,
221                      int * outputStatus, int * outputIterations,
222                      bool stopOnFirstInfeasible=true,
223                      bool alwaysFinish=false);
224  //@}
225
226  /**@name Needed for functionality of OsiSimplexInterface */
227  //@{
228  /** Pivot in a variable and out a variable.  Returns 0 if okay,
229      1 if inaccuracy forced re-factorization, -1 if would be singular.
230      Also updates primal/dual infeasibilities.
231      Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
232  */
233  int pivot();
234
235  /** Pivot in a variable and choose an outgoing one.  Assumes primal
236      feasible - will not go through a bound.  Returns step length in theta
237      Returns ray in ray_ (or NULL if no pivot)
238      Return codes as before but -1 means no acceptable pivot
239  */
240  int primalPivotResult();
241 
242  /** Pivot out a variable and choose an incoing one.  Assumes dual
243      feasible - will not go through a reduced cost. 
244      Returns step length in theta
245      Returns ray in ray_ (or NULL if no pivot)
246      Return codes as before but -1 means no acceptable pivot
247  */
248  int dualPivotResult();
249
250  /** Common bits of coding for dual and primal.  Return s0 if okay,
251      1 if bad matrix, 2 if very bad factorization
252  */
253  int startup(int ifValuesPass);
254  void finish();
255 
256  /** Factorizes and returns true if optimal.  Used by user */
257  bool statusOfProblem();
258  //@}
259
260  /**@name most useful gets and sets */
261  //@{
262  /// If problem is primal feasible
263  inline bool primalFeasible() const
264         { return (numberPrimalInfeasibilities_==0);};
265  /// If problem is dual feasible
266  inline bool dualFeasible() const
267         { return (numberDualInfeasibilities_==0);};
268  /// factorization
269  inline ClpFactorization * factorization() const 
270          { return factorization_;};
271  /// Sparsity on or off
272  bool sparseFactorization() const;
273  void setSparseFactorization(bool value);
274  /// Factorization frequency
275  int factorizationFrequency() const;
276  void setFactorizationFrequency(int value);
277  /// Dual bound
278  inline double dualBound() const
279          { return dualBound_;};
280  void setDualBound(double value);
281  /// Infeasibility cost
282  inline double infeasibilityCost() const
283          { return infeasibilityCost_;};
284  void setInfeasibilityCost(double value);
285  /** Amount of print out:
286      0 - none
287      1 - just final
288      2 - just factorizations
289      3 - as 2 plus a bit more
290      4 - verbose
291      above that 8,16,32 etc just for selective debug
292  */
293  /** Perturbation:
294      -50 to +50 - perturb by this power of ten (-6 sounds good)
295      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
296      101 - we are perturbed
297      102 - don't try perturbing again
298      default is 100
299  */
300  inline int perturbation() const
301    { return perturbation_;};
302  void setPerturbation(int value);
303  /// Current (or last) algorithm
304  inline int algorithm() const 
305  {return algorithm_; } ;
306  /// Set algorithm
307  inline void setAlgorithm(int value)
308  {algorithm_=value; } ;
309  /// Sum of dual infeasibilities
310  inline double sumDualInfeasibilities() const 
311          { return sumDualInfeasibilities_;} ;
312  /// Number of dual infeasibilities
313  inline int numberDualInfeasibilities() const 
314          { return numberDualInfeasibilities_;} ;
315  /// Sum of primal infeasibilities
316  inline double sumPrimalInfeasibilities() const 
317          { return sumPrimalInfeasibilities_;} ;
318  /// Number of primal infeasibilities
319  inline int numberPrimalInfeasibilities() const 
320          { return numberPrimalInfeasibilities_;} ;
321  /** Save model to file, returns 0 if success.  This is designed for
322      use outside algorithms so does not save iterating arrays etc.
323  It does not save any messaging information.
324  Does not save scaling values.
325  It does not know about all types of virtual functions.
326  */
327  int saveModel(const char * fileName);
328  /** Restore model from file, returns 0 if success,
329      deletes current model */
330  int restoreModel(const char * fileName);
331 
332  /** Just check solution (for external use) - sets sum of
333      infeasibilities etc */
334  void checkSolution();
335  /// Useful row length arrays (0,1,2,3,4,5)
336  inline CoinIndexedVector * rowArray(int index) const
337  { return rowArray_[index];};
338  /// Useful column length arrays (0,1,2,3,4,5)
339  inline CoinIndexedVector * columnArray(int index) const
340  { return columnArray_[index];};
341  //@}
342
343  /******************** End of most useful part **************/
344  /**@name Functions less likely to be useful to casual user */
345  //@{
346  /** Given an existing factorization computes and checks
347      primal and dual solutions.  Uses input arrays for variables at
348      bounds.  Returns feasibility states */
349  int getSolution (  const double * rowActivities,
350                     const double * columnActivities);
351  /** Given an existing factorization computes and checks
352      primal and dual solutions.  Uses current problem arrays for
353      bounds.  Returns feasibility states */
354  int getSolution ();
355  /** Constructs a non linear cost from list of non-linearities (columns only)
356      First lower of each column is taken as real lower
357      Last lower is taken as real upper and cost ignored
358
359      Returns nonzero if bad data e.g. lowers not monotonic
360  */
361  int createPiecewiseLinearCosts(const int * starts,
362                   const double * lower, const double * gradient);
363  /** Return model - updates any scalars */
364  void returnModel(ClpSimplex & otherModel);
365  /** Factorizes using current basis. 
366      solveType - 1 iterating, 0 initial, -1 external
367      If 10 added then in primal values pass
368      Return codes are as from ClpFactorization unless initial factorization
369      when total number of singularities is returned
370  */
371    /// Save data
372    ClpDataSave saveData() ;
373    /// Restore data
374    void restoreData(ClpDataSave saved);
375  int internalFactorize(int solveType);
376  /// Clean up status
377  void cleanStatus();
378  /// Factorizes using current basis. For external use
379  int factorize();
380  /** Computes duals from scratch. If givenDjs then
381      allows for nonzero basic djs */
382  void computeDuals(double * givenDjs);
383  /// Computes primals from scratch
384  void computePrimals (  const double * rowActivities,
385                     const double * columnActivities);
386  /**
387     Unpacks one column of the matrix into indexed array
388     Uses sequenceIn_
389     Also applies scaling if needed
390  */
391  void unpack(CoinIndexedVector * rowArray) const ;
392  /**
393     Unpacks one column of the matrix into indexed array
394     Slack if sequence>= numberColumns
395     Also applies scaling if needed
396  */
397  void unpack(CoinIndexedVector * rowArray,int sequence) const;
398  /**
399     Unpacks one column of the matrix into indexed array
400     ** as packed vector
401     Uses sequenceIn_
402     Also applies scaling if needed
403  */
404  void unpackPacked(CoinIndexedVector * rowArray) ;
405  /**
406     Unpacks one column of the matrix into indexed array
407     ** as packed vector
408     Slack if sequence>= numberColumns
409     Also applies scaling if needed
410  */
411  void unpackPacked(CoinIndexedVector * rowArray,int sequence);
412 
413  /**
414      This does basis housekeeping and does values for in/out variables.
415      Can also decide to re-factorize
416  */
417  int housekeeping(double objectiveChange);
418  /** This sets largest infeasibility and most infeasible and sum
419      and number of infeasibilities (Primal) */
420  void checkPrimalSolution(const double * rowActivities=NULL,
421                           const double * columnActivies=NULL);
422  /** This sets largest infeasibility and most infeasible and sum
423      and number of infeasibilities (Dual) */
424  void checkDualSolution();
425  /** For advanced use.  When doing iterative solves things can get
426      nasty so on values pass if incoming solution has largest
427      infeasibility < incomingInfeasibility throw out variables
428      from basis until largest infeasibility < allowedInfeasibility
429      or incoming largest infeasibility.
430      If allowedInfeasibility>= incomingInfeasibility this is
431      always possible altough you may end up with an all slack basis.
432
433      Defaults are 1.0,10.0
434  */
435  void setValuesPassAction(float incomingInfeasibility,
436                           float allowedInfeasibility);
437  //@}
438  /**@name Matrix times vector methods
439     They can be faster if scalar is +- 1
440     These are covers so user need not worry about scaling
441     Also for simplex I am not using basic/non-basic split */
442  //@{
443    /** Return <code>y + A * x * scalar</code> in <code>y</code>.
444        @pre <code>x</code> must be of size <code>numColumns()</code>
445        @pre <code>y</code> must be of size <code>numRows()</code> */
446   void times(double scalar,
447                       const double * x, double * y) const;
448    /** Return <code>y + x * scalar * A</code> in <code>y</code>.
449        @pre <code>x</code> must be of size <code>numRows()</code>
450        @pre <code>y</code> must be of size <code>numColumns()</code> */
451    void transposeTimes(double scalar,
452                                const double * x, double * y) const ;
453  //@}
454
455  /**@name most useful gets and sets */
456  //@{
457  /// Worst column primal infeasibility
458  inline double columnPrimalInfeasibility() const 
459          { return columnPrimalInfeasibility_;} ;
460  /// Sequence of worst (-1 if feasible)
461  inline int columnPrimalSequence() const 
462          { return columnPrimalSequence_;} ;
463  /// Worst row primal infeasibility
464  inline double rowPrimalInfeasibility() const 
465          { return rowPrimalInfeasibility_;} ;
466  /// Sequence of worst (-1 if feasible)
467  inline int rowPrimalSequence() const 
468          { return rowPrimalSequence_;} ;
469  /** Worst column dual infeasibility (note - these may not be as meaningful
470      if the problem is primal infeasible */
471  inline double columnDualInfeasibility() const 
472          { return columnDualInfeasibility_;} ;
473  /// Sequence of worst (-1 if feasible)
474  inline int columnDualSequence() const 
475          { return columnDualSequence_;} ;
476  /// Worst row dual infeasibility
477  inline double rowDualInfeasibility() const 
478          { return rowDualInfeasibility_;} ;
479  /// Sequence of worst (-1 if feasible)
480  inline int rowDualSequence() const 
481          { return rowDualSequence_;} ;
482  /// Primal tolerance needed to make dual feasible (<largeTolerance)
483  inline double primalToleranceToGetOptimal() const 
484          { return primalToleranceToGetOptimal_;} ;
485  /// Remaining largest dual infeasibility
486  inline double remainingDualInfeasibility() const 
487          { return remainingDualInfeasibility_;} ;
488  /// Large bound value (for complementarity etc)
489  inline double largeValue() const 
490          { return largeValue_;} ;
491  void setLargeValue( double value) ;
492  /// Largest error on Ax-b
493  inline double largestPrimalError() const
494          { return largestPrimalError_;} ;
495  /// Largest error on basic duals
496  inline double largestDualError() const
497          { return largestDualError_;} ;
498  /// Largest difference between input primal solution and computed
499  inline double largestSolutionError() const
500          { return largestSolutionError_;} ;
501  /// Basic variables pivoting on which rows
502  inline int * pivotVariable() const
503          { return pivotVariable_;};
504  /// Scaling of objective 12345.0 (auto), 0.0 (off), other user
505  inline double objectiveScale() const 
506          { return objectiveScale_;} ;
507  inline void setObjectiveScale(double value)
508          { objectiveScale_ = value;} ;
509  /// Current dual tolerance
510  inline double currentDualTolerance() const 
511          { return dualTolerance_;} ;
512  inline void setCurrentDualTolerance(double value)
513          { dualTolerance_ = value;} ;
514  /// Current primal tolerance
515  inline double currentPrimalTolerance() const 
516          { return primalTolerance_;} ;
517  inline void setCurrentPrimalTolerance(double value)
518          { primalTolerance_ = value;} ;
519  /// How many iterative refinements to do
520  inline int numberRefinements() const 
521          { return numberRefinements_;} ;
522  void setNumberRefinements( int value) ;
523  /// Alpha (pivot element) for use by classes e.g. steepestedge
524  inline double alpha() const { return alpha_;};
525  /// Reduced cost of last incoming for use by classes e.g. steepestedge
526  inline double dualIn() const { return dualIn_;};
527  /// Pivot Row for use by classes e.g. steepestedge
528  inline int pivotRow() const{ return pivotRow_;};
529  /// value of incoming variable (in Dual)
530  double valueIncomingDual() const;
531  //@}
532
533  protected:
534  /**@name protected methods */
535  //@{
536  /** May change basis and then returns number changed.
537      Computation of solutions may be overriden by given pi and solution
538  */
539  int gutsOfSolution ( double * givenDuals,
540                       const double * givenPrimals,
541                       bool valuesPass=false);
542  /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
543  void gutsOfDelete(int type);
544  /// Does most of copying
545  void gutsOfCopy(const ClpSimplex & rhs);
546  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
547      1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
548      8 bit does solution scaling in
549      16 bit does rowArray and columnArray indexed vectors
550      and makes row copy if wanted, also sets columnStart_ etc
551      Also creates scaling arrays if needed.  It does scaling if needed.
552      16 also moves solutions etc in to work arrays
553      On 16 returns false if problem "bad" i.e. matrix or bounds bad
554  */
555  bool createRim(int what,bool makeRowCopy=false);
556  /** releases above arrays and does solution scaling out.  May also
557      get rid of factorization data -
558      0 get rid of nothing, 1 get rid of arrays, 2 also factorization
559  */
560  void deleteRim(int getRidOfFactorizationData=2);
561  /// Sanity check on input rim data (after scaling) - returns true if okay
562  bool sanityCheck();
563  //@}
564  public:
565  /**@name public methods */
566  //@{
567  /** Return row or column sections - not as much needed as it
568      once was.  These just map into single arrays */
569  inline double * solutionRegion(int section)
570  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
571  inline double * djRegion(int section)
572  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
573  inline double * lowerRegion(int section)
574  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
575  inline double * upperRegion(int section)
576  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
577  inline double * costRegion(int section)
578  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
579  /// Return region as single array
580  inline double * solutionRegion()
581  { return solution_;};
582  inline double * djRegion()
583  { return dj_;};
584  inline double * lowerRegion()
585  { return lower_;};
586  inline double * upperRegion()
587  { return upper_;};
588  inline double * costRegion()
589  { return cost_;};
590  inline Status getStatus(int sequence) const
591  {return static_cast<Status> (status_[sequence]&7);};
592  inline void setStatus(int sequence, Status status)
593  {
594    unsigned char & st_byte = status_[sequence];
595    st_byte &= ~7;
596    st_byte |= status;
597  };
598  /** Normally the first factorization does sparse coding because
599      the factorization could be singular.  This allows initial dense
600      factorization when it is known to be safe
601  */
602  void setInitialDenseFactorization(bool onOff);
603  bool  initialDenseFactorization() const;
604  /** Return sequence In or Out */
605  inline int sequenceIn() const
606  {return sequenceIn_;};
607  inline int sequenceOut() const
608  {return sequenceOut_;};
609  /** Set sequenceIn or Out */
610  inline void  setSequenceIn(int sequence)
611  { sequenceIn_=sequence;};
612  inline void  setSequenceOut(int sequence)
613  { sequenceOut_=sequence;};
614  /** Return direction In or Out */
615  inline int directionIn() const
616  {return directionIn_;};
617  inline int directionOut() const
618  {return directionOut_;};
619  /** Set directionIn or Out */
620  inline void  setDirectionIn(int direction)
621  { directionIn_=direction;};
622  inline void  setDirectionOut(int direction)
623  { directionOut_=direction;};
624  /// Returns 1 if sequence indicates column
625  inline int isColumn(int sequence) const
626  { return sequence<numberColumns_ ? 1 : 0;};
627  /// Returns sequence number within section
628  inline int sequenceWithin(int sequence) const
629  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
630  /// Return row or column values
631  inline double solution(int sequence)
632  { return solution_[sequence];};
633  /// Return address of row or column values
634  inline double & solutionAddress(int sequence)
635  { return solution_[sequence];};
636  inline double reducedCost(int sequence)
637   { return dj_[sequence];};
638  inline double & reducedCostAddress(int sequence)
639   { return dj_[sequence];};
640  inline double lower(int sequence)
641  { return lower_[sequence];};
642  /// Return address of row or column lower bound
643  inline double & lowerAddress(int sequence)
644  { return lower_[sequence];};
645  inline double upper(int sequence)
646  { return upper_[sequence];};
647  /// Return address of row or column upper bound
648  inline double & upperAddress(int sequence)
649  { return upper_[sequence];};
650  inline double cost(int sequence)
651  { return cost_[sequence];};
652  /// Return address of row or column cost
653  inline double & costAddress(int sequence)
654  { return cost_[sequence];};
655  /// Return original lower bound
656  inline double originalLower(int iSequence) const
657  { if (iSequence<numberColumns_) return columnLower_[iSequence]; else
658    return rowLower_[iSequence-numberColumns_];};
659  /// Return original lower bound
660  inline double originalUpper(int iSequence) const
661  { if (iSequence<numberColumns_) return columnUpper_[iSequence]; else
662    return rowUpper_[iSequence-numberColumns_];};
663  /// Theta (pivot change)
664  inline double theta() const
665  { return theta_;};
666  /// Scaling
667  const double * rowScale() const {return rowScale_;};
668  const double * columnScale() const {return columnScale_;};
669  void setRowScale(double * scale) { rowScale_ = scale;};
670  void setColumnScale(double * scale) { columnScale_ = scale;};
671  /// Return pointer to details of costs
672  inline ClpNonLinearCost * nonLinearCost() const
673  { return nonLinearCost_;};
674  //@}
675  /**@name status methods */
676  //@{
677  inline void setFakeBound(int sequence, FakeBound fakeBound)
678  {
679    unsigned char & st_byte = status_[sequence];
680    st_byte &= ~24;
681    st_byte |= fakeBound<<3;
682  };
683  inline FakeBound getFakeBound(int sequence) const
684  {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
685  inline void setRowStatus(int sequence, Status status)
686  {
687    unsigned char & st_byte = status_[sequence+numberColumns_];
688    st_byte &= ~7;
689    st_byte |= status;
690  };
691  inline Status getRowStatus(int sequence) const
692  {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
693  inline void setColumnStatus(int sequence, Status status)
694  {
695    unsigned char & st_byte = status_[sequence];
696    st_byte &= ~7;
697    st_byte |= status;
698  };
699  inline Status getColumnStatus(int sequence) const
700  {return static_cast<Status> (status_[sequence]&7);};
701  inline void setPivoted( int sequence)
702  { status_[sequence] |= 32;};
703  inline void clearPivoted( int sequence)
704  { status_[sequence] &= ~32; };
705  inline bool pivoted(int sequence) const
706  {return (((status_[sequence]>>5)&1)!=0);};
707  inline void setFlagged( int sequence)
708  {
709    status_[sequence] |= 64;
710  };
711  inline void clearFlagged( int sequence)
712  {
713    status_[sequence] &= ~64;
714  };
715  inline bool flagged(int sequence) const
716  {return (((status_[sequence]>>6)&1)!=0);};
717  /** Set up status array (can be used by OsiClp).
718      Also can be used to set up all slack basis */
719  void createStatus() ;
720  inline void allSlackBasis()
721  { createStatus();};
722   
723  /// So we know when to be cautious
724  inline int lastBadIteration() const
725  {return lastBadIteration_;};
726  /// Progress flag - at present 0 bit says artificials out
727  inline int progressFlag() const
728  {return progressFlag_;};
729  /// Force re-factorization early
730  inline void forceFactorization(int value)
731  { forceFactorization_ = value;};
732  /// Raw objective value (so always minimize in primal)
733  inline double rawObjectiveValue() const
734  { return objectiveValue_;};
735  //@}
736
737////////////////// data //////////////////
738protected:
739
740  /**@name data.  Many arrays have a row part and a column part.
741   There is a single array with both - columns then rows and
742   then normally two arrays pointing to rows and columns.  The
743   single array is the owner of memory
744  */
745  //@{
746  /// Worst column primal infeasibility
747  double columnPrimalInfeasibility_;
748  /// Worst row primal infeasibility
749  double rowPrimalInfeasibility_;
750  /// Sequence of worst (-1 if feasible)
751  int columnPrimalSequence_;
752  /// Sequence of worst (-1 if feasible)
753  int rowPrimalSequence_;
754  /// Worst column dual infeasibility
755  double columnDualInfeasibility_;
756  /// Worst row dual infeasibility
757  double rowDualInfeasibility_;
758  /// Sequence of worst (-1 if feasible)
759  int columnDualSequence_;
760  /// Sequence of worst (-1 if feasible)
761  int rowDualSequence_;
762  /// Primal tolerance needed to make dual feasible (<largeTolerance)
763  double primalToleranceToGetOptimal_;
764  /// Remaining largest dual infeasibility
765  double remainingDualInfeasibility_;
766  /// Large bound value (for complementarity etc)
767  double largeValue_;
768  /// Scaling of objective
769  double objectiveScale_;
770  /// Largest error on Ax-b
771  double largestPrimalError_;
772  /// Largest error on basic duals
773  double largestDualError_;
774  /// Largest difference between input primal solution and computed
775  double largestSolutionError_;
776  /// Dual bound
777  double dualBound_;
778  /// Alpha (pivot element)
779  double alpha_;
780  /// Theta (pivot change)
781  double theta_;
782  /// Lower Bound on In variable
783  double lowerIn_;
784  /// Value of In variable
785  double valueIn_;
786  /// Upper Bound on In variable
787  double upperIn_;
788  /// Reduced cost of In variable
789  double dualIn_;
790  /// Lower Bound on Out variable
791  double lowerOut_;
792  /// Value of Out variable
793  double valueOut_;
794  /// Upper Bound on Out variable
795  double upperOut_;
796  /// Infeasibility (dual) or ? (primal) of Out variable
797  double dualOut_;
798  /// Current dual tolerance for algorithm
799  double dualTolerance_;
800  /// Current primal tolerance for algorithm
801  double primalTolerance_;
802  /// Sum of dual infeasibilities
803  double sumDualInfeasibilities_;
804  /// Sum of primal infeasibilities
805  double sumPrimalInfeasibilities_;
806  /// Weight assigned to being infeasible in primal
807  double infeasibilityCost_;
808  /// Sum of Dual infeasibilities using tolerance based on error in duals
809  double sumOfRelaxedDualInfeasibilities_;
810  /// Sum of Primal infeasibilities using tolerance based on error in primals
811  double sumOfRelaxedPrimalInfeasibilities_;
812  /// Working copy of lower bounds (Owner of arrays below)
813  double * lower_;
814  /// Row lower bounds - working copy
815  double * rowLowerWork_;
816  /// Column lower bounds - working copy
817  double * columnLowerWork_;
818  /// Working copy of upper bounds (Owner of arrays below)
819  double * upper_;
820  /// Row upper bounds - working copy
821  double * rowUpperWork_;
822  /// Column upper bounds - working copy
823  double * columnUpperWork_;
824  /// Working copy of objective (Owner of arrays below)
825  double * cost_;
826  /// Row objective - working copy
827  double * rowObjectiveWork_;
828  /// Column objective - working copy
829  double * objectiveWork_;
830  /// Useful row length arrays
831  CoinIndexedVector * rowArray_[6];
832  /// Useful column length arrays
833  CoinIndexedVector * columnArray_[6];
834  /// Sequence of In variable
835  int sequenceIn_;
836  /// Direction of In, 1 going up, -1 going down, 0 not a clude
837  int directionIn_;
838  /// Sequence of Out variable
839  int sequenceOut_;
840  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
841  int directionOut_;
842  /// Pivot Row
843  int pivotRow_;
844  /// Last good iteration (immediately after a re-factorization)
845  int lastGoodIteration_;
846  /// Working copy of reduced costs (Owner of arrays below)
847  double * dj_;
848  /// Reduced costs of slacks not same as duals (or - duals)
849  double * rowReducedCost_;
850  /// Possible scaled reduced costs
851  double * reducedCostWork_;
852  /// Working copy of primal solution (Owner of arrays below)
853  double * solution_;
854  /// Row activities - working copy
855  double * rowActivityWork_;
856  /// Column activities - working copy
857  double * columnActivityWork_;
858  /// Number of dual infeasibilities
859  int numberDualInfeasibilities_;
860  /// Number of dual infeasibilities (without free)
861  int numberDualInfeasibilitiesWithoutFree_;
862  /// Number of primal infeasibilities
863  int numberPrimalInfeasibilities_;
864  /// How many iterative refinements to do
865  int numberRefinements_;
866  /// dual row pivot choice
867  ClpDualRowPivot * dualRowPivot_;
868  /// primal column pivot choice
869  ClpPrimalColumnPivot * primalColumnPivot_;
870  /// Basic variables pivoting on which rows
871  int * pivotVariable_;
872  /// factorization
873  ClpFactorization * factorization_;
874  /// Row scale factors for matrix
875  // ****** get working simply then make coding more efficient
876  // on full matrix operations
877  double * rowScale_;
878  /// Saved version of solution
879  double * savedSolution_;
880  /// Column scale factors
881  double * columnScale_;
882  /// Scale flag, 0 none, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic
883  int scalingFlag_;
884  /// Number of times code has tentatively thought optimal
885  int numberTimesOptimal_;
886  /// If change has been made (first attempt at stopping looping)
887  int changeMade_;
888  /// Algorithm >0 == Primal, <0 == Dual
889  int algorithm_;
890  /** Now for some reliability aids
891      This forces re-factorization early */
892  int forceFactorization_;
893  /** Perturbation:
894      -50 to +50 - perturb by this power of ten (-6 sounds good)
895      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
896      101 - we are perturbed
897      102 - don't try perturbing again
898      default is 100
899  */
900  int perturbation_;
901  /// Saved status regions
902  unsigned char * saveStatus_;
903  /** Very wasteful way of dealing with infeasibilities in primal.
904      However it will allow non-linearities and use of dual
905      analysis.  If it doesn't work it can easily be replaced.
906  */
907  ClpNonLinearCost * nonLinearCost_;
908  /** For advanced options
909      1 - Don't keep changing infeasibility weight
910      2 - Keep nonLinearCost round solves
911      4 - Force outgoing variables to exact bound (primal)
912      8 - Safe to use dense initial factorization
913  */
914  unsigned int specialOptions_;
915  /// So we know when to be cautious
916  int lastBadIteration_;
917  /// Can be used for count of fake bounds (dual) or fake costs (primal)
918  int numberFake_;
919  /// Progress flag - at present 0 bit says artificials out, 1 free in
920  int progressFlag_;
921  /// First free/super-basic variable (-1 if none)
922  int firstFree_;
923  /** For advanced use.  When doing iterative solves things can get
924      nasty so on values pass if incoming solution has largest
925      infeasibility < incomingInfeasibility throw out variables
926      from basis until largest infeasibility < allowedInfeasibility.
927      if allowedInfeasibility>= incomingInfeasibility this is
928      always possible altough you may end up with an all slack basis.
929
930      Defaults are 1.0,10.0
931  */
932  float incomingInfeasibility_;
933  float allowedInfeasibility_;
934  /// For dealing with all issues of cycling etc
935  ClpSimplexProgress * progress_;
936  //@}
937};
938//#############################################################################
939/** A function that tests the methods in the ClpSimplex class. The
940    only reason for it not to be a member method is that this way it doesn't
941    have to be compiled into the library. And that's a gain, because the
942    library should be compiled with optimization on, but this method should be
943    compiled with debugging.
944
945    It also does some testing of ClpFactorization class
946 */
947void
948ClpSimplexUnitTest(const std::string & mpsDir,
949                   const std::string & netlibDir);
950
951
952/// For saving extra information to see if looping.
953class ClpSimplexProgress {
954
955public:
956
957
958  /**@name Constructors and destructor and copy */
959  //@{
960  /// Default constructor
961    ClpSimplexProgress (  );
962
963  /// Constructor from model
964    ClpSimplexProgress ( ClpSimplex * model );
965
966  /// Copy constructor.
967  ClpSimplexProgress(const ClpSimplexProgress &);
968
969  /// Assignment operator. This copies the data
970    ClpSimplexProgress & operator=(const ClpSimplexProgress & rhs);
971  /// Destructor
972   ~ClpSimplexProgress (  );
973  //@}
974
975  /**@name Check progress */
976  //@{
977  /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
978      >=0 if give up and use as problem status
979  */
980    int looping (  );
981  /// Start check at beginning of whileIterating
982  void startCheck();
983  /// Returns cycle length in whileIterating
984  int cycle(int in, int out,int wayIn,int wayOut); 
985
986  /// Returns previous objective (if -1) - current if (0)
987  double lastObjective(int back=1) const;
988  /// Modify objective e.g. if dual infeasible in dual
989  void modifyObjective(double value);
990  /// Returns previous iteration number (if -1) - current if (0)
991  int lastIterationNumber(int back=1) const;
992
993  //@}
994  /**@name Data  */
995#define CLP_PROGRESS 5
996  //@{
997  /// Objective values
998  double objective_[CLP_PROGRESS];
999  /// Sum of infeasibilities for algorithm
1000  double infeasibility_[CLP_PROGRESS];
1001  /// Pointer back to model so we can get information
1002  ClpSimplex * model_;
1003  /// Number of infeasibilities
1004  int numberInfeasibilities_[CLP_PROGRESS];
1005  /// Iteration number at which occurred
1006  int iterationNumber_[CLP_PROGRESS];
1007  /// Number of times checked (so won't stop too early)
1008  int numberTimes_;
1009  /// Number of times it looked like loop
1010  int numberBadTimes_;
1011#define CLP_CYCLE 12
1012  /// For cycle checking
1013  int in_[CLP_CYCLE];
1014  int out_[CLP_CYCLE];
1015  char way_[CLP_CYCLE];
1016  //@}
1017};
1018#endif
Note: See TracBrowser for help on using the repository browser.