source: trunk/Clp/src/ClpSimplex.hpp @ 2071

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

for build pthreads for use in Cbc analyze

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