source: branches/devel/Clp/src/ClpSimplex.hpp @ 918

Last change on this file since 918 was 918, checked in by forrest, 13 years ago

for messages

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.6 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;
24class CoinModel;
25class OsiClpSolverInterface;
26class CoinWarmStartBasis;
27class ClpDisasterHandler;
28
29/** This solves LPs using the simplex method
30
31    It inherits from ClpModel and all its arrays are created at
32    algorithm time. Originally I tried to work with model arrays
33    but for simplicity of coding I changed to single arrays with
34    structural variables then row variables.  Some coding is still
35    based on old style and needs cleaning up.
36
37    For a description of algorithms:
38
39    for dual see ClpSimplexDual.hpp and at top of ClpSimplexDual.cpp
40    for primal see ClpSimplexPrimal.hpp and at top of ClpSimplexPrimal.cpp
41
42    There is an algorithm data member.  + for primal variations
43    and - for dual variations
44
45    This file also includes (at end) a very simple class ClpSimplexProgress
46    which is where anti-looping stuff should migrate to
47
48*/
49
50class ClpSimplex : public ClpModel {
51   friend void ClpSimplexUnitTest(const std::string & mpsDir,
52                                  const std::string & netlibDir);
53
54public:
55  /** enums for status of various sorts.
56      First 4 match CoinWarmStartBasis,
57      isFixed means fixed at lower bound and out of basis
58  */
59  enum Status {
60    isFree = 0x00,
61    basic = 0x01,
62    atUpperBound = 0x02,
63    atLowerBound = 0x03,
64    superBasic = 0x04,
65    isFixed = 0x05
66  };
67  // For Dual
68  enum FakeBound {
69    noFake = 0x00,
70    bothFake = 0x01,
71    upperFake = 0x02,
72    lowerFake = 0x03
73  };
74
75  /**@name Constructors and destructor and copy */
76  //@{
77  /// Default constructor
78    ClpSimplex (bool emptyMessages = false  );
79
80  /** Copy constructor. May scale depending on mode
81      -1 leave mode as is
82      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
83  */
84  ClpSimplex(const ClpSimplex & rhs, int scalingMode =-1);
85  /** Copy constructor from model. May scale depending on mode
86      -1 leave mode as is
87      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
88  */
89  ClpSimplex(const ClpModel & rhs, int scalingMode=-1);
90  /** Subproblem constructor.  A subset of whole model is created from the
91      row and column lists given.  The new order is given by list order and
92      duplicates are allowed.  Name and integer information can be dropped
93      Can optionally modify rhs to take into account variables NOT in list
94      in this case duplicates are not allowed (also see getbackSolution)
95  */
96  ClpSimplex (const ClpModel * wholeModel,
97              int numberRows, const int * whichRows,
98              int numberColumns, const int * whichColumns,
99              bool dropNames=true, bool dropIntegers=true,
100              bool fixOthers=false);
101  /** This constructor modifies original ClpSimplex and stores
102      original stuff in created ClpSimplex.  It is only to be used in
103      conjunction with originalModel */
104  ClpSimplex (ClpSimplex * wholeModel,
105              int numberColumns, const int * whichColumns);
106  /** This copies back stuff from miniModel and then deletes miniModel.
107      Only to be used with mini constructor */
108  void originalModel(ClpSimplex * miniModel);
109  /** Array persistence flag
110      If 0 then as now (delete/new)
111      1 then only do arrays if bigger needed
112      2 as 1 but give a bit extra if bigger needed
113  */
114  void setPersistenceFlag(int value);
115  /**
116     If you are re-using the same matrix again and again then the setup time
117     to do scaling may be significant.  Also you may not want to initialize all values
118     or return all values (especially if infeasible).  While an auxiliary model exists
119     it will be faster.  If options -1 then model is switched off.  Otherwise switched on
120     with following options.
121     1 - rhs is constant
122     2 - bounds are constant
123     4 - objective is constant
124     8 - solution in by basis and no djs etc in
125     16 - no duals out (but reduced costs)
126     32 - no output if infeasible
127  */
128  void auxiliaryModel(int options);
129  /// Switch off e.g. if people using presolve
130  void deleteAuxiliaryModel();
131  /// See if we have auxiliary model
132  inline bool usingAuxiliaryModel() const
133  { return auxiliaryModel_!=NULL;};
134  /// Assignment operator. This copies the data
135    ClpSimplex & operator=(const ClpSimplex & rhs);
136  /// Destructor
137   ~ClpSimplex (  );
138  // Ones below are just ClpModel with some changes
139  /** Loads a problem (the constraints on the
140        rows are given by lower and upper bounds). If a pointer is 0 then the
141        following values are the default:
142        <ul>
143          <li> <code>colub</code>: all columns have upper bound infinity
144          <li> <code>collb</code>: all columns have lower bound 0
145          <li> <code>rowub</code>: all rows have upper bound infinity
146          <li> <code>rowlb</code>: all rows have lower bound -infinity
147          <li> <code>obj</code>: all variables have 0 objective coefficient
148        </ul>
149    */
150  void loadProblem (  const ClpMatrixBase& matrix,
151                     const double* collb, const double* colub,   
152                     const double* obj,
153                     const double* rowlb, const double* rowub,
154                      const double * rowObjective=NULL);
155  void loadProblem (  const CoinPackedMatrix& matrix,
156                     const double* collb, const double* colub,   
157                     const double* obj,
158                     const double* rowlb, const double* rowub,
159                      const double * rowObjective=NULL);
160
161  /** Just like the other loadProblem() method except that the matrix is
162        given in a standard column major ordered format (without gaps). */
163  void loadProblem (  const int numcols, const int numrows,
164                     const CoinBigIndex* start, const int* index,
165                     const double* value,
166                     const double* collb, const double* colub,   
167                     const double* obj,
168                      const double* rowlb, const double* rowub,
169                      const double * rowObjective=NULL);
170  /// This one is for after presolve to save memory
171  void loadProblem (  const int numcols, const int numrows,
172                     const CoinBigIndex* start, const int* index,
173                      const double* value,const int * length,
174                     const double* collb, const double* colub,   
175                     const double* obj,
176                      const double* rowlb, const double* rowub,
177                      const double * rowObjective=NULL);
178  /** This loads a model from a coinModel object - returns number of errors.
179      If keepSolution true and size is same as current then
180      keeps current status and solution
181  */
182  int loadProblem (  CoinModel & modelObject,bool keepSolution=false);
183  /// Read an mps file from the given filename
184  int readMps(const char *filename,
185              bool keepNames=false,
186              bool ignoreErrors = false);
187  /// Read GMPL files from the given filenames
188  int readGMPL(const char *filename,const char * dataName,
189               bool keepNames=false);
190  /** Borrow model.  This is so we dont have to copy large amounts
191      of data around.  It assumes a derived class wants to overwrite
192      an empty model with a real one - while it does an algorithm.
193      This is same as ClpModel one, but sets scaling on etc. */
194  void borrowModel(ClpModel & otherModel);
195  void borrowModel(ClpSimplex & otherModel);
196   /// Pass in Event handler (cloned and deleted at end)
197   void passInEventHandler(const ClpEventHandler * eventHandler);
198  /// Puts solution back into small model
199  void getbackSolution(const ClpSimplex & smallModel,const int * whichRow, const int * whichColumn);
200  //@}
201
202  /**@name Functions most useful to user */
203  //@{
204  /** General solve algorithm which can do presolve.
205      See  ClpSolve.hpp for options
206   */
207  int initialSolve(ClpSolve & options);
208  /// Default initial solve
209  int initialSolve();
210  /// Dual initial solve
211  int initialDualSolve();
212  /// Primal initial solve
213  int initialPrimalSolve();
214 /// Barrier initial solve
215  int initialBarrierSolve();
216  /// Barrier initial solve, not to be followed by crossover
217  int initialBarrierNoCrossSolve();
218  /** Dual algorithm - see ClpSimplexDual.hpp for method.
219      ifValuesPass==2 just does values pass and then stops.
220
221      startFinishOptions - bits
222      1 - do not delete work areas and factorization at end
223      2 - use old factorization if same number of rows
224      4 - skip as much initialization of work areas as possible
225          (based on whatsChanged in clpmodel.hpp) ** work in progress
226      maybe other bits later
227  */
228  int dual(int ifValuesPass=0, int startFinishOptions=0);
229  // If using Debug
230  int dualDebug(int ifValuesPass=0, int startFinishOptions=0);
231  /** Primal algorithm - see ClpSimplexPrimal.hpp for method.
232      ifValuesPass==2 just does values pass and then stops.
233
234      startFinishOptions - bits
235      1 - do not delete work areas and factorization at end
236      2 - use old factorization if same number of rows
237      4 - skip as much initialization of work areas as possible
238          (based on whatsChanged in clpmodel.hpp) ** work in progress
239      maybe other bits later
240  */
241  int primal(int ifValuesPass=0, int startFinishOptions=0);
242  /** Solves nonlinear problem using SLP - may be used as crash
243      for other algorithms when number of iterations small.
244      Also exits if all problematical variables are changing
245      less than deltaTolerance
246  */
247  int nonlinearSLP(int numberPasses,double deltaTolerance);
248  /** Solves using barrier (assumes you have good cholesky factor code).
249      Does crossover to simplex if asked*/
250  int barrier(bool crossover=true);
251  /** Solves non-linear using reduced gradient.  Phase = 0 get feasible,
252      =1 use solution */
253  int reducedGradient(int phase=0);
254  /**
255     When scaling is on it is possible that the scaled problem
256     is feasible but the unscaled is not.  Clp returns a secondary
257     status code to that effect.  This option allows for a cleanup.
258     If you use it I would suggest 1.
259     This only affects actions when scaled optimal
260     0 - no action
261     1 - clean up using dual if primal infeasibility
262     2 - clean up using dual if dual infeasibility
263     3 - clean up using dual if primal or dual infeasibility
264     11,12,13 - as 1,2,3 but use primal
265
266     return code as dual/primal
267  */
268  int cleanup(int cleanupScaling);
269  /** Dual ranging.
270      This computes increase/decrease in cost for each given variable and corresponding
271      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
272      and numberColumns.. for artificials/slacks.
273      For non-basic variables the information is trivial to compute and the change in cost is just minus the
274      reduced cost and the sequence number will be that of the non-basic variables.
275      For basic variables a ratio test is between the reduced costs for non-basic variables
276      and the row of the tableau corresponding to the basic variable.
277      The increase/decrease value is always >= 0.0
278
279      Up to user to provide correct length arrays where each array is of length numberCheck.
280      which contains list of variables for which information is desired.  All other
281      arrays will be filled in by function.  If fifth entry in which is variable 7 then fifth entry in output arrays
282      will be information for variable 7.
283
284      If valueIncrease/Decrease not NULL (both must be NULL or both non NULL) then these are filled with
285      the value of variable if such a change in cost were made (the existing bounds are ignored)
286
287      Returns non-zero if infeasible unbounded etc
288  */
289  int dualRanging(int numberCheck,const int * which,
290                  double * costIncrease, int * sequenceIncrease,
291                  double * costDecrease, int * sequenceDecrease,
292                  double * valueIncrease=NULL, double * valueDecrease=NULL);
293  /** Primal ranging.
294      This computes increase/decrease in value for each given variable and corresponding
295      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
296      and numberColumns.. for artificials/slacks.
297      This should only be used for non-basic variabls as otherwise information is pretty useless
298      For basic variables the sequence number will be that of the basic variables.
299
300      Up to user to provide correct length arrays where each array is of length numberCheck.
301      which contains list of variables for which information is desired.  All other
302      arrays will be filled in by function.  If fifth entry in which is variable 7 then fifth entry in output arrays
303      will be information for variable 7.
304
305      Returns non-zero if infeasible unbounded etc
306  */
307  int primalRanging(int numberCheck,const int * which,
308                  double * valueIncrease, int * sequenceIncrease,
309                  double * valueDecrease, int * sequenceDecrease);
310  /** Write the basis in MPS format to the specified file.
311      If writeValues true writes values of structurals
312      (and adds VALUES to end of NAME card)
313     
314      Row and column names may be null.
315      formatType is
316      <ul>
317      <li> 0 - normal
318      <li> 1 - extra accuracy
319      <li> 2 - IEEE hex (later)
320      </ul>
321     
322      Returns non-zero on I/O error
323  */
324  int writeBasis(const char *filename,
325                 bool writeValues=false,
326                 int formatType=0) const;
327  /** Read a basis from the given filename,
328      returns -1 on file error, 0 if no values, 1 if values */
329  int readBasis(const char *filename);
330  /// Returns a basis (to be deleted by user)
331  CoinWarmStartBasis * getBasis() const;
332  /// Passes in factorization
333  void setFactorization( ClpFactorization & factorization);
334  /** Tightens primal bounds to make dual faster.  Unless
335      fixed or doTight>10, bounds are slightly looser than they could be.
336      This is to make dual go faster and is probably not needed
337      with a presolve.  Returns non-zero if problem infeasible.
338
339      Fudge for branch and bound - put bounds on columns of factor *
340      largest value (at continuous) - should improve stability
341      in branch and bound on infeasible branches (0.0 is off)
342  */
343  int tightenPrimalBounds(double factor=0.0,int doTight=0);
344  /** Crash - at present just aimed at dual, returns
345      -2 if dual preferred and crash basis created
346      -1 if dual preferred and all slack basis preferred
347       0 if basis going in was not all slack
348       1 if primal preferred and all slack basis preferred
349       2 if primal preferred and crash basis created.
350       
351       if gap between bounds <="gap" variables can be flipped
352       ( If pivot -1 then can be made super basic!)
353
354       If "pivot" is
355       -1 No pivoting - always primal
356       0 No pivoting (so will just be choice of algorithm)
357       1 Simple pivoting e.g. gub
358       2 Mini iterations
359  */
360  int crash(double gap,int pivot);
361  /// Sets row pivot choice algorithm in dual
362  void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
363  /// Sets column pivot choice algorithm in primal
364  void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
365  /** For strong branching.  On input lower and upper are new bounds
366      while on output they are change in objective function values
367      (>1.0e50 infeasible).
368      Return code is 0 if nothing interesting, -1 if infeasible both
369      ways and +1 if infeasible one way (check values to see which one(s))
370      Solutions are filled in as well - even down, odd up - also
371      status and number of iterations
372  */
373  int strongBranching(int numberVariables,const int * variables,
374                      double * newLower, double * newUpper,
375                      double ** outputSolution,
376                      int * outputStatus, int * outputIterations,
377                      bool stopOnFirstInfeasible=true,
378                      bool alwaysFinish=false,
379                      int startFinishOptions=0);
380  //@}
381
382  /**@name Needed for functionality of OsiSimplexInterface */
383  //@{
384  /** Pivot in a variable and out a variable.  Returns 0 if okay,
385      1 if inaccuracy forced re-factorization, -1 if would be singular.
386      Also updates primal/dual infeasibilities.
387      Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
388  */
389  int pivot();
390
391  /** Pivot in a variable and choose an outgoing one.  Assumes primal
392      feasible - will not go through a bound.  Returns step length in theta
393      Returns ray in ray_ (or NULL if no pivot)
394      Return codes as before but -1 means no acceptable pivot
395  */
396  int primalPivotResult();
397 
398  /** Pivot out a variable and choose an incoing one.  Assumes dual
399      feasible - will not go through a reduced cost. 
400      Returns step length in theta
401      Returns ray in ray_ (or NULL if no pivot)
402      Return codes as before but -1 means no acceptable pivot
403  */
404  int dualPivotResult();
405
406  /** Common bits of coding for dual and primal.  Return 0 if okay,
407      1 if bad matrix, 2 if very bad factorization
408
409      startFinishOptions - bits
410      1 - do not delete work areas and factorization at end
411      2 - use old factorization if same number of rows
412      4 - skip as much initialization of work areas as possible
413          (based on whatsChanged in clpmodel.hpp) ** work in progress
414      maybe other bits later
415     
416  */
417  int startup(int ifValuesPass,int startFinishOptions=0);
418  void finish(int startFinishOptions=0);
419 
420  /** Factorizes and returns true if optimal.  Used by user */
421  bool statusOfProblem(bool initial=false);
422  /// If user left factorization frequency then compute
423  void defaultFactorizationFrequency();
424  //@}
425
426  /**@name most useful gets and sets */
427  //@{
428  /// If problem is primal feasible
429  inline bool primalFeasible() const
430         { return (numberPrimalInfeasibilities_==0);};
431  /// If problem is dual feasible
432  inline bool dualFeasible() const
433         { return (numberDualInfeasibilities_==0);};
434  /// factorization
435  inline ClpFactorization * factorization() const 
436          { return factorization_;};
437  /// Sparsity on or off
438  bool sparseFactorization() const;
439  void setSparseFactorization(bool value);
440  /// Factorization frequency
441  int factorizationFrequency() const;
442  void setFactorizationFrequency(int value);
443  /// Dual bound
444  inline double dualBound() const
445          { return dualBound_;};
446  void setDualBound(double value);
447  /// Infeasibility cost
448  inline double infeasibilityCost() const
449          { return infeasibilityCost_;};
450  void setInfeasibilityCost(double value);
451  /** Amount of print out:
452      0 - none
453      1 - just final
454      2 - just factorizations
455      3 - as 2 plus a bit more
456      4 - verbose
457      above that 8,16,32 etc just for selective debug
458  */
459  /** Perturbation:
460      50  - switch on perturbation
461      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
462      101 - we are perturbed
463      102 - don't try perturbing again
464      default is 100
465      others are for playing
466  */
467  inline int perturbation() const
468    { return perturbation_;};
469  void setPerturbation(int value);
470  /// Current (or last) algorithm
471  inline int algorithm() const 
472  {return algorithm_; } ;
473  /// Set algorithm
474  inline void setAlgorithm(int value)
475  {algorithm_=value; } ;
476  /// Sum of dual infeasibilities
477  inline double sumDualInfeasibilities() const 
478          { return sumDualInfeasibilities_;} ;
479  inline void setSumDualInfeasibilities(double value)
480          { sumDualInfeasibilities_=value;} ;
481  /// Sum of relaxed dual infeasibilities
482  inline double sumOfRelaxedDualInfeasibilities() const 
483          { return sumOfRelaxedDualInfeasibilities_;} ;
484  inline void setSumOfRelaxedDualInfeasibilities(double value)
485          { sumOfRelaxedDualInfeasibilities_=value;} ;
486  /// Number of dual infeasibilities
487  inline int numberDualInfeasibilities() const 
488          { return numberDualInfeasibilities_;} ;
489  inline void setNumberDualInfeasibilities(int value)
490          { numberDualInfeasibilities_=value;} ;
491  /// Number of dual infeasibilities (without free)
492  inline int numberDualInfeasibilitiesWithoutFree() const 
493          { return numberDualInfeasibilitiesWithoutFree_;} ;
494  /// Sum of primal infeasibilities
495  inline double sumPrimalInfeasibilities() const 
496          { return sumPrimalInfeasibilities_;} ;
497  inline void setSumPrimalInfeasibilities(double value)
498          { sumPrimalInfeasibilities_=value;} ;
499  /// Sum of relaxed primal infeasibilities
500  inline double sumOfRelaxedPrimalInfeasibilities() const 
501          { return sumOfRelaxedPrimalInfeasibilities_;} ;
502  inline void setSumOfRelaxedPrimalInfeasibilities(double value)
503          { sumOfRelaxedPrimalInfeasibilities_=value;} ;
504  /// Number of primal infeasibilities
505  inline int numberPrimalInfeasibilities() const 
506          { return numberPrimalInfeasibilities_;} ;
507  inline void setNumberPrimalInfeasibilities(int value)
508          { numberPrimalInfeasibilities_=value;} ;
509  /** Save model to file, returns 0 if success.  This is designed for
510      use outside algorithms so does not save iterating arrays etc.
511  It does not save any messaging information.
512  Does not save scaling values.
513  It does not know about all types of virtual functions.
514  */
515  int saveModel(const char * fileName);
516  /** Restore model from file, returns 0 if success,
517      deletes current model */
518  int restoreModel(const char * fileName);
519 
520  /** Just check solution (for external use) - sets sum of
521      infeasibilities etc.
522      If setToBounds 0 then primal column values not changed
523      and used to compute primal row activity values.  If 1 or 2
524      then status used - so all nonbasic variables set to
525      indicated bound and if any values changed (or ==2)  basic values re-computed.
526  */
527  void checkSolution(int setToBounds=false);
528  /** Just check solution (for internal use) - sets sum of
529      infeasibilities etc. */
530  void checkSolutionInternal();
531  /// Useful row length arrays (0,1,2,3,4,5)
532  inline CoinIndexedVector * rowArray(int index) const
533  { return rowArray_[index];};
534  /// Useful column length arrays (0,1,2,3,4,5)
535  inline CoinIndexedVector * columnArray(int index) const
536  { return columnArray_[index];};
537  //@}
538
539  /******************** End of most useful part **************/
540  /**@name Functions less likely to be useful to casual user */
541  //@{
542  /** Given an existing factorization computes and checks
543      primal and dual solutions.  Uses input arrays for variables at
544      bounds.  Returns feasibility states */
545  int getSolution (  const double * rowActivities,
546                     const double * columnActivities);
547  /** Given an existing factorization computes and checks
548      primal and dual solutions.  Uses current problem arrays for
549      bounds.  Returns feasibility states */
550  int getSolution ();
551  /** Constructs a non linear cost from list of non-linearities (columns only)
552      First lower of each column is taken as real lower
553      Last lower is taken as real upper and cost ignored
554
555      Returns nonzero if bad data e.g. lowers not monotonic
556  */
557  int createPiecewiseLinearCosts(const int * starts,
558                   const double * lower, const double * gradient);
559  /** Return model - updates any scalars */
560  void returnModel(ClpSimplex & otherModel);
561  /** Factorizes using current basis. 
562      solveType - 1 iterating, 0 initial, -1 external
563      If 10 added then in primal values pass
564      Return codes are as from ClpFactorization unless initial factorization
565      when total number of singularities is returned.
566      Special case is numberRows_+1 -> all slack basis.
567  */
568  int internalFactorize(int solveType);
569  /// Save data
570  ClpDataSave saveData() ;
571  /// Restore data
572  void restoreData(ClpDataSave saved);
573  /// Clean up status
574  void cleanStatus();
575  /// Factorizes using current basis. For external use
576  int factorize();
577  /** Computes duals from scratch. If givenDjs then
578      allows for nonzero basic djs */
579  void computeDuals(double * givenDjs);
580  /// Computes primals from scratch
581  void computePrimals (  const double * rowActivities,
582                     const double * columnActivities);
583  /** Adds multiple of a column into an array */
584  void add(double * array,
585                   int column, double multiplier) const;
586  /**
587     Unpacks one column of the matrix into indexed array
588     Uses sequenceIn_
589     Also applies scaling if needed
590  */
591  void unpack(CoinIndexedVector * rowArray) const ;
592  /**
593     Unpacks one column of the matrix into indexed array
594     Slack if sequence>= numberColumns
595     Also applies scaling if needed
596  */
597  void unpack(CoinIndexedVector * rowArray,int sequence) const;
598  /**
599     Unpacks one column of the matrix into indexed array
600     ** as packed vector
601     Uses sequenceIn_
602     Also applies scaling if needed
603  */
604  void unpackPacked(CoinIndexedVector * rowArray) ;
605  /**
606     Unpacks one column of the matrix into indexed array
607     ** as packed vector
608     Slack if sequence>= numberColumns
609     Also applies scaling if needed
610  */
611  void unpackPacked(CoinIndexedVector * rowArray,int sequence);
612protected: 
613  /**
614      This does basis housekeeping and does values for in/out variables.
615      Can also decide to re-factorize
616  */
617  int housekeeping(double objectiveChange);
618  /** This sets largest infeasibility and most infeasible and sum
619      and number of infeasibilities (Primal) */
620  void checkPrimalSolution(const double * rowActivities=NULL,
621                           const double * columnActivies=NULL);
622  /** This sets largest infeasibility and most infeasible and sum
623      and number of infeasibilities (Dual) */
624  void checkDualSolution();
625  /** This sets sum and number of infeasibilities (Dual and Primal) */
626  void checkBothSolutions();
627public:
628  /** For advanced use.  When doing iterative solves things can get
629      nasty so on values pass if incoming solution has largest
630      infeasibility < incomingInfeasibility throw out variables
631      from basis until largest infeasibility < allowedInfeasibility
632      or incoming largest infeasibility.
633      If allowedInfeasibility>= incomingInfeasibility this is
634      always possible altough you may end up with an all slack basis.
635
636      Defaults are 1.0,10.0
637  */
638  void setValuesPassAction(float incomingInfeasibility,
639                           float allowedInfeasibility);
640  //@}
641  /**@name most useful gets and sets */
642  //@{
643  // On reflection I doubt whether anyone uses so test
644private:
645  /// Worst column primal infeasibility
646  inline double columnPrimalInfeasibility() const 
647          { return columnPrimalInfeasibility_;} ;
648  /// Sequence of worst (-1 if feasible)
649  inline int columnPrimalSequence() const 
650          { return columnPrimalSequence_;} ;
651  /// Worst row primal infeasibility
652  inline double rowPrimalInfeasibility() const 
653          { return rowPrimalInfeasibility_;} ;
654  /// Sequence of worst (-1 if feasible)
655  inline int rowPrimalSequence() const 
656          { return rowPrimalSequence_;} ;
657  /** Worst column dual infeasibility (note - these may not be as meaningful
658      if the problem is primal infeasible */
659  inline double columnDualInfeasibility() const 
660          { return columnDualInfeasibility_;} ;
661  /// Sequence of worst (-1 if feasible)
662  inline int columnDualSequence() const 
663          { return columnDualSequence_;} ;
664  /// Worst row dual infeasibility
665  inline double rowDualInfeasibility() const 
666          { return rowDualInfeasibility_;} ;
667  /// Sequence of worst (-1 if feasible)
668  inline int rowDualSequence() const 
669          { return rowDualSequence_;} ;
670  /// Primal tolerance needed to make dual feasible (<largeTolerance)
671  inline double primalToleranceToGetOptimal() const 
672          { return primalToleranceToGetOptimal_;} ;
673  /// Remaining largest dual infeasibility
674  inline double remainingDualInfeasibility() const 
675          { return remainingDualInfeasibility_;} ;
676  /// Largest difference between input primal solution and computed
677  inline double largestSolutionError() const
678          { return largestSolutionError_;} ;
679public:
680  /// Disaster handler
681  inline void setDisasterHandler(ClpDisasterHandler * handler)
682  { disasterArea_= handler;};
683  /// Large bound value (for complementarity etc)
684  inline double largeValue() const 
685          { return largeValue_;} ;
686  void setLargeValue( double value) ;
687  /// Largest error on Ax-b
688  inline double largestPrimalError() const
689          { return largestPrimalError_;} ;
690  /// Largest error on basic duals
691  inline double largestDualError() const
692          { return largestDualError_;} ;
693  /// Largest error on Ax-b
694  inline void setLargestPrimalError(double value)
695          { largestPrimalError_=value;} ;
696  /// Largest error on basic duals
697  inline void setLargestDualError(double value)
698          { largestDualError_=value;} ;
699  /// Basic variables pivoting on which rows
700  inline int * pivotVariable() const
701          { return pivotVariable_;};
702  /// If automatic scaling on
703  inline bool automaticScaling() const
704  { return automaticScale_!=0;};
705  inline void setAutomaticScaling(bool onOff)
706  { automaticScale_ = onOff ? 1: 0;}; 
707  /// Current dual tolerance
708  inline double currentDualTolerance() const 
709          { return dualTolerance_;} ;
710  inline void setCurrentDualTolerance(double value)
711          { dualTolerance_ = value;} ;
712  /// Current primal tolerance
713  inline double currentPrimalTolerance() const 
714          { return primalTolerance_;} ;
715  inline void setCurrentPrimalTolerance(double value)
716          { primalTolerance_ = value;} ;
717  /// How many iterative refinements to do
718  inline int numberRefinements() const 
719          { return numberRefinements_;} ;
720  void setNumberRefinements( int value) ;
721  /// Alpha (pivot element) for use by classes e.g. steepestedge
722  inline double alpha() const { return alpha_;};
723  inline void setAlpha(double value) { alpha_ = value;};
724  /// Reduced cost of last incoming for use by classes e.g. steepestedge
725  inline double dualIn() const { return dualIn_;};
726  /// Pivot Row for use by classes e.g. steepestedge
727  inline int pivotRow() const{ return pivotRow_;};
728  inline void setPivotRow(int value) { pivotRow_=value;};
729  /// value of incoming variable (in Dual)
730  double valueIncomingDual() const;
731  //@}
732
733  protected:
734  /**@name protected methods */
735  //@{
736  /** May change basis and then returns number changed.
737      Computation of solutions may be overriden by given pi and solution
738  */
739  int gutsOfSolution ( double * givenDuals,
740                       const double * givenPrimals,
741                       bool valuesPass=false);
742  /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
743  void gutsOfDelete(int type);
744  /// Does most of copying
745  void gutsOfCopy(const ClpSimplex & rhs);
746  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
747      1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
748      8 bit does solution scaling in
749      16 bit does rowArray and columnArray indexed vectors
750      and makes row copy if wanted, also sets columnStart_ etc
751      Also creates scaling arrays if needed.  It does scaling if needed.
752      16 also moves solutions etc in to work arrays
753      On 16 returns false if problem "bad" i.e. matrix or bounds bad
754      If startFinishOptions is -1 then called by user in getSolution
755      so do arrays but keep pivotVariable_
756  */
757  bool createRim(int what,bool makeRowCopy=false,int startFinishOptions=0);
758  /** releases above arrays and does solution scaling out.  May also
759      get rid of factorization data -
760      0 get rid of nothing, 1 get rid of arrays, 2 also factorization
761  */
762  void deleteRim(int getRidOfFactorizationData=2);
763  /// Sanity check on input rim data (after scaling) - returns true if okay
764  bool sanityCheck();
765  //@}
766  public:
767  /**@name public methods */
768  //@{
769  /** Return row or column sections - not as much needed as it
770      once was.  These just map into single arrays */
771  inline double * solutionRegion(int section) const
772  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
773  inline double * djRegion(int section) const
774  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
775  inline double * lowerRegion(int section) const
776  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
777  inline double * upperRegion(int section) const
778  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
779  inline double * costRegion(int section) const
780  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
781  /// Return region as single array
782  inline double * solutionRegion() const
783  { return solution_;};
784  inline double * djRegion() const
785  { return dj_;};
786  inline double * lowerRegion() const
787  { return lower_;};
788  inline double * upperRegion() const
789  { return upper_;};
790  inline double * costRegion() const
791  { return cost_;};
792  inline Status getStatus(int sequence) const
793  {return static_cast<Status> (status_[sequence]&7);};
794  inline void setStatus(int sequence, Status status)
795  {
796    unsigned char & st_byte = status_[sequence];
797    st_byte &= ~7;
798    st_byte |= status;
799  };
800  /** Normally the first factorization does sparse coding because
801      the factorization could be singular.  This allows initial dense
802      factorization when it is known to be safe
803  */
804  void setInitialDenseFactorization(bool onOff);
805  bool  initialDenseFactorization() const;
806  /** Return sequence In or Out */
807  inline int sequenceIn() const
808  {return sequenceIn_;};
809  inline int sequenceOut() const
810  {return sequenceOut_;};
811  /** Set sequenceIn or Out */
812  inline void  setSequenceIn(int sequence)
813  { sequenceIn_=sequence;};
814  inline void  setSequenceOut(int sequence)
815  { sequenceOut_=sequence;};
816  /** Return direction In or Out */
817  inline int directionIn() const
818  {return directionIn_;};
819  inline int directionOut() const
820  {return directionOut_;};
821  /** Set directionIn or Out */
822  inline void  setDirectionIn(int direction)
823  { directionIn_=direction;};
824  inline void  setDirectionOut(int direction)
825  { directionOut_=direction;};
826  /// Value of Out variable
827  inline double valueOut() const
828  { return valueOut_;};
829  /// Returns 1 if sequence indicates column
830  inline int isColumn(int sequence) const
831  { return sequence<numberColumns_ ? 1 : 0;};
832  /// Returns sequence number within section
833  inline int sequenceWithin(int sequence) const
834  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
835  /// Return row or column values
836  inline double solution(int sequence)
837  { return solution_[sequence];};
838  /// Return address of row or column values
839  inline double & solutionAddress(int sequence)
840  { return solution_[sequence];};
841  inline double reducedCost(int sequence)
842   { return dj_[sequence];};
843  inline double & reducedCostAddress(int sequence)
844   { return dj_[sequence];};
845  inline double lower(int sequence)
846  { return lower_[sequence];};
847  /// Return address of row or column lower bound
848  inline double & lowerAddress(int sequence)
849  { return lower_[sequence];};
850  inline double upper(int sequence)
851  { return upper_[sequence];};
852  /// Return address of row or column upper bound
853  inline double & upperAddress(int sequence)
854  { return upper_[sequence];};
855  inline double cost(int sequence)
856  { return cost_[sequence];};
857  /// Return address of row or column cost
858  inline double & costAddress(int sequence)
859  { return cost_[sequence];};
860  /// Return original lower bound
861  inline double originalLower(int iSequence) const
862  { if (iSequence<numberColumns_) return columnLower_[iSequence]; else
863    return rowLower_[iSequence-numberColumns_];};
864  /// Return original lower bound
865  inline double originalUpper(int iSequence) const
866  { if (iSequence<numberColumns_) return columnUpper_[iSequence]; else
867    return rowUpper_[iSequence-numberColumns_];};
868  /// Theta (pivot change)
869  inline double theta() const
870  { return theta_;};
871  /// Return pointer to details of costs
872  inline ClpNonLinearCost * nonLinearCost() const
873  { return nonLinearCost_;};
874  //@}
875  /**@name status methods */
876  //@{
877  inline void setFakeBound(int sequence, FakeBound fakeBound)
878  {
879    unsigned char & st_byte = status_[sequence];
880    st_byte &= ~24;
881    st_byte |= fakeBound<<3;
882  };
883  inline FakeBound getFakeBound(int sequence) const
884  {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
885  inline void setRowStatus(int sequence, Status status)
886  {
887    unsigned char & st_byte = status_[sequence+numberColumns_];
888    st_byte &= ~7;
889    st_byte |= status;
890  };
891  inline Status getRowStatus(int sequence) const
892  {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
893  inline void setColumnStatus(int sequence, Status status)
894  {
895    unsigned char & st_byte = status_[sequence];
896    st_byte &= ~7;
897    st_byte |= status;
898  };
899  inline Status getColumnStatus(int sequence) const
900  {return static_cast<Status> (status_[sequence]&7);};
901  inline void setPivoted( int sequence)
902  { status_[sequence] |= 32;};
903  inline void clearPivoted( int sequence)
904  { status_[sequence] &= ~32; };
905  inline bool pivoted(int sequence) const
906  {return (((status_[sequence]>>5)&1)!=0);};
907  /// To flag a variable (not inline to allow for column generation)
908  void setFlagged( int sequence);
909  inline void clearFlagged( int sequence)
910  {
911    status_[sequence] &= ~64;
912  };
913  inline bool flagged(int sequence) const
914  {return ((status_[sequence]&64)!=0);};
915  /// To say row active in primal pivot row choice
916  inline void setActive( int iRow)
917  {
918    status_[iRow] |= 128;
919  };
920  inline void clearActive( int iRow)
921  {
922    status_[iRow] &= ~128;
923  };
924  inline bool active(int iRow) const
925  {return ((status_[iRow]&128)!=0);};
926  /** Set up status array (can be used by OsiClp).
927      Also can be used to set up all slack basis */
928  void createStatus() ;
929  /** Sets up all slack basis and resets solution to
930      as it was after initial load or readMps */
931  void allSlackBasis(bool resetSolution=false);
932   
933  /// So we know when to be cautious
934  inline int lastBadIteration() const
935  {return lastBadIteration_;};
936  /// Progress flag - at present 0 bit says artificials out
937  inline int progressFlag() const
938  {return progressFlag_;};
939  /// Force re-factorization early
940  inline void forceFactorization(int value)
941  { forceFactorization_ = value;};
942  /// Raw objective value (so always minimize in primal)
943  inline double rawObjectiveValue() const
944  { return objectiveValue_;};
945   /// Compute objective value from solution and put in objectiveValue_
946  void computeObjectiveValue(bool useWorkingSolution=false);
947  /// Compute minimization objective value from internal solution without perturbation
948  double computeInternalObjectiveValue();
949  /** Number of extra rows.  These are ones which will be dynamically created
950      each iteration.  This is for GUB but may have other uses.
951  */
952  inline int numberExtraRows() const
953  { return numberExtraRows_;};
954  /** Maximum number of basic variables - can be more than number of rows if GUB
955  */
956  inline int maximumBasic() const
957  { return maximumBasic_;};
958  /// Create C++ lines to get to current state
959  void generateCpp( FILE * fp,bool defaultFactor=false);
960  /** For advanced options
961      1 - Don't keep changing infeasibility weight
962      2 - Keep nonLinearCost round solves
963      4 - Force outgoing variables to exact bound (primal)
964      8 - Safe to use dense initial factorization
965      16 -Just use basic variables for operation if column generation
966      32 -Clean up with primal before strong branching
967      64 -Treat problem as feasible until last minute (i.e. minimize infeasibilities)
968      128 - Switch off all matrix sanity checks
969      256 - No row copy
970      512 - If not in values pass, solution guaranteed, skip as much as possible
971      1024 - In branch and bound
972      2048 - Don't bother to re-factorize if < 20 iterations
973      4096 - Skip some optimality checks
974      8192 - Do Primal when cleaning up primal
975      16384 - In fast dual (so we can switch off things)
976      32678 - called from Osi
977      65356 - keep arrays around as much as possible
978      NOTE - many applications can call Clp but there may be some short cuts
979             which are taken which are not guaranteed safe from all applications.
980             Vetted applications will have a bit set and the code may test this
981             At present I expect a few such applications - if too many I will
982             have to re-think.  It is up to application owner to change the code
983             if she/he needs these short cuts.  I will not debug unless in Coin
984             repository.  See COIN_CLP_VETTED comments.
985      0x01000000 is Cbc (and in branch and bound)
986      0x02000000 is in a different branch and bound
987  */
988#define COIN_CBC_USING_CLP 0x01000000
989  inline unsigned int specialOptions() const
990  { return specialOptions_;};
991  void setSpecialOptions(unsigned int value);
992  /// Gets clean and emptyish factorization
993  ClpFactorization * getEmptyFactorization();
994  /// May delete or may make clean and emptyish factorization
995  void setEmptyFactorization();
996  //@}
997
998  ///@name Basis handling
999  // These are only to be used using startFinishOptions (ClpSimplexDual, ClpSimplexPrimal)
1000  // *** At present only without scaling
1001  // *** Slacks havve -1.0 element (so == row activity) - take care
1002  ///Get a row of the tableau (slack part in slack if not NULL)
1003  void getBInvARow(int row, double* z, double * slack=NULL);
1004 
1005  ///Get a row of the basis inverse
1006  void getBInvRow(int row, double* z);
1007 
1008  ///Get a column of the tableau
1009  void getBInvACol(int col, double* vec);
1010 
1011  ///Get a column of the basis inverse
1012  void getBInvCol(int col, double* vec);
1013 
1014  /** Get basic indices (order of indices corresponds to the
1015      order of elements in a vector retured by getBInvACol() and
1016      getBInvCol()).
1017  */
1018  void getBasics(int* index);
1019 
1020  //@}
1021    //-------------------------------------------------------------------------
1022    /**@name Changing bounds on variables and constraints */
1023    //@{
1024       /** Set an objective function coefficient */
1025       void setObjectiveCoefficient( int elementIndex, double elementValue );
1026       /** Set an objective function coefficient */
1027       inline void setObjCoeff( int elementIndex, double elementValue )
1028       { setObjectiveCoefficient( elementIndex, elementValue);};
1029
1030      /** Set a single column lower bound<br>
1031          Use -DBL_MAX for -infinity. */
1032       void setColumnLower( int elementIndex, double elementValue );
1033     
1034      /** Set a single column upper bound<br>
1035          Use DBL_MAX for infinity. */
1036       void setColumnUpper( int elementIndex, double elementValue );
1037
1038      /** Set a single column lower and upper bound */
1039      void setColumnBounds( int elementIndex,
1040        double lower, double upper );
1041
1042      /** Set the bounds on a number of columns simultaneously<br>
1043          The default implementation just invokes setColLower() and
1044          setColUpper() over and over again.
1045          @param indexFirst,indexLast pointers to the beginning and after the
1046                 end of the array of the indices of the variables whose
1047                 <em>either</em> bound changes
1048          @param boundList the new lower/upper bound pairs for the variables
1049      */
1050      void setColumnSetBounds(const int* indexFirst,
1051                                   const int* indexLast,
1052                                   const double* boundList);
1053     
1054      /** Set a single column lower bound<br>
1055          Use -DBL_MAX for -infinity. */
1056       inline void setColLower( int elementIndex, double elementValue )
1057       { setColumnLower(elementIndex, elementValue);};
1058      /** Set a single column upper bound<br>
1059          Use DBL_MAX for infinity. */
1060       inline void setColUpper( int elementIndex, double elementValue )
1061       { setColumnUpper(elementIndex, elementValue);};
1062
1063      /** Set a single column lower and upper bound */
1064      inline void setColBounds( int elementIndex,
1065        double lower, double upper )
1066       { setColumnBounds(elementIndex, lower, upper);};
1067
1068      /** Set the bounds on a number of columns simultaneously<br>
1069          @param indexFirst,indexLast pointers to the beginning and after the
1070                 end of the array of the indices of the variables whose
1071                 <em>either</em> bound changes
1072          @param boundList the new lower/upper bound pairs for the variables
1073      */
1074      inline void setColSetBounds(const int* indexFirst,
1075                                   const int* indexLast,
1076                                   const double* boundList)
1077      { setColumnSetBounds(indexFirst, indexLast, boundList);};
1078     
1079      /** Set a single row lower bound<br>
1080          Use -DBL_MAX for -infinity. */
1081      void setRowLower( int elementIndex, double elementValue );
1082     
1083      /** Set a single row upper bound<br>
1084          Use DBL_MAX for infinity. */
1085      void setRowUpper( int elementIndex, double elementValue ) ;
1086   
1087      /** Set a single row lower and upper bound */
1088      void setRowBounds( int elementIndex,
1089                                 double lower, double upper ) ;
1090   
1091      /** Set the bounds on a number of rows simultaneously<br>
1092          @param indexFirst,indexLast pointers to the beginning and after the
1093                 end of the array of the indices of the constraints whose
1094                 <em>either</em> bound changes
1095          @param boundList the new lower/upper bound pairs for the constraints
1096      */
1097      void setRowSetBounds(const int* indexFirst,
1098                                   const int* indexLast,
1099                                   const double* boundList);
1100   
1101    //@}
1102
1103////////////////// data //////////////////
1104protected:
1105
1106  /**@name data.  Many arrays have a row part and a column part.
1107   There is a single array with both - columns then rows and
1108   then normally two arrays pointing to rows and columns.  The
1109   single array is the owner of memory
1110  */
1111  //@{
1112  /// Worst column primal infeasibility
1113  double columnPrimalInfeasibility_;
1114  /// Worst row primal infeasibility
1115  double rowPrimalInfeasibility_;
1116  /// Sequence of worst (-1 if feasible)
1117  int columnPrimalSequence_;
1118  /// Sequence of worst (-1 if feasible)
1119  int rowPrimalSequence_;
1120  /// Worst column dual infeasibility
1121  double columnDualInfeasibility_;
1122  /// Worst row dual infeasibility
1123  double rowDualInfeasibility_;
1124  /// Sequence of worst (-1 if feasible)
1125  int columnDualSequence_;
1126  /// Sequence of worst (-1 if feasible)
1127  int rowDualSequence_;
1128  /// Primal tolerance needed to make dual feasible (<largeTolerance)
1129  double primalToleranceToGetOptimal_;
1130  /// Remaining largest dual infeasibility
1131  double remainingDualInfeasibility_;
1132  /// Large bound value (for complementarity etc)
1133  double largeValue_;
1134  /// Largest error on Ax-b
1135  double largestPrimalError_;
1136  /// Largest error on basic duals
1137  double largestDualError_;
1138  /// Largest difference between input primal solution and computed
1139  double largestSolutionError_;
1140  /// Dual bound
1141  double dualBound_;
1142  /// Alpha (pivot element)
1143  double alpha_;
1144  /// Theta (pivot change)
1145  double theta_;
1146  /// Lower Bound on In variable
1147  double lowerIn_;
1148  /// Value of In variable
1149  double valueIn_;
1150  /// Upper Bound on In variable
1151  double upperIn_;
1152  /// Reduced cost of In variable
1153  double dualIn_;
1154  /// Lower Bound on Out variable
1155  double lowerOut_;
1156  /// Value of Out variable
1157  double valueOut_;
1158  /// Upper Bound on Out variable
1159  double upperOut_;
1160  /// Infeasibility (dual) or ? (primal) of Out variable
1161  double dualOut_;
1162  /// Current dual tolerance for algorithm
1163  double dualTolerance_;
1164  /// Current primal tolerance for algorithm
1165  double primalTolerance_;
1166  /// Sum of dual infeasibilities
1167  double sumDualInfeasibilities_;
1168  /// Sum of primal infeasibilities
1169  double sumPrimalInfeasibilities_;
1170  /// Weight assigned to being infeasible in primal
1171  double infeasibilityCost_;
1172  /// Sum of Dual infeasibilities using tolerance based on error in duals
1173  double sumOfRelaxedDualInfeasibilities_;
1174  /// Sum of Primal infeasibilities using tolerance based on error in primals
1175  double sumOfRelaxedPrimalInfeasibilities_;
1176  /// Acceptable pivot value just after factorization
1177  double acceptablePivot_;
1178  /// Working copy of lower bounds (Owner of arrays below)
1179  double * lower_;
1180  /// Row lower bounds - working copy
1181  double * rowLowerWork_;
1182  /// Column lower bounds - working copy
1183  double * columnLowerWork_;
1184  /// Working copy of upper bounds (Owner of arrays below)
1185  double * upper_;
1186  /// Row upper bounds - working copy
1187  double * rowUpperWork_;
1188  /// Column upper bounds - working copy
1189  double * columnUpperWork_;
1190  /// Working copy of objective (Owner of arrays below)
1191  double * cost_;
1192  /// Row objective - working copy
1193  double * rowObjectiveWork_;
1194  /// Column objective - working copy
1195  double * objectiveWork_;
1196  /// Useful row length arrays
1197  CoinIndexedVector * rowArray_[6];
1198  /// Useful column length arrays
1199  CoinIndexedVector * columnArray_[6];
1200  /// Sequence of In variable
1201  int sequenceIn_;
1202  /// Direction of In, 1 going up, -1 going down, 0 not a clude
1203  int directionIn_;
1204  /// Sequence of Out variable
1205  int sequenceOut_;
1206  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
1207  int directionOut_;
1208  /// Pivot Row
1209  int pivotRow_;
1210  /// Last good iteration (immediately after a re-factorization)
1211  int lastGoodIteration_;
1212  /// Working copy of reduced costs (Owner of arrays below)
1213  double * dj_;
1214  /// Reduced costs of slacks not same as duals (or - duals)
1215  double * rowReducedCost_;
1216  /// Possible scaled reduced costs
1217  double * reducedCostWork_;
1218  /// Working copy of primal solution (Owner of arrays below)
1219  double * solution_;
1220  /// Row activities - working copy
1221  double * rowActivityWork_;
1222  /// Column activities - working copy
1223  double * columnActivityWork_;
1224  /// Auxiliary model
1225  ClpSimplex * auxiliaryModel_;
1226  /// Number of dual infeasibilities
1227  int numberDualInfeasibilities_;
1228  /// Number of dual infeasibilities (without free)
1229  int numberDualInfeasibilitiesWithoutFree_;
1230  /// Number of primal infeasibilities
1231  int numberPrimalInfeasibilities_;
1232  /// How many iterative refinements to do
1233  int numberRefinements_;
1234  /// dual row pivot choice
1235  ClpDualRowPivot * dualRowPivot_;
1236  /// primal column pivot choice
1237  ClpPrimalColumnPivot * primalColumnPivot_;
1238  /// Basic variables pivoting on which rows
1239  int * pivotVariable_;
1240  /// factorization
1241  ClpFactorization * factorization_;
1242  /// Saved version of solution
1243  double * savedSolution_;
1244  /// Number of times code has tentatively thought optimal
1245  int numberTimesOptimal_;
1246  /// Disaster handler
1247  ClpDisasterHandler * disasterArea_;
1248  /// If change has been made (first attempt at stopping looping)
1249  int changeMade_;
1250  /// Algorithm >0 == Primal, <0 == Dual
1251  int algorithm_;
1252  /** Now for some reliability aids
1253      This forces re-factorization early */
1254  int forceFactorization_;
1255  /** Perturbation:
1256      -50 to +50 - perturb by this power of ten (-6 sounds good)
1257      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
1258      101 - we are perturbed
1259      102 - don't try perturbing again
1260      default is 100
1261  */
1262  int perturbation_;
1263  /// Saved status regions
1264  unsigned char * saveStatus_;
1265  /** Very wasteful way of dealing with infeasibilities in primal.
1266      However it will allow non-linearities and use of dual
1267      analysis.  If it doesn't work it can easily be replaced.
1268  */
1269  ClpNonLinearCost * nonLinearCost_;
1270  /** For advanced options
1271      See get and set for meaning
1272  */
1273  unsigned int specialOptions_;
1274  /// So we know when to be cautious
1275  int lastBadIteration_;
1276  /// So we know when to open up again
1277  int lastFlaggedIteration_;
1278  /// Can be used for count of fake bounds (dual) or fake costs (primal)
1279  int numberFake_;
1280  /// Can be used for count of changed costs (dual) or changed bounds (primal)
1281  int numberChanged_;
1282  /// Progress flag - at present 0 bit says artificials out, 1 free in
1283  int progressFlag_;
1284  /// First free/super-basic variable (-1 if none)
1285  int firstFree_;
1286  /** Number of extra rows.  These are ones which will be dynamically created
1287      each iteration.  This is for GUB but may have other uses.
1288  */
1289  int numberExtraRows_;
1290  /** Maximum number of basic variables - can be more than number of rows if GUB
1291  */
1292  int maximumBasic_;
1293  /** For advanced use.  When doing iterative solves things can get
1294      nasty so on values pass if incoming solution has largest
1295      infeasibility < incomingInfeasibility throw out variables
1296      from basis until largest infeasibility < allowedInfeasibility.
1297      if allowedInfeasibility>= incomingInfeasibility this is
1298      always possible altough you may end up with an all slack basis.
1299
1300      Defaults are 1.0,10.0
1301  */
1302  float incomingInfeasibility_;
1303  float allowedInfeasibility_;
1304  /// Automatic scaling of objective and rhs and bounds
1305  int automaticScale_;
1306  /// For dealing with all issues of cycling etc
1307  ClpSimplexProgress * progress_;
1308public:
1309  /// Spare int array for passing information [0]!=0 switches on
1310  mutable int spareIntArray_[4];
1311  /// Spare double array for passing information [0]!=0 switches on
1312  mutable double spareDoubleArray_[4];
1313protected:
1314  /// Allow OsiClp certain perks
1315  friend class OsiClpSolverInterface;
1316  //@}
1317};
1318//#############################################################################
1319/** A function that tests the methods in the ClpSimplex class. The
1320    only reason for it not to be a member method is that this way it doesn't
1321    have to be compiled into the library. And that's a gain, because the
1322    library should be compiled with optimization on, but this method should be
1323    compiled with debugging.
1324
1325    It also does some testing of ClpFactorization class
1326 */
1327void
1328ClpSimplexUnitTest(const std::string & mpsDir,
1329                   const std::string & netlibDir);
1330
1331
1332/// For saving extra information to see if looping.
1333class ClpSimplexProgress {
1334
1335public:
1336
1337
1338  /**@name Constructors and destructor and copy */
1339  //@{
1340  /// Default constructor
1341    ClpSimplexProgress (  );
1342
1343  /// Constructor from model
1344    ClpSimplexProgress ( ClpSimplex * model );
1345
1346  /// Copy constructor.
1347  ClpSimplexProgress(const ClpSimplexProgress &);
1348
1349  /// Assignment operator. This copies the data
1350    ClpSimplexProgress & operator=(const ClpSimplexProgress & rhs);
1351  /// Destructor
1352   ~ClpSimplexProgress (  );
1353  //@}
1354
1355  /**@name Check progress */
1356  //@{
1357  /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
1358      >=0 if give up and use as problem status
1359  */
1360    int looping (  );
1361  /// Start check at beginning of whileIterating
1362  void startCheck();
1363  /// Returns cycle length in whileIterating
1364  int cycle(int in, int out,int wayIn,int wayOut); 
1365
1366  /// Returns previous objective (if -1) - current if (0)
1367  double lastObjective(int back=1) const;
1368  /// Set real primal infeasibility and move back
1369  void setInfeasibility(double value);
1370  /// Returns real primal infeasibility (if -1) - current if (0)
1371  double lastInfeasibility(int back=1) const;
1372  /// Modify objective e.g. if dual infeasible in dual
1373  void modifyObjective(double value);
1374  /// Returns previous iteration number (if -1) - current if (0)
1375  int lastIterationNumber(int back=1) const;
1376  /// clears all iteration numbers (to switch off panic)
1377  void clearIterationNumbers();
1378  /// Odd state
1379  inline void newOddState()
1380  { oddState_= - oddState_-1;};
1381  inline void endOddState()
1382  { oddState_=abs(oddState_);};
1383  inline void clearOddState() 
1384  { oddState_=0;};
1385  inline int oddState() const
1386  { return oddState_;};
1387  /// number of bad times
1388  inline int badTimes() const
1389  { return numberBadTimes_;};
1390  inline void clearBadTimes()
1391  { numberBadTimes_=0;};
1392
1393  //@}
1394  /**@name Data  */
1395#define CLP_PROGRESS 5
1396  //#define CLP_PROGRESS_WEIGHT 10
1397  //@{
1398  /// Objective values
1399  double objective_[CLP_PROGRESS];
1400  /// Sum of infeasibilities for algorithm
1401  double infeasibility_[CLP_PROGRESS];
1402  /// Sum of real primal infeasibilities for primal
1403  double realInfeasibility_[CLP_PROGRESS];
1404#ifdef CLP_PROGRESS_WEIGHT
1405  /// Objective values for weights
1406  double objectiveWeight_[CLP_PROGRESS_WEIGHT];
1407  /// Sum of infeasibilities for algorithm for weights
1408  double infeasibilityWeight_[CLP_PROGRESS_WEIGHT];
1409  /// Sum of real primal infeasibilities for primal for weights
1410  double realInfeasibilityWeight_[CLP_PROGRESS_WEIGHT];
1411  /// Drop  for weights
1412  double drop_;
1413  /// Best? for weights
1414  double best_;
1415#endif
1416  /// Initial weight for weights
1417  double initialWeight_;
1418#define CLP_CYCLE 12
1419  /// For cycle checking
1420  //double obj_[CLP_CYCLE];
1421  int in_[CLP_CYCLE];
1422  int out_[CLP_CYCLE];
1423  char way_[CLP_CYCLE];
1424  /// Pointer back to model so we can get information
1425  ClpSimplex * model_;
1426  /// Number of infeasibilities
1427  int numberInfeasibilities_[CLP_PROGRESS];
1428  /// Iteration number at which occurred
1429  int iterationNumber_[CLP_PROGRESS];
1430#ifdef CLP_PROGRESS_WEIGHT
1431  /// Number of infeasibilities for weights
1432  int numberInfeasibilitiesWeight_[CLP_PROGRESS_WEIGHT];
1433  /// Iteration number at which occurred for weights
1434  int iterationNumberWeight_[CLP_PROGRESS_WEIGHT];
1435#endif
1436  /// Number of times checked (so won't stop too early)
1437  int numberTimes_;
1438  /// Number of times it looked like loop
1439  int numberBadTimes_;
1440  /// If things are in an odd state
1441  int oddState_;
1442  //@}
1443};
1444// For Devex stuff
1445#define DEVEX_TRY_NORM 1.0e-4
1446#define DEVEX_ADD_ONE 1.0
1447#endif
Note: See TracBrowser for help on using the repository browser.