summaryrefslogtreecommitdiff
path: root/build/Bonmin/include/coin/ClpSimplex.hpp
diff options
context:
space:
mode:
authorHarpreet2016-08-04 15:25:44 +0530
committerHarpreet2016-08-04 15:25:44 +0530
commit9fd2976931c088dc523974afb901e96bad20f73c (patch)
tree22502de6e6988d5cd595290d11266f8432ad825b /build/Bonmin/include/coin/ClpSimplex.hpp
downloadFOSSEE-Optim-toolbox-development-9fd2976931c088dc523974afb901e96bad20f73c.tar.gz
FOSSEE-Optim-toolbox-development-9fd2976931c088dc523974afb901e96bad20f73c.tar.bz2
FOSSEE-Optim-toolbox-development-9fd2976931c088dc523974afb901e96bad20f73c.zip
initial add
Diffstat (limited to 'build/Bonmin/include/coin/ClpSimplex.hpp')
-rw-r--r--build/Bonmin/include/coin/ClpSimplex.hpp1797
1 files changed, 1797 insertions, 0 deletions
diff --git a/build/Bonmin/include/coin/ClpSimplex.hpp b/build/Bonmin/include/coin/ClpSimplex.hpp
new file mode 100644
index 0000000..bab4506
--- /dev/null
+++ b/build/Bonmin/include/coin/ClpSimplex.hpp
@@ -0,0 +1,1797 @@
+/* $Id: ClpSimplex.hpp 2114 2015-02-10 12:12:46Z forrest $ */
+// Copyright (C) 2002, International Business Machines
+// Corporation and others. All Rights Reserved.
+// This code is licensed under the terms of the Eclipse Public License (EPL).
+/*
+ Authors
+
+ John Forrest
+
+ */
+#ifndef ClpSimplex_H
+#define ClpSimplex_H
+
+#include <iostream>
+#include <cfloat>
+#include "ClpModel.hpp"
+#include "ClpMatrixBase.hpp"
+#include "ClpSolve.hpp"
+#include "ClpConfig.h"
+class ClpDualRowPivot;
+class ClpPrimalColumnPivot;
+class ClpFactorization;
+class CoinIndexedVector;
+class ClpNonLinearCost;
+class ClpNodeStuff;
+class CoinStructuredModel;
+class OsiClpSolverInterface;
+class CoinWarmStartBasis;
+class ClpDisasterHandler;
+class ClpConstraint;
+/*
+ May want to use Clp defaults so that with ABC defined but not used
+ it behaves as Clp (and ABC used will be different than if not defined)
+ */
+#ifdef ABC_INHERIT
+#ifndef CLP_INHERIT_MODE
+#define CLP_INHERIT_MODE 1
+#endif
+#ifndef ABC_CLP_DEFAULTS
+#define ABC_CLP_DEFAULTS 0
+#endif
+#else
+#undef ABC_CLP_DEFAULTS
+#define ABC_CLP_DEFAULTS 1
+#endif
+#ifdef CLP_HAS_ABC
+#include "AbcCommon.hpp"
+class AbcTolerancesEtc;
+class AbcSimplex;
+#include "CoinAbcCommon.hpp"
+#endif
+/** This solves LPs using the simplex method
+
+ It inherits from ClpModel and all its arrays are created at
+ algorithm time. Originally I tried to work with model arrays
+ but for simplicity of coding I changed to single arrays with
+ structural variables then row variables. Some coding is still
+ based on old style and needs cleaning up.
+
+ For a description of algorithms:
+
+ for dual see ClpSimplexDual.hpp and at top of ClpSimplexDual.cpp
+ for primal see ClpSimplexPrimal.hpp and at top of ClpSimplexPrimal.cpp
+
+ There is an algorithm data member. + for primal variations
+ and - for dual variations
+
+*/
+
+class ClpSimplex : public ClpModel {
+ friend void ClpSimplexUnitTest(const std::string & mpsDir);
+
+public:
+ /** enums for status of various sorts.
+ First 4 match CoinWarmStartBasis,
+ isFixed means fixed at lower bound and out of basis
+ */
+ enum Status {
+ isFree = 0x00,
+ basic = 0x01,
+ atUpperBound = 0x02,
+ atLowerBound = 0x03,
+ superBasic = 0x04,
+ isFixed = 0x05
+ };
+ // For Dual
+ enum FakeBound {
+ noFake = 0x00,
+ lowerFake = 0x01,
+ upperFake = 0x02,
+ bothFake = 0x03
+ };
+
+ /**@name Constructors and destructor and copy */
+ //@{
+ /// Default constructor
+ ClpSimplex (bool emptyMessages = false );
+
+ /** Copy constructor. May scale depending on mode
+ -1 leave mode as is
+ 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
+ */
+ ClpSimplex(const ClpSimplex & rhs, int scalingMode = -1);
+ /** Copy constructor from model. May scale depending on mode
+ -1 leave mode as is
+ 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
+ */
+ ClpSimplex(const ClpModel & rhs, int scalingMode = -1);
+ /** Subproblem constructor. A subset of whole model is created from the
+ row and column lists given. The new order is given by list order and
+ duplicates are allowed. Name and integer information can be dropped
+ Can optionally modify rhs to take into account variables NOT in list
+ in this case duplicates are not allowed (also see getbackSolution)
+ */
+ ClpSimplex (const ClpModel * wholeModel,
+ int numberRows, const int * whichRows,
+ int numberColumns, const int * whichColumns,
+ bool dropNames = true, bool dropIntegers = true,
+ bool fixOthers = false);
+ /** Subproblem constructor. A subset of whole model is created from the
+ row and column lists given. The new order is given by list order and
+ duplicates are allowed. Name and integer information can be dropped
+ Can optionally modify rhs to take into account variables NOT in list
+ in this case duplicates are not allowed (also see getbackSolution)
+ */
+ ClpSimplex (const ClpSimplex * wholeModel,
+ int numberRows, const int * whichRows,
+ int numberColumns, const int * whichColumns,
+ bool dropNames = true, bool dropIntegers = true,
+ bool fixOthers = false);
+ /** This constructor modifies original ClpSimplex and stores
+ original stuff in created ClpSimplex. It is only to be used in
+ conjunction with originalModel */
+ ClpSimplex (ClpSimplex * wholeModel,
+ int numberColumns, const int * whichColumns);
+ /** This copies back stuff from miniModel and then deletes miniModel.
+ Only to be used with mini constructor */
+ void originalModel(ClpSimplex * miniModel);
+ inline int abcState() const
+ { return abcState_;}
+ inline void setAbcState(int state)
+ { abcState_=state;}
+#ifdef ABC_INHERIT
+ inline AbcSimplex * abcSimplex() const
+ { return abcSimplex_;}
+ inline void setAbcSimplex(AbcSimplex * simplex)
+ { abcSimplex_=simplex;}
+ /// Returns 0 if dual can be skipped
+ int doAbcDual();
+ /// Returns 0 if primal can be skipped
+ int doAbcPrimal(int ifValuesPass);
+#endif
+ /** Array persistence flag
+ If 0 then as now (delete/new)
+ 1 then only do arrays if bigger needed
+ 2 as 1 but give a bit extra if bigger needed
+ */
+ void setPersistenceFlag(int value);
+ /// Save a copy of model with certain state - normally without cuts
+ void makeBaseModel();
+ /// Switch off base model
+ void deleteBaseModel();
+ /// See if we have base model
+ inline ClpSimplex * baseModel() const {
+ return baseModel_;
+ }
+ /** Reset to base model (just size and arrays needed)
+ If model NULL use internal copy
+ */
+ void setToBaseModel(ClpSimplex * model = NULL);
+ /// Assignment operator. This copies the data
+ ClpSimplex & operator=(const ClpSimplex & rhs);
+ /// Destructor
+ ~ClpSimplex ( );
+ // Ones below are just ClpModel with some changes
+ /** Loads a problem (the constraints on the
+ rows are given by lower and upper bounds). If a pointer is 0 then the
+ following values are the default:
+ <ul>
+ <li> <code>colub</code>: all columns have upper bound infinity
+ <li> <code>collb</code>: all columns have lower bound 0
+ <li> <code>rowub</code>: all rows have upper bound infinity
+ <li> <code>rowlb</code>: all rows have lower bound -infinity
+ <li> <code>obj</code>: all variables have 0 objective coefficient
+ </ul>
+ */
+ void loadProblem ( const ClpMatrixBase& matrix,
+ const double* collb, const double* colub,
+ const double* obj,
+ const double* rowlb, const double* rowub,
+ const double * rowObjective = NULL);
+ void loadProblem ( const CoinPackedMatrix& matrix,
+ const double* collb, const double* colub,
+ const double* obj,
+ const double* rowlb, const double* rowub,
+ const double * rowObjective = NULL);
+
+ /** Just like the other loadProblem() method except that the matrix is
+ given in a standard column major ordered format (without gaps). */
+ void loadProblem ( const int numcols, const int numrows,
+ const CoinBigIndex* start, const int* index,
+ const double* value,
+ const double* collb, const double* colub,
+ const double* obj,
+ const double* rowlb, const double* rowub,
+ const double * rowObjective = NULL);
+ /// This one is for after presolve to save memory
+ void loadProblem ( const int numcols, const int numrows,
+ const CoinBigIndex* start, const int* index,
+ const double* value, const int * length,
+ const double* collb, const double* colub,
+ const double* obj,
+ const double* rowlb, const double* rowub,
+ const double * rowObjective = NULL);
+ /** This loads a model from a coinModel object - returns number of errors.
+ If keepSolution true and size is same as current then
+ keeps current status and solution
+ */
+ int loadProblem ( CoinModel & modelObject, bool keepSolution = false);
+ /// Read an mps file from the given filename
+ int readMps(const char *filename,
+ bool keepNames = false,
+ bool ignoreErrors = false);
+ /// Read GMPL files from the given filenames
+ int readGMPL(const char *filename, const char * dataName,
+ bool keepNames = false);
+ /// Read file in LP format from file with name filename.
+ /// See class CoinLpIO for description of this format.
+ int readLp(const char *filename, const double epsilon = 1e-5);
+ /** Borrow model. This is so we dont have to copy large amounts
+ of data around. It assumes a derived class wants to overwrite
+ an empty model with a real one - while it does an algorithm.
+ This is same as ClpModel one, but sets scaling on etc. */
+ void borrowModel(ClpModel & otherModel);
+ void borrowModel(ClpSimplex & otherModel);
+ /// Pass in Event handler (cloned and deleted at end)
+ void passInEventHandler(const ClpEventHandler * eventHandler);
+ /// Puts solution back into small model
+ void getbackSolution(const ClpSimplex & smallModel, const int * whichRow, const int * whichColumn);
+ /** Load nonlinear part of problem from AMPL info
+ Returns 0 if linear
+ 1 if quadratic objective
+ 2 if quadratic constraints
+ 3 if nonlinear objective
+ 4 if nonlinear constraints
+ -1 on failure
+ */
+ int loadNonLinear(void * info, int & numberConstraints,
+ ClpConstraint ** & constraints);
+#ifdef ABC_INHERIT
+ /// Loads tolerances etc
+ void loadTolerancesEtc(const AbcTolerancesEtc & data);
+ /// Unloads tolerances etc
+ void unloadTolerancesEtc(AbcTolerancesEtc & data);
+#endif
+ //@}
+
+ /**@name Functions most useful to user */
+ //@{
+ /** General solve algorithm which can do presolve.
+ See ClpSolve.hpp for options
+ */
+ int initialSolve(ClpSolve & options);
+ /// Default initial solve
+ int initialSolve();
+ /// Dual initial solve
+ int initialDualSolve();
+ /// Primal initial solve
+ int initialPrimalSolve();
+ /// Barrier initial solve
+ int initialBarrierSolve();
+ /// Barrier initial solve, not to be followed by crossover
+ int initialBarrierNoCrossSolve();
+ /** Dual algorithm - see ClpSimplexDual.hpp for method.
+ ifValuesPass==2 just does values pass and then stops.
+
+ startFinishOptions - bits
+ 1 - do not delete work areas and factorization at end
+ 2 - use old factorization if same number of rows
+ 4 - skip as much initialization of work areas as possible
+ (based on whatsChanged in clpmodel.hpp) ** work in progress
+ maybe other bits later
+ */
+ int dual(int ifValuesPass = 0, int startFinishOptions = 0);
+ // If using Debug
+ int dualDebug(int ifValuesPass = 0, int startFinishOptions = 0);
+ /** Primal algorithm - see ClpSimplexPrimal.hpp for method.
+ ifValuesPass==2 just does values pass and then stops.
+
+ startFinishOptions - bits
+ 1 - do not delete work areas and factorization at end
+ 2 - use old factorization if same number of rows
+ 4 - skip as much initialization of work areas as possible
+ (based on whatsChanged in clpmodel.hpp) ** work in progress
+ maybe other bits later
+ */
+ int primal(int ifValuesPass = 0, int startFinishOptions = 0);
+ /** Solves nonlinear problem using SLP - may be used as crash
+ for other algorithms when number of iterations small.
+ Also exits if all problematical variables are changing
+ less than deltaTolerance
+ */
+ int nonlinearSLP(int numberPasses, double deltaTolerance);
+ /** Solves problem with nonlinear constraints using SLP - may be used as crash
+ for other algorithms when number of iterations small.
+ Also exits if all problematical variables are changing
+ less than deltaTolerance
+ */
+ int nonlinearSLP(int numberConstraints, ClpConstraint ** constraints,
+ int numberPasses, double deltaTolerance);
+ /** Solves using barrier (assumes you have good cholesky factor code).
+ Does crossover to simplex if asked*/
+ int barrier(bool crossover = true);
+ /** Solves non-linear using reduced gradient. Phase = 0 get feasible,
+ =1 use solution */
+ int reducedGradient(int phase = 0);
+ /// Solve using structure of model and maybe in parallel
+ int solve(CoinStructuredModel * model);
+#ifdef ABC_INHERIT
+ /** solvetype 0 for dual, 1 for primal
+ startup 1 for values pass
+ interrupt whether to pass across interrupt handler
+ add 10 to return AbcSimplex
+ */
+ AbcSimplex * dealWithAbc(int solveType,int startUp,bool interrupt=false);
+ //void dealWithAbc(int solveType,int startUp,bool interrupt=false);
+#endif
+ /** This loads a model from a CoinStructuredModel object - returns number of errors.
+ If originalOrder then keep to order stored in blocks,
+ otherwise first column/rows correspond to first block - etc.
+ If keepSolution true and size is same as current then
+ keeps current status and solution
+ */
+ int loadProblem ( CoinStructuredModel & modelObject,
+ bool originalOrder = true, bool keepSolution = false);
+ /**
+ When scaling is on it is possible that the scaled problem
+ is feasible but the unscaled is not. Clp returns a secondary
+ status code to that effect. This option allows for a cleanup.
+ If you use it I would suggest 1.
+ This only affects actions when scaled optimal
+ 0 - no action
+ 1 - clean up using dual if primal infeasibility
+ 2 - clean up using dual if dual infeasibility
+ 3 - clean up using dual if primal or dual infeasibility
+ 11,12,13 - as 1,2,3 but use primal
+
+ return code as dual/primal
+ */
+ int cleanup(int cleanupScaling);
+ /** Dual ranging.
+ This computes increase/decrease in cost for each given variable and corresponding
+ sequence numbers which would change basis. Sequence numbers are 0..numberColumns
+ and numberColumns.. for artificials/slacks.
+ For non-basic variables the information is trivial to compute and the change in cost is just minus the
+ reduced cost and the sequence number will be that of the non-basic variables.
+ For basic variables a ratio test is between the reduced costs for non-basic variables
+ and the row of the tableau corresponding to the basic variable.
+ The increase/decrease value is always >= 0.0
+
+ Up to user to provide correct length arrays where each array is of length numberCheck.
+ which contains list of variables for which information is desired. All other
+ arrays will be filled in by function. If fifth entry in which is variable 7 then fifth entry in output arrays
+ will be information for variable 7.
+
+ If valueIncrease/Decrease not NULL (both must be NULL or both non NULL) then these are filled with
+ the value of variable if such a change in cost were made (the existing bounds are ignored)
+
+ Returns non-zero if infeasible unbounded etc
+ */
+ int dualRanging(int numberCheck, const int * which,
+ double * costIncrease, int * sequenceIncrease,
+ double * costDecrease, int * sequenceDecrease,
+ double * valueIncrease = NULL, double * valueDecrease = NULL);
+ /** Primal ranging.
+ This computes increase/decrease in value for each given variable and corresponding
+ sequence numbers which would change basis. Sequence numbers are 0..numberColumns
+ and numberColumns.. for artificials/slacks.
+ This should only be used for non-basic variabls as otherwise information is pretty useless
+ For basic variables the sequence number will be that of the basic variables.
+
+ Up to user to provide correct length arrays where each array is of length numberCheck.
+ which contains list of variables for which information is desired. All other
+ arrays will be filled in by function. If fifth entry in which is variable 7 then fifth entry in output arrays
+ will be information for variable 7.
+
+ Returns non-zero if infeasible unbounded etc
+ */
+ int primalRanging(int numberCheck, const int * which,
+ double * valueIncrease, int * sequenceIncrease,
+ double * valueDecrease, int * sequenceDecrease);
+ /**
+ Modifies coefficients etc and if necessary pivots in and out.
+ All at same status will be done (basis may go singular).
+ User can tell which others have been done (i.e. if status matches).
+ If called from outside will change status and return 0.
+ If called from event handler returns non-zero if user has to take action.
+ indices>=numberColumns are slacks (obviously no coefficients)
+ status array is (char) Status enum
+ */
+ int modifyCoefficientsAndPivot(int number,
+ const int * which,
+ const CoinBigIndex * start,
+ const int * row,
+ const double * newCoefficient,
+ const unsigned char * newStatus=NULL,
+ const double * newLower=NULL,
+ const double * newUpper=NULL,
+ const double * newObjective=NULL);
+ /** Take out duplicate rows (includes scaled rows and intersections).
+ On exit whichRows has rows to delete - return code is number can be deleted
+ or -1 if would be infeasible.
+ If tolerance is -1.0 use primalTolerance for equality rows and infeasibility
+ If cleanUp not zero then spend more time trying to leave more stable row
+ and make row bounds exact multiple of cleanUp if close enough
+ */
+ int outDuplicateRows(int numberLook,int * whichRows, bool noOverlaps=false, double tolerance=-1.0,
+ double cleanUp=0.0);
+ /** Try simple crash like techniques to get closer to primal feasibility
+ returns final sum of infeasibilities */
+ double moveTowardsPrimalFeasible();
+ /** Try simple crash like techniques to remove super basic slacks
+ but only if > threshold */
+ void removeSuperBasicSlacks(int threshold=0);
+ /** Mini presolve (faster)
+ Char arrays must be numberRows and numberColumns long
+ on entry second part must be filled in as follows -
+ 0 - possible
+ >0 - take out and do something (depending on value - TBD)
+ -1 row/column can't vanish but can have entries removed/changed
+ -2 don't touch at all
+ on exit <=0 ones will be in presolved problem
+ struct will be created and will be long enough
+ (information on length etc in first entry)
+ user must delete struct
+ */
+ ClpSimplex * miniPresolve(char * rowType, char * columnType,void ** info);
+ /// After mini presolve
+ void miniPostsolve(const ClpSimplex * presolvedModel,void * info);
+ /// mini presolve and solve
+ void miniSolve(char * rowType, char *columnType,int algorithm, int startUp);
+ /** Write the basis in MPS format to the specified file.
+ If writeValues true writes values of structurals
+ (and adds VALUES to end of NAME card)
+
+ Row and column names may be null.
+ formatType is
+ <ul>
+ <li> 0 - normal
+ <li> 1 - extra accuracy
+ <li> 2 - IEEE hex (later)
+ </ul>
+
+ Returns non-zero on I/O error
+ */
+ int writeBasis(const char *filename,
+ bool writeValues = false,
+ int formatType = 0) const;
+ /** Read a basis from the given filename,
+ returns -1 on file error, 0 if no values, 1 if values */
+ int readBasis(const char *filename);
+ /// Returns a basis (to be deleted by user)
+ CoinWarmStartBasis * getBasis() const;
+ /// Passes in factorization
+ void setFactorization( ClpFactorization & factorization);
+ // Swaps factorization
+ ClpFactorization * swapFactorization( ClpFactorization * factorization);
+ /// Copies in factorization to existing one
+ void copyFactorization( ClpFactorization & factorization);
+ /** Tightens primal bounds to make dual faster. Unless
+ fixed or doTight>10, bounds are slightly looser than they could be.
+ This is to make dual go faster and is probably not needed
+ with a presolve. Returns non-zero if problem infeasible.
+
+ Fudge for branch and bound - put bounds on columns of factor *
+ largest value (at continuous) - should improve stability
+ in branch and bound on infeasible branches (0.0 is off)
+ */
+ int tightenPrimalBounds(double factor = 0.0, int doTight = 0, bool tightIntegers = false);
+ /** Crash - at present just aimed at dual, returns
+ -2 if dual preferred and crash basis created
+ -1 if dual preferred and all slack basis preferred
+ 0 if basis going in was not all slack
+ 1 if primal preferred and all slack basis preferred
+ 2 if primal preferred and crash basis created.
+
+ if gap between bounds <="gap" variables can be flipped
+ ( If pivot -1 then can be made super basic!)
+
+ If "pivot" is
+ -1 No pivoting - always primal
+ 0 No pivoting (so will just be choice of algorithm)
+ 1 Simple pivoting e.g. gub
+ 2 Mini iterations
+ */
+ int crash(double gap, int pivot);
+ /// Sets row pivot choice algorithm in dual
+ void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
+ /// Sets column pivot choice algorithm in primal
+ void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
+ /// Create a hotstart point of the optimization process
+ void markHotStart(void * & saveStuff);
+ /// Optimize starting from the hotstart
+ void solveFromHotStart(void * saveStuff);
+ /// Delete the snapshot
+ void unmarkHotStart(void * saveStuff);
+ /** For strong branching. On input lower and upper are new bounds
+ while on output they are change in objective function values
+ (>1.0e50 infeasible).
+ Return code is 0 if nothing interesting, -1 if infeasible both
+ ways and +1 if infeasible one way (check values to see which one(s))
+ Solutions are filled in as well - even down, odd up - also
+ status and number of iterations
+ */
+ int strongBranching(int numberVariables, const int * variables,
+ double * newLower, double * newUpper,
+ double ** outputSolution,
+ int * outputStatus, int * outputIterations,
+ bool stopOnFirstInfeasible = true,
+ bool alwaysFinish = false,
+ int startFinishOptions = 0);
+ /// Fathom - 1 if solution
+ int fathom(void * stuff);
+ /** Do up to N deep - returns
+ -1 - no solution nNodes_ valid nodes
+ >= if solution and that node gives solution
+ ClpNode array is 2**N long. Values for N and
+ array are in stuff (nNodes_ also in stuff) */
+ int fathomMany(void * stuff);
+ /// Double checks OK
+ double doubleCheck();
+ /// Starts Fast dual2
+ int startFastDual2(ClpNodeStuff * stuff);
+ /// Like Fast dual
+ int fastDual2(ClpNodeStuff * stuff);
+ /// Stops Fast dual2
+ void stopFastDual2(ClpNodeStuff * stuff);
+ /** Deals with crunch aspects
+ mode 0 - in
+ 1 - out with solution
+ 2 - out without solution
+ returns small model or NULL
+ */
+ ClpSimplex * fastCrunch(ClpNodeStuff * stuff, int mode);
+ //@}
+
+ /**@name Needed for functionality of OsiSimplexInterface */
+ //@{
+ /** Pivot in a variable and out a variable. Returns 0 if okay,
+ 1 if inaccuracy forced re-factorization, -1 if would be singular.
+ Also updates primal/dual infeasibilities.
+ Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
+ */
+ int pivot();
+
+ /** Pivot in a variable and choose an outgoing one. Assumes primal
+ feasible - will not go through a bound. Returns step length in theta
+ Returns ray in ray_ (or NULL if no pivot)
+ Return codes as before but -1 means no acceptable pivot
+ */
+ int primalPivotResult();
+
+ /** Pivot out a variable and choose an incoing one. Assumes dual
+ feasible - will not go through a reduced cost.
+ Returns step length in theta
+ Return codes as before but -1 means no acceptable pivot
+ */
+ int dualPivotResultPart1();
+ /** Do actual pivot
+ state is 0 if need tableau column, 1 if in rowArray_[1]
+ */
+ int pivotResultPart2(int algorithm,int state);
+
+ /** Common bits of coding for dual and primal. Return 0 if okay,
+ 1 if bad matrix, 2 if very bad factorization
+
+ startFinishOptions - bits
+ 1 - do not delete work areas and factorization at end
+ 2 - use old factorization if same number of rows
+ 4 - skip as much initialization of work areas as possible
+ (based on whatsChanged in clpmodel.hpp) ** work in progress
+ maybe other bits later
+
+ */
+ int startup(int ifValuesPass, int startFinishOptions = 0);
+ void finish(int startFinishOptions = 0);
+
+ /** Factorizes and returns true if optimal. Used by user */
+ bool statusOfProblem(bool initial = false);
+ /// If user left factorization frequency then compute
+ void defaultFactorizationFrequency();
+ /// Copy across enabled stuff from one solver to another
+ void copyEnabledStuff(const ClpSimplex * rhs);
+ //@}
+
+ /**@name most useful gets and sets */
+ //@{
+ /// If problem is primal feasible
+ inline bool primalFeasible() const {
+ return (numberPrimalInfeasibilities_ == 0);
+ }
+ /// If problem is dual feasible
+ inline bool dualFeasible() const {
+ return (numberDualInfeasibilities_ == 0);
+ }
+ /// factorization
+ inline ClpFactorization * factorization() const {
+ return factorization_;
+ }
+ /// Sparsity on or off
+ bool sparseFactorization() const;
+ void setSparseFactorization(bool value);
+ /// Factorization frequency
+ int factorizationFrequency() const;
+ void setFactorizationFrequency(int value);
+ /// Dual bound
+ inline double dualBound() const {
+ return dualBound_;
+ }
+ void setDualBound(double value);
+ /// Infeasibility cost
+ inline double infeasibilityCost() const {
+ return infeasibilityCost_;
+ }
+ void setInfeasibilityCost(double value);
+ /** Amount of print out:
+ 0 - none
+ 1 - just final
+ 2 - just factorizations
+ 3 - as 2 plus a bit more
+ 4 - verbose
+ above that 8,16,32 etc just for selective debug
+ */
+ /** Perturbation:
+ 50 - switch on perturbation
+ 100 - auto perturb if takes too long (1.0e-6 largest nonzero)
+ 101 - we are perturbed
+ 102 - don't try perturbing again
+ default is 100
+ others are for playing
+ */
+ inline int perturbation() const {
+ return perturbation_;
+ }
+ void setPerturbation(int value);
+ /// Current (or last) algorithm
+ inline int algorithm() const {
+ return algorithm_;
+ }
+ /// Set algorithm
+ inline void setAlgorithm(int value) {
+ algorithm_ = value;
+ }
+ /// Return true if the objective limit test can be relied upon
+ bool isObjectiveLimitTestValid() const ;
+ /// Sum of dual infeasibilities
+ inline double sumDualInfeasibilities() const {
+ return sumDualInfeasibilities_;
+ }
+ inline void setSumDualInfeasibilities(double value) {
+ sumDualInfeasibilities_ = value;
+ }
+ /// Sum of relaxed dual infeasibilities
+ inline double sumOfRelaxedDualInfeasibilities() const {
+ return sumOfRelaxedDualInfeasibilities_;
+ }
+ inline void setSumOfRelaxedDualInfeasibilities(double value) {
+ sumOfRelaxedDualInfeasibilities_ = value;
+ }
+ /// Number of dual infeasibilities
+ inline int numberDualInfeasibilities() const {
+ return numberDualInfeasibilities_;
+ }
+ inline void setNumberDualInfeasibilities(int value) {
+ numberDualInfeasibilities_ = value;
+ }
+ /// Number of dual infeasibilities (without free)
+ inline int numberDualInfeasibilitiesWithoutFree() const {
+ return numberDualInfeasibilitiesWithoutFree_;
+ }
+ /// Sum of primal infeasibilities
+ inline double sumPrimalInfeasibilities() const {
+ return sumPrimalInfeasibilities_;
+ }
+ inline void setSumPrimalInfeasibilities(double value) {
+ sumPrimalInfeasibilities_ = value;
+ }
+ /// Sum of relaxed primal infeasibilities
+ inline double sumOfRelaxedPrimalInfeasibilities() const {
+ return sumOfRelaxedPrimalInfeasibilities_;
+ }
+ inline void setSumOfRelaxedPrimalInfeasibilities(double value) {
+ sumOfRelaxedPrimalInfeasibilities_ = value;
+ }
+ /// Number of primal infeasibilities
+ inline int numberPrimalInfeasibilities() const {
+ return numberPrimalInfeasibilities_;
+ }
+ inline void setNumberPrimalInfeasibilities(int value) {
+ numberPrimalInfeasibilities_ = value;
+ }
+ /** Save model to file, returns 0 if success. This is designed for
+ use outside algorithms so does not save iterating arrays etc.
+ It does not save any messaging information.
+ Does not save scaling values.
+ It does not know about all types of virtual functions.
+ */
+ int saveModel(const char * fileName);
+ /** Restore model from file, returns 0 if success,
+ deletes current model */
+ int restoreModel(const char * fileName);
+
+ /** Just check solution (for external use) - sets sum of
+ infeasibilities etc.
+ If setToBounds 0 then primal column values not changed
+ and used to compute primal row activity values. If 1 or 2
+ then status used - so all nonbasic variables set to
+ indicated bound and if any values changed (or ==2) basic values re-computed.
+ */
+ void checkSolution(int setToBounds = 0);
+ /** Just check solution (for internal use) - sets sum of
+ infeasibilities etc. */
+ void checkSolutionInternal();
+ /// Check unscaled primal solution but allow for rounding error
+ void checkUnscaledSolution();
+ /// Useful row length arrays (0,1,2,3,4,5)
+ inline CoinIndexedVector * rowArray(int index) const {
+ return rowArray_[index];
+ }
+ /// Useful column length arrays (0,1,2,3,4,5)
+ inline CoinIndexedVector * columnArray(int index) const {
+ return columnArray_[index];
+ }
+ //@}
+
+ /******************** End of most useful part **************/
+ /**@name Functions less likely to be useful to casual user */
+ //@{
+ /** Given an existing factorization computes and checks
+ primal and dual solutions. Uses input arrays for variables at
+ bounds. Returns feasibility states */
+ int getSolution ( const double * rowActivities,
+ const double * columnActivities);
+ /** Given an existing factorization computes and checks
+ primal and dual solutions. Uses current problem arrays for
+ bounds. Returns feasibility states */
+ int getSolution ();
+ /** Constructs a non linear cost from list of non-linearities (columns only)
+ First lower of each column is taken as real lower
+ Last lower is taken as real upper and cost ignored
+
+ Returns nonzero if bad data e.g. lowers not monotonic
+ */
+ int createPiecewiseLinearCosts(const int * starts,
+ const double * lower, const double * gradient);
+ /// dual row pivot choice
+ inline ClpDualRowPivot * dualRowPivot() const {
+ return dualRowPivot_;
+ }
+ /// primal column pivot choice
+ inline ClpPrimalColumnPivot * primalColumnPivot() const {
+ return primalColumnPivot_;
+ }
+ /// Returns true if model looks OK
+ inline bool goodAccuracy() const {
+ return (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7);
+ }
+ /** Return model - updates any scalars */
+ void returnModel(ClpSimplex & otherModel);
+ /** Factorizes using current basis.
+ solveType - 1 iterating, 0 initial, -1 external
+ If 10 added then in primal values pass
+ Return codes are as from ClpFactorization unless initial factorization
+ when total number of singularities is returned.
+ Special case is numberRows_+1 -> all slack basis.
+ */
+ int internalFactorize(int solveType);
+ /// Save data
+ ClpDataSave saveData() ;
+ /// Restore data
+ void restoreData(ClpDataSave saved);
+ /// Clean up status
+ void cleanStatus();
+ /// Factorizes using current basis. For external use
+ int factorize();
+ /** Computes duals from scratch. If givenDjs then
+ allows for nonzero basic djs */
+ void computeDuals(double * givenDjs);
+ /// Computes primals from scratch
+ void computePrimals ( const double * rowActivities,
+ const double * columnActivities);
+ /** Adds multiple of a column into an array */
+ void add(double * array,
+ int column, double multiplier) const;
+ /**
+ Unpacks one column of the matrix into indexed array
+ Uses sequenceIn_
+ Also applies scaling if needed
+ */
+ void unpack(CoinIndexedVector * rowArray) const ;
+ /**
+ Unpacks one column of the matrix into indexed array
+ Slack if sequence>= numberColumns
+ Also applies scaling if needed
+ */
+ void unpack(CoinIndexedVector * rowArray, int sequence) const;
+ /**
+ Unpacks one column of the matrix into indexed array
+ ** as packed vector
+ Uses sequenceIn_
+ Also applies scaling if needed
+ */
+ void unpackPacked(CoinIndexedVector * rowArray) ;
+ /**
+ Unpacks one column of the matrix into indexed array
+ ** as packed vector
+ Slack if sequence>= numberColumns
+ Also applies scaling if needed
+ */
+ void unpackPacked(CoinIndexedVector * rowArray, int sequence);
+#ifndef CLP_USER_DRIVEN
+protected:
+#endif
+ /**
+ This does basis housekeeping and does values for in/out variables.
+ Can also decide to re-factorize
+ */
+ int housekeeping(double objectiveChange);
+ /** This sets largest infeasibility and most infeasible and sum
+ and number of infeasibilities (Primal) */
+ void checkPrimalSolution(const double * rowActivities = NULL,
+ const double * columnActivies = NULL);
+ /** This sets largest infeasibility and most infeasible and sum
+ and number of infeasibilities (Dual) */
+ void checkDualSolution();
+ /** This sets sum and number of infeasibilities (Dual and Primal) */
+ void checkBothSolutions();
+ /** If input negative scales objective so maximum <= -value
+ and returns scale factor used. If positive unscales and also
+ redoes dual stuff
+ */
+ double scaleObjective(double value);
+ /// Solve using Dantzig-Wolfe decomposition and maybe in parallel
+ int solveDW(CoinStructuredModel * model, ClpSolve & options);
+ /// Solve using Benders decomposition and maybe in parallel
+ int solveBenders(CoinStructuredModel * model, ClpSolve & options);
+public:
+ /** For advanced use. When doing iterative solves things can get
+ nasty so on values pass if incoming solution has largest
+ infeasibility < incomingInfeasibility throw out variables
+ from basis until largest infeasibility < allowedInfeasibility
+ or incoming largest infeasibility.
+ If allowedInfeasibility>= incomingInfeasibility this is
+ always possible altough you may end up with an all slack basis.
+
+ Defaults are 1.0,10.0
+ */
+ void setValuesPassAction(double incomingInfeasibility,
+ double allowedInfeasibility);
+ /** Get a clean factorization - i.e. throw out singularities
+ may do more later */
+ int cleanFactorization(int ifValuesPass);
+ //@}
+ /**@name most useful gets and sets */
+ //@{
+public:
+ /// Initial value for alpha accuracy calculation (-1.0 off)
+ inline double alphaAccuracy() const {
+ return alphaAccuracy_;
+ }
+ inline void setAlphaAccuracy(double value) {
+ alphaAccuracy_ = value;
+ }
+public:
+ /// Objective value
+ //inline double objectiveValue() const {
+ //return (objectiveValue_-bestPossibleImprovement_)*optimizationDirection_ - dblParam_[ClpObjOffset];
+ //}
+ /// Set disaster handler
+ inline void setDisasterHandler(ClpDisasterHandler * handler) {
+ disasterArea_ = handler;
+ }
+ /// Get disaster handler
+ inline ClpDisasterHandler * disasterHandler() const {
+ return disasterArea_;
+ }
+ /// Large bound value (for complementarity etc)
+ inline double largeValue() const {
+ return largeValue_;
+ }
+ void setLargeValue( double value) ;
+ /// Largest error on Ax-b
+ inline double largestPrimalError() const {
+ return largestPrimalError_;
+ }
+ /// Largest error on basic duals
+ inline double largestDualError() const {
+ return largestDualError_;
+ }
+ /// Largest error on Ax-b
+ inline void setLargestPrimalError(double value) {
+ largestPrimalError_ = value;
+ }
+ /// Largest error on basic duals
+ inline void setLargestDualError(double value) {
+ largestDualError_ = value;
+ }
+ /// Get zero tolerance
+ inline double zeroTolerance() const {
+ return zeroTolerance_;/*factorization_->zeroTolerance();*/
+ }
+ /// Set zero tolerance
+ inline void setZeroTolerance( double value) {
+ zeroTolerance_ = value;
+ }
+ /// Basic variables pivoting on which rows
+ inline int * pivotVariable() const {
+ return pivotVariable_;
+ }
+ /// If automatic scaling on
+ inline bool automaticScaling() const {
+ return automaticScale_ != 0;
+ }
+ inline void setAutomaticScaling(bool onOff) {
+ automaticScale_ = onOff ? 1 : 0;
+ }
+ /// Current dual tolerance
+ inline double currentDualTolerance() const {
+ return dualTolerance_;
+ }
+ inline void setCurrentDualTolerance(double value) {
+ dualTolerance_ = value;
+ }
+ /// Current primal tolerance
+ inline double currentPrimalTolerance() const {
+ return primalTolerance_;
+ }
+ inline void setCurrentPrimalTolerance(double value) {
+ primalTolerance_ = value;
+ }
+ /// How many iterative refinements to do
+ inline int numberRefinements() const {
+ return numberRefinements_;
+ }
+ void setNumberRefinements( int value) ;
+ /// Alpha (pivot element) for use by classes e.g. steepestedge
+ inline double alpha() const {
+ return alpha_;
+ }
+ inline void setAlpha(double value) {
+ alpha_ = value;
+ }
+ /// Reduced cost of last incoming for use by classes e.g. steepestedge
+ inline double dualIn() const {
+ return dualIn_;
+ }
+ /// Set reduced cost of last incoming to force error
+ inline void setDualIn(double value) {
+ dualIn_ = value;
+ }
+ /// Pivot Row for use by classes e.g. steepestedge
+ inline int pivotRow() const {
+ return pivotRow_;
+ }
+ inline void setPivotRow(int value) {
+ pivotRow_ = value;
+ }
+ /// value of incoming variable (in Dual)
+ double valueIncomingDual() const;
+ //@}
+
+#ifndef CLP_USER_DRIVEN
+protected:
+#endif
+ /**@name protected methods */
+ //@{
+ /** May change basis and then returns number changed.
+ Computation of solutions may be overriden by given pi and solution
+ */
+ int gutsOfSolution ( double * givenDuals,
+ const double * givenPrimals,
+ bool valuesPass = false);
+ /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
+ void gutsOfDelete(int type);
+ /// Does most of copying
+ void gutsOfCopy(const ClpSimplex & rhs);
+ /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
+ 1 bit does rows (now and columns), (2 bit does column bounds), 4 bit does objective(s).
+ 8 bit does solution scaling in
+ 16 bit does rowArray and columnArray indexed vectors
+ and makes row copy if wanted, also sets columnStart_ etc
+ Also creates scaling arrays if needed. It does scaling if needed.
+ 16 also moves solutions etc in to work arrays
+ On 16 returns false if problem "bad" i.e. matrix or bounds bad
+ If startFinishOptions is -1 then called by user in getSolution
+ so do arrays but keep pivotVariable_
+ */
+ bool createRim(int what, bool makeRowCopy = false, int startFinishOptions = 0);
+ /// Does rows and columns
+ void createRim1(bool initial);
+ /// Does objective
+ void createRim4(bool initial);
+ /// Does rows and columns and objective
+ void createRim5(bool initial);
+ /** releases above arrays and does solution scaling out. May also
+ get rid of factorization data -
+ 0 get rid of nothing, 1 get rid of arrays, 2 also factorization
+ */
+ void deleteRim(int getRidOfFactorizationData = 2);
+ /// Sanity check on input rim data (after scaling) - returns true if okay
+ bool sanityCheck();
+ //@}
+public:
+ /**@name public methods */
+ //@{
+ /** Return row or column sections - not as much needed as it
+ once was. These just map into single arrays */
+ inline double * solutionRegion(int section) const {
+ if (!section) return rowActivityWork_;
+ else return columnActivityWork_;
+ }
+ inline double * djRegion(int section) const {
+ if (!section) return rowReducedCost_;
+ else return reducedCostWork_;
+ }
+ inline double * lowerRegion(int section) const {
+ if (!section) return rowLowerWork_;
+ else return columnLowerWork_;
+ }
+ inline double * upperRegion(int section) const {
+ if (!section) return rowUpperWork_;
+ else return columnUpperWork_;
+ }
+ inline double * costRegion(int section) const {
+ if (!section) return rowObjectiveWork_;
+ else return objectiveWork_;
+ }
+ /// Return region as single array
+ inline double * solutionRegion() const {
+ return solution_;
+ }
+ inline double * djRegion() const {
+ return dj_;
+ }
+ inline double * lowerRegion() const {
+ return lower_;
+ }
+ inline double * upperRegion() const {
+ return upper_;
+ }
+ inline double * costRegion() const {
+ return cost_;
+ }
+ inline Status getStatus(int sequence) const {
+ return static_cast<Status> (status_[sequence] & 7);
+ }
+ inline void setStatus(int sequence, Status newstatus) {
+ unsigned char & st_byte = status_[sequence];
+ st_byte = static_cast<unsigned char>(st_byte & ~7);
+ st_byte = static_cast<unsigned char>(st_byte | newstatus);
+ }
+ /// Start or reset using maximumRows_ and Columns_ - true if change
+ bool startPermanentArrays();
+ /** Normally the first factorization does sparse coding because
+ the factorization could be singular. This allows initial dense
+ factorization when it is known to be safe
+ */
+ void setInitialDenseFactorization(bool onOff);
+ bool initialDenseFactorization() const;
+ /** Return sequence In or Out */
+ inline int sequenceIn() const {
+ return sequenceIn_;
+ }
+ inline int sequenceOut() const {
+ return sequenceOut_;
+ }
+ /** Set sequenceIn or Out */
+ inline void setSequenceIn(int sequence) {
+ sequenceIn_ = sequence;
+ }
+ inline void setSequenceOut(int sequence) {
+ sequenceOut_ = sequence;
+ }
+ /** Return direction In or Out */
+ inline int directionIn() const {
+ return directionIn_;
+ }
+ inline int directionOut() const {
+ return directionOut_;
+ }
+ /** Set directionIn or Out */
+ inline void setDirectionIn(int direction) {
+ directionIn_ = direction;
+ }
+ inline void setDirectionOut(int direction) {
+ directionOut_ = direction;
+ }
+ /// Value of Out variable
+ inline double valueOut() const {
+ return valueOut_;
+ }
+ /// Set value of out variable
+ inline void setValueOut(double value) {
+ valueOut_ = value;
+ }
+ /// Dual value of Out variable
+ inline double dualOut() const {
+ return dualOut_;
+ }
+ /// Set dual value of out variable
+ inline void setDualOut(double value) {
+ dualOut_ = value;
+ }
+ /// Set lower of out variable
+ inline void setLowerOut(double value) {
+ lowerOut_ = value;
+ }
+ /// Set upper of out variable
+ inline void setUpperOut(double value) {
+ upperOut_ = value;
+ }
+ /// Set theta of out variable
+ inline void setTheta(double value) {
+ theta_ = value;
+ }
+ /// Returns 1 if sequence indicates column
+ inline int isColumn(int sequence) const {
+ return sequence < numberColumns_ ? 1 : 0;
+ }
+ /// Returns sequence number within section
+ inline int sequenceWithin(int sequence) const {
+ return sequence < numberColumns_ ? sequence : sequence - numberColumns_;
+ }
+ /// Return row or column values
+ inline double solution(int sequence) {
+ return solution_[sequence];
+ }
+ /// Return address of row or column values
+ inline double & solutionAddress(int sequence) {
+ return solution_[sequence];
+ }
+ inline double reducedCost(int sequence) {
+ return dj_[sequence];
+ }
+ inline double & reducedCostAddress(int sequence) {
+ return dj_[sequence];
+ }
+ inline double lower(int sequence) {
+ return lower_[sequence];
+ }
+ /// Return address of row or column lower bound
+ inline double & lowerAddress(int sequence) {
+ return lower_[sequence];
+ }
+ inline double upper(int sequence) {
+ return upper_[sequence];
+ }
+ /// Return address of row or column upper bound
+ inline double & upperAddress(int sequence) {
+ return upper_[sequence];
+ }
+ inline double cost(int sequence) {
+ return cost_[sequence];
+ }
+ /// Return address of row or column cost
+ inline double & costAddress(int sequence) {
+ return cost_[sequence];
+ }
+ /// Return original lower bound
+ inline double originalLower(int iSequence) const {
+ if (iSequence < numberColumns_) return columnLower_[iSequence];
+ else
+ return rowLower_[iSequence-numberColumns_];
+ }
+ /// Return original lower bound
+ inline double originalUpper(int iSequence) const {
+ if (iSequence < numberColumns_) return columnUpper_[iSequence];
+ else
+ return rowUpper_[iSequence-numberColumns_];
+ }
+ /// Theta (pivot change)
+ inline double theta() const {
+ return theta_;
+ }
+ /** Best possible improvement using djs (primal) or
+ obj change by flipping bounds to make dual feasible (dual) */
+ inline double bestPossibleImprovement() const {
+ return bestPossibleImprovement_;
+ }
+ /// Return pointer to details of costs
+ inline ClpNonLinearCost * nonLinearCost() const {
+ return nonLinearCost_;
+ }
+ /** Return more special options
+ 1 bit - if presolve says infeasible in ClpSolve return
+ 2 bit - if presolved problem infeasible return
+ 4 bit - keep arrays like upper_ around
+ 8 bit - if factorization kept can still declare optimal at once
+ 16 bit - if checking replaceColumn accuracy before updating
+ 32 bit - say optimal if primal feasible!
+ 64 bit - give up easily in dual (and say infeasible)
+ 128 bit - no objective, 0-1 and in B&B
+ 256 bit - in primal from dual or vice versa
+ 512 bit - alternative use of solveType_
+ 1024 bit - don't do row copy of factorization
+ 2048 bit - perturb in complete fathoming
+ 4096 bit - try more for complete fathoming
+ 8192 bit - don't even think of using primal if user asks for dual (and vv)
+ 16384 bit - in initialSolve so be more flexible
+ 32768 bit - don't swap algorithms from dual if small infeasibility
+ 65536 bit - perturb in postsolve cleanup (even if < 10000 rows)
+ 131072 bit (*3) initial stateDualColumn
+ 524288 bit - stop when primal feasible
+ */
+ inline int moreSpecialOptions() const {
+ return moreSpecialOptions_;
+ }
+ /** Set more special options
+ 1 bit - if presolve says infeasible in ClpSolve return
+ 2 bit - if presolved problem infeasible return
+ 4 bit - keep arrays like upper_ around
+ 8 bit - no free or superBasic variables
+ 16 bit - if checking replaceColumn accuracy before updating
+ 32 bit - say optimal if primal feasible!
+ 64 bit - give up easily in dual (and say infeasible)
+ 128 bit - no objective, 0-1 and in B&B
+ 256 bit - in primal from dual or vice versa
+ 512 bit - alternative use of solveType_
+ 1024 bit - don't do row copy of factorization
+ 2048 bit - perturb in complete fathoming
+ 4096 bit - try more for complete fathoming
+ 8192 bit - don't even think of using primal if user asks for dual (and vv)
+ 16384 bit - in initialSolve so be more flexible
+ 32768 bit - don't swap algorithms from dual if small infeasibility
+ 65536 bit - perturb in postsolve cleanup (even if < 10000 rows)
+ 131072 bit (*3) initial stateDualColumn
+ 524288 bit - stop when primal feasible
+ 1048576 bit - don't perturb even if long time
+ 2097152 bit - no primal in fastDual2 if feasible
+ 4194304 bit - tolerances have been changed by code
+ 8388608 bit - tolerances are dynamic (at first)
+ */
+ inline void setMoreSpecialOptions(int value) {
+ moreSpecialOptions_ = value;
+ }
+ //@}
+ /**@name status methods */
+ //@{
+ inline void setFakeBound(int sequence, FakeBound fakeBound) {
+ unsigned char & st_byte = status_[sequence];
+ st_byte = static_cast<unsigned char>(st_byte & ~24);
+ st_byte = static_cast<unsigned char>(st_byte | (fakeBound << 3));
+ }
+ inline FakeBound getFakeBound(int sequence) const {
+ return static_cast<FakeBound> ((status_[sequence] >> 3) & 3);
+ }
+ inline void setRowStatus(int sequence, Status newstatus) {
+ unsigned char & st_byte = status_[sequence+numberColumns_];
+ st_byte = static_cast<unsigned char>(st_byte & ~7);
+ st_byte = static_cast<unsigned char>(st_byte | newstatus);
+ }
+ inline Status getRowStatus(int sequence) const {
+ return static_cast<Status> (status_[sequence+numberColumns_] & 7);
+ }
+ inline void setColumnStatus(int sequence, Status newstatus) {
+ unsigned char & st_byte = status_[sequence];
+ st_byte = static_cast<unsigned char>(st_byte & ~7);
+ st_byte = static_cast<unsigned char>(st_byte | newstatus);
+ }
+ inline Status getColumnStatus(int sequence) const {
+ return static_cast<Status> (status_[sequence] & 7);
+ }
+ inline void setPivoted( int sequence) {
+ status_[sequence] = static_cast<unsigned char>(status_[sequence] | 32);
+ }
+ inline void clearPivoted( int sequence) {
+ status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~32);
+ }
+ inline bool pivoted(int sequence) const {
+ return (((status_[sequence] >> 5) & 1) != 0);
+ }
+ /// To flag a variable (not inline to allow for column generation)
+ void setFlagged( int sequence);
+ inline void clearFlagged( int sequence) {
+ status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~64);
+ }
+ inline bool flagged(int sequence) const {
+ return ((status_[sequence] & 64) != 0);
+ }
+ /// To say row active in primal pivot row choice
+ inline void setActive( int iRow) {
+ status_[iRow] = static_cast<unsigned char>(status_[iRow] | 128);
+ }
+ inline void clearActive( int iRow) {
+ status_[iRow] = static_cast<unsigned char>(status_[iRow] & ~128);
+ }
+ inline bool active(int iRow) const {
+ return ((status_[iRow] & 128) != 0);
+ }
+ /// To say perturbed
+ inline void setPerturbed( int iSequence) {
+ status_[iSequence] = static_cast<unsigned char>(status_[iSequence] | 128);
+ }
+ inline void clearPerturbed( int iSequence) {
+ status_[iSequence] = static_cast<unsigned char>(status_[iSequence] & ~128);
+ }
+ inline bool perturbed(int iSequence) const {
+ return ((status_[iSequence] & 128) != 0);
+ }
+ /** Set up status array (can be used by OsiClp).
+ Also can be used to set up all slack basis */
+ void createStatus() ;
+ /** Sets up all slack basis and resets solution to
+ as it was after initial load or readMps */
+ void allSlackBasis(bool resetSolution = false);
+
+ /// So we know when to be cautious
+ inline int lastBadIteration() const {
+ return lastBadIteration_;
+ }
+ /// Set so we know when to be cautious
+ inline void setLastBadIteration(int value) {
+ lastBadIteration_=value;
+ }
+ /// Progress flag - at present 0 bit says artificials out
+ inline int progressFlag() const {
+ return (progressFlag_ & 3);
+ }
+ /// For dealing with all issues of cycling etc
+ inline ClpSimplexProgress * progress()
+ { return &progress_;}
+ /// Force re-factorization early value
+ inline int forceFactorization() const {
+ return forceFactorization_ ;
+ }
+ /// Force re-factorization early
+ inline void forceFactorization(int value) {
+ forceFactorization_ = value;
+ }
+ /// Raw objective value (so always minimize in primal)
+ inline double rawObjectiveValue() const {
+ return objectiveValue_;
+ }
+ /// Compute objective value from solution and put in objectiveValue_
+ void computeObjectiveValue(bool useWorkingSolution = false);
+ /// Compute minimization objective value from internal solution without perturbation
+ double computeInternalObjectiveValue();
+ /** Infeasibility/unbounded ray (NULL returned if none/wrong)
+ Up to user to use delete [] on these arrays. */
+ double * infeasibilityRay(bool fullRay=false) const;
+ /** Number of extra rows. These are ones which will be dynamically created
+ each iteration. This is for GUB but may have other uses.
+ */
+ inline int numberExtraRows() const {
+ return numberExtraRows_;
+ }
+ /** Maximum number of basic variables - can be more than number of rows if GUB
+ */
+ inline int maximumBasic() const {
+ return maximumBasic_;
+ }
+ /// Iteration when we entered dual or primal
+ inline int baseIteration() const {
+ return baseIteration_;
+ }
+ /// Create C++ lines to get to current state
+ void generateCpp( FILE * fp, bool defaultFactor = false);
+ /// Gets clean and emptyish factorization
+ ClpFactorization * getEmptyFactorization();
+ /// May delete or may make clean and emptyish factorization
+ void setEmptyFactorization();
+ /// Move status and solution across
+ void moveInfo(const ClpSimplex & rhs, bool justStatus = false);
+ //@}
+
+ ///@name Basis handling
+ // These are only to be used using startFinishOptions (ClpSimplexDual, ClpSimplexPrimal)
+ // *** At present only without scaling
+ // *** Slacks havve -1.0 element (so == row activity) - take care
+ ///Get a row of the tableau (slack part in slack if not NULL)
+ void getBInvARow(int row, double* z, double * slack = NULL);
+
+ ///Get a row of the basis inverse
+ void getBInvRow(int row, double* z);
+
+ ///Get a column of the tableau
+ void getBInvACol(int col, double* vec);
+
+ ///Get a column of the basis inverse
+ void getBInvCol(int col, double* vec);
+
+ /** Get basic indices (order of indices corresponds to the
+ order of elements in a vector retured by getBInvACol() and
+ getBInvCol()).
+ */
+ void getBasics(int* index);
+
+ //@}
+ //-------------------------------------------------------------------------
+ /**@name Changing bounds on variables and constraints */
+ //@{
+ /** Set an objective function coefficient */
+ void setObjectiveCoefficient( int elementIndex, double elementValue );
+ /** Set an objective function coefficient */
+ inline void setObjCoeff( int elementIndex, double elementValue ) {
+ setObjectiveCoefficient( elementIndex, elementValue);
+ }
+
+ /** Set a single column lower bound<br>
+ Use -DBL_MAX for -infinity. */
+ void setColumnLower( int elementIndex, double elementValue );
+
+ /** Set a single column upper bound<br>
+ Use DBL_MAX for infinity. */
+ void setColumnUpper( int elementIndex, double elementValue );
+
+ /** Set a single column lower and upper bound */
+ void setColumnBounds( int elementIndex,
+ double lower, double upper );
+
+ /** Set the bounds on a number of columns simultaneously<br>
+ The default implementation just invokes setColLower() and
+ setColUpper() over and over again.
+ @param indexFirst,indexLast pointers to the beginning and after the
+ end of the array of the indices of the variables whose
+ <em>either</em> bound changes
+ @param boundList the new lower/upper bound pairs for the variables
+ */
+ void setColumnSetBounds(const int* indexFirst,
+ const int* indexLast,
+ const double* boundList);
+
+ /** Set a single column lower bound<br>
+ Use -DBL_MAX for -infinity. */
+ inline void setColLower( int elementIndex, double elementValue ) {
+ setColumnLower(elementIndex, elementValue);
+ }
+ /** Set a single column upper bound<br>
+ Use DBL_MAX for infinity. */
+ inline void setColUpper( int elementIndex, double elementValue ) {
+ setColumnUpper(elementIndex, elementValue);
+ }
+
+ /** Set a single column lower and upper bound */
+ inline void setColBounds( int elementIndex,
+ double newlower, double newupper ) {
+ setColumnBounds(elementIndex, newlower, newupper);
+ }
+
+ /** Set the bounds on a number of columns simultaneously<br>
+ @param indexFirst,indexLast pointers to the beginning and after the
+ end of the array of the indices of the variables whose
+ <em>either</em> bound changes
+ @param boundList the new lower/upper bound pairs for the variables
+ */
+ inline void setColSetBounds(const int* indexFirst,
+ const int* indexLast,
+ const double* boundList) {
+ setColumnSetBounds(indexFirst, indexLast, boundList);
+ }
+
+ /** Set a single row lower bound<br>
+ Use -DBL_MAX for -infinity. */
+ void setRowLower( int elementIndex, double elementValue );
+
+ /** Set a single row upper bound<br>
+ Use DBL_MAX for infinity. */
+ void setRowUpper( int elementIndex, double elementValue ) ;
+
+ /** Set a single row lower and upper bound */
+ void setRowBounds( int elementIndex,
+ double lower, double upper ) ;
+
+ /** Set the bounds on a number of rows simultaneously<br>
+ @param indexFirst,indexLast pointers to the beginning and after the
+ end of the array of the indices of the constraints whose
+ <em>either</em> bound changes
+ @param boundList the new lower/upper bound pairs for the constraints
+ */
+ void setRowSetBounds(const int* indexFirst,
+ const int* indexLast,
+ const double* boundList);
+ /// Resizes rim part of model
+ void resize (int newNumberRows, int newNumberColumns);
+
+ //@}
+
+////////////////// data //////////////////
+protected:
+
+ /**@name data. Many arrays have a row part and a column part.
+ There is a single array with both - columns then rows and
+ then normally two arrays pointing to rows and columns. The
+ single array is the owner of memory
+ */
+ //@{
+ /** Best possible improvement using djs (primal) or
+ obj change by flipping bounds to make dual feasible (dual) */
+ double bestPossibleImprovement_;
+ /// Zero tolerance
+ double zeroTolerance_;
+ /// Sequence of worst (-1 if feasible)
+ int columnPrimalSequence_;
+ /// Sequence of worst (-1 if feasible)
+ int rowPrimalSequence_;
+ /// "Best" objective value
+ double bestObjectiveValue_;
+ /// More special options - see set for details
+ int moreSpecialOptions_;
+ /// Iteration when we entered dual or primal
+ int baseIteration_;
+ /// Primal tolerance needed to make dual feasible (<largeTolerance)
+ double primalToleranceToGetOptimal_;
+ /// Large bound value (for complementarity etc)
+ double largeValue_;
+ /// Largest error on Ax-b
+ double largestPrimalError_;
+ /// Largest error on basic duals
+ double largestDualError_;
+ /// For computing whether to re-factorize
+ double alphaAccuracy_;
+ /// Dual bound
+ double dualBound_;
+ /// Alpha (pivot element)
+ double alpha_;
+ /// Theta (pivot change)
+ double theta_;
+ /// Lower Bound on In variable
+ double lowerIn_;
+ /// Value of In variable
+ double valueIn_;
+ /// Upper Bound on In variable
+ double upperIn_;
+ /// Reduced cost of In variable
+ double dualIn_;
+ /// Lower Bound on Out variable
+ double lowerOut_;
+ /// Value of Out variable
+ double valueOut_;
+ /// Upper Bound on Out variable
+ double upperOut_;
+ /// Infeasibility (dual) or ? (primal) of Out variable
+ double dualOut_;
+ /// Current dual tolerance for algorithm
+ double dualTolerance_;
+ /// Current primal tolerance for algorithm
+ double primalTolerance_;
+ /// Sum of dual infeasibilities
+ double sumDualInfeasibilities_;
+ /// Sum of primal infeasibilities
+ double sumPrimalInfeasibilities_;
+ /// Weight assigned to being infeasible in primal
+ double infeasibilityCost_;
+ /// Sum of Dual infeasibilities using tolerance based on error in duals
+ double sumOfRelaxedDualInfeasibilities_;
+ /// Sum of Primal infeasibilities using tolerance based on error in primals
+ double sumOfRelaxedPrimalInfeasibilities_;
+ /// Acceptable pivot value just after factorization
+ double acceptablePivot_;
+ /// Minimum primal tolerance
+ double minimumPrimalTolerance_;
+ /// Last few infeasibilities
+#define CLP_INFEAS_SAVE 5
+ double averageInfeasibility_[CLP_INFEAS_SAVE];
+ /// Working copy of lower bounds (Owner of arrays below)
+ double * lower_;
+ /// Row lower bounds - working copy
+ double * rowLowerWork_;
+ /// Column lower bounds - working copy
+ double * columnLowerWork_;
+ /// Working copy of upper bounds (Owner of arrays below)
+ double * upper_;
+ /// Row upper bounds - working copy
+ double * rowUpperWork_;
+ /// Column upper bounds - working copy
+ double * columnUpperWork_;
+ /// Working copy of objective (Owner of arrays below)
+ double * cost_;
+ /// Row objective - working copy
+ double * rowObjectiveWork_;
+ /// Column objective - working copy
+ double * objectiveWork_;
+ /// Useful row length arrays
+ CoinIndexedVector * rowArray_[6];
+ /// Useful column length arrays
+ CoinIndexedVector * columnArray_[6];
+ /// Sequence of In variable
+ int sequenceIn_;
+ /// Direction of In, 1 going up, -1 going down, 0 not a clude
+ int directionIn_;
+ /// Sequence of Out variable
+ int sequenceOut_;
+ /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
+ int directionOut_;
+ /// Pivot Row
+ int pivotRow_;
+ /// Last good iteration (immediately after a re-factorization)
+ int lastGoodIteration_;
+ /// Working copy of reduced costs (Owner of arrays below)
+ double * dj_;
+ /// Reduced costs of slacks not same as duals (or - duals)
+ double * rowReducedCost_;
+ /// Possible scaled reduced costs
+ double * reducedCostWork_;
+ /// Working copy of primal solution (Owner of arrays below)
+ double * solution_;
+ /// Row activities - working copy
+ double * rowActivityWork_;
+ /// Column activities - working copy
+ double * columnActivityWork_;
+ /// Number of dual infeasibilities
+ int numberDualInfeasibilities_;
+ /// Number of dual infeasibilities (without free)
+ int numberDualInfeasibilitiesWithoutFree_;
+ /// Number of primal infeasibilities
+ int numberPrimalInfeasibilities_;
+ /// How many iterative refinements to do
+ int numberRefinements_;
+ /// dual row pivot choice
+ ClpDualRowPivot * dualRowPivot_;
+ /// primal column pivot choice
+ ClpPrimalColumnPivot * primalColumnPivot_;
+ /// Basic variables pivoting on which rows
+ int * pivotVariable_;
+ /// factorization
+ ClpFactorization * factorization_;
+ /// Saved version of solution
+ double * savedSolution_;
+ /// Number of times code has tentatively thought optimal
+ int numberTimesOptimal_;
+ /// Disaster handler
+ ClpDisasterHandler * disasterArea_;
+ /// If change has been made (first attempt at stopping looping)
+ int changeMade_;
+ /// Algorithm >0 == Primal, <0 == Dual
+ int algorithm_;
+ /** Now for some reliability aids
+ This forces re-factorization early */
+ int forceFactorization_;
+ /** Perturbation:
+ -50 to +50 - perturb by this power of ten (-6 sounds good)
+ 100 - auto perturb if takes too long (1.0e-6 largest nonzero)
+ 101 - we are perturbed
+ 102 - don't try perturbing again
+ default is 100
+ */
+ int perturbation_;
+ /// Saved status regions
+ unsigned char * saveStatus_;
+ /** Very wasteful way of dealing with infeasibilities in primal.
+ However it will allow non-linearities and use of dual
+ analysis. If it doesn't work it can easily be replaced.
+ */
+ ClpNonLinearCost * nonLinearCost_;
+ /// So we know when to be cautious
+ int lastBadIteration_;
+ /// So we know when to open up again
+ int lastFlaggedIteration_;
+ /// Can be used for count of fake bounds (dual) or fake costs (primal)
+ int numberFake_;
+ /// Can be used for count of changed costs (dual) or changed bounds (primal)
+ int numberChanged_;
+ /// Progress flag - at present 0 bit says artificials out, 1 free in
+ int progressFlag_;
+ /// First free/super-basic variable (-1 if none)
+ int firstFree_;
+ /** Number of extra rows. These are ones which will be dynamically created
+ each iteration. This is for GUB but may have other uses.
+ */
+ int numberExtraRows_;
+ /** Maximum number of basic variables - can be more than number of rows if GUB
+ */
+ int maximumBasic_;
+ /// If may skip final factorize then allow up to this pivots (default 20)
+ int dontFactorizePivots_;
+ /** For advanced use. When doing iterative solves things can get
+ nasty so on values pass if incoming solution has largest
+ infeasibility < incomingInfeasibility throw out variables
+ from basis until largest infeasibility < allowedInfeasibility.
+ if allowedInfeasibility>= incomingInfeasibility this is
+ always possible altough you may end up with an all slack basis.
+
+ Defaults are 1.0,10.0
+ */
+ double incomingInfeasibility_;
+ double allowedInfeasibility_;
+ /// Automatic scaling of objective and rhs and bounds
+ int automaticScale_;
+ /// Maximum perturbation array size (take out when code rewritten)
+ int maximumPerturbationSize_;
+ /// Perturbation array (maximumPerturbationSize_)
+ double * perturbationArray_;
+ /// A copy of model with certain state - normally without cuts
+ ClpSimplex * baseModel_;
+ /// For dealing with all issues of cycling etc
+ ClpSimplexProgress progress_;
+#ifdef ABC_INHERIT
+ AbcSimplex * abcSimplex_;
+#define CLP_ABC_WANTED 1
+#define CLP_ABC_WANTED_PARALLEL 2
+#define CLP_ABC_FULL_DONE 8
+ // bits 256,512,1024 for crash
+#endif
+#define CLP_ABC_BEEN_FEASIBLE 65536
+ int abcState_;
+ /// Number of degenerate pivots since last perturbed
+ int numberDegeneratePivots_;
+public:
+ /// Spare int array for passing information [0]!=0 switches on
+ mutable int spareIntArray_[4];
+ /// Spare double array for passing information [0]!=0 switches on
+ mutable double spareDoubleArray_[4];
+protected:
+ /// Allow OsiClp certain perks
+ friend class OsiClpSolverInterface;
+ /// And OsiCLP
+ friend class OsiCLPSolverInterface;
+ //@}
+};
+//#############################################################################
+/** A function that tests the methods in the ClpSimplex class. The
+ only reason for it not to be a member method is that this way it doesn't
+ have to be compiled into the library. And that's a gain, because the
+ library should be compiled with optimization on, but this method should be
+ compiled with debugging.
+
+ It also does some testing of ClpFactorization class
+ */
+void
+ClpSimplexUnitTest(const std::string & mpsDir);
+
+// For Devex stuff
+#define DEVEX_TRY_NORM 1.0e-4
+#define DEVEX_ADD_ONE 1.0
+#if defined(ABC_INHERIT) || defined(CBC_THREAD) || defined(THREADS_IN_ANALYZE)
+// Use pthreads
+#include <pthread.h>
+typedef struct {
+ double result;
+ //const CoinIndexedVector * constVector; // can get rid of
+ //CoinIndexedVector * vectors[2]; // can get rid of
+ void * extraInfo;
+ void * extraInfo2;
+ int status;
+ int stuff[4];
+} CoinThreadInfo;
+class CoinPthreadStuff {
+public:
+ /**@name Constructors and destructor and copy */
+ //@{
+ /** Main constructor
+ */
+ CoinPthreadStuff (int numberThreads=0,
+ void * parallelManager(void * stuff)=NULL);
+ /// Assignment operator. This copies the data
+ CoinPthreadStuff & operator=(const CoinPthreadStuff & rhs);
+ /// Destructor
+ ~CoinPthreadStuff ( );
+ /// set stop start
+ inline void setStopStart(int value)
+ { stopStart_=value;}
+#ifndef NUMBER_THREADS
+#define NUMBER_THREADS 8
+#endif
+ // For waking up thread
+ inline pthread_mutex_t * mutexPointer(int which,int thread=0)
+ { return mutex_+which+3*thread;}
+#ifdef PTHREAD_BARRIER_SERIAL_THREAD
+ inline pthread_barrier_t * barrierPointer()
+ { return &barrier_;}
+#endif
+ inline int whichLocked(int thread=0) const
+ { return locked_[thread];}
+ inline CoinThreadInfo * threadInfoPointer(int thread=0)
+ { return threadInfo_+thread;}
+ void startParallelTask(int type,int iThread,void * info=NULL);
+ int waitParallelTask(int type, int & iThread,bool allowIdle);
+ void waitAllTasks();
+ /// so thread can find out which one it is
+ int whichThread() const;
+ void sayIdle(int iThread);
+ //void startThreads(int numberThreads);
+ //void stopThreads();
+ // For waking up thread
+ pthread_mutex_t mutex_[3*(NUMBER_THREADS+1)];
+#ifdef PTHREAD_BARRIER_SERIAL_THREAD
+ pthread_barrier_t barrier_;
+#endif
+ CoinThreadInfo threadInfo_[NUMBER_THREADS+1];
+ pthread_t abcThread_[NUMBER_THREADS+1];
+ int locked_[NUMBER_THREADS+1];
+ int stopStart_;
+ int numberThreads_;
+};
+void * clp_parallelManager(void * stuff);
+#endif
+#endif