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

Last change on this file since 223 was 223, checked in by forrest, 17 years ago

Stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.5 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  //@}
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  /// To flag a variable
708  inline void setFlagged( int sequence)
709  {
710    status_[sequence] |= 64;
711  };
712  inline void clearFlagged( int sequence)
713  {
714    status_[sequence] &= ~64;
715  };
716  inline bool flagged(int sequence) const
717  {return ((status_[sequence]&64)!=0);};
718  /// To say variable active in primal pivot row choice
719  inline void setActive( int sequence)
720  {
721    status_[sequence] |= 128;
722  };
723  inline void clearActive( int sequence)
724  {
725    status_[sequence] &= ~128;
726  };
727  inline bool active(int sequence) const
728  {return ((status_[sequence]&128)!=0);};
729  /** Set up status array (can be used by OsiClp).
730      Also can be used to set up all slack basis */
731  void createStatus() ;
732  inline void allSlackBasis()
733  { createStatus();};
734   
735  /// So we know when to be cautious
736  inline int lastBadIteration() const
737  {return lastBadIteration_;};
738  /// Progress flag - at present 0 bit says artificials out
739  inline int progressFlag() const
740  {return progressFlag_;};
741  /// Force re-factorization early
742  inline void forceFactorization(int value)
743  { forceFactorization_ = value;};
744  /// Raw objective value (so always minimize in primal)
745  inline double rawObjectiveValue() const
746  { return objectiveValue_;};
747  //@}
748
749////////////////// data //////////////////
750protected:
751
752  /**@name data.  Many arrays have a row part and a column part.
753   There is a single array with both - columns then rows and
754   then normally two arrays pointing to rows and columns.  The
755   single array is the owner of memory
756  */
757  //@{
758  /// Worst column primal infeasibility
759  double columnPrimalInfeasibility_;
760  /// Worst row primal infeasibility
761  double rowPrimalInfeasibility_;
762  /// Sequence of worst (-1 if feasible)
763  int columnPrimalSequence_;
764  /// Sequence of worst (-1 if feasible)
765  int rowPrimalSequence_;
766  /// Worst column dual infeasibility
767  double columnDualInfeasibility_;
768  /// Worst row dual infeasibility
769  double rowDualInfeasibility_;
770  /// Sequence of worst (-1 if feasible)
771  int columnDualSequence_;
772  /// Sequence of worst (-1 if feasible)
773  int rowDualSequence_;
774  /// Primal tolerance needed to make dual feasible (<largeTolerance)
775  double primalToleranceToGetOptimal_;
776  /// Remaining largest dual infeasibility
777  double remainingDualInfeasibility_;
778  /// Large bound value (for complementarity etc)
779  double largeValue_;
780  /// Scaling of objective
781  double objectiveScale_;
782  /// Largest error on Ax-b
783  double largestPrimalError_;
784  /// Largest error on basic duals
785  double largestDualError_;
786  /// Largest difference between input primal solution and computed
787  double largestSolutionError_;
788  /// Dual bound
789  double dualBound_;
790  /// Alpha (pivot element)
791  double alpha_;
792  /// Theta (pivot change)
793  double theta_;
794  /// Lower Bound on In variable
795  double lowerIn_;
796  /// Value of In variable
797  double valueIn_;
798  /// Upper Bound on In variable
799  double upperIn_;
800  /// Reduced cost of In variable
801  double dualIn_;
802  /// Lower Bound on Out variable
803  double lowerOut_;
804  /// Value of Out variable
805  double valueOut_;
806  /// Upper Bound on Out variable
807  double upperOut_;
808  /// Infeasibility (dual) or ? (primal) of Out variable
809  double dualOut_;
810  /// Current dual tolerance for algorithm
811  double dualTolerance_;
812  /// Current primal tolerance for algorithm
813  double primalTolerance_;
814  /// Sum of dual infeasibilities
815  double sumDualInfeasibilities_;
816  /// Sum of primal infeasibilities
817  double sumPrimalInfeasibilities_;
818  /// Weight assigned to being infeasible in primal
819  double infeasibilityCost_;
820  /// Sum of Dual infeasibilities using tolerance based on error in duals
821  double sumOfRelaxedDualInfeasibilities_;
822  /// Sum of Primal infeasibilities using tolerance based on error in primals
823  double sumOfRelaxedPrimalInfeasibilities_;
824  /// Working copy of lower bounds (Owner of arrays below)
825  double * lower_;
826  /// Row lower bounds - working copy
827  double * rowLowerWork_;
828  /// Column lower bounds - working copy
829  double * columnLowerWork_;
830  /// Working copy of upper bounds (Owner of arrays below)
831  double * upper_;
832  /// Row upper bounds - working copy
833  double * rowUpperWork_;
834  /// Column upper bounds - working copy
835  double * columnUpperWork_;
836  /// Working copy of objective (Owner of arrays below)
837  double * cost_;
838  /// Row objective - working copy
839  double * rowObjectiveWork_;
840  /// Column objective - working copy
841  double * objectiveWork_;
842  /// Useful row length arrays
843  CoinIndexedVector * rowArray_[6];
844  /// Useful column length arrays
845  CoinIndexedVector * columnArray_[6];
846  /// Sequence of In variable
847  int sequenceIn_;
848  /// Direction of In, 1 going up, -1 going down, 0 not a clude
849  int directionIn_;
850  /// Sequence of Out variable
851  int sequenceOut_;
852  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
853  int directionOut_;
854  /// Pivot Row
855  int pivotRow_;
856  /// Last good iteration (immediately after a re-factorization)
857  int lastGoodIteration_;
858  /// Working copy of reduced costs (Owner of arrays below)
859  double * dj_;
860  /// Reduced costs of slacks not same as duals (or - duals)
861  double * rowReducedCost_;
862  /// Possible scaled reduced costs
863  double * reducedCostWork_;
864  /// Working copy of primal solution (Owner of arrays below)
865  double * solution_;
866  /// Row activities - working copy
867  double * rowActivityWork_;
868  /// Column activities - working copy
869  double * columnActivityWork_;
870  /// Number of dual infeasibilities
871  int numberDualInfeasibilities_;
872  /// Number of dual infeasibilities (without free)
873  int numberDualInfeasibilitiesWithoutFree_;
874  /// Number of primal infeasibilities
875  int numberPrimalInfeasibilities_;
876  /// How many iterative refinements to do
877  int numberRefinements_;
878  /// dual row pivot choice
879  ClpDualRowPivot * dualRowPivot_;
880  /// primal column pivot choice
881  ClpPrimalColumnPivot * primalColumnPivot_;
882  /// Basic variables pivoting on which rows
883  int * pivotVariable_;
884  /// factorization
885  ClpFactorization * factorization_;
886  /// Row scale factors for matrix
887  // ****** get working simply then make coding more efficient
888  // on full matrix operations
889  double * rowScale_;
890  /// Saved version of solution
891  double * savedSolution_;
892  /// Column scale factors
893  double * columnScale_;
894  /// Scale flag, 0 none, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic
895  int scalingFlag_;
896  /// Number of times code has tentatively thought optimal
897  int numberTimesOptimal_;
898  /// If change has been made (first attempt at stopping looping)
899  int changeMade_;
900  /// Algorithm >0 == Primal, <0 == Dual
901  int algorithm_;
902  /** Now for some reliability aids
903      This forces re-factorization early */
904  int forceFactorization_;
905  /** Perturbation:
906      -50 to +50 - perturb by this power of ten (-6 sounds good)
907      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
908      101 - we are perturbed
909      102 - don't try perturbing again
910      default is 100
911  */
912  int perturbation_;
913  /// Saved status regions
914  unsigned char * saveStatus_;
915  /** Very wasteful way of dealing with infeasibilities in primal.
916      However it will allow non-linearities and use of dual
917      analysis.  If it doesn't work it can easily be replaced.
918  */
919  ClpNonLinearCost * nonLinearCost_;
920  /** For advanced options
921      1 - Don't keep changing infeasibility weight
922      2 - Keep nonLinearCost round solves
923      4 - Force outgoing variables to exact bound (primal)
924      8 - Safe to use dense initial factorization
925  */
926  unsigned int specialOptions_;
927  /// So we know when to be cautious
928  int lastBadIteration_;
929  /// Can be used for count of fake bounds (dual) or fake costs (primal)
930  int numberFake_;
931  /// Progress flag - at present 0 bit says artificials out, 1 free in
932  int progressFlag_;
933  /// First free/super-basic variable (-1 if none)
934  int firstFree_;
935  /** For advanced use.  When doing iterative solves things can get
936      nasty so on values pass if incoming solution has largest
937      infeasibility < incomingInfeasibility throw out variables
938      from basis until largest infeasibility < allowedInfeasibility.
939      if allowedInfeasibility>= incomingInfeasibility this is
940      always possible altough you may end up with an all slack basis.
941
942      Defaults are 1.0,10.0
943  */
944  float incomingInfeasibility_;
945  float allowedInfeasibility_;
946  /// For dealing with all issues of cycling etc
947  ClpSimplexProgress * progress_;
948  //@}
949};
950//#############################################################################
951/** A function that tests the methods in the ClpSimplex class. The
952    only reason for it not to be a member method is that this way it doesn't
953    have to be compiled into the library. And that's a gain, because the
954    library should be compiled with optimization on, but this method should be
955    compiled with debugging.
956
957    It also does some testing of ClpFactorization class
958 */
959void
960ClpSimplexUnitTest(const std::string & mpsDir,
961                   const std::string & netlibDir);
962
963
964/// For saving extra information to see if looping.
965class ClpSimplexProgress {
966
967public:
968
969
970  /**@name Constructors and destructor and copy */
971  //@{
972  /// Default constructor
973    ClpSimplexProgress (  );
974
975  /// Constructor from model
976    ClpSimplexProgress ( ClpSimplex * model );
977
978  /// Copy constructor.
979  ClpSimplexProgress(const ClpSimplexProgress &);
980
981  /// Assignment operator. This copies the data
982    ClpSimplexProgress & operator=(const ClpSimplexProgress & rhs);
983  /// Destructor
984   ~ClpSimplexProgress (  );
985  //@}
986
987  /**@name Check progress */
988  //@{
989  /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
990      >=0 if give up and use as problem status
991  */
992    int looping (  );
993  /// Start check at beginning of whileIterating
994  void startCheck();
995  /// Returns cycle length in whileIterating
996  int cycle(int in, int out,int wayIn,int wayOut); 
997
998  /// Returns previous objective (if -1) - current if (0)
999  double lastObjective(int back=1) const;
1000  /// Modify objective e.g. if dual infeasible in dual
1001  void modifyObjective(double value);
1002  /// Returns previous iteration number (if -1) - current if (0)
1003  int lastIterationNumber(int back=1) const;
1004
1005  //@}
1006  /**@name Data  */
1007#define CLP_PROGRESS 5
1008  //@{
1009  /// Objective values
1010  double objective_[CLP_PROGRESS];
1011  /// Sum of infeasibilities for algorithm
1012  double infeasibility_[CLP_PROGRESS];
1013  /// Pointer back to model so we can get information
1014  ClpSimplex * model_;
1015  /// Number of infeasibilities
1016  int numberInfeasibilities_[CLP_PROGRESS];
1017  /// Iteration number at which occurred
1018  int iterationNumber_[CLP_PROGRESS];
1019  /// Number of times checked (so won't stop too early)
1020  int numberTimes_;
1021  /// Number of times it looked like loop
1022  int numberBadTimes_;
1023#define CLP_CYCLE 12
1024  /// For cycle checking
1025  int in_[CLP_CYCLE];
1026  int out_[CLP_CYCLE];
1027  char way_[CLP_CYCLE];
1028  //@}
1029};
1030#endif
Note: See TracBrowser for help on using the repository browser.