source: stable/1.15/Clp/src/ClpSimplex.hpp @ 1995

Last change on this file since 1995 was 1989, checked in by forrest, 6 years ago

Minor stability fixes and an option to allow perturbation after presolve

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 67.1 KB
Line 
1/* $Id: ClpSimplex.hpp 1989 2013-11-11 17:27:32Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5/*
6   Authors
7
8   John Forrest
9
10 */
11#ifndef ClpSimplex_H
12#define ClpSimplex_H
13
14#include <iostream>
15#include <cfloat>
16#include "ClpModel.hpp"
17#include "ClpMatrixBase.hpp"
18#include "ClpSolve.hpp"
19#include "ClpConfig.h"
20class ClpDualRowPivot;
21class ClpPrimalColumnPivot;
22class ClpFactorization;
23class CoinIndexedVector;
24class ClpNonLinearCost;
25class ClpNodeStuff;
26class CoinStructuredModel;
27class OsiClpSolverInterface;
28class CoinWarmStartBasis;
29class ClpDisasterHandler;
30class ClpConstraint;
31#ifdef CLP_HAS_ABC
32#include "AbcCommon.hpp"
33class AbcTolerancesEtc;
34class AbcSimplex;
35#include "CoinAbcCommon.hpp"
36#endif
37/** This solves LPs using the simplex method
38
39    It inherits from ClpModel and all its arrays are created at
40    algorithm time. Originally I tried to work with model arrays
41    but for simplicity of coding I changed to single arrays with
42    structural variables then row variables.  Some coding is still
43    based on old style and needs cleaning up.
44
45    For a description of algorithms:
46
47    for dual see ClpSimplexDual.hpp and at top of ClpSimplexDual.cpp
48    for primal see ClpSimplexPrimal.hpp and at top of ClpSimplexPrimal.cpp
49
50    There is an algorithm data member.  + for primal variations
51    and - for dual variations
52
53*/
54
55class ClpSimplex : public ClpModel {
56     friend void ClpSimplexUnitTest(const std::string & mpsDir);
57
58public:
59     /** enums for status of various sorts.
60         First 4 match CoinWarmStartBasis,
61         isFixed means fixed at lower bound and out of basis
62     */
63     enum Status {
64          isFree = 0x00,
65          basic = 0x01,
66          atUpperBound = 0x02,
67          atLowerBound = 0x03,
68          superBasic = 0x04,
69          isFixed = 0x05
70     };
71     // For Dual
72     enum FakeBound {
73          noFake = 0x00,
74          lowerFake = 0x01,
75          upperFake = 0x02,
76          bothFake = 0x03
77     };
78
79     /**@name Constructors and destructor and copy */
80     //@{
81     /// Default constructor
82     ClpSimplex (bool emptyMessages = false  );
83
84     /** Copy constructor. May scale depending on mode
85         -1 leave mode as is
86         0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
87     */
88     ClpSimplex(const ClpSimplex & rhs, int scalingMode = -1);
89     /** Copy constructor from model. May scale depending on mode
90         -1 leave mode as is
91         0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
92     */
93     ClpSimplex(const ClpModel & rhs, int scalingMode = -1);
94     /** Subproblem constructor.  A subset of whole model is created from the
95         row and column lists given.  The new order is given by list order and
96         duplicates are allowed.  Name and integer information can be dropped
97         Can optionally modify rhs to take into account variables NOT in list
98         in this case duplicates are not allowed (also see getbackSolution)
99     */
100     ClpSimplex (const ClpModel * wholeModel,
101                 int numberRows, const int * whichRows,
102                 int numberColumns, const int * whichColumns,
103                 bool dropNames = true, bool dropIntegers = true,
104                 bool fixOthers = false);
105     /** Subproblem constructor.  A subset of whole model is created from the
106         row and column lists given.  The new order is given by list order and
107         duplicates are allowed.  Name and integer information can be dropped
108         Can optionally modify rhs to take into account variables NOT in list
109         in this case duplicates are not allowed (also see getbackSolution)
110     */
111     ClpSimplex (const ClpSimplex * wholeModel,
112                 int numberRows, const int * whichRows,
113                 int numberColumns, const int * whichColumns,
114                 bool dropNames = true, bool dropIntegers = true,
115                 bool fixOthers = false);
116     /** This constructor modifies original ClpSimplex and stores
117         original stuff in created ClpSimplex.  It is only to be used in
118         conjunction with originalModel */
119     ClpSimplex (ClpSimplex * wholeModel,
120                 int numberColumns, const int * whichColumns);
121     /** This copies back stuff from miniModel and then deletes miniModel.
122         Only to be used with mini constructor */
123     void originalModel(ClpSimplex * miniModel);
124  inline int abcState() const
125  { return abcState_;}
126  inline void setAbcState(int state)
127  { abcState_=state;}
128#ifdef ABC_INHERIT
129  inline AbcSimplex * abcSimplex() const
130  { return abcSimplex_;}
131  inline void setAbcSimplex(AbcSimplex * simplex)
132  { abcSimplex_=simplex;}
133  /// Returns 0 if dual can be skipped
134  int doAbcDual();
135  /// Returns 0 if primal can be skipped
136  int doAbcPrimal(int ifValuesPass);
137#endif
138     /** Array persistence flag
139         If 0 then as now (delete/new)
140         1 then only do arrays if bigger needed
141         2 as 1 but give a bit extra if bigger needed
142     */
143     void setPersistenceFlag(int value);
144     /// Save a copy of model with certain state - normally without cuts
145     void makeBaseModel();
146     /// Switch off base model
147     void deleteBaseModel();
148     /// See if we have base model
149     inline ClpSimplex *  baseModel() const {
150          return baseModel_;
151     }
152     /** Reset to base model (just size and arrays needed)
153         If model NULL use internal copy
154     */
155     void setToBaseModel(ClpSimplex * model = NULL);
156     /// Assignment operator. This copies the data
157     ClpSimplex & operator=(const ClpSimplex & rhs);
158     /// Destructor
159     ~ClpSimplex (  );
160     // Ones below are just ClpModel with some changes
161     /** Loads a problem (the constraints on the
162           rows are given by lower and upper bounds). If a pointer is 0 then the
163           following values are the default:
164           <ul>
165             <li> <code>colub</code>: all columns have upper bound infinity
166             <li> <code>collb</code>: all columns have lower bound 0
167             <li> <code>rowub</code>: all rows have upper bound infinity
168             <li> <code>rowlb</code>: all rows have lower bound -infinity
169         <li> <code>obj</code>: all variables have 0 objective coefficient
170           </ul>
171       */
172     void loadProblem (  const ClpMatrixBase& matrix,
173                         const double* collb, const double* colub,
174                         const double* obj,
175                         const double* rowlb, const double* rowub,
176                         const double * rowObjective = NULL);
177     void loadProblem (  const CoinPackedMatrix& matrix,
178                         const double* collb, const double* colub,
179                         const double* obj,
180                         const double* rowlb, const double* rowub,
181                         const double * rowObjective = NULL);
182
183     /** Just like the other loadProblem() method except that the matrix is
184       given in a standard column major ordered format (without gaps). */
185     void loadProblem (  const int numcols, const int numrows,
186                         const CoinBigIndex* start, const int* index,
187                         const double* value,
188                         const double* collb, const double* colub,
189                         const double* obj,
190                         const double* rowlb, const double* rowub,
191                         const double * rowObjective = NULL);
192     /// This one is for after presolve to save memory
193     void loadProblem (  const int numcols, const int numrows,
194                         const CoinBigIndex* start, const int* index,
195                         const double* value, const int * length,
196                         const double* collb, const double* colub,
197                         const double* obj,
198                         const double* rowlb, const double* rowub,
199                         const double * rowObjective = NULL);
200     /** This loads a model from a coinModel object - returns number of errors.
201         If keepSolution true and size is same as current then
202         keeps current status and solution
203     */
204     int loadProblem (  CoinModel & modelObject, bool keepSolution = false);
205     /// Read an mps file from the given filename
206     int readMps(const char *filename,
207                 bool keepNames = false,
208                 bool ignoreErrors = false);
209     /// Read GMPL files from the given filenames
210     int readGMPL(const char *filename, const char * dataName,
211                  bool keepNames = false);
212     /// Read file in LP format from file with name filename.
213     /// See class CoinLpIO for description of this format.
214     int readLp(const char *filename, const double epsilon = 1e-5);
215     /** Borrow model.  This is so we dont have to copy large amounts
216         of data around.  It assumes a derived class wants to overwrite
217         an empty model with a real one - while it does an algorithm.
218         This is same as ClpModel one, but sets scaling on etc. */
219     void borrowModel(ClpModel & otherModel);
220     void borrowModel(ClpSimplex & otherModel);
221     /// Pass in Event handler (cloned and deleted at end)
222     void passInEventHandler(const ClpEventHandler * eventHandler);
223     /// Puts solution back into small model
224     void getbackSolution(const ClpSimplex & smallModel, const int * whichRow, const int * whichColumn);
225     /** Load nonlinear part of problem from AMPL info
226         Returns 0 if linear
227         1 if quadratic objective
228         2 if quadratic constraints
229         3 if nonlinear objective
230         4 if nonlinear constraints
231         -1 on failure
232     */
233     int loadNonLinear(void * info, int & numberConstraints,
234                       ClpConstraint ** & constraints);
235#ifdef ABC_INHERIT
236  /// Loads tolerances etc
237  void loadTolerancesEtc(const AbcTolerancesEtc & data);
238  /// Unloads tolerances etc
239  void unloadTolerancesEtc(AbcTolerancesEtc & data);
240#endif
241     //@}
242
243     /**@name Functions most useful to user */
244     //@{
245     /** General solve algorithm which can do presolve.
246         See  ClpSolve.hpp for options
247      */
248     int initialSolve(ClpSolve & options);
249     /// Default initial solve
250     int initialSolve();
251     /// Dual initial solve
252     int initialDualSolve();
253     /// Primal initial solve
254     int initialPrimalSolve();
255/// Barrier initial solve
256     int initialBarrierSolve();
257     /// Barrier initial solve, not to be followed by crossover
258     int initialBarrierNoCrossSolve();
259     /** Dual algorithm - see ClpSimplexDual.hpp for method.
260         ifValuesPass==2 just does values pass and then stops.
261
262         startFinishOptions - bits
263         1 - do not delete work areas and factorization at end
264         2 - use old factorization if same number of rows
265         4 - skip as much initialization of work areas as possible
266             (based on whatsChanged in clpmodel.hpp) ** work in progress
267         maybe other bits later
268     */
269     int dual(int ifValuesPass = 0, int startFinishOptions = 0);
270     // If using Debug
271     int dualDebug(int ifValuesPass = 0, int startFinishOptions = 0);
272     /** Primal algorithm - see ClpSimplexPrimal.hpp for method.
273         ifValuesPass==2 just does values pass and then stops.
274
275         startFinishOptions - bits
276         1 - do not delete work areas and factorization at end
277         2 - use old factorization if same number of rows
278         4 - skip as much initialization of work areas as possible
279             (based on whatsChanged in clpmodel.hpp) ** work in progress
280         maybe other bits later
281     */
282     int primal(int ifValuesPass = 0, int startFinishOptions = 0);
283     /** Solves nonlinear problem using SLP - may be used as crash
284         for other algorithms when number of iterations small.
285         Also exits if all problematical variables are changing
286         less than deltaTolerance
287     */
288     int nonlinearSLP(int numberPasses, double deltaTolerance);
289     /** Solves problem with nonlinear constraints using SLP - may be used as crash
290         for other algorithms when number of iterations small.
291         Also exits if all problematical variables are changing
292         less than deltaTolerance
293     */
294     int nonlinearSLP(int numberConstraints, ClpConstraint ** constraints,
295                      int numberPasses, double deltaTolerance);
296     /** Solves using barrier (assumes you have good cholesky factor code).
297         Does crossover to simplex if asked*/
298     int barrier(bool crossover = true);
299     /** Solves non-linear using reduced gradient.  Phase = 0 get feasible,
300         =1 use solution */
301     int reducedGradient(int phase = 0);
302     /// Solve using structure of model and maybe in parallel
303     int solve(CoinStructuredModel * model);
304#ifdef ABC_INHERIT
305  /** solvetype 0 for dual, 1 for primal
306      startup 1 for values pass
307      interrupt whether to pass across interrupt handler
308  */
309  void dealWithAbc(int solveType,int startUp,bool interrupt=false);
310#endif
311     /** This loads a model from a CoinStructuredModel object - returns number of errors.
312         If originalOrder then keep to order stored in blocks,
313         otherwise first column/rows correspond to first block - etc.
314         If keepSolution true and size is same as current then
315         keeps current status and solution
316     */
317     int loadProblem (  CoinStructuredModel & modelObject,
318                        bool originalOrder = true, bool keepSolution = false);
319     /**
320        When scaling is on it is possible that the scaled problem
321        is feasible but the unscaled is not.  Clp returns a secondary
322        status code to that effect.  This option allows for a cleanup.
323        If you use it I would suggest 1.
324        This only affects actions when scaled optimal
325        0 - no action
326        1 - clean up using dual if primal infeasibility
327        2 - clean up using dual if dual infeasibility
328        3 - clean up using dual if primal or dual infeasibility
329        11,12,13 - as 1,2,3 but use primal
330
331        return code as dual/primal
332     */
333     int cleanup(int cleanupScaling);
334     /** Dual ranging.
335         This computes increase/decrease in cost for each given variable and corresponding
336         sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
337         and numberColumns.. for artificials/slacks.
338         For non-basic variables the information is trivial to compute and the change in cost is just minus the
339         reduced cost and the sequence number will be that of the non-basic variables.
340         For basic variables a ratio test is between the reduced costs for non-basic variables
341         and the row of the tableau corresponding to the basic variable.
342         The increase/decrease value is always >= 0.0
343
344         Up to user to provide correct length arrays where each array is of length numberCheck.
345         which contains list of variables for which information is desired.  All other
346         arrays will be filled in by function.  If fifth entry in which is variable 7 then fifth entry in output arrays
347         will be information for variable 7.
348
349         If valueIncrease/Decrease not NULL (both must be NULL or both non NULL) then these are filled with
350         the value of variable if such a change in cost were made (the existing bounds are ignored)
351
352         Returns non-zero if infeasible unbounded etc
353     */
354     int dualRanging(int numberCheck, const int * which,
355                     double * costIncrease, int * sequenceIncrease,
356                     double * costDecrease, int * sequenceDecrease,
357                     double * valueIncrease = NULL, double * valueDecrease = NULL);
358     /** Primal ranging.
359         This computes increase/decrease in value for each given variable and corresponding
360         sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
361         and numberColumns.. for artificials/slacks.
362         This should only be used for non-basic variabls as otherwise information is pretty useless
363         For basic variables the sequence number will be that of the basic variables.
364
365         Up to user to provide correct length arrays where each array is of length numberCheck.
366         which contains list of variables for which information is desired.  All other
367         arrays will be filled in by function.  If fifth entry in which is variable 7 then fifth entry in output arrays
368         will be information for variable 7.
369
370         Returns non-zero if infeasible unbounded etc
371     */
372     int primalRanging(int numberCheck, const int * which,
373                       double * valueIncrease, int * sequenceIncrease,
374                       double * valueDecrease, int * sequenceDecrease);
375     /**
376        Modifies coefficients etc and if necessary pivots in and out.
377        All at same status will be done (basis may go singular).
378        User can tell which others have been done (i.e. if status matches).
379        If called from outside will change status and return 0.
380        If called from event handler returns non-zero if user has to take action.
381        indices>=numberColumns are slacks (obviously no coefficients)
382        status array is (char) Status enum
383     */
384     int modifyCoefficientsAndPivot(int number,
385                                 const int * which,
386                                 const CoinBigIndex * start,
387                                 const int * row,
388                                 const double * newCoefficient,
389                                 const unsigned char * newStatus=NULL,
390                                 const double * newLower=NULL,
391                                 const double * newUpper=NULL,
392                                 const double * newObjective=NULL);
393     /** Take out duplicate rows (includes scaled rows and intersections).
394         On exit whichRows has rows to delete - return code is number can be deleted
395         or -1 if would be infeasible.
396         If tolerance is -1.0 use primalTolerance for equality rows and infeasibility
397         If cleanUp not zero then spend more time trying to leave more stable row
398         and make row bounds exact multiple of cleanUp if close enough
399     */
400     int outDuplicateRows(int numberLook,int * whichRows, bool noOverlaps=false, double tolerance=-1.0,
401                          double cleanUp=0.0);
402     /** Try simple crash like techniques to get closer to primal feasibility
403         returns final sum of infeasibilities */
404     double moveTowardsPrimalFeasible();
405     /** Try simple crash like techniques to remove super basic slacks
406         but only if > threshold */
407     void removeSuperBasicSlacks(int threshold=0);
408     /** Mini presolve (faster)
409         Char arrays must be numberRows and numberColumns long
410         on entry second part must be filled in as follows -
411         0 - possible
412         >0 - take out and do something (depending on value - TBD)
413         -1 row/column can't vanish but can have entries removed/changed
414         -2 don't touch at all
415         on exit <=0 ones will be in presolved problem
416         struct will be created and will be long enough
417         (information on length etc in first entry)
418         user must delete struct
419     */
420     ClpSimplex * miniPresolve(char * rowType, char * columnType,void ** info);
421     /// After mini presolve
422     void miniPostsolve(const ClpSimplex * presolvedModel,void * info);
423     /** Write the basis in MPS format to the specified file.
424         If writeValues true writes values of structurals
425         (and adds VALUES to end of NAME card)
426
427         Row and column names may be null.
428         formatType is
429         <ul>
430         <li> 0 - normal
431         <li> 1 - extra accuracy
432         <li> 2 - IEEE hex (later)
433         </ul>
434
435         Returns non-zero on I/O error
436     */
437     int writeBasis(const char *filename,
438                    bool writeValues = false,
439                    int formatType = 0) const;
440     /** Read a basis from the given filename,
441         returns -1 on file error, 0 if no values, 1 if values */
442     int readBasis(const char *filename);
443     /// Returns a basis (to be deleted by user)
444     CoinWarmStartBasis * getBasis() const;
445     /// Passes in factorization
446     void setFactorization( ClpFactorization & factorization);
447     // Swaps factorization
448     ClpFactorization * swapFactorization( ClpFactorization * factorization);
449     /// Copies in factorization to existing one
450     void copyFactorization( ClpFactorization & factorization);
451     /** Tightens primal bounds to make dual faster.  Unless
452         fixed or doTight>10, bounds are slightly looser than they could be.
453         This is to make dual go faster and is probably not needed
454         with a presolve.  Returns non-zero if problem infeasible.
455
456         Fudge for branch and bound - put bounds on columns of factor *
457         largest value (at continuous) - should improve stability
458         in branch and bound on infeasible branches (0.0 is off)
459     */
460     int tightenPrimalBounds(double factor = 0.0, int doTight = 0, bool tightIntegers = false);
461     /** Crash - at present just aimed at dual, returns
462         -2 if dual preferred and crash basis created
463         -1 if dual preferred and all slack basis preferred
464          0 if basis going in was not all slack
465          1 if primal preferred and all slack basis preferred
466          2 if primal preferred and crash basis created.
467
468          if gap between bounds <="gap" variables can be flipped
469          ( If pivot -1 then can be made super basic!)
470
471          If "pivot" is
472          -1 No pivoting - always primal
473          0 No pivoting (so will just be choice of algorithm)
474          1 Simple pivoting e.g. gub
475          2 Mini iterations
476     */
477     int crash(double gap, int pivot);
478     /// Sets row pivot choice algorithm in dual
479     void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
480     /// Sets column pivot choice algorithm in primal
481     void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
482     /** For strong branching.  On input lower and upper are new bounds
483         while on output they are change in objective function values
484         (>1.0e50 infeasible).
485         Return code is 0 if nothing interesting, -1 if infeasible both
486         ways and +1 if infeasible one way (check values to see which one(s))
487         Solutions are filled in as well - even down, odd up - also
488         status and number of iterations
489     */
490     int strongBranching(int numberVariables, const int * variables,
491                         double * newLower, double * newUpper,
492                         double ** outputSolution,
493                         int * outputStatus, int * outputIterations,
494                         bool stopOnFirstInfeasible = true,
495                         bool alwaysFinish = false,
496                         int startFinishOptions = 0);
497     /// Fathom - 1 if solution
498     int fathom(void * stuff);
499     /** Do up to N deep - returns
500         -1 - no solution nNodes_ valid nodes
501         >= if solution and that node gives solution
502         ClpNode array is 2**N long.  Values for N and
503         array are in stuff (nNodes_ also in stuff) */
504     int fathomMany(void * stuff);
505     /// Double checks OK
506     double doubleCheck();
507     /// Starts Fast dual2
508     int startFastDual2(ClpNodeStuff * stuff);
509     /// Like Fast dual
510     int fastDual2(ClpNodeStuff * stuff);
511     /// Stops Fast dual2
512     void stopFastDual2(ClpNodeStuff * stuff);
513     /** Deals with crunch aspects
514         mode 0 - in
515              1 - out with solution
516          2 - out without solution
517         returns small model or NULL
518     */
519     ClpSimplex * fastCrunch(ClpNodeStuff * stuff, int mode);
520     //@}
521
522     /**@name Needed for functionality of OsiSimplexInterface */
523     //@{
524     /** Pivot in a variable and out a variable.  Returns 0 if okay,
525         1 if inaccuracy forced re-factorization, -1 if would be singular.
526         Also updates primal/dual infeasibilities.
527         Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
528     */
529     int pivot();
530
531     /** Pivot in a variable and choose an outgoing one.  Assumes primal
532         feasible - will not go through a bound.  Returns step length in theta
533         Returns ray in ray_ (or NULL if no pivot)
534         Return codes as before but -1 means no acceptable pivot
535     */
536     int primalPivotResult();
537
538     /** Pivot out a variable and choose an incoing one.  Assumes dual
539         feasible - will not go through a reduced cost.
540         Returns step length in theta
541         Return codes as before but -1 means no acceptable pivot
542     */
543     int dualPivotResultPart1();
544     /** Do actual pivot
545         state is 0 if need tableau column, 1 if in rowArray_[1]
546     */
547  int pivotResultPart2(int algorithm,int state);
548
549     /** Common bits of coding for dual and primal.  Return 0 if okay,
550         1 if bad matrix, 2 if very bad factorization
551
552         startFinishOptions - bits
553         1 - do not delete work areas and factorization at end
554         2 - use old factorization if same number of rows
555         4 - skip as much initialization of work areas as possible
556             (based on whatsChanged in clpmodel.hpp) ** work in progress
557         maybe other bits later
558
559     */
560     int startup(int ifValuesPass, int startFinishOptions = 0);
561     void finish(int startFinishOptions = 0);
562
563     /** Factorizes and returns true if optimal.  Used by user */
564     bool statusOfProblem(bool initial = false);
565     /// If user left factorization frequency then compute
566     void defaultFactorizationFrequency();
567     /// Copy across enabled stuff from one solver to another
568     void copyEnabledStuff(const ClpSimplex * rhs);
569     //@}
570
571     /**@name most useful gets and sets */
572     //@{
573     /// If problem is primal feasible
574     inline bool primalFeasible() const {
575          return (numberPrimalInfeasibilities_ == 0);
576     }
577     /// If problem is dual feasible
578     inline bool dualFeasible() const {
579          return (numberDualInfeasibilities_ == 0);
580     }
581     /// factorization
582     inline ClpFactorization * factorization() const {
583          return factorization_;
584     }
585     /// Sparsity on or off
586     bool sparseFactorization() const;
587     void setSparseFactorization(bool value);
588     /// Factorization frequency
589     int factorizationFrequency() const;
590     void setFactorizationFrequency(int value);
591     /// Dual bound
592     inline double dualBound() const {
593          return dualBound_;
594     }
595     void setDualBound(double value);
596     /// Infeasibility cost
597     inline double infeasibilityCost() const {
598          return infeasibilityCost_;
599     }
600     void setInfeasibilityCost(double value);
601     /** Amount of print out:
602         0 - none
603         1 - just final
604         2 - just factorizations
605         3 - as 2 plus a bit more
606         4 - verbose
607         above that 8,16,32 etc just for selective debug
608     */
609     /** Perturbation:
610         50  - switch on perturbation
611         100 - auto perturb if takes too long (1.0e-6 largest nonzero)
612         101 - we are perturbed
613         102 - don't try perturbing again
614         default is 100
615         others are for playing
616     */
617     inline int perturbation() const {
618          return perturbation_;
619     }
620     void setPerturbation(int value);
621     /// Current (or last) algorithm
622     inline int algorithm() const {
623          return algorithm_;
624     }
625     /// Set algorithm
626     inline void setAlgorithm(int value) {
627          algorithm_ = value;
628     }
629     /// Return true if the objective limit test can be relied upon
630     bool isObjectiveLimitTestValid() const ;
631     /// Sum of dual infeasibilities
632     inline double sumDualInfeasibilities() const {
633          return sumDualInfeasibilities_;
634     }
635     inline void setSumDualInfeasibilities(double value) {
636          sumDualInfeasibilities_ = value;
637     }
638     /// Sum of relaxed dual infeasibilities
639     inline double sumOfRelaxedDualInfeasibilities() const {
640          return sumOfRelaxedDualInfeasibilities_;
641     }
642     inline void setSumOfRelaxedDualInfeasibilities(double value) {
643          sumOfRelaxedDualInfeasibilities_ = value;
644     }
645     /// Number of dual infeasibilities
646     inline int numberDualInfeasibilities() const {
647          return numberDualInfeasibilities_;
648     }
649     inline void setNumberDualInfeasibilities(int value) {
650          numberDualInfeasibilities_ = value;
651     }
652     /// Number of dual infeasibilities (without free)
653     inline int numberDualInfeasibilitiesWithoutFree() const {
654          return numberDualInfeasibilitiesWithoutFree_;
655     }
656     /// Sum of primal infeasibilities
657     inline double sumPrimalInfeasibilities() const {
658          return sumPrimalInfeasibilities_;
659     }
660     inline void setSumPrimalInfeasibilities(double value) {
661          sumPrimalInfeasibilities_ = value;
662     }
663     /// Sum of relaxed primal infeasibilities
664     inline double sumOfRelaxedPrimalInfeasibilities() const {
665          return sumOfRelaxedPrimalInfeasibilities_;
666     }
667     inline void setSumOfRelaxedPrimalInfeasibilities(double value) {
668          sumOfRelaxedPrimalInfeasibilities_ = value;
669     }
670     /// Number of primal infeasibilities
671     inline int numberPrimalInfeasibilities() const {
672          return numberPrimalInfeasibilities_;
673     }
674     inline void setNumberPrimalInfeasibilities(int value) {
675          numberPrimalInfeasibilities_ = value;
676     }
677     /** Save model to file, returns 0 if success.  This is designed for
678         use outside algorithms so does not save iterating arrays etc.
679     It does not save any messaging information.
680     Does not save scaling values.
681     It does not know about all types of virtual functions.
682     */
683     int saveModel(const char * fileName);
684     /** Restore model from file, returns 0 if success,
685         deletes current model */
686     int restoreModel(const char * fileName);
687
688     /** Just check solution (for external use) - sets sum of
689         infeasibilities etc.
690         If setToBounds 0 then primal column values not changed
691         and used to compute primal row activity values.  If 1 or 2
692         then status used - so all nonbasic variables set to
693         indicated bound and if any values changed (or ==2)  basic values re-computed.
694     */
695     void checkSolution(int setToBounds = 0);
696     /** Just check solution (for internal use) - sets sum of
697         infeasibilities etc. */
698     void checkSolutionInternal();
699     /// Check unscaled primal solution but allow for rounding error
700     void checkUnscaledSolution();
701     /// Useful row length arrays (0,1,2,3,4,5)
702     inline CoinIndexedVector * rowArray(int index) const {
703          return rowArray_[index];
704     }
705     /// Useful column length arrays (0,1,2,3,4,5)
706     inline CoinIndexedVector * columnArray(int index) const {
707          return columnArray_[index];
708     }
709     //@}
710
711     /******************** End of most useful part **************/
712     /**@name Functions less likely to be useful to casual user */
713     //@{
714     /** Given an existing factorization computes and checks
715         primal and dual solutions.  Uses input arrays for variables at
716         bounds.  Returns feasibility states */
717     int getSolution (  const double * rowActivities,
718                        const double * columnActivities);
719     /** Given an existing factorization computes and checks
720         primal and dual solutions.  Uses current problem arrays for
721         bounds.  Returns feasibility states */
722     int getSolution ();
723     /** Constructs a non linear cost from list of non-linearities (columns only)
724         First lower of each column is taken as real lower
725         Last lower is taken as real upper and cost ignored
726
727         Returns nonzero if bad data e.g. lowers not monotonic
728     */
729     int createPiecewiseLinearCosts(const int * starts,
730                                    const double * lower, const double * gradient);
731     /// dual row pivot choice
732     inline ClpDualRowPivot * dualRowPivot() const {
733          return dualRowPivot_;
734     }
735     /// primal column pivot choice
736     inline ClpPrimalColumnPivot * primalColumnPivot() const {
737          return primalColumnPivot_;
738     }
739     /// Returns true if model looks OK
740     inline bool goodAccuracy() const {
741          return (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7);
742     }
743     /** Return model - updates any scalars */
744     void returnModel(ClpSimplex & otherModel);
745     /** Factorizes using current basis.
746         solveType - 1 iterating, 0 initial, -1 external
747         If 10 added then in primal values pass
748         Return codes are as from ClpFactorization unless initial factorization
749         when total number of singularities is returned.
750         Special case is numberRows_+1 -> all slack basis.
751     */
752     int internalFactorize(int solveType);
753     /// Save data
754     ClpDataSave saveData() ;
755     /// Restore data
756     void restoreData(ClpDataSave saved);
757     /// Clean up status
758     void cleanStatus();
759     /// Factorizes using current basis. For external use
760     int factorize();
761     /** Computes duals from scratch. If givenDjs then
762         allows for nonzero basic djs */
763     void computeDuals(double * givenDjs);
764     /// Computes primals from scratch
765     void computePrimals (  const double * rowActivities,
766                            const double * columnActivities);
767     /** Adds multiple of a column into an array */
768     void add(double * array,
769              int column, double multiplier) const;
770     /**
771        Unpacks one column of the matrix into indexed array
772        Uses sequenceIn_
773        Also applies scaling if needed
774     */
775     void unpack(CoinIndexedVector * rowArray) const ;
776     /**
777        Unpacks one column of the matrix into indexed array
778        Slack if sequence>= numberColumns
779        Also applies scaling if needed
780     */
781     void unpack(CoinIndexedVector * rowArray, int sequence) const;
782     /**
783        Unpacks one column of the matrix into indexed array
784        ** as packed vector
785        Uses sequenceIn_
786        Also applies scaling if needed
787     */
788     void unpackPacked(CoinIndexedVector * rowArray) ;
789     /**
790        Unpacks one column of the matrix into indexed array
791        ** as packed vector
792        Slack if sequence>= numberColumns
793        Also applies scaling if needed
794     */
795     void unpackPacked(CoinIndexedVector * rowArray, int sequence);
796#ifndef CLP_USER_DRIVEN
797protected:
798#endif
799     /**
800         This does basis housekeeping and does values for in/out variables.
801         Can also decide to re-factorize
802     */
803     int housekeeping(double objectiveChange);
804     /** This sets largest infeasibility and most infeasible and sum
805         and number of infeasibilities (Primal) */
806     void checkPrimalSolution(const double * rowActivities = NULL,
807                              const double * columnActivies = NULL);
808     /** This sets largest infeasibility and most infeasible and sum
809         and number of infeasibilities (Dual) */
810     void checkDualSolution();
811     /** This sets sum and number of infeasibilities (Dual and Primal) */
812     void checkBothSolutions();
813     /**  If input negative scales objective so maximum <= -value
814          and returns scale factor used.  If positive unscales and also
815          redoes dual stuff
816     */
817     double scaleObjective(double value);
818     /// Solve using Dantzig-Wolfe decomposition and maybe in parallel
819     int solveDW(CoinStructuredModel * model);
820     /// Solve using Benders decomposition and maybe in parallel
821     int solveBenders(CoinStructuredModel * model);
822public:
823     /** For advanced use.  When doing iterative solves things can get
824         nasty so on values pass if incoming solution has largest
825         infeasibility < incomingInfeasibility throw out variables
826         from basis until largest infeasibility < allowedInfeasibility
827         or incoming largest infeasibility.
828         If allowedInfeasibility>= incomingInfeasibility this is
829         always possible altough you may end up with an all slack basis.
830
831         Defaults are 1.0,10.0
832     */
833     void setValuesPassAction(double incomingInfeasibility,
834                              double allowedInfeasibility);
835     /** Get a clean factorization - i.e. throw out singularities
836         may do more later */
837     int cleanFactorization(int ifValuesPass);
838     //@}
839     /**@name most useful gets and sets */
840     //@{
841public:
842     /// Initial value for alpha accuracy calculation (-1.0 off)
843     inline double alphaAccuracy() const {
844          return alphaAccuracy_;
845     }
846     inline void setAlphaAccuracy(double value) {
847          alphaAccuracy_ = value;
848     }
849public:
850     /// Objective value
851     //inline double objectiveValue() const {
852     //return (objectiveValue_-bestPossibleImprovement_)*optimizationDirection_ - dblParam_[ClpObjOffset];
853     //}
854     /// Set disaster handler
855     inline void setDisasterHandler(ClpDisasterHandler * handler) {
856          disasterArea_ = handler;
857     }
858     /// Get disaster handler
859     inline ClpDisasterHandler * disasterHandler() const {
860          return disasterArea_;
861     }
862     /// Large bound value (for complementarity etc)
863     inline double largeValue() const {
864          return largeValue_;
865     }
866     void setLargeValue( double value) ;
867     /// Largest error on Ax-b
868     inline double largestPrimalError() const {
869          return largestPrimalError_;
870     }
871     /// Largest error on basic duals
872     inline double largestDualError() const {
873          return largestDualError_;
874     }
875     /// Largest error on Ax-b
876     inline void setLargestPrimalError(double value) {
877          largestPrimalError_ = value;
878     }
879     /// Largest error on basic duals
880     inline void setLargestDualError(double value) {
881          largestDualError_ = value;
882     }
883     /// Get zero tolerance
884     inline double zeroTolerance() const {
885          return zeroTolerance_;/*factorization_->zeroTolerance();*/
886     }
887     /// Set zero tolerance
888     inline void setZeroTolerance( double value) {
889          zeroTolerance_ = value;
890     }
891     /// Basic variables pivoting on which rows
892     inline int * pivotVariable() const {
893          return pivotVariable_;
894     }
895     /// If automatic scaling on
896     inline bool automaticScaling() const {
897          return automaticScale_ != 0;
898     }
899     inline void setAutomaticScaling(bool onOff) {
900          automaticScale_ = onOff ? 1 : 0;
901     }
902     /// Current dual tolerance
903     inline double currentDualTolerance() const {
904          return dualTolerance_;
905     }
906     inline void setCurrentDualTolerance(double value) {
907          dualTolerance_ = value;
908     }
909     /// Current primal tolerance
910     inline double currentPrimalTolerance() const {
911          return primalTolerance_;
912     }
913     inline void setCurrentPrimalTolerance(double value) {
914          primalTolerance_ = value;
915     }
916     /// How many iterative refinements to do
917     inline int numberRefinements() const {
918          return numberRefinements_;
919     }
920     void setNumberRefinements( int value) ;
921     /// Alpha (pivot element) for use by classes e.g. steepestedge
922     inline double alpha() const {
923          return alpha_;
924     }
925     inline void setAlpha(double value) {
926          alpha_ = value;
927     }
928     /// Reduced cost of last incoming for use by classes e.g. steepestedge
929     inline double dualIn() const {
930          return dualIn_;
931     }
932     /// Set reduced cost of last incoming to force error
933     inline void setDualIn(double value) {
934          dualIn_ = value;
935     }
936     /// Pivot Row for use by classes e.g. steepestedge
937     inline int pivotRow() const {
938          return pivotRow_;
939     }
940     inline void setPivotRow(int value) {
941          pivotRow_ = value;
942     }
943     /// value of incoming variable (in Dual)
944     double valueIncomingDual() const;
945     //@}
946
947#ifndef CLP_USER_DRIVEN
948protected:
949#endif
950     /**@name protected methods */
951     //@{
952     /** May change basis and then returns number changed.
953         Computation of solutions may be overriden by given pi and solution
954     */
955     int gutsOfSolution ( double * givenDuals,
956                          const double * givenPrimals,
957                          bool valuesPass = false);
958     /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
959     void gutsOfDelete(int type);
960     /// Does most of copying
961     void gutsOfCopy(const ClpSimplex & rhs);
962     /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
963         1 bit does rows (now and columns), (2 bit does column bounds), 4 bit does objective(s).
964         8 bit does solution scaling in
965         16 bit does rowArray and columnArray indexed vectors
966         and makes row copy if wanted, also sets columnStart_ etc
967         Also creates scaling arrays if needed.  It does scaling if needed.
968         16 also moves solutions etc in to work arrays
969         On 16 returns false if problem "bad" i.e. matrix or bounds bad
970         If startFinishOptions is -1 then called by user in getSolution
971         so do arrays but keep pivotVariable_
972     */
973     bool createRim(int what, bool makeRowCopy = false, int startFinishOptions = 0);
974     /// Does rows and columns
975     void createRim1(bool initial);
976     /// Does objective
977     void createRim4(bool initial);
978     /// Does rows and columns and objective
979     void createRim5(bool initial);
980     /** releases above arrays and does solution scaling out.  May also
981         get rid of factorization data -
982         0 get rid of nothing, 1 get rid of arrays, 2 also factorization
983     */
984     void deleteRim(int getRidOfFactorizationData = 2);
985     /// Sanity check on input rim data (after scaling) - returns true if okay
986     bool sanityCheck();
987     //@}
988public:
989     /**@name public methods */
990     //@{
991     /** Return row or column sections - not as much needed as it
992         once was.  These just map into single arrays */
993     inline double * solutionRegion(int section) const {
994          if (!section) return rowActivityWork_;
995          else return columnActivityWork_;
996     }
997     inline double * djRegion(int section) const {
998          if (!section) return rowReducedCost_;
999          else return reducedCostWork_;
1000     }
1001     inline double * lowerRegion(int section) const {
1002          if (!section) return rowLowerWork_;
1003          else return columnLowerWork_;
1004     }
1005     inline double * upperRegion(int section) const {
1006          if (!section) return rowUpperWork_;
1007          else return columnUpperWork_;
1008     }
1009     inline double * costRegion(int section) const {
1010          if (!section) return rowObjectiveWork_;
1011          else return objectiveWork_;
1012     }
1013     /// Return region as single array
1014     inline double * solutionRegion() const {
1015          return solution_;
1016     }
1017     inline double * djRegion() const {
1018          return dj_;
1019     }
1020     inline double * lowerRegion() const {
1021          return lower_;
1022     }
1023     inline double * upperRegion() const {
1024          return upper_;
1025     }
1026     inline double * costRegion() const {
1027          return cost_;
1028     }
1029     inline Status getStatus(int sequence) const {
1030          return static_cast<Status> (status_[sequence] & 7);
1031     }
1032     inline void setStatus(int sequence, Status newstatus) {
1033          unsigned char & st_byte = status_[sequence];
1034          st_byte = static_cast<unsigned char>(st_byte & ~7);
1035          st_byte = static_cast<unsigned char>(st_byte | newstatus);
1036     }
1037     /// Start or reset using maximumRows_ and Columns_ - true if change
1038     bool startPermanentArrays();
1039     /** Normally the first factorization does sparse coding because
1040         the factorization could be singular.  This allows initial dense
1041         factorization when it is known to be safe
1042     */
1043     void setInitialDenseFactorization(bool onOff);
1044     bool  initialDenseFactorization() const;
1045     /** Return sequence In or Out */
1046     inline int sequenceIn() const {
1047          return sequenceIn_;
1048     }
1049     inline int sequenceOut() const {
1050          return sequenceOut_;
1051     }
1052     /** Set sequenceIn or Out */
1053     inline void  setSequenceIn(int sequence) {
1054          sequenceIn_ = sequence;
1055     }
1056     inline void  setSequenceOut(int sequence) {
1057          sequenceOut_ = sequence;
1058     }
1059     /** Return direction In or Out */
1060     inline int directionIn() const {
1061          return directionIn_;
1062     }
1063     inline int directionOut() const {
1064          return directionOut_;
1065     }
1066     /** Set directionIn or Out */
1067     inline void  setDirectionIn(int direction) {
1068          directionIn_ = direction;
1069     }
1070     inline void  setDirectionOut(int direction) {
1071          directionOut_ = direction;
1072     }
1073     /// Value of Out variable
1074     inline double valueOut() const {
1075          return valueOut_;
1076     }
1077     /// Set value of out variable
1078     inline void setValueOut(double value) {
1079          valueOut_ = value;
1080     }
1081     /// Dual value of Out variable
1082     inline double dualOut() const {
1083          return dualOut_;
1084     }
1085     /// Set dual value of out variable
1086     inline void setDualOut(double value) {
1087          dualOut_ = value;
1088     }
1089     /// Set lower of out variable
1090     inline void setLowerOut(double value) {
1091          lowerOut_ = value;
1092     }
1093     /// Set upper of out variable
1094     inline void setUpperOut(double value) {
1095          upperOut_ = value;
1096     }
1097     /// Set theta of out variable
1098     inline void setTheta(double value) {
1099          theta_ = value;
1100     }
1101     /// Returns 1 if sequence indicates column
1102     inline int isColumn(int sequence) const {
1103          return sequence < numberColumns_ ? 1 : 0;
1104     }
1105     /// Returns sequence number within section
1106     inline int sequenceWithin(int sequence) const {
1107          return sequence < numberColumns_ ? sequence : sequence - numberColumns_;
1108     }
1109     /// Return row or column values
1110     inline double solution(int sequence) {
1111          return solution_[sequence];
1112     }
1113     /// Return address of row or column values
1114     inline double & solutionAddress(int sequence) {
1115          return solution_[sequence];
1116     }
1117     inline double reducedCost(int sequence) {
1118          return dj_[sequence];
1119     }
1120     inline double & reducedCostAddress(int sequence) {
1121          return dj_[sequence];
1122     }
1123     inline double lower(int sequence) {
1124          return lower_[sequence];
1125     }
1126     /// Return address of row or column lower bound
1127     inline double & lowerAddress(int sequence) {
1128          return lower_[sequence];
1129     }
1130     inline double upper(int sequence) {
1131          return upper_[sequence];
1132     }
1133     /// Return address of row or column upper bound
1134     inline double & upperAddress(int sequence) {
1135          return upper_[sequence];
1136     }
1137     inline double cost(int sequence) {
1138          return cost_[sequence];
1139     }
1140     /// Return address of row or column cost
1141     inline double & costAddress(int sequence) {
1142          return cost_[sequence];
1143     }
1144     /// Return original lower bound
1145     inline double originalLower(int iSequence) const {
1146          if (iSequence < numberColumns_) return columnLower_[iSequence];
1147          else
1148               return rowLower_[iSequence-numberColumns_];
1149     }
1150     /// Return original lower bound
1151     inline double originalUpper(int iSequence) const {
1152          if (iSequence < numberColumns_) return columnUpper_[iSequence];
1153          else
1154               return rowUpper_[iSequence-numberColumns_];
1155     }
1156     /// Theta (pivot change)
1157     inline double theta() const {
1158          return theta_;
1159     }
1160     /** Best possible improvement using djs (primal) or
1161         obj change by flipping bounds to make dual feasible (dual) */
1162     inline double bestPossibleImprovement() const {
1163          return bestPossibleImprovement_;
1164     }
1165     /// Return pointer to details of costs
1166     inline ClpNonLinearCost * nonLinearCost() const {
1167          return nonLinearCost_;
1168     }
1169     /** Return more special options
1170         1 bit - if presolve says infeasible in ClpSolve return
1171         2 bit - if presolved problem infeasible return
1172         4 bit - keep arrays like upper_ around
1173         8 bit - if factorization kept can still declare optimal at once
1174         16 bit - if checking replaceColumn accuracy before updating
1175         32 bit - say optimal if primal feasible!
1176         64 bit - give up easily in dual (and say infeasible)
1177         128 bit - no objective, 0-1 and in B&B
1178         256 bit - in primal from dual or vice versa
1179         512 bit - alternative use of solveType_
1180         1024 bit - don't do row copy of factorization
1181         2048 bit - perturb in complete fathoming
1182         4096 bit - try more for complete fathoming
1183         8192 bit - don't even think of using primal if user asks for dual (and vv)
1184         16384 bit - in initialSolve so be more flexible
1185         debug
1186         32768 bit - do dual in netlibd
1187         65536 (*3) initial stateDualColumn
1188     */
1189     inline int moreSpecialOptions() const {
1190          return moreSpecialOptions_;
1191     }
1192     /** Set more special options
1193         1 bit - if presolve says infeasible in ClpSolve return
1194         2 bit - if presolved problem infeasible return
1195         4 bit - keep arrays like upper_ around
1196         8 bit - no free or superBasic variables
1197         16 bit - if checking replaceColumn accuracy before updating
1198         32 bit - say optimal if primal feasible!
1199         64 bit - give up easily in dual (and say infeasible)
1200         128 bit - no objective, 0-1 and in B&B
1201         256 bit - in primal from dual or vice versa
1202         512 bit - alternative use of solveType_
1203         1024 bit - don't do row copy of factorization
1204         2048 bit - perturb in complete fathoming
1205         4096 bit - try more for complete fathoming
1206         8192 bit - don't even think of using primal if user asks for dual (and vv)
1207         16384 bit - in initialSolve so be more flexible
1208         32768 bit - don't swap algorithms from dual if small infeasibility
1209         65536 bit - perturb in postsolve cleanup (even if < 10000 rows)
1210     */
1211     inline void setMoreSpecialOptions(int value) {
1212          moreSpecialOptions_ = value;
1213     }
1214     //@}
1215     /**@name status methods */
1216     //@{
1217     inline void setFakeBound(int sequence, FakeBound fakeBound) {
1218          unsigned char & st_byte = status_[sequence];
1219          st_byte = static_cast<unsigned char>(st_byte & ~24);
1220          st_byte = static_cast<unsigned char>(st_byte | (fakeBound << 3));
1221     }
1222     inline FakeBound getFakeBound(int sequence) const {
1223          return static_cast<FakeBound> ((status_[sequence] >> 3) & 3);
1224     }
1225     inline void setRowStatus(int sequence, Status newstatus) {
1226          unsigned char & st_byte = status_[sequence+numberColumns_];
1227          st_byte = static_cast<unsigned char>(st_byte & ~7);
1228          st_byte = static_cast<unsigned char>(st_byte | newstatus);
1229     }
1230     inline Status getRowStatus(int sequence) const {
1231          return static_cast<Status> (status_[sequence+numberColumns_] & 7);
1232     }
1233     inline void setColumnStatus(int sequence, Status newstatus) {
1234          unsigned char & st_byte = status_[sequence];
1235          st_byte = static_cast<unsigned char>(st_byte & ~7);
1236          st_byte = static_cast<unsigned char>(st_byte | newstatus);
1237     }
1238     inline Status getColumnStatus(int sequence) const {
1239          return static_cast<Status> (status_[sequence] & 7);
1240     }
1241     inline void setPivoted( int sequence) {
1242          status_[sequence] = static_cast<unsigned char>(status_[sequence] | 32);
1243     }
1244     inline void clearPivoted( int sequence) {
1245          status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~32);
1246     }
1247     inline bool pivoted(int sequence) const {
1248          return (((status_[sequence] >> 5) & 1) != 0);
1249     }
1250     /// To flag a variable (not inline to allow for column generation)
1251     void setFlagged( int sequence);
1252     inline void clearFlagged( int sequence) {
1253          status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~64);
1254     }
1255     inline bool flagged(int sequence) const {
1256          return ((status_[sequence] & 64) != 0);
1257     }
1258     /// To say row active in primal pivot row choice
1259     inline void setActive( int iRow) {
1260          status_[iRow] = static_cast<unsigned char>(status_[iRow] | 128);
1261     }
1262     inline void clearActive( int iRow) {
1263          status_[iRow] = static_cast<unsigned char>(status_[iRow] & ~128);
1264     }
1265     inline bool active(int iRow) const {
1266          return ((status_[iRow] & 128) != 0);
1267     }
1268     /** Set up status array (can be used by OsiClp).
1269         Also can be used to set up all slack basis */
1270     void createStatus() ;
1271     /** Sets up all slack basis and resets solution to
1272         as it was after initial load or readMps */
1273     void allSlackBasis(bool resetSolution = false);
1274
1275     /// So we know when to be cautious
1276     inline int lastBadIteration() const {
1277          return lastBadIteration_;
1278     }
1279     /// Set so we know when to be cautious
1280     inline void setLastBadIteration(int value) {
1281          lastBadIteration_=value;
1282     }
1283     /// Progress flag - at present 0 bit says artificials out
1284     inline int progressFlag() const {
1285          return (progressFlag_ & 3);
1286     }
1287     /// For dealing with all issues of cycling etc
1288     inline ClpSimplexProgress * progress()
1289     { return &progress_;}
1290     /// Force re-factorization early value
1291     inline int forceFactorization() const {
1292          return forceFactorization_ ;
1293     }
1294     /// Force re-factorization early
1295     inline void forceFactorization(int value) {
1296          forceFactorization_ = value;
1297     }
1298     /// Raw objective value (so always minimize in primal)
1299     inline double rawObjectiveValue() const {
1300          return objectiveValue_;
1301     }
1302     /// Compute objective value from solution and put in objectiveValue_
1303     void computeObjectiveValue(bool useWorkingSolution = false);
1304     /// Compute minimization objective value from internal solution without perturbation
1305     double computeInternalObjectiveValue();
1306     /** Infeasibility/unbounded ray (NULL returned if none/wrong)
1307         Up to user to use delete [] on these arrays.  */
1308     double * infeasibilityRay(bool fullRay=false) const;
1309     /** Number of extra rows.  These are ones which will be dynamically created
1310         each iteration.  This is for GUB but may have other uses.
1311     */
1312     inline int numberExtraRows() const {
1313          return numberExtraRows_;
1314     }
1315     /** Maximum number of basic variables - can be more than number of rows if GUB
1316     */
1317     inline int maximumBasic() const {
1318          return maximumBasic_;
1319     }
1320     /// Iteration when we entered dual or primal
1321     inline int baseIteration() const {
1322          return baseIteration_;
1323     }
1324     /// Create C++ lines to get to current state
1325     void generateCpp( FILE * fp, bool defaultFactor = false);
1326     /// Gets clean and emptyish factorization
1327     ClpFactorization * getEmptyFactorization();
1328     /// May delete or may make clean and emptyish factorization
1329     void setEmptyFactorization();
1330     /// Move status and solution across
1331     void moveInfo(const ClpSimplex & rhs, bool justStatus = false);
1332     //@}
1333
1334     ///@name Basis handling
1335     // These are only to be used using startFinishOptions (ClpSimplexDual, ClpSimplexPrimal)
1336     // *** At present only without scaling
1337     // *** Slacks havve -1.0 element (so == row activity) - take care
1338     ///Get a row of the tableau (slack part in slack if not NULL)
1339     void getBInvARow(int row, double* z, double * slack = NULL);
1340
1341     ///Get a row of the basis inverse
1342     void getBInvRow(int row, double* z);
1343
1344     ///Get a column of the tableau
1345     void getBInvACol(int col, double* vec);
1346
1347     ///Get a column of the basis inverse
1348     void getBInvCol(int col, double* vec);
1349
1350     /** Get basic indices (order of indices corresponds to the
1351         order of elements in a vector retured by getBInvACol() and
1352         getBInvCol()).
1353     */
1354     void getBasics(int* index);
1355
1356     //@}
1357     //-------------------------------------------------------------------------
1358     /**@name Changing bounds on variables and constraints */
1359     //@{
1360     /** Set an objective function coefficient */
1361     void setObjectiveCoefficient( int elementIndex, double elementValue );
1362     /** Set an objective function coefficient */
1363     inline void setObjCoeff( int elementIndex, double elementValue ) {
1364          setObjectiveCoefficient( elementIndex, elementValue);
1365     }
1366
1367     /** Set a single column lower bound<br>
1368         Use -DBL_MAX for -infinity. */
1369     void setColumnLower( int elementIndex, double elementValue );
1370
1371     /** Set a single column upper bound<br>
1372         Use DBL_MAX for infinity. */
1373     void setColumnUpper( int elementIndex, double elementValue );
1374
1375     /** Set a single column lower and upper bound */
1376     void setColumnBounds( int elementIndex,
1377                           double lower, double upper );
1378
1379     /** Set the bounds on a number of columns simultaneously<br>
1380         The default implementation just invokes setColLower() and
1381         setColUpper() over and over again.
1382         @param indexFirst,indexLast pointers to the beginning and after the
1383            end of the array of the indices of the variables whose
1384        <em>either</em> bound changes
1385         @param boundList the new lower/upper bound pairs for the variables
1386     */
1387     void setColumnSetBounds(const int* indexFirst,
1388                             const int* indexLast,
1389                             const double* boundList);
1390
1391     /** Set a single column lower bound<br>
1392         Use -DBL_MAX for -infinity. */
1393     inline void setColLower( int elementIndex, double elementValue ) {
1394          setColumnLower(elementIndex, elementValue);
1395     }
1396     /** Set a single column upper bound<br>
1397         Use DBL_MAX for infinity. */
1398     inline void setColUpper( int elementIndex, double elementValue ) {
1399          setColumnUpper(elementIndex, elementValue);
1400     }
1401
1402     /** Set a single column lower and upper bound */
1403     inline void setColBounds( int elementIndex,
1404                               double newlower, double newupper ) {
1405          setColumnBounds(elementIndex, newlower, newupper);
1406     }
1407
1408     /** Set the bounds on a number of columns simultaneously<br>
1409         @param indexFirst,indexLast pointers to the beginning and after the
1410            end of the array of the indices of the variables whose
1411        <em>either</em> bound changes
1412         @param boundList the new lower/upper bound pairs for the variables
1413     */
1414     inline void setColSetBounds(const int* indexFirst,
1415                                 const int* indexLast,
1416                                 const double* boundList) {
1417          setColumnSetBounds(indexFirst, indexLast, boundList);
1418     }
1419
1420     /** Set a single row lower bound<br>
1421         Use -DBL_MAX for -infinity. */
1422     void setRowLower( int elementIndex, double elementValue );
1423
1424     /** Set a single row upper bound<br>
1425         Use DBL_MAX for infinity. */
1426     void setRowUpper( int elementIndex, double elementValue ) ;
1427
1428     /** Set a single row lower and upper bound */
1429     void setRowBounds( int elementIndex,
1430                        double lower, double upper ) ;
1431
1432     /** Set the bounds on a number of rows simultaneously<br>
1433         @param indexFirst,indexLast pointers to the beginning and after the
1434            end of the array of the indices of the constraints whose
1435        <em>either</em> bound changes
1436         @param boundList the new lower/upper bound pairs for the constraints
1437     */
1438     void setRowSetBounds(const int* indexFirst,
1439                          const int* indexLast,
1440                          const double* boundList);
1441     /// Resizes rim part of model
1442     void resize (int newNumberRows, int newNumberColumns);
1443
1444     //@}
1445
1446////////////////// data //////////////////
1447protected:
1448
1449     /**@name data.  Many arrays have a row part and a column part.
1450      There is a single array with both - columns then rows and
1451      then normally two arrays pointing to rows and columns.  The
1452      single array is the owner of memory
1453     */
1454     //@{
1455     /** Best possible improvement using djs (primal) or
1456         obj change by flipping bounds to make dual feasible (dual) */
1457     double bestPossibleImprovement_;
1458     /// Zero tolerance
1459     double zeroTolerance_;
1460     /// Sequence of worst (-1 if feasible)
1461     int columnPrimalSequence_;
1462     /// Sequence of worst (-1 if feasible)
1463     int rowPrimalSequence_;
1464     /// "Best" objective value
1465     double bestObjectiveValue_;
1466     /// More special options - see set for details
1467     int moreSpecialOptions_;
1468     /// Iteration when we entered dual or primal
1469     int baseIteration_;
1470     /// Primal tolerance needed to make dual feasible (<largeTolerance)
1471     double primalToleranceToGetOptimal_;
1472     /// Large bound value (for complementarity etc)
1473     double largeValue_;
1474     /// Largest error on Ax-b
1475     double largestPrimalError_;
1476     /// Largest error on basic duals
1477     double largestDualError_;
1478     /// For computing whether to re-factorize
1479     double alphaAccuracy_;
1480     /// Dual bound
1481     double dualBound_;
1482     /// Alpha (pivot element)
1483     double alpha_;
1484     /// Theta (pivot change)
1485     double theta_;
1486     /// Lower Bound on In variable
1487     double lowerIn_;
1488     /// Value of In variable
1489     double valueIn_;
1490     /// Upper Bound on In variable
1491     double upperIn_;
1492     /// Reduced cost of In variable
1493     double dualIn_;
1494     /// Lower Bound on Out variable
1495     double lowerOut_;
1496     /// Value of Out variable
1497     double valueOut_;
1498     /// Upper Bound on Out variable
1499     double upperOut_;
1500     /// Infeasibility (dual) or ? (primal) of Out variable
1501     double dualOut_;
1502     /// Current dual tolerance for algorithm
1503     double dualTolerance_;
1504     /// Current primal tolerance for algorithm
1505     double primalTolerance_;
1506     /// Sum of dual infeasibilities
1507     double sumDualInfeasibilities_;
1508     /// Sum of primal infeasibilities
1509     double sumPrimalInfeasibilities_;
1510     /// Weight assigned to being infeasible in primal
1511     double infeasibilityCost_;
1512     /// Sum of Dual infeasibilities using tolerance based on error in duals
1513     double sumOfRelaxedDualInfeasibilities_;
1514     /// Sum of Primal infeasibilities using tolerance based on error in primals
1515     double sumOfRelaxedPrimalInfeasibilities_;
1516     /// Acceptable pivot value just after factorization
1517     double acceptablePivot_;
1518     /// Working copy of lower bounds (Owner of arrays below)
1519     double * lower_;
1520     /// Row lower bounds - working copy
1521     double * rowLowerWork_;
1522     /// Column lower bounds - working copy
1523     double * columnLowerWork_;
1524     /// Working copy of upper bounds (Owner of arrays below)
1525     double * upper_;
1526     /// Row upper bounds - working copy
1527     double * rowUpperWork_;
1528     /// Column upper bounds - working copy
1529     double * columnUpperWork_;
1530     /// Working copy of objective (Owner of arrays below)
1531     double * cost_;
1532     /// Row objective - working copy
1533     double * rowObjectiveWork_;
1534     /// Column objective - working copy
1535     double * objectiveWork_;
1536     /// Useful row length arrays
1537     CoinIndexedVector * rowArray_[6];
1538     /// Useful column length arrays
1539     CoinIndexedVector * columnArray_[6];
1540     /// Sequence of In variable
1541     int sequenceIn_;
1542     /// Direction of In, 1 going up, -1 going down, 0 not a clude
1543     int directionIn_;
1544     /// Sequence of Out variable
1545     int sequenceOut_;
1546     /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
1547     int directionOut_;
1548     /// Pivot Row
1549     int pivotRow_;
1550     /// Last good iteration (immediately after a re-factorization)
1551     int lastGoodIteration_;
1552     /// Working copy of reduced costs (Owner of arrays below)
1553     double * dj_;
1554     /// Reduced costs of slacks not same as duals (or - duals)
1555     double * rowReducedCost_;
1556     /// Possible scaled reduced costs
1557     double * reducedCostWork_;
1558     /// Working copy of primal solution (Owner of arrays below)
1559     double * solution_;
1560     /// Row activities - working copy
1561     double * rowActivityWork_;
1562     /// Column activities - working copy
1563     double * columnActivityWork_;
1564     /// Number of dual infeasibilities
1565     int numberDualInfeasibilities_;
1566     /// Number of dual infeasibilities (without free)
1567     int numberDualInfeasibilitiesWithoutFree_;
1568     /// Number of primal infeasibilities
1569     int numberPrimalInfeasibilities_;
1570     /// How many iterative refinements to do
1571     int numberRefinements_;
1572     /// dual row pivot choice
1573     ClpDualRowPivot * dualRowPivot_;
1574     /// primal column pivot choice
1575     ClpPrimalColumnPivot * primalColumnPivot_;
1576     /// Basic variables pivoting on which rows
1577     int * pivotVariable_;
1578     /// factorization
1579     ClpFactorization * factorization_;
1580     /// Saved version of solution
1581     double * savedSolution_;
1582     /// Number of times code has tentatively thought optimal
1583     int numberTimesOptimal_;
1584     /// Disaster handler
1585     ClpDisasterHandler * disasterArea_;
1586     /// If change has been made (first attempt at stopping looping)
1587     int changeMade_;
1588     /// Algorithm >0 == Primal, <0 == Dual
1589     int algorithm_;
1590     /** Now for some reliability aids
1591         This forces re-factorization early */
1592     int forceFactorization_;
1593     /** Perturbation:
1594         -50 to +50 - perturb by this power of ten (-6 sounds good)
1595         100 - auto perturb if takes too long (1.0e-6 largest nonzero)
1596         101 - we are perturbed
1597         102 - don't try perturbing again
1598         default is 100
1599     */
1600     int perturbation_;
1601     /// Saved status regions
1602     unsigned char * saveStatus_;
1603     /** Very wasteful way of dealing with infeasibilities in primal.
1604         However it will allow non-linearities and use of dual
1605         analysis.  If it doesn't work it can easily be replaced.
1606     */
1607     ClpNonLinearCost * nonLinearCost_;
1608     /// So we know when to be cautious
1609     int lastBadIteration_;
1610     /// So we know when to open up again
1611     int lastFlaggedIteration_;
1612     /// Can be used for count of fake bounds (dual) or fake costs (primal)
1613     int numberFake_;
1614     /// Can be used for count of changed costs (dual) or changed bounds (primal)
1615     int numberChanged_;
1616     /// Progress flag - at present 0 bit says artificials out, 1 free in
1617     int progressFlag_;
1618     /// First free/super-basic variable (-1 if none)
1619     int firstFree_;
1620     /** Number of extra rows.  These are ones which will be dynamically created
1621         each iteration.  This is for GUB but may have other uses.
1622     */
1623     int numberExtraRows_;
1624     /** Maximum number of basic variables - can be more than number of rows if GUB
1625     */
1626     int maximumBasic_;
1627     /// If may skip final factorize then allow up to this pivots (default 20)
1628     int dontFactorizePivots_;
1629     /** For advanced use.  When doing iterative solves things can get
1630         nasty so on values pass if incoming solution has largest
1631         infeasibility < incomingInfeasibility throw out variables
1632         from basis until largest infeasibility < allowedInfeasibility.
1633         if allowedInfeasibility>= incomingInfeasibility this is
1634         always possible altough you may end up with an all slack basis.
1635
1636         Defaults are 1.0,10.0
1637     */
1638     double incomingInfeasibility_;
1639     double allowedInfeasibility_;
1640     /// Automatic scaling of objective and rhs and bounds
1641     int automaticScale_;
1642     /// Maximum perturbation array size (take out when code rewritten)
1643     int maximumPerturbationSize_;
1644     /// Perturbation array (maximumPerturbationSize_)
1645     double * perturbationArray_;
1646     /// A copy of model with certain state - normally without cuts
1647     ClpSimplex * baseModel_;
1648     /// For dealing with all issues of cycling etc
1649     ClpSimplexProgress progress_;
1650#ifdef ABC_INHERIT
1651  AbcSimplex * abcSimplex_;
1652#define CLP_ABC_WANTED 1
1653#define CLP_ABC_WANTED_PARALLEL 2
1654#define CLP_ABC_FULL_DONE 8
1655  // bits 256,512,1024 for crash
1656#endif
1657#define CLP_ABC_BEEN_FEASIBLE 65536
1658  int abcState_;
1659public:
1660     /// Spare int array for passing information [0]!=0 switches on
1661     mutable int spareIntArray_[4];
1662     /// Spare double array for passing information [0]!=0 switches on
1663     mutable double spareDoubleArray_[4];
1664protected:
1665     /// Allow OsiClp certain perks
1666     friend class OsiClpSolverInterface;
1667     //@}
1668};
1669//#############################################################################
1670/** A function that tests the methods in the ClpSimplex class. The
1671    only reason for it not to be a member method is that this way it doesn't
1672    have to be compiled into the library. And that's a gain, because the
1673    library should be compiled with optimization on, but this method should be
1674    compiled with debugging.
1675
1676    It also does some testing of ClpFactorization class
1677 */
1678void
1679ClpSimplexUnitTest(const std::string & mpsDir);
1680
1681// For Devex stuff
1682#define DEVEX_TRY_NORM 1.0e-4
1683#define DEVEX_ADD_ONE 1.0
1684#endif
Note: See TracBrowser for help on using the repository browser.