source: stable/1.16/Clp/src/ClpSimplex.hpp @ 2114

Last change on this file since 2114 was 2114, checked in by forrest, 5 years ago

was not totally threadsafe!

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