ClpSimplex Class Reference

#include <ClpSimplex.hpp>

Inheritance diagram for ClpSimplex:

List of all members.

Public Types
enum	Status {   isFree = 0x00, basic = 0x01, atUpperBound = 0x02, atLowerBound = 0x03,   superBasic = 0x04, isFixed = 0x05 }
enum	FakeBound { noFake = 0x00, bothFake = 0x01, upperFake = 0x02, lowerFake = 0x03 }
Public Member Functions
Constructors and destructor and copy
	ClpSimplex ()
	Default constructor.
	ClpSimplex (const ClpSimplex &)
	Copy constructor.
	ClpSimplex (const ClpModel &)
	Copy constructor from model.
	ClpSimplex (const ClpModel *wholeModel, int numberRows, const int *whichRows, int numberColumns, const int *whichColumns, bool dropNames=true, bool dropIntegers=true)
	ClpSimplex (ClpSimplex *wholeModel, int numberColumns, const int *whichColumns)
void	originalModel (ClpSimplex *miniModel)
ClpSimplex &	operator= (const ClpSimplex &rhs)
	Assignment operator. This copies the data.
	~ClpSimplex ()
	Destructor.
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)
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)
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 one is for after presolve to save memory.
int	readMps (const char *filename, bool keepNames=false, bool ignoreErrors=false)
	Read an mps file from the given filename.
void	borrowModel (ClpModel &otherModel)
Functions most useful to user
int	initialSolve (ClpSolve &options)
int	initialSolve ()
	Default initial solve.
int	initialDualSolve ()
	Dual initial solve.
int	initialPrimalSolve ()
	Primal initial solve.
int	dual (int ifValuesPass=0)
int	primal (int ifValuesPass=0)
int	quadraticSLP (int numberPasses, double deltaTolerance)
int	quadraticPrimal (int phase=2)
	Solves quadratic using Dantzig's algorithm - primal.
void	setFactorization (ClpFactorization &factorization)
	Passes in factorization.
void	scaling (int mode=1)
	Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later).
int	scalingFlag () const
	Gets scalingFlag.
int	tightenPrimalBounds (double factor=0.0)
int	crash (double gap, int pivot)
void	setDualRowPivotAlgorithm (ClpDualRowPivot &choice)
	Sets row pivot choice algorithm in dual.
void	setPrimalColumnPivotAlgorithm (ClpPrimalColumnPivot &choice)
	Sets column pivot choice algorithm in primal.
int	strongBranching (int numberVariables, const int *variables, double *newLower, double *newUpper, double **outputSolution, int *outputStatus, int *outputIterations, bool stopOnFirstInfeasible=true, bool alwaysFinish=false)
Needed for functionality of OsiSimplexInterface
int	pivot ()
int	primalPivotResult ()
int	dualPivotResult ()
int	startup (int ifValuesPass)
void	finish ()
bool	statusOfProblem ()
most useful gets and sets
bool	primalFeasible () const
	If problem is primal feasible.
bool	dualFeasible () const
	If problem is dual feasible.
ClpFactorization *	factorization () const
	factorization
bool	sparseFactorization () const
	Sparsity on or off.
void	setSparseFactorization (bool value)
int	factorizationFrequency () const
	Factorization frequency.
void	setFactorizationFrequency (int value)
double	dualBound () const
	Dual bound.
void	setDualBound (double value)
double	infeasibilityCost () const
	Infeasibility cost.
void	setInfeasibilityCost (double value)
int	perturbation () const
void	setPerturbation (int value)
int	algorithm () const
	Current (or last) algorithm.
void	setAlgorithm (int value)
	Set algorithm.
double	sumDualInfeasibilities () const
	Sum of dual infeasibilities.
int	numberDualInfeasibilities () const
	Number of dual infeasibilities.
double	sumPrimalInfeasibilities () const
	Sum of primal infeasibilities.
void	setSumPrimalInfeasibilities (double value)
int	numberPrimalInfeasibilities () const
	Number of primal infeasibilities.
int	saveModel (const char *fileName)
int	restoreModel (const char *fileName)
void	checkSolution ()
CoinIndexedVector *	rowArray (int index) const
	Useful row length arrays (0,1,2,3,4,5).
CoinIndexedVector *	columnArray (int index) const
	Useful column length arrays (0,1,2,3,4,5).
Functions less likely to be useful to casual user
int	getSolution (const double *rowActivities, const double *columnActivities)
int	getSolution ()
int	createPiecewiseLinearCosts (const int *starts, const double *lower, const double *gradient)
void	returnModel (ClpSimplex &otherModel)
ClpDataSave	saveData ()
	Save data.
void	restoreData (ClpDataSave saved)
	Restore data.
int	internalFactorize (int solveType)
void	cleanStatus ()
	Clean up status.
int	factorize ()
	Factorizes using current basis. For external use.
void	computeDuals (double *givenDjs)
void	computePrimals (const double *rowActivities, const double *columnActivities)
	Computes primals from scratch.
void	unpack (CoinIndexedVector *rowArray) const
void	unpack (CoinIndexedVector *rowArray, int sequence) const
void	unpackPacked (CoinIndexedVector *rowArray)
void	unpackPacked (CoinIndexedVector *rowArray, int sequence)
int	housekeeping (double objectiveChange)
void	checkPrimalSolution (const double *rowActivities=NULL, const double *columnActivies=NULL)
void	checkDualSolution ()
void	setValuesPassAction (float incomingInfeasibility, float allowedInfeasibility)
Matrix times vector methods
They can be faster if scalar is +- 1 These are covers so user need not worry about scaling Also for simplex I am not using basic/non-basic split
void	times (double scalar, const double *x, double *y) const
void	transposeTimes (double scalar, const double *x, double *y) const
most useful gets and sets
double	columnPrimalInfeasibility () const
	Worst column primal infeasibility.
int	columnPrimalSequence () const
	Sequence of worst (-1 if feasible).
double	rowPrimalInfeasibility () const
	Worst row primal infeasibility.
int	rowPrimalSequence () const
	Sequence of worst (-1 if feasible).
double	columnDualInfeasibility () const
int	columnDualSequence () const
	Sequence of worst (-1 if feasible).
double	rowDualInfeasibility () const
	Worst row dual infeasibility.
int	rowDualSequence () const
	Sequence of worst (-1 if feasible).
double	primalToleranceToGetOptimal () const
	Primal tolerance needed to make dual feasible (<largeTolerance).
double	remainingDualInfeasibility () const
	Remaining largest dual infeasibility.
double	largeValue () const
	Large bound value (for complementarity etc).
void	setLargeValue (double value)
double	largestPrimalError () const
	Largest error on Ax-b.
double	largestDualError () const
	Largest error on basic duals.
double	largestSolutionError () const
	Largest difference between input primal solution and computed.
int *	pivotVariable () const
	Basic variables pivoting on which rows.
double	objectiveScale () const
	Scaling of objective 12345.0 (auto), 0.0 (off), other user.
void	setObjectiveScale (double value)
double	currentDualTolerance () const
	Current dual tolerance.
void	setCurrentDualTolerance (double value)
double	currentPrimalTolerance () const
	Current primal tolerance.
void	setCurrentPrimalTolerance (double value)
int	numberRefinements () const
	How many iterative refinements to do.
void	setNumberRefinements (int value)
double	alpha () const
	Alpha (pivot element) for use by classes e.g. steepestedge.
double	dualIn () const
	Reduced cost of last incoming for use by classes e.g. steepestedge.
int	pivotRow () const
	Pivot Row for use by classes e.g. steepestedge.
double	valueIncomingDual () const
	value of incoming variable (in Dual)
public methods
double *	solutionRegion (int section)
double *	djRegion (int section)
double *	lowerRegion (int section)
double *	upperRegion (int section)
double *	costRegion (int section)
double *	solutionRegion ()
	Return region as single array.
double *	djRegion ()
double *	lowerRegion ()
double *	upperRegion ()
double *	costRegion ()
Status	getStatus (int sequence) const
void	setStatus (int sequence, Status status)
void	setInitialDenseFactorization (bool onOff)
bool	initialDenseFactorization () const
int	sequenceIn () const
int	sequenceOut () const
void	setSequenceIn (int sequence)
void	setSequenceOut (int sequence)
int	directionIn () const
int	directionOut () const
void	setDirectionIn (int direction)
void	setDirectionOut (int direction)
int	isColumn (int sequence) const
	Returns 1 if sequence indicates column.
int	sequenceWithin (int sequence) const
	Returns sequence number within section.
double	solution (int sequence)
	Return row or column values.
double &	solutionAddress (int sequence)
	Return address of row or column values.
double	reducedCost (int sequence)
double &	reducedCostAddress (int sequence)
double	lower (int sequence)
double &	lowerAddress (int sequence)
	Return address of row or column lower bound.
double	upper (int sequence)
double &	upperAddress (int sequence)
	Return address of row or column upper bound.
double	cost (int sequence)
double &	costAddress (int sequence)
	Return address of row or column cost.
double	originalLower (int iSequence) const
	Return original lower bound.
double	originalUpper (int iSequence) const
	Return original lower bound.
double	theta () const
	Theta (pivot change).
const double *	rowScale () const
	Scaling.
const double *	columnScale () const
void	setRowScale (double *scale)
void	setColumnScale (double *scale)
ClpNonLinearCost *	nonLinearCost () const
	Return pointer to details of costs.
status methods
void	setFakeBound (int sequence, FakeBound fakeBound)
FakeBound	getFakeBound (int sequence) const
void	setRowStatus (int sequence, Status status)
Status	getRowStatus (int sequence) const
void	setColumnStatus (int sequence, Status status)
Status	getColumnStatus (int sequence) const
void	setPivoted (int sequence)
void	clearPivoted (int sequence)
bool	pivoted (int sequence) const
void	setFlagged (int sequence)
	To flag a variable.
void	clearFlagged (int sequence)
bool	flagged (int sequence) const
void	setActive (int iRow)
	To say row active in primal pivot row choice.
void	clearActive (int iRow)
bool	active (int iRow) const
void	createStatus ()
void	allSlackBasis ()
int	lastBadIteration () const
	So we know when to be cautious.
int	progressFlag () const
	Progress flag - at present 0 bit says artificials out.
void	forceFactorization (int value)
	Force re-factorization early.
double	rawObjectiveValue () const
	Raw objective value (so always minimize in primal).
Protected Member Functions
protected methods
int	gutsOfSolution (double *givenDuals, const double *givenPrimals, bool valuesPass=false)
void	gutsOfDelete (int type)
	Does most of deletion (0 = all, 1 = most, 2 most + factorization).
void	gutsOfCopy (const ClpSimplex &rhs)
	Does most of copying.
bool	createRim (int what, bool makeRowCopy=false)
void	deleteRim (int getRidOfFactorizationData=2)
bool	sanityCheck ()
	Sanity check on input rim data (after scaling) - returns true if okay.
Protected Attributes
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
double	columnPrimalInfeasibility_
	Worst column primal infeasibility.
double	rowPrimalInfeasibility_
	Worst row primal infeasibility.
int	columnPrimalSequence_
	Sequence of worst (-1 if feasible).
int	rowPrimalSequence_
	Sequence of worst (-1 if feasible).
double	columnDualInfeasibility_
	Worst column dual infeasibility.
double	rowDualInfeasibility_
	Worst row dual infeasibility.
int	columnDualSequence_
	Sequence of worst (-1 if feasible).
int	rowDualSequence_
	Sequence of worst (-1 if feasible).
double	primalToleranceToGetOptimal_
	Primal tolerance needed to make dual feasible (<largeTolerance).
double	remainingDualInfeasibility_
	Remaining largest dual infeasibility.
double	largeValue_
	Large bound value (for complementarity etc).
double	objectiveScale_
	Scaling of objective.
double	largestPrimalError_
	Largest error on Ax-b.
double	largestDualError_
	Largest error on basic duals.
double	largestSolutionError_
	Largest difference between input primal solution and computed.
double	dualBound_
	Dual bound.
double	alpha_
	Alpha (pivot element).
double	theta_
	Theta (pivot change).
double	lowerIn_
	Lower Bound on In variable.
double	valueIn_
	Value of In variable.
double	upperIn_
	Upper Bound on In variable.
double	dualIn_
	Reduced cost of In variable.
double	lowerOut_
	Lower Bound on Out variable.
double	valueOut_
	Value of Out variable.
double	upperOut_
	Upper Bound on Out variable.
double	dualOut_
	Infeasibility (dual) or ? (primal) of Out variable.
double	dualTolerance_
	Current dual tolerance for algorithm.
double	primalTolerance_
	Current primal tolerance for algorithm.
double	sumDualInfeasibilities_
	Sum of dual infeasibilities.
double	sumPrimalInfeasibilities_
	Sum of primal infeasibilities.
double	infeasibilityCost_
	Weight assigned to being infeasible in primal.
double	sumOfRelaxedDualInfeasibilities_
	Sum of Dual infeasibilities using tolerance based on error in duals.
double	sumOfRelaxedPrimalInfeasibilities_
	Sum of Primal infeasibilities using tolerance based on error in primals.
double *	lower_
	Working copy of lower bounds (Owner of arrays below).
double *	rowLowerWork_
	Row lower bounds - working copy.
double *	columnLowerWork_
	Column lower bounds - working copy.
double *	upper_
	Working copy of upper bounds (Owner of arrays below).
double *	rowUpperWork_
	Row upper bounds - working copy.
double *	columnUpperWork_
	Column upper bounds - working copy.
double *	cost_
	Working copy of objective (Owner of arrays below).
double *	rowObjectiveWork_
	Row objective - working copy.
double *	objectiveWork_
	Column objective - working copy.
CoinIndexedVector *	rowArray_ [6]
	Useful row length arrays.
CoinIndexedVector *	columnArray_ [6]
	Useful column length arrays.
int	sequenceIn_
	Sequence of In variable.
int	directionIn_
	Direction of In, 1 going up, -1 going down, 0 not a clude.
int	sequenceOut_
	Sequence of Out variable.
int	directionOut_
	Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic.
int	pivotRow_
	Pivot Row.
int	lastGoodIteration_
	Last good iteration (immediately after a re-factorization).
double *	dj_
	Working copy of reduced costs (Owner of arrays below).
double *	rowReducedCost_
	Reduced costs of slacks not same as duals (or - duals).
double *	reducedCostWork_
	Possible scaled reduced costs.
double *	solution_
	Working copy of primal solution (Owner of arrays below).
double *	rowActivityWork_
	Row activities - working copy.
double *	columnActivityWork_
	Column activities - working copy.
int	numberDualInfeasibilities_
	Number of dual infeasibilities.
int	numberDualInfeasibilitiesWithoutFree_
	Number of dual infeasibilities (without free).
int	numberPrimalInfeasibilities_
	Number of primal infeasibilities.
int	numberRefinements_
	How many iterative refinements to do.
ClpDualRowPivot *	dualRowPivot_
	dual row pivot choice
ClpPrimalColumnPivot *	primalColumnPivot_
	primal column pivot choice
int *	pivotVariable_
	Basic variables pivoting on which rows.
ClpFactorization *	factorization_
	factorization
double *	rowScale_
	Row scale factors for matrix.
double *	savedSolution_
	Saved version of solution.
double *	columnScale_
	Column scale factors.
int	scalingFlag_
	Scale flag, 0 none, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic.
int	numberTimesOptimal_
	Number of times code has tentatively thought optimal.
int	changeMade_
	If change has been made (first attempt at stopping looping).
int	algorithm_
	Algorithm >0 == Primal, <0 == Dual.
int	forceFactorization_
int	perturbation_
unsigned char *	saveStatus_
	Saved status regions.
ClpNonLinearCost *	nonLinearCost_
unsigned int	specialOptions_
int	lastBadIteration_
	So we know when to be cautious.
int	numberFake_
	Can be used for count of fake bounds (dual) or fake costs (primal).
int	progressFlag_
	Progress flag - at present 0 bit says artificials out, 1 free in.
int	firstFree_
	First free/super-basic variable (-1 if none).
float	incomingInfeasibility_
float	allowedInfeasibility_
ClpSimplexProgress *	progress_
	For dealing with all issues of cycling etc.
Friends
void	ClpSimplexUnitTest (const std::string &mpsDir, const std::string &netlibDir)

---

Detailed Description

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

This file also includes (at end) a very simple class ClpSimplexProgress which is where anti-looping stuff should migrate to

Definition at line 46 of file ClpSimplex.hpp.

---

Member Enumeration Documentation

enum ClpSimplex::Status

enums for status of various sorts. First 4 match CoinWarmStartBasis, isFixed means fixed at lower bound and out of basis Definition at line 55 of file ClpSimplex.hpp. { isFree = 0x00, basic = 0x01, atUpperBound = 0x02, atLowerBound = 0x03, superBasic = 0x04, isFixed = 0x05 };

---

Constructor & Destructor Documentation

ClpSimplex::ClpSimplex (  const ClpModel *  wholeModel, int  numberRows, const int *  whichRows, int  numberColumns, const int *  whichColumns, bool  dropNames = true, bool  dropIntegers = true )

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 Definition at line 130 of file ClpSimplex.cpp. References columnArray_, dualRowPivot_, factorization_, primalColumnPivot_, rowArray_, and saveStatus_. : ClpModel(rhs, numberRows, whichRow, numberColumns,whichColumn,dropNames,dropIntegers), columnPrimalInfeasibility_(0.0), rowPrimalInfeasibility_(0.0), columnPrimalSequence_(-2), rowPrimalSequence_(-2), columnDualInfeasibility_(0.0), rowDualInfeasibility_(0.0), columnDualSequence_(-2), rowDualSequence_(-2), primalToleranceToGetOptimal_(-1.0), remainingDualInfeasibility_(0.0), largeValue_(1.0e15), largestPrimalError_(0.0), largestDualError_(0.0), largestSolutionError_(0.0), dualBound_(1.0e10), alpha_(0.0), theta_(0.0), lowerIn_(0.0), valueIn_(0.0), upperIn_(0.0), dualIn_(0.0), lowerOut_(-1), valueOut_(-1), upperOut_(-1), dualOut_(-1), dualTolerance_(0.0), primalTolerance_(0.0), sumDualInfeasibilities_(0.0), sumPrimalInfeasibilities_(0.0), infeasibilityCost_(1.0e10), sumOfRelaxedDualInfeasibilities_(0.0), sumOfRelaxedPrimalInfeasibilities_(0.0), lower_(NULL), rowLowerWork_(NULL), columnLowerWork_(NULL), upper_(NULL), rowUpperWork_(NULL), columnUpperWork_(NULL), cost_(NULL), rowObjectiveWork_(NULL), objectiveWork_(NULL), sequenceIn_(-1), directionIn_(-1), sequenceOut_(-1), directionOut_(-1), pivotRow_(-1), lastGoodIteration_(-100), dj_(NULL), rowReducedCost_(NULL), reducedCostWork_(NULL), solution_(NULL), rowActivityWork_(NULL), columnActivityWork_(NULL), numberDualInfeasibilities_(0), numberDualInfeasibilitiesWithoutFree_(0), numberPrimalInfeasibilities_(100), numberRefinements_(0), pivotVariable_(NULL), factorization_(NULL), rowScale_(NULL), savedSolution_(NULL), columnScale_(NULL), scalingFlag_(3), numberTimesOptimal_(0), changeMade_(1), algorithm_(0), forceFactorization_(-1), perturbation_(100), nonLinearCost_(NULL), specialOptions_(0), lastBadIteration_(-999999), numberFake_(0), progressFlag_(0), firstFree_(-1), incomingInfeasibility_(1.0), allowedInfeasibility_(10.0), progress_(NULL) { int i; for (i=0;i<6;i++) { rowArray_[i]=NULL; columnArray_[i]=NULL; } saveStatus_=NULL; // get an empty factorization so we can set tolerances etc factorization_ = new ClpFactorization(); // say Steepest pricing dualRowPivot_ = new ClpDualRowSteepest(); // say Steepest pricing primalColumnPivot_ = new ClpPrimalColumnSteepest(); solveType_=1; // say simplex based life form }

ClpSimplex::ClpSimplex (  ClpSimplex *  wholeModel, int  numberColumns, const int *  whichColumns )

This constructor modifies original ClpSimplex and stores original stuff in created ClpSimplex. It is only to be used in conjunction with originalModel Definition at line 4944 of file ClpSimplex.cpp. References ClpNonLinearCost::checkInfeasibilities(), columnActivityWork_, columnLowerWork_, columnScale_, columnUpperWork_, cost_, createRim(), dj_, ClpModel::getDblParam(), getStatus(), lower_, ClpModel::matrix_, nonLinearCost_, ClpModel::numberColumns_, ClpNonLinearCost::numberInfeasibilities(), ClpModel::numberRows_, objectiveWork_, pivotVariable_, primalColumnPivot_, reducedCostWork_, rowActivityWork_, ClpModel::rowCopy_, rowLowerWork_, rowObjectiveWork_, rowReducedCost_, rowScale_, rowUpperWork_, savedSolution_, saveStatus_, ClpPrimalColumnPivot::saveWeights(), ClpModel::setDblParam(), solution_, ClpModel::status_, ClpMatrixBase::subsetClone(), ClpNonLinearCost::sumInfeasibilities(), ClpMatrixBase::times(), and upper_. { // Set up dummy row selection list numberRows_ = wholeModel->numberRows_; int * whichRow = new int [numberRows_]; int iRow; for (iRow=0;iRow<numberRows_;iRow++) whichRow[iRow]=iRow; // ClpModel stuff (apart from numberColumns_) matrix_ = wholeModel->matrix_; rowCopy_ = wholeModel->rowCopy_; if (wholeModel->rowCopy_) { // note reversal of order wholeModel->rowCopy_ = wholeModel->rowCopy_->subsetClone(numberRows_,whichRow, numberColumns,whichColumns); } else { wholeModel->rowCopy_=NULL; } assert (wholeModel->matrix_); wholeModel->matrix_ = wholeModel->matrix_->subsetClone(numberRows_,whichRow, numberColumns,whichColumns); delete [] whichRow; numberColumns_ = wholeModel->numberColumns_; // Now ClpSimplex stuff and status_ ClpPrimalColumnSteepest * steep = dynamic_cast< ClpPrimalColumnSteepest*>(wholeModel->primalColumnPivot_); #ifdef NDEBUG if (!steep) abort(); #else assert (steep); #endif delete wholeModel->primalColumnPivot_; wholeModel->primalColumnPivot_ = new ClpPrimalColumnSteepest(0); nonLinearCost_ = wholeModel->nonLinearCost_; // Now main arrays int iColumn; int numberTotal = numberRows_+numberColumns; printf("%d %d %d\n",numberTotal,numberRows_,numberColumns); // mapping int * mapping = new int[numberRows_+numberColumns_]; for (iColumn=0;iColumn<numberColumns_;iColumn++) mapping[iColumn]=-1; for (iRow=0;iRow<numberRows_;iRow++) mapping[iRow+numberColumns_] = iRow+numberColumns; // Redo costs and bounds of whole model wholeModel->createRim(5,false); lower_ = wholeModel->lower_; wholeModel->lower_ = new double [numberTotal]; memcpy(wholeModel->lower_+numberColumns,lower_+numberColumns_,numberRows_*sizeof(double)); for (iColumn=0;iColumn<numberColumns;iColumn++) { int jColumn = whichColumns[iColumn]; wholeModel->lower_[iColumn]=lower_[jColumn]; // and pointer back mapping[jColumn]=iColumn; } #ifdef CLP_DEBUG for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) printf("mapx %d %d\n",iColumn,mapping[iColumn]); #endif // Re-define pivotVariable_ for (iRow=0;iRow<numberRows_;iRow++) { int iPivot = wholeModel->pivotVariable_[iRow]; wholeModel->pivotVariable_[iRow]=mapping[iPivot]; #ifdef CLP_DEBUG printf("p row %d, pivot %d -> %d\n",iRow,iPivot,mapping[iPivot]); #endif assert (wholeModel->pivotVariable_[iRow]>=0); } // Reverse mapping (so extended version of whichColumns) for (iColumn=0;iColumn<numberColumns;iColumn++) mapping[iColumn]=whichColumns[iColumn]; for (;iColumn<numberRows_+numberColumns;iColumn++) mapping[iColumn] = iColumn + (numberColumns_-numberColumns); #ifdef CLP_DEBUG for (iColumn=0;iColumn<numberRows_+numberColumns;iColumn++) printf("map %d %d\n",iColumn,mapping[iColumn]); #endif // Save mapping somewhere - doesn't matter rowUpper_ = (double *) mapping; upper_ = wholeModel->upper_; wholeModel->upper_ = new double [numberTotal]; for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; wholeModel->upper_[iColumn]=upper_[jColumn]; } cost_ = wholeModel->cost_; wholeModel->cost_ = new double [numberTotal]; for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; wholeModel->cost_[iColumn]=cost_[jColumn]; } dj_ = wholeModel->dj_; wholeModel->dj_ = new double [numberTotal]; for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; wholeModel->dj_[iColumn]=dj_[jColumn]; } solution_ = wholeModel->solution_; wholeModel->solution_ = new double [numberTotal]; for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; wholeModel->solution_[iColumn]=solution_[jColumn]; } // now see what variables left out do to row solution double * rowSolution = wholeModel->solution_+numberColumns; double * fullSolution = solution_; double * sumFixed = new double[numberRows_]; memset (sumFixed,0,numberRows_*sizeof(double)); // zero out ones in small problem for (iColumn=0;iColumn<numberColumns;iColumn++) { int jColumn = mapping[iColumn]; fullSolution[jColumn]=0.0; } // Get objective offset double originalOffset; wholeModel->getDblParam(ClpObjOffset,originalOffset); double offset=0.0; const double * cost = cost_; for (iColumn=0;iColumn<numberColumns_;iColumn++) offset += fullSolution[iColumn]*cost[iColumn]; wholeModel->setDblParam(ClpObjOffset,originalOffset-offset); setDblParam(ClpObjOffset,originalOffset); matrix_->times(1.0,fullSolution,sumFixed,wholeModel->rowScale_,wholeModel->columnScale_); double * lower = lower_+numberColumns; double * upper = upper_+numberColumns; double fixed=0.0; for (iRow=0;iRow<numberRows_;iRow++) { fixed += fabs(sumFixed[iRow]); if (lower[iRow]>-1.0e50) lower[iRow] -= sumFixed[iRow]; if (upper[iRow]<1.0e50) upper[iRow] -= sumFixed[iRow]; rowSolution[iRow] -= sumFixed[iRow]; } printf("offset %g sumfixed %g\n",offset,fixed); delete [] sumFixed; columnScale_ = wholeModel->columnScale_; if (columnScale_) { wholeModel->columnScale_ = new double [numberTotal]; for (iColumn=0;iColumn<numberColumns;iColumn++) { int jColumn = mapping[iColumn]; wholeModel->columnScale_[iColumn]=columnScale_[jColumn]; } } status_ = wholeModel->status_; wholeModel->status_ = new unsigned char [numberTotal]; for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; wholeModel->status_[iColumn]=status_[jColumn]; } savedSolution_ = wholeModel->savedSolution_; if (savedSolution_) { wholeModel->savedSolution_ = new double [numberTotal]; for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; wholeModel->savedSolution_[iColumn]=savedSolution_[jColumn]; } } saveStatus_ = wholeModel->saveStatus_; if (saveStatus_) { wholeModel->saveStatus_ = new unsigned char [numberTotal]; for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; wholeModel->saveStatus_[iColumn]=saveStatus_[jColumn]; } } wholeModel->numberColumns_ = numberColumns; // Initialize weights wholeModel->primalColumnPivot_->saveWeights(wholeModel,2); // Costs wholeModel->nonLinearCost_ = new ClpNonLinearCost(wholeModel); wholeModel->nonLinearCost_->checkInfeasibilities(); printf("after contraction %d infeasibilities summing to %g\n", nonLinearCost_->numberInfeasibilities(),nonLinearCost_->sumInfeasibilities()); // Redo some stuff wholeModel->reducedCostWork_ = wholeModel->dj_; wholeModel->rowReducedCost_ = wholeModel->dj_+wholeModel->numberColumns_; wholeModel->columnActivityWork_ = wholeModel->solution_; wholeModel->rowActivityWork_ = wholeModel->solution_+wholeModel->numberColumns_; wholeModel->objectiveWork_ = wholeModel->cost_; wholeModel->rowObjectiveWork_ = wholeModel->cost_+wholeModel->numberColumns_; wholeModel->rowLowerWork_ = wholeModel->lower_+wholeModel->numberColumns_; wholeModel->columnLowerWork_ = wholeModel->lower_; wholeModel->rowUpperWork_ = wholeModel->upper_+wholeModel->numberColumns_; wholeModel->columnUpperWork_ = wholeModel->upper_; #ifndef NDEBUG // Check status ClpSimplex * xxxx = wholeModel; int nBasic=0; for (iColumn=0;iColumn<xxxx->numberRows_+xxxx->numberColumns_;iColumn++) if (xxxx->getStatus(iColumn)==basic) nBasic++; assert (nBasic==xxxx->numberRows_); for (iRow=0;iRow<xxxx->numberRows_;iRow++) { int iPivot=xxxx->pivotVariable_[iRow]; assert (xxxx->getStatus(iPivot)==basic); } #endif }

---

Member Function Documentation

void ClpSimplex::borrowModel (  ClpModel &  otherModel  )

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. Reimplemented from ClpModel. Definition at line 2957 of file ClpSimplex.cpp. References ClpModel::borrowModel(), createStatus(), setDualRowPivotAlgorithm(), and setPrimalColumnPivotAlgorithm(). { ClpModel::borrowModel(otherModel); createStatus(); ClpDualRowSteepest steep1; setDualRowPivotAlgorithm(steep1); ClpPrimalColumnSteepest steep2; setPrimalColumnPivotAlgorithm(steep2); }

void ClpSimplex::checkDualSolution (   )

This sets largest infeasibility and most infeasible and sum and number of infeasibilities (Dual) Definition at line 1718 of file ClpSimplex.cpp. References columnActivityWork_, columnDualInfeasibility_, columnDualSequence_, columnLowerWork_, columnPrimalInfeasibility_, columnUpperWork_, dualTolerance_, firstFree_, largestDualError_, largeValue_, ClpSimplexProgress::lastIterationNumber(), numberDualInfeasibilities_, numberDualInfeasibilitiesWithoutFree_, primalTolerance_, primalToleranceToGetOptimal_, progress_, reducedCostWork_, remainingDualInfeasibility_, rowActivityWork_, rowDualInfeasibility_, rowDualSequence_, rowLowerWork_, rowPrimalInfeasibility_, rowReducedCost_, rowUpperWork_, sumDualInfeasibilities_, and sumOfRelaxedDualInfeasibilities_. Referenced by checkSolution(), gutsOfSolution(), and ClpSimplexDual::statusOfProblemInDual(). { int iRow,iColumn; sumDualInfeasibilities_=0.0; numberDualInfeasibilities_=0; numberDualInfeasibilitiesWithoutFree_=0; columnDualInfeasibility_=0.0; columnDualSequence_=-1; rowDualInfeasibility_=0.0; rowDualSequence_=-1; int firstFreePrimal = -1; int firstFreeDual = -1; int numberSuperBasicWithDj=0; primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_, columnPrimalInfeasibility_); remainingDualInfeasibility_=0.0; double relaxedTolerance=dualTolerance_; // we can't really trust infeasibilities if there is dual error double error = min(1.0e-3,largestDualError_); // allow tolerance at least slightly bigger than standard relaxedTolerance = relaxedTolerance + error; sumOfRelaxedDualInfeasibilities_ = 0.0; for (iColumn=0;iColumn<numberColumns_;iColumn++) { if (getColumnStatus(iColumn) != basic&&!flagged(iColumn)) { // not basic double distanceUp = columnUpperWork_[iColumn]- columnActivityWork_[iColumn]; double distanceDown = columnActivityWork_[iColumn] - columnLowerWork_[iColumn]; if (distanceUp>primalTolerance_) { double value = reducedCostWork_[iColumn]; // Check if "free" if (distanceDown>primalTolerance_) { if (fabs(value)>1.0e2*relaxedTolerance) { numberSuperBasicWithDj++; if (firstFreeDual<0) firstFreeDual = iColumn; } if (firstFreePrimal<0) firstFreePrimal = iColumn; } // should not be negative if (value<0.0) { value = - value; if (value>columnDualInfeasibility_) { columnDualInfeasibility_=value; columnDualSequence_=iColumn; } if (value>dualTolerance_) { if (getColumnStatus(iColumn) != isFree) { numberDualInfeasibilitiesWithoutFree_ ++; sumDualInfeasibilities_ += value-dualTolerance_; if (value>relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance; numberDualInfeasibilities_ ++; } else { // free so relax a lot value *= 0.01; if (value>dualTolerance_) { sumDualInfeasibilities_ += value-dualTolerance_; if (value>relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance; numberDualInfeasibilities_ ++; } } // maybe we can make feasible by increasing tolerance if (distanceUp<largeValue_) { if (distanceUp>primalToleranceToGetOptimal_) primalToleranceToGetOptimal_=distanceUp; } else { //gap too big for any tolerance remainingDualInfeasibility_= max(remainingDualInfeasibility_,value); } } } } if (distanceDown>primalTolerance_) { double value = reducedCostWork_[iColumn]; // should not be positive if (value>0.0) { if (value>columnDualInfeasibility_) { columnDualInfeasibility_=value; columnDualSequence_=iColumn; } if (value>dualTolerance_) { sumDualInfeasibilities_ += value-dualTolerance_; if (value>relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance; numberDualInfeasibilities_ ++; if (getColumnStatus(iColumn) != isFree) numberDualInfeasibilitiesWithoutFree_ ++; // maybe we can make feasible by increasing tolerance if (distanceDown<largeValue_&& distanceDown>primalToleranceToGetOptimal_) primalToleranceToGetOptimal_=distanceDown; } } } } } for (iRow=0;iRow<numberRows_;iRow++) { if (getRowStatus(iRow) != basic&&!flagged(iRow+numberColumns_)) { // not basic double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow]; double distanceDown = rowActivityWork_[iRow] -rowLowerWork_[iRow]; if (distanceUp>primalTolerance_) { double value = rowReducedCost_[iRow]; // Check if "free" if (distanceDown>primalTolerance_) { if (fabs(value)>1.0e2*relaxedTolerance) { numberSuperBasicWithDj++; if (firstFreeDual<0) firstFreeDual = iRow+numberColumns_; } if (firstFreePrimal<0) firstFreePrimal = iRow+numberColumns_; } // should not be negative if (value<0.0) { value = - value; if (value>rowDualInfeasibility_) { rowDualInfeasibility_=value; rowDualSequence_=iRow; } if (value>dualTolerance_) { sumDualInfeasibilities_ += value-dualTolerance_; if (value>relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance; numberDualInfeasibilities_ ++; if (getRowStatus(iRow) != isFree) numberDualInfeasibilitiesWithoutFree_ ++; // maybe we can make feasible by increasing tolerance if (distanceUp<largeValue_) { if (distanceUp>primalToleranceToGetOptimal_) primalToleranceToGetOptimal_=distanceUp; } else { //gap too big for any tolerance remainingDualInfeasibility_= max(remainingDualInfeasibility_,value); } } } } if (distanceDown>primalTolerance_) { double value = rowReducedCost_[iRow]; // should not be positive if (value>0.0) { if (value>rowDualInfeasibility_) { rowDualInfeasibility_=value; rowDualSequence_=iRow; } if (value>dualTolerance_) { sumDualInfeasibilities_ += value-dualTolerance_; if (value>relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance; numberDualInfeasibilities_ ++; if (getRowStatus(iRow) != isFree) numberDualInfeasibilitiesWithoutFree_ ++; // maybe we can make feasible by increasing tolerance if (distanceDown<largeValue_&& distanceDown>primalToleranceToGetOptimal_) primalToleranceToGetOptimal_=distanceDown; } } } } } if (algorithm_<0&&firstFreeDual>=0) { // dual firstFree_ = firstFreeDual; } else if (numberSuperBasicWithDj|| (progress_&&progress_->lastIterationNumber(0)<=0)) { firstFree_=firstFreePrimal; } }

void ClpSimplex::checkPrimalSolution (  const double *  rowActivities = NULL, const double *  columnActivies = NULL )

This sets largest infeasibility and most infeasible and sum and number of infeasibilities (Primal) Definition at line 1645 of file ClpSimplex.cpp. References columnActivityWork_, columnLowerWork_, columnPrimalInfeasibility_, columnPrimalSequence_, columnUpperWork_, largestPrimalError_, largestSolutionError_, numberPrimalInfeasibilities_, objectiveWork_, primalTolerance_, rowActivityWork_, rowLowerWork_, rowObjectiveWork_, rowPrimalInfeasibility_, rowPrimalSequence_, rowUpperWork_, solution(), sumOfRelaxedPrimalInfeasibilities_, and sumPrimalInfeasibilities_. Referenced by checkSolution(), gutsOfSolution(), and ClpSimplexDual::strongBranching(). { double * solution; int iRow,iColumn; objectiveValue_ = 0.0; // now look at primal solution columnPrimalInfeasibility_=0.0; columnPrimalSequence_=-1; rowPrimalInfeasibility_=0.0; rowPrimalSequence_=-1; largestSolutionError_=0.0; solution = rowActivityWork_; sumPrimalInfeasibilities_=0.0; numberPrimalInfeasibilities_=0; double primalTolerance = primalTolerance_; double relaxedTolerance=primalTolerance_; // we can't really trust infeasibilities if there is primal error double error = min(1.0e-3,largestPrimalError_); // allow tolerance at least slightly bigger than standard relaxedTolerance = relaxedTolerance + error; sumOfRelaxedPrimalInfeasibilities_ = 0.0; for (iRow=0;iRow<numberRows_;iRow++) { //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic); double infeasibility=0.0; objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow]; if (solution[iRow]>rowUpperWork_[iRow]) { infeasibility=solution[iRow]-rowUpperWork_[iRow]; } else if (solution[iRow]<rowLowerWork_[iRow]) { infeasibility=rowLowerWork_[iRow]-solution[iRow]; } if (infeasibility>primalTolerance) { sumPrimalInfeasibilities_ += infeasibility-primalTolerance_; if (infeasibility>relaxedTolerance) sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance; numberPrimalInfeasibilities_ ++; } if (infeasibility>rowPrimalInfeasibility_) { rowPrimalInfeasibility_=infeasibility; rowPrimalSequence_=iRow; } infeasibility = fabs(rowActivities[iRow]-solution[iRow]); if (infeasibility>largestSolutionError_) largestSolutionError_=infeasibility; } solution = columnActivityWork_; for (iColumn=0;iColumn<numberColumns_;iColumn++) { //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic); double infeasibility=0.0; objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn]; if (solution[iColumn]>columnUpperWork_[iColumn]) { infeasibility=solution[iColumn]-columnUpperWork_[iColumn]; } else if (solution[iColumn]<columnLowerWork_[iColumn]) { infeasibility=columnLowerWork_[iColumn]-solution[iColumn]; } if (infeasibility>columnPrimalInfeasibility_) { columnPrimalInfeasibility_=infeasibility; columnPrimalSequence_=iColumn; } if (infeasibility>primalTolerance) { sumPrimalInfeasibilities_ += infeasibility-primalTolerance_; if (infeasibility>relaxedTolerance) sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance; numberPrimalInfeasibilities_ ++; } infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]); if (infeasibility>largestSolutionError_) largestSolutionError_=infeasibility; } }

void ClpSimplex::checkSolution (   )

Just check solution (for external use) - sets sum of infeasibilities etc Definition at line 3673 of file ClpSimplex.cpp. References checkDualSolution(), checkPrimalSolution(), columnActivityWork_, createRim(), deleteRim(), dj_, dualTolerance_, numberDualInfeasibilities_, numberPrimalInfeasibilities_, ClpModel::objectiveValue(), primalTolerance_, rowActivityWork_, and solution_. Referenced by initialSolve(). { // put in standard form createRim(7+8+16); dualTolerance_=dblParam_[ClpDualTolerance]; primalTolerance_=dblParam_[ClpPrimalTolerance]; checkPrimalSolution( rowActivityWork_, columnActivityWork_); checkDualSolution(); if (!numberDualInfeasibilities_&& !numberPrimalInfeasibilities_) problemStatus_=0; else problemStatus_=-1; #ifdef CLP_DEBUG int i; double value=0.0; for (i=0;i<numberRows_+numberColumns_;i++) value += dj_[i]*solution_[i]; printf("dual value %g, primal %g\n",value,objectiveValue()); #endif // release extra memory deleteRim(0); }

double ClpSimplex::columnDualInfeasibility (   )  const [inline]

Worst column dual infeasibility (note - these may not be as meaningful if the problem is primal infeasible Definition at line 474 of file ClpSimplex.hpp. References columnDualInfeasibility_. { return columnDualInfeasibility_;} ;

void ClpSimplex::computeDuals (  double *  givenDjs  )

Computes duals from scratch. If givenDjs then allows for nonzero basic djs Definition at line 492 of file ClpSimplex.cpp. References CoinIndexedVector::checkClear(), CoinIndexedVector::clear(), cost_, CoinIndexedVector::denseVector(), dj_, factorization_, CoinIndexedVector::getIndices(), largestDualError_, numberRefinements_, objectiveWork_, pivotVariable_, reducedCostWork_, CoinIndexedVector::reserve(), rowArray_, rowObjectiveWork_, rowReducedCost_, CoinIndexedVector::setNumElements(), and transposeTimes(). Referenced by gutsOfSolution(), ClpSimplexPrimal::primal(), and ClpSimplexDual::statusOfProblemInDual(). { //work space CoinIndexedVector * workSpace = rowArray_[0]; CoinIndexedVector arrayVector; arrayVector.reserve(numberRows_+1); CoinIndexedVector previousVector; previousVector.reserve(numberRows_+1); int iRow; #ifdef CLP_DEBUG workSpace->checkClear(); #endif double * array = arrayVector.denseVector(); int * index = arrayVector.getIndices(); int number=0; if (!givenDjs) { for (iRow=0;iRow<numberRows_;iRow++) { int iPivot=pivotVariable_[iRow]; double value = cost_[iPivot]; if (value) { array[iRow]=value; index[number++]=iRow; } } } else { // dual values pass - djs may not be zero for (iRow=0;iRow<numberRows_;iRow++) { int iPivot=pivotVariable_[iRow]; // make sure zero if done if (!pivoted(iPivot)) givenDjs[iPivot]=0.0; double value =cost_[iPivot]-givenDjs[iPivot]; if (value) { array[iRow]=value; index[number++]=iRow; } } } arrayVector.setNumElements(number); // Btran basic costs and get as accurate as possible double lastError=COIN_DBL_MAX; int iRefine; double * work = workSpace->denseVector(); CoinIndexedVector * thisVector = &arrayVector; CoinIndexedVector * lastVector = &previousVector; factorization_->updateColumnTranspose(workSpace,thisVector); for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) { // check basic reduced costs zero largestDualError_=0.0; // would be faster to do just for basic but this reduces code ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_); transposeTimes(-1.0,array,reducedCostWork_); if (!givenDjs) { for (iRow=0;iRow<numberRows_;iRow++) { int iPivot=pivotVariable_[iRow]; double value; if (iPivot>=numberColumns_) { // slack value = rowObjectiveWork_[iPivot-numberColumns_] + array[iPivot-numberColumns_]; } else { // column value = reducedCostWork_[iPivot]; } work[iRow]=value; if (fabs(value)>largestDualError_) { largestDualError_=fabs(value); } } } else { for (iRow=0;iRow<numberRows_;iRow++) { int iPivot=pivotVariable_[iRow]; if (iPivot>=numberColumns_) { // slack work[iRow] = rowObjectiveWork_[iPivot-numberColumns_] + array[iPivot-numberColumns_]-givenDjs[iPivot]; } else { // column work[iRow] = reducedCostWork_[iPivot]- givenDjs[iPivot]; } if (fabs(work[iRow])>largestDualError_) { largestDualError_=fabs(work[iRow]); //assert (largestDualError_<1.0e-7); //if (largestDualError_>1.0e-7) //printf("large dual error %g\n",largestDualError_); } } } if (largestDualError_>=lastError) { // restore CoinIndexedVector * temp = thisVector; thisVector = lastVector; lastVector=temp; break; } if (iRefine<numberRefinements_&&largestDualError_>1.0e-10 &&!givenDjs) { // try and make better // save this CoinIndexedVector * temp = thisVector; thisVector = lastVector; lastVector=temp; int * indexOut = thisVector->getIndices(); int number=0; array = thisVector->denseVector(); thisVector->clear(); double multiplier = 131072.0; for (iRow=0;iRow<numberRows_;iRow++) { double value = multiplier*work[iRow]; if (value) { array[iRow]=value; indexOut[number++]=iRow; work[iRow]=0.0; } work[iRow]=0.0; } thisVector->setNumElements(number); lastError=largestDualError_; factorization_->updateColumnTranspose(workSpace,thisVector); multiplier = 1.0/multiplier; double * previous = lastVector->denseVector(); number=0; for (iRow=0;iRow<numberRows_;iRow++) { double value = previous[iRow] + multiplier*array[iRow]; if (value) { array[iRow]=value; indexOut[number++]=iRow; } else { array[iRow]=0.0; } } thisVector->setNumElements(number); } else { break; } } ClpFillN(work,numberRows_,0.0); // now look at dual solution array = thisVector->denseVector(); for (iRow=0;iRow<numberRows_;iRow++) { // slack double value = array[iRow]; dual_[iRow]=value; value += rowObjectiveWork_[iRow]; rowReducedCost_[iRow]=value; } ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_); transposeTimes(-1.0,dual_,reducedCostWork_); // If necessary - override results if (givenDjs) { // restore accurate duals memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double)); } }

int ClpSimplex::crash (  double  gap, int  pivot )

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" is 0 No pivoting (so will just be choice of algorithm) 1 Simple pivoting e.g. gub 2 Mini iterations Definition at line 3711 of file ClpSimplex.cpp. References createStatus(), ClpModel::matrix(), CoinMessageHandler::message(), ClpModel::objective(), and solution(). { assert(!rowObjective_); // not coded int iColumn; int numberBad=0; int numberBasic=0; double dualTolerance=dblParam_[ClpDualTolerance]; //double primalTolerance=dblParam_[ClpPrimalTolerance]; int returnCode=0; // If no basis then make all slack one if (!status_) createStatus(); for (iColumn=0;iColumn<numberColumns_;iColumn++) { if (getColumnStatus(iColumn)==basic) numberBasic++; } if (!numberBasic) { // all slack double * dj = new double [numberColumns_]; double * solution = columnActivity_; const double * linearObjective = objective(); //double objectiveValue=0.0; int iColumn; double direction = optimizationDirection_; // direction is actually scale out not scale in if (direction) direction = 1.0/direction; for (iColumn=0;iColumn<numberColumns_;iColumn++) dj[iColumn] = direction*linearObjective[iColumn]; for (iColumn=0;iColumn<numberColumns_;iColumn++) { // assume natural place is closest to zero double lowerBound = columnLower_[iColumn]; double upperBound = columnUpper_[iColumn]; if (lowerBound>-1.0e20||upperBound<1.0e20) { bool atLower; if (fabs(upperBound)<fabs(lowerBound)) { atLower=false; setColumnStatus(iColumn,atUpperBound); solution[iColumn]=upperBound; } else { atLower=true; setColumnStatus(iColumn,atLowerBound); solution[iColumn]=lowerBound; } if (dj[iColumn]<0.0) { // should be at upper bound if (atLower) { // can we flip if (upperBound-lowerBound<=gap) { columnActivity_[iColumn]=upperBound; setColumnStatus(iColumn,atUpperBound); } else if (dj[iColumn]<-dualTolerance) { numberBad++; } } } else if (dj[iColumn]>0.0) { // should be at lower bound if (!atLower) { // can we flip if (upperBound-lowerBound<=gap) { columnActivity_[iColumn]=lowerBound; setColumnStatus(iColumn,atLowerBound); } else if (dj[iColumn]>dualTolerance) { numberBad++; } } } } else { // free setColumnStatus(iColumn,isFree); if (fabs(dj[iColumn])>dualTolerance) numberBad++; } } if (numberBad||pivot) { if (!pivot) { delete [] dj; returnCode = 1; } else { // see if can be made dual feasible with gubs etc double * pi = new double[numberRows_]; memset (pi,0,numberRows_*sizeof(double)); int * way = new int[numberColumns_]; int numberIn = 0; // Get column copy CoinPackedMatrix * columnCopy = matrix(); // Get a row copy in standard format CoinPackedMatrix copy; copy.reverseOrderedCopyOf(*columnCopy); // get matrix data pointers const int * column = copy.getIndices(); const CoinBigIndex * rowStart = copy.getVectorStarts(); const int * rowLength = copy.getVectorLengths(); const double * elementByRow = copy.getElements(); //const int * row = columnCopy->getIndices(); //const CoinBigIndex * columnStart = columnCopy->getVectorStarts(); //const int * columnLength = columnCopy->getVectorLengths(); //const double * element = columnCopy->getElements(); // if equality row and bounds mean artificial in basis bad // then do anyway for (iColumn=0;iColumn<numberColumns_;iColumn++) { // - if we want to reduce dj, + if we want to increase int thisWay = 100; double lowerBound = columnLower_[iColumn]; double upperBound = columnUpper_[iColumn]; if (upperBound>lowerBound) { switch(getColumnStatus(iColumn)) { case basic: thisWay=0; case ClpSimplex::isFixed: break; case isFree: case superBasic: if (dj[iColumn]<-dualTolerance) thisWay = 1; else if (dj[iColumn]>dualTolerance) thisWay = -1; else thisWay =0; break; case atUpperBound: if (dj[iColumn]>dualTolerance) thisWay = -1; else if (dj[iColumn]<-dualTolerance) thisWay = -3; else thisWay = -2; break; case atLowerBound: if (dj[iColumn]<-dualTolerance) thisWay = 1; else if (dj[iColumn]>dualTolerance) thisWay = 3; else thisWay = 2; break; } } way[iColumn] = thisWay; } /*if (!numberBad) printf("Was dual feasible before passes - rows %d\n", numberRows_);*/ int lastNumberIn = -100000; int numberPasses=5; while (numberIn>lastNumberIn+numberRows_/100) { lastNumberIn = numberIn; // we need to maximize chance of doing good int iRow; for (iRow=0;iRow<numberRows_;iRow++) { double lowerBound = rowLower_[iRow]; double upperBound = rowUpper_[iRow]; if (getRowStatus(iRow)==basic) { // see if we can find a column to pivot on int j; // down is amount pi can go down double maximumDown = COIN_DBL_MAX; double maximumUp = COIN_DBL_MAX; double minimumDown =0.0; double minimumUp =0.0; int iUp=-1; int iDown=-1; int iUpB=-1; int iDownB=-1; if (lowerBound<-1.0e20) maximumUp = -1.0; if (upperBound>1.0e20) maximumDown = -1.0; for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) { int iColumn = column[j]; double value = elementByRow[j]; double djValue = dj[iColumn]; /* way - -3 - okay at upper bound with negative dj -2 - marginal at upper bound with zero dj - can only decrease -1 - bad at upper bound 0 - we can never pivot on this row 1 - bad at lower bound 2 - marginal at lower bound with zero dj - can only increase 3 - okay at lower bound with positive dj 100 - fine we can just ignore */ if (way[iColumn]!=100) { switch(way[iColumn]) { case -3: if (value>0.0) { if (maximumDown*value>-djValue) { maximumDown = -djValue/value; iDown = iColumn; } } else { if (-maximumUp*value>-djValue) { maximumUp = djValue/value; iUp = iColumn; } } break; case -2: if (value>0.0) { maximumDown = 0.0; } else { maximumUp = 0.0; } break; case -1: // see if could be satisfied // dj value > 0 if (value>0.0) { maximumDown=0.0; if (maximumUp*value<djValue-dualTolerance) { maximumUp = 0.0; // would improve but not enough } else { if (minimumUp*value<djValue) { minimumUp = djValue/value; iUpB = iColumn; } } } else { maximumUp=0.0; if (-maximumDown*value<djValue-dualTolerance) { maximumDown = 0.0; // would improve but not enough } else { if (-minimumDown*value<djValue) { minimumDown = -djValue/value; iDownB = iColumn; } } } break; case 0: maximumDown = -1.0; maximumUp=-1.0; break; case 1: // see if could be satisfied // dj value < 0 if (value>0.0) { maximumUp=0.0; if (maximumDown*value<-djValue-dualTolerance) { maximumDown = 0.0; // would improve but not enough } else { if (minimumDown*value<-djValue) { minimumDown = -djValue/value; iDownB = iColumn; } } } else { maximumDown=0.0; if (-maximumUp*value<-djValue-dualTolerance) { maximumUp = 0.0; // would improve but not enough } else { if (-minimumUp*value<-djValue) { minimumUp = djValue/value; iUpB = iColumn; } } } break; case 2: if (value>0.0) { maximumUp = 0.0; } else { maximumDown = 0.0; } break; case 3: if (value>0.0) { if (maximumUp*value>djValue) { maximumUp = djValue/value; iUp = iColumn; } } else { if (-maximumDown*value>djValue) { maximumDown = -djValue/value; iDown = iColumn; } } break; default: break; } } } if (iUpB>=0) iUp=iUpB; if (maximumUp<=dualTolerance||maximumUp<minimumUp) iUp=-1; if (iDownB>=0) iDown=iDownB; if (maximumDown<=dualTolerance||maximumDown<minimumDown) iDown=-1; if (iUp>=0||iDown>=0) { // do something if (iUp>=0&&iDown>=0) { if (maximumDown>maximumUp) iUp=-1; } double change; int kColumn; if (iUp>=0) { kColumn=iUp; change=maximumUp; // just do minimum if was dual infeasible // ? only if maximum large? if (minimumUp>0.0) change=minimumUp; setRowStatus(iRow,atUpperBound); } else { kColumn=iDown; change=-maximumDown; // just do minimum if was dual infeasible // ? only if maximum large? if (minimumDown>0.0) change=-minimumDown; setRowStatus(iRow,atLowerBound); } assert (fabs(change)<1.0e20); setColumnStatus(kColumn,basic); numberIn++; pi[iRow]=change; for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) { int iColumn = column[j]; double value = elementByRow[j]; double djValue = dj[iColumn]-change*value; dj[iColumn]=djValue; if (abs(way[iColumn])==1) { numberBad--; /*if (!numberBad) printf("Became dual feasible at row %d out of %d\n", iRow, numberRows_);*/ lastNumberIn=-1000000; } int thisWay = 100; double lowerBound = columnLower_[iColumn]; double upperBound = columnUpper_[iColumn]; if (upperBound>lowerBound) { switch(getColumnStatus(iColumn)) { case basic: thisWay=0; case isFixed: break; case isFree: case superBasic: if (djValue<-dualTolerance) thisWay = 1; else if (djValue>dualTolerance) thisWay = -1; else { thisWay =0; abort();} break; case atUpperBound: if (djValue>dualTolerance) { thisWay =-1; abort();} else if (djValue<-dualTolerance) thisWay = -3; else thisWay = -2; break; case atLowerBound: if (djValue<-dualTolerance) { thisWay =1; abort();} else if (djValue>dualTolerance) thisWay = 3; else thisWay = 2; break; } } way[iColumn] = thisWay; } } } } if (numberIn==lastNumberIn||numberBad||pivot<2) break; if (!(--numberPasses)) break; //printf("%d put in so far\n",numberIn); } // last attempt to flip for (iColumn=0;iColumn<numberColumns_;iColumn++) { double lowerBound = columnLower_[iColumn]; double upperBound = columnUpper_[iColumn]; if (upperBound-lowerBound<=gap&&upperBound>lowerBound) { double djValue=dj[iColumn]; switch(getColumnStatus(iColumn)) { case basic: case ClpSimplex::isFixed: break; case isFree: case superBasic: break; case atUpperBound: if (djValue>dualTolerance) { setColumnStatus(iColumn,atUpperBound); solution[iColumn]=upperBound; } break; case atLowerBound: if (djValue<-dualTolerance) { setColumnStatus(iColumn,atUpperBound); solution[iColumn]=upperBound; } break; } } } delete [] pi; delete [] dj; delete [] way; handler_->message(CLP_CRASH,messages_) <<numberIn <<numberBad <<CoinMessageEol; returnCode = -1; } } else { delete [] dj; returnCode = -1; } //cleanStatus(); } return returnCode; }

int ClpSimplex::createPiecewiseLinearCosts (  const int *  starts, const double *  lower, const double *  gradient )

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 Definition at line 4577 of file ClpSimplex.cpp. References nonLinearCost_, and specialOptions_. { delete nonLinearCost_; // Set up feasible bounds and check monotonicity int iColumn; int returnCode=0; for (iColumn=0;iColumn<numberColumns_;iColumn++) { int iIndex = starts[iColumn]; int end = starts[iColumn+1]-1; columnLower_[iColumn] = lower[iIndex]; columnUpper_[iColumn] = lower[end]; double value = columnLower_[iColumn]; iIndex++; for (;iIndex<end;iIndex++) { if (lower[iIndex]<value) returnCode++; // not monotonic value=lower[iIndex]; } } nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient); specialOptions_ |= 2; // say keep return returnCode; }

bool ClpSimplex::createRim (  int  what, bool  makeRowCopy = false )  [protected]

puts in format I like (rowLower,rowUpper) also see StandardMatrix 1 bit does rows, 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 Definition at line 1961 of file ClpSimplex.cpp. References ClpMatrixBase::allElementsInRange(), columnActivityWork_, columnArray_, columnLowerWork_, columnScale_, columnUpperWork_, cost_, dj_, factorization_, ClpPackedMatrix::getElements(), ClpPackedMatrix::getIndices(), ClpPackedMatrix::getVectorStarts(), CoinMessageHandler::logLevel(), lower_, ClpModel::objective(), objectiveWork_, pivotVariable_, reducedCostWork_, CoinIndexedVector::reserve(), ClpMatrixBase::reverseOrderedCopy(), rowActivityWork_, rowArray_, rowLowerWork_, rowObjectiveWork_, rowReducedCost_, rowScale_, rowUpperWork_, sanityCheck(), ClpMatrixBase::scale(), scalingFlag_, CoinMessageHandler::setLogLevel(), solution_, and upper_. Referenced by ClpSimplexDual::changeBounds(), checkSolution(), ClpSimplex(), factorize(), getSolution(), ClpSimplexPrimal::primal(), startup(), statusOfProblem(), ClpSimplexDual::statusOfProblemInDual(), ClpSimplexPrimalQuadratic::statusOfProblemInPrimal(), ClpSimplexPrimal::statusOfProblemInPrimal(), ClpSimplexDual::strongBranching(), ClpSimplexPrimal::unPerturb(), and ClpSimplexPrimal::whileIterating(). { bool goodMatrix=true; int saveLevel=handler_->logLevel(); if (problemStatus_==10) { handler_->setLogLevel(0); // switch off messages } else if (factorization_) { // match up factorization messages if (handler_->logLevel()<3) factorization_->messageLevel(0); else factorization_->messageLevel(3); } int i; if ((what&(16+32))!=0) { // move information to work arrays double direction = optimizationDirection_; // direction is actually scale out not scale in if (direction) direction = 1.0/direction; if (direction!=1.0) { // reverse all dual signs for (i=0;i<numberColumns_;i++) reducedCost_[i] *= direction; for (i=0;i<numberRows_;i++) dual_[i] *= direction; } // row reduced costs if (!dj_) { dj_ = new double[numberRows_+numberColumns_]; reducedCostWork_ = dj_; rowReducedCost_ = dj_+numberColumns_; memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double)); memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double)); } if (!solution_||(what&32)!=0) { if (!solution_) solution_ = new double[numberRows_+numberColumns_]; columnActivityWork_ = solution_; rowActivityWork_ = solution_+numberColumns_; memcpy(columnActivityWork_,columnActivity_, numberColumns_*sizeof(double)); memcpy(rowActivityWork_,rowActivity_, numberRows_*sizeof(double)); } } if ((what&16)!=0) { //check matrix if (!matrix_) matrix_=new ClpPackedMatrix(); if (!matrix_->allElementsInRange(this,smallElement_,1.0e20)) { problemStatus_=4; goodMatrix= false; } if (makeRowCopy) { delete rowCopy_; // may return NULL if can't give row copy rowCopy_ = matrix_->reverseOrderedCopy(); #ifdef TAKEOUT { ClpPackedMatrix* rowCopy = dynamic_cast< ClpPackedMatrix*>(rowCopy_); const int * column = rowCopy->getIndices(); const CoinBigIndex * rowStart = rowCopy->getVectorStarts(); const double * element = rowCopy->getElements(); int i; for (i=133;i<numberRows_;i++) { if (rowStart[i+1]-rowStart[i]==10||rowStart[i+1]-rowStart[i]==15) printf("Row %d has %d elements\n",i,rowStart[i+1]-rowStart[i]); } } #endif } } if ((what&4)!=0) { delete [] cost_; // extra copy with original costs int nTotal = numberRows_+numberColumns_; //cost_ = new double[2*nTotal]; cost_ = new double[nTotal]; objectiveWork_ = cost_; rowObjectiveWork_ = cost_+numberColumns_; memcpy(objectiveWork_,objective(),numberColumns_*sizeof(double)); if (rowObjective_) memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double)); else memset(rowObjectiveWork_,0,numberRows_*sizeof(double)); // and initialize changes to zero //memset(cost_+nTotal,0,nTotal*sizeof(double)); } if ((what&1)!=0) { delete [] lower_; delete [] upper_; lower_ = new double[numberColumns_+numberRows_]; upper_ = new double[numberColumns_+numberRows_]; rowLowerWork_ = lower_+numberColumns_; columnLowerWork_ = lower_; rowUpperWork_ = upper_+numberColumns_; columnUpperWork_ = upper_; memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double)); memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double)); memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double)); memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double)); // clean up any mismatches on infinity for (i=0;i<numberColumns_;i++) { if (columnLowerWork_[i]<-1.0e30) columnLowerWork_[i] = -COIN_DBL_MAX; if (columnUpperWork_[i]>1.0e30) columnUpperWork_[i] = COIN_DBL_MAX; } // clean up any mismatches on infinity for (i=0;i<numberRows_;i++) { if (rowLowerWork_[i]<-1.0e30) rowLowerWork_[i] = -COIN_DBL_MAX; if (rowUpperWork_[i]>1.0e30) rowUpperWork_[i] = COIN_DBL_MAX; } } // do scaling if needed if (scalingFlag_>0&&!rowScale_) { if (matrix_->scale(this)) scalingFlag_=-scalingFlag_; // not scaled after all } if ((what&4)!=0) { double direction = optimizationDirection_; // direction is actually scale out not scale in if (direction) direction = 1.0/direction; // but also scale by scale factors // not really sure about row scaling if (!rowScale_) { if (direction!=1.0) { for (i=0;i<numberRows_;i++) rowObjectiveWork_[i] *= direction; for (i=0;i<numberColumns_;i++) objectiveWork_[i] *= direction; } } else { for (i=0;i<numberRows_;i++) rowObjectiveWork_[i] *= direction/rowScale_[i]; for (i=0;i<numberColumns_;i++) objectiveWork_[i] *= direction*columnScale_[i]; } } if ((what&(1+32))!=0&&rowScale_) { for (i=0;i<numberColumns_;i++) { double multiplier = 1.0/columnScale_[i]; if (columnLowerWork_[i]>-1.0e50) columnLowerWork_[i] *= multiplier; if (columnUpperWork_[i]<1.0e50) columnUpperWork_[i] *= multiplier; } for (i=0;i<numberRows_;i++) { if (rowLowerWork_[i]>-1.0e50) rowLowerWork_[i] *= rowScale_[i]; if (rowUpperWork_[i]<1.0e50) rowUpperWork_[i] *= rowScale_[i]; } } if ((what&(8+32))!=0&&rowScale_) { // on entry for (i=0;i<numberColumns_;i++) { columnActivityWork_[i] /= columnScale_[i]; reducedCostWork_[i] *= columnScale_[i]; } for (i=0;i<numberRows_;i++) { rowActivityWork_[i] *= rowScale_[i]; dual_[i] /= rowScale_[i]; rowReducedCost_[i] = dual_[i]; } } if ((what&16)!=0) { // check rim of problem okay if (!sanityCheck()) goodMatrix=false; } // we need to treat matrix as if each element by rowScaleIn and columnScaleout?? // maybe we need to move scales to SimplexModel for factorization? if ((what&8)!=0&&!pivotVariable_) { pivotVariable_=new int[numberRows_]; } if ((what&16)!=0&&!rowArray_[2]) { // get some arrays int iRow,iColumn; // these are "indexed" arrays so we always know where nonzeros are /********************************************************** rowArray_[3] is long enough for rows+columns *********************************************************/ for (iRow=0;iRow<4;iRow++) { if (!rowArray_[iRow]) { rowArray_[iRow]=new CoinIndexedVector(); int length =numberRows_+factorization_->maximumPivots(); if (iRow==3) length += numberColumns_; rowArray_[iRow]->reserve(length); } } for (iColumn=0;iColumn<2;iColumn++) { if (!columnArray_[iColumn]) { columnArray_[iColumn]=new CoinIndexedVector(); if (!iColumn) columnArray_[iColumn]->reserve(numberColumns_); else columnArray_[iColumn]->reserve(max(numberRows_,numberColumns_)); } } } double primalTolerance=dblParam_[ClpPrimalTolerance]; if ((what&1)!=0) { // fix any variables with tiny gaps for (i=0;i<numberColumns_;i++) { if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) { if (columnLowerWork_[i]>=0.0) { columnUpperWork_[i] = columnLowerWork_[i]; } else if (columnUpperWork_[i]<=0.0) { columnLowerWork_[i] = columnUpperWork_[i]; } else { columnUpperWork_[i] = 0.0; columnLowerWork_[i] = 0.0; } } } for (i=0;i<numberRows_;i++) { if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) { if (rowLowerWork_[i]>=0.0) { rowUpperWork_[i] = rowLowerWork_[i]; } else if (rowUpperWork_[i]<=0.0) { rowLowerWork_[i] = rowUpperWork_[i]; } else { rowUpperWork_[i] = 0.0; rowLowerWork_[i] = 0.0; } } } } if (problemStatus_==10) { problemStatus_=-1; handler_->setLogLevel(saveLevel); // switch back messages } return goodMatrix; }

void ClpSimplex::createStatus (   )

Set up status array (can be used by OsiClp). Also can be used to set up all slack basis Definition at line 3568 of file ClpSimplex.cpp. Referenced by borrowModel(), cleanStatus(), crash(), loadProblem(), and readMps(). { if(!status_) status_ = new unsigned char [numberColumns_+numberRows_]; memset(status_,0,(numberColumns_+numberRows_)*sizeof(char)); int i; // set column status to one nearest zero for (i=0;i<numberColumns_;i++) { #if 0 if (columnLower_[i]>=0.0) { setColumnStatus(i,atLowerBound); } else if (columnUpper_[i]<=0.0) { setColumnStatus(i,atUpperBound); } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) { // free setColumnStatus(i,isFree); } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) { setColumnStatus(i,atLowerBound); } else { setColumnStatus(i,atUpperBound); } #else setColumnStatus(i,atLowerBound); #endif } for (i=0;i<numberRows_;i++) { setRowStatus(i,basic); } }

void ClpSimplex::deleteRim (  int  getRidOfFactorizationData = 2  )  [protected]

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 Definition at line 2208 of file ClpSimplex.cpp. References columnActivityWork_, columnLowerWork_, columnScale_, columnUpperWork_, dualTolerance_, ClpModel::gutsOfDelete(), primalTolerance_, reducedCostWork_, rowActivityWork_, rowLowerWork_, rowObjectiveWork_, rowScale_, rowUpperWork_, and scalingFlag_. Referenced by checkSolution(), factorize(), getSolution(), statusOfProblem(), and ClpSimplexDual::strongBranching(). { int i; if (problemStatus_!=1&&problemStatus_!=2) { delete [] ray_; ray_=NULL; } // ray may be null if in branch and bound if (rowScale_) { // Collect infeasibilities int numberPrimalScaled=0; int numberPrimalUnscaled=0; int numberDualScaled=0; int numberDualUnscaled=0; for (i=0;i<numberColumns_;i++) { double scaleFactor = columnScale_[i]; double valueScaled = columnActivityWork_[i]; if (valueScaled<columnLowerWork_[i]-primalTolerance_) numberPrimalScaled++; else if (valueScaled>columnUpperWork_[i]+primalTolerance_) numberPrimalScaled++; columnActivity_[i] = valueScaled*scaleFactor; double value = columnActivity_[i]; if (value<columnLower_[i]-primalTolerance_) numberPrimalUnscaled++; else if (value>columnUpper_[i]+primalTolerance_) numberPrimalUnscaled++; double valueScaledDual = reducedCostWork_[i]; if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_) numberDualScaled++; if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_) numberDualScaled++; reducedCost_[i] = valueScaledDual/scaleFactor; double valueDual = reducedCost_[i]; if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_) numberDualUnscaled++; if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_) numberDualUnscaled++; } for (i=0;i<numberRows_;i++) { double scaleFactor = rowScale_[i]; double valueScaled = rowActivityWork_[i]; if (valueScaled<rowLowerWork_[i]-primalTolerance_) numberPrimalScaled++; else if (valueScaled>rowUpperWork_[i]+primalTolerance_) numberPrimalScaled++; rowActivity_[i] = valueScaled/scaleFactor; double value = rowActivity_[i]; if (value<rowLower_[i]-primalTolerance_) numberPrimalUnscaled++; else if (value>rowUpper_[i]+primalTolerance_) numberPrimalUnscaled++; double valueScaledDual = dual_[i]+rowObjectiveWork_[i];; if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_) numberDualScaled++; if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_) numberDualScaled++; dual_[i] = valueScaledDual*scaleFactor; double valueDual = dual_[i]; if (rowObjective_) valueDual += rowObjective_[i]; if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_) numberDualUnscaled++; if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_) numberDualUnscaled++; } if (!problemStatus_&&!secondaryStatus_) { // See if we need to set secondary status if (numberPrimalUnscaled) { if (numberDualUnscaled) secondaryStatus_=4; else secondaryStatus_=2; } else { if (numberDualUnscaled) secondaryStatus_=3; } } if (problemStatus_==2) { for (i=0;i<numberColumns_;i++) { ray_[i] *= columnScale_[i]; } } else if (problemStatus_==1&&ray_) { for (i=0;i<numberRows_;i++) { ray_[i] *= rowScale_[i]; } } } else { for (i=0;i<numberColumns_;i++) { columnActivity_[i] = columnActivityWork_[i]; reducedCost_[i] = reducedCostWork_[i]; } for (i=0;i<numberRows_;i++) { rowActivity_[i] = rowActivityWork_[i]; } } if (optimizationDirection_!=1.0) { // and modify all dual signs for (i=0;i<numberColumns_;i++) reducedCost_[i] *= optimizationDirection_; for (i=0;i<numberRows_;i++) dual_[i] *= optimizationDirection_; } // scaling may have been turned off scalingFlag_ = abs(scalingFlag_); if(getRidOfFactorizationData>=0) gutsOfDelete(getRidOfFactorizationData+1); }

int ClpSimplex::directionIn (   )  const [inline]

Return direction In or Out Definition at line 618 of file ClpSimplex.hpp. References directionIn_. {return directionIn_;};

int ClpSimplex::dual (  int  ifValuesPass = 0  )

Dual algorithm - see ClpSimplexDual.hpp for method Reimplemented in ClpSimplexDual. Definition at line 2880 of file ClpSimplex.cpp. References perturbation_, and setInitialDenseFactorization(). { assert (ifValuesPass>=0&&ifValuesPass<2); int returnCode = ((ClpSimplexDual *) this)->dual(ifValuesPass); if (problemStatus_==10) { //printf("Cleaning up with primal\n"); int savePerturbation = perturbation_; perturbation_=100; bool denseFactorization = initialDenseFactorization(); // It will be safe to allow dense setInitialDenseFactorization(true); returnCode = ((ClpSimplexPrimal *) this)->primal(1); setInitialDenseFactorization(denseFactorization); perturbation_=savePerturbation; if (problemStatus_==10) problemStatus_=0; } return returnCode; }

int ClpSimplex::dualPivotResult (   )

Pivot out a variable and choose an incoing one. Assumes dual feasible - will not go through a reduced cost. Returns step length in theta Returns ray in ray_ (or NULL if no pivot) Return codes as before but -1 means no acceptable pivot Definition at line 4348 of file ClpSimplex.cpp. { return ((ClpSimplexDual *) this)->pivotResult(); }

int ClpSimplex::getSolution (   )

Given an existing factorization computes and checks primal and dual solutions. Uses current problem arrays for bounds. Returns feasibility states Definition at line 671 of file ClpSimplex.cpp. References columnActivityWork_, and rowActivityWork_. { double * rowActivities = new double[numberRows_]; double * columnActivities = new double[numberColumns_]; ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities); ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities); int status = getSolution( rowActivities, columnActivities); delete [] rowActivities; delete [] columnActivities; return status; }

int ClpSimplex::getSolution (  const double *  rowActivities, const double *  columnActivities )

Given an existing factorization computes and checks primal and dual solutions. Uses input arrays for variables at bounds. Returns feasibility states Definition at line 655 of file ClpSimplex.cpp. References createRim(), deleteRim(), factorization_, and gutsOfSolution(). { if (!factorization_->status()) { // put in standard form createRim(7+8+16+32); // do work gutsOfSolution ( NULL,NULL); // release extra memory deleteRim(0); } return factorization_->status(); }

int ClpSimplex::gutsOfSolution (  double *  givenDuals, const double *  givenPrimals, bool  valuesPass = false )  [protected]

May change basis and then returns number changed. Computation of solutions may be overriden by given pi and solution Definition at line 244 of file ClpSimplex.cpp. References algorithm_, ClpNonLinearCost::changeInCost(), checkDualSolution(), ClpNonLinearCost::checkInfeasibilities(), checkPrimalSolution(), columnActivityWork_, computeDuals(), computePrimals(), firstFree_, incomingInfeasibility_, largestDualError_, ClpNonLinearCost::largestInfeasibility(), largestPrimalError_, CoinMessageHandler::logLevel(), CoinMessageHandler::message(), nonLinearCost_, ClpNonLinearCost::numberInfeasibilities(), pivotVariable_, primalTolerance_, rowActivityWork_, solution_, specialOptions_, and times(). Referenced by ClpSimplexDual::dual(), ClpSimplexDual::fastDual(), getSolution(), pivot(), ClpSimplexPrimal::primal(), startup(), statusOfProblem(), ClpSimplexDual::statusOfProblemInDual(), ClpSimplexPrimalQuadratic::statusOfProblemInPrimal(), ClpSimplexPrimal::statusOfProblemInPrimal(), ClpSimplexPrimal::unPerturb(), ClpSimplexPrimal::whileIterating(), and ClpSimplexDual::whileIterating(). { // if values pass, save values of basic variables double * save = NULL; double oldValue=0.0; if (valuesPass) { assert(algorithm_>0); // only primal at present assert(nonLinearCost_); int iRow; checkPrimalSolution( rowActivityWork_, columnActivityWork_); // get correct bounds on all variables nonLinearCost_->checkInfeasibilities(primalTolerance_); oldValue = nonLinearCost_->largestInfeasibility(); save = new double[numberRows_]; for (iRow=0;iRow<numberRows_;iRow++) { int iPivot=pivotVariable_[iRow]; save[iRow] = solution_[iPivot]; } } // do work computePrimals(rowActivityWork_, columnActivityWork_); // If necessary - override results if (givenPrimals) { memcpy(columnActivityWork_,givenPrimals,numberColumns_*sizeof(double)); memset(rowActivityWork_,0,numberRows_*sizeof(double)); times(-1.0,columnActivityWork_,rowActivityWork_); } double objectiveModification = 0.0; if (algorithm_>0&&nonLinearCost_!=NULL) { // primal algorithm // get correct bounds on all variables // If 4 bit set - Force outgoing variables to exact bound (primal) if ((specialOptions_&4)==0) nonLinearCost_->checkInfeasibilities(primalTolerance_); else nonLinearCost_->checkInfeasibilities(0.0); objectiveModification += nonLinearCost_->changeInCost(); if (nonLinearCost_->numberInfeasibilities()) handler_->message(CLP_SIMPLEX_NONLINEAR,messages_) <<nonLinearCost_->changeInCost() <<nonLinearCost_->numberInfeasibilities() <<CoinMessageEol; } if (valuesPass) { #ifdef CLP_DEBUG std::cout<<"Largest given infeasibility "<<oldValue <<" now "<<nonLinearCost_->largestInfeasibility()<<std::endl; #endif int numberOut=0; if (oldValue<incomingInfeasibility_ &&nonLinearCost_->largestInfeasibility()> max(incomingInfeasibility_,allowedInfeasibility_)|| largestPrimalError_>1.0e-3) { printf("Original largest infeas %g, now %g, primalError %g\n", oldValue,nonLinearCost_->largestInfeasibility(), largestPrimalError_); // throw out up to 1000 structurals int iRow; int * sort = new int[numberRows_]; // first put back solution and store difference for (iRow=0;iRow<numberRows_;iRow++) { int iPivot=pivotVariable_[iRow]; double difference = fabs(solution_[iPivot]-save[iRow]); solution_[iPivot]=save[iRow]; save[iRow]=difference; } for (iRow=0;iRow<numberRows_;iRow++) { int iPivot=pivotVariable_[iRow]; if (iPivot<numberColumns_) { // column double difference= save[iRow]; if (difference>1.0e-4) { sort[numberOut]=iPivot; save[numberOut++]=difference; } } } CoinSort_2(save, save + numberOut, sort, CoinFirstGreater_2<double, int>()); numberOut = min(1000,numberOut); for (iRow=0;iRow<numberOut;iRow++) { int iColumn=sort[iRow]; setColumnStatus(iColumn,superBasic); } delete [] sort; } delete [] save; if (numberOut) return numberOut; } computeDuals(givenDuals); // now check solutions checkPrimalSolution( rowActivityWork_, columnActivityWork_); objectiveValue_ += objectiveModification; checkDualSolution(); if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2|| largestDualError_>1.0e-2)) handler_->message(CLP_SIMPLEX_ACCURACY,messages_) <<largestPrimalError_ <<largestDualError_ <<CoinMessageEol; // Switch off false values pass indicator if (!valuesPass&&algorithm_>0) firstFree_ = -1; return 0; }

int ClpSimplex::housekeeping (  double  objectiveChange  )

This does basis housekeeping and does values for in/out variables. Can also decide to re-factorize Definition at line 1142 of file ClpSimplex.cpp. References algorithm_, alpha_, changeMade_, ClpSimplexProgress::cycle(), directionIn_, directionOut_, dualIn_, dualOut_, factorization_, forceFactorization_, ClpModel::hitMaximumIterations(), isColumn(), lower_, CoinMessageHandler::message(), ClpModel::objectiveValue(), pivotRow_, pivotVariable_, CoinMessageHandler::printing(), progress_, progressFlag_, sequenceIn(), sequenceIn_, sequenceOut_, sequenceWithin(), setFlagged(), solution_, ClpSimplexProgress::startCheck(), theta_, upper_, valueIn_, and valueOut_. Referenced by pivot(), ClpSimplexPrimal::pivotResult(), ClpSimplexPrimalQuadratic::whileIterating(), and ClpSimplexDual::whileIterating(). { numberIterations_++; changeMade_++; // something has happened // incoming variable handler_->message(CLP_SIMPLEX_HOUSE1,messages_) <<directionOut_ <<directionIn_<<theta_ <<dualOut_<<dualIn_<<alpha_ <<CoinMessageEol; if (getStatus(sequenceIn_)==isFree) { handler_->message(CLP_SIMPLEX_FREEIN,messages_) <<sequenceIn_ <<CoinMessageEol; } // change of incoming char rowcol[]={'R','C'}; if (pivotRow_>=0) pivotVariable_[pivotRow_]=sequenceIn(); if (upper_[sequenceIn_]>1.0e20&&lower_[sequenceIn_]<-1.0e20) progressFlag_ |= 2; // making real progress solution_[sequenceIn_]=valueIn_; if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12) progressFlag_ |= 1; // making real progress if (sequenceIn_!=sequenceOut_) { //assert( getStatus(sequenceOut_)== basic); setStatus(sequenceIn_,basic); if (upper_[sequenceOut_]-lower_[sequenceOut_]>0) { // As Nonlinear costs may have moved bounds (to more feasible) // Redo using value if (fabs(valueOut_-lower_[sequenceOut_])<fabs(valueOut_-upper_[sequenceOut_])) { // going to lower setStatus(sequenceOut_,atLowerBound); } else { // going to upper setStatus(sequenceOut_,atUpperBound); } } else { // fixed setStatus(sequenceOut_,isFixed); } solution_[sequenceOut_]=valueOut_; } else { // flip from bound to bound // As Nonlinear costs may have moved bounds (to more feasible) // Redo using value if (fabs(valueIn_-lower_[sequenceIn_])<fabs(valueIn_-upper_[sequenceIn_])) { // as if from upper bound setStatus(sequenceIn_, atLowerBound); } else { // as if from lower bound setStatus(sequenceIn_, atUpperBound); } } objectiveValue_ += objectiveChange; handler_->message(CLP_SIMPLEX_HOUSE2,messages_) <<numberIterations_<<objectiveValue() <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_) <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_); handler_->printing(algorithm_<0)<<theta_<<dualOut_; handler_->printing(algorithm_>0)<<dualIn_<<theta_; handler_->message()<<CoinMessageEol; if (hitMaximumIterations()) return 2; #if 0 // check for small cycles int cycle=progress_->cycle(sequenceIn_,sequenceOut_, directionIn_,directionOut_); if (cycle>0) { printf("Cycle of %d\n",cycle); // reset progress_->startCheck(); if (factorization_->pivots()>cycle) { forceFactorization_=cycle-1; } else { // need to reject something int iSequence; if (algorithm_<0) iSequence = sequenceIn_; else iSequence = sequenceOut_; char x = isColumn(iSequence) ? 'C' :'R'; handler_->message(CLP_SIMPLEX_FLAG,messages_) <<x<<sequenceWithin(iSequence) <<CoinMessageEol; setFlagged(iSequence); printf("flagging %d\n",iSequence); } return 1; } #endif // only time to re-factorize if one before real time // this is so user won't be surprised that maximumPivots has exact meaning if (factorization_->pivots()==factorization_->maximumPivots()) { return 1; } else { if (forceFactorization_>0&& factorization_->pivots()==forceFactorization_) { // relax forceFactorization_ = (3+5*forceFactorization_)/4; if (forceFactorization_>factorization_->maximumPivots()) forceFactorization_ = -1; //off return 1; } else { // carry on iterating return 0; } } }

int ClpSimplex::initialSolve (  ClpSolve &  options  )

General solve algorithm which can do presolve. See ClpSolve.hpp for options Definition at line 53 of file ClpSolve.cpp. References checkSolution(), Idiot::crash(), ClpModel::deleteColumns(), ClpModel::dualColumnSolution(), ClpModel::dualRowSolution(), dualTolerance_, Idiot::getDropEnoughFeasibility(), ClpSolve::getExtraInfo(), ClpPlusMinusOneMatrix::getIndices(), ClpModel::getNumElements(), ClpSolve::getPresolvePasses(), ClpSolve::getPresolveType(), ClpSolve::getSolveType(), ClpSolve::getSpecialOption(), Idiot::getStrategy(), ClpPackedMatrix::matrix(), ClpModel::maximumIterations(), CoinMessageHandler::message(), ClpModel::numberColumns_, ClpModel::numberIterations(), ClpModel::numberRows_, ClpModel::objectiveValue(), perturbation(), perturbation_, primal(), ClpModel::primalColumnSolution(), ClpInterior::primalDual(), ClpModel::primalRowSolution(), primalTolerance_, CoinMessageHandler::printing(), ClpModel::rowLower(), ClpModel::rowUpper(), scalingFlag_, ClpModel::setDblParam(), Idiot::setDropEnoughFeasibility(), Idiot::setDropEnoughWeighted(), Idiot::setLightweight(), ClpModel::setLogLevel(), setPerturbation(), Idiot::setReduceIterations(), Idiot::setStartingWeight(), Idiot::setStrategy(), solution(), and ClpModel::status(). { ClpSolve::SolveType method=options.getSolveType(); ClpSolve::PresolveType presolve = options.getPresolveType(); int saveMaxIterations = maximumIterations(); int finalStatus=-1; int numberIterations=0; double time1 = CoinCpuTime(); double timeX = time1; double time2; ClpMatrixBase * saveMatrix=NULL; ClpSimplex * model2 = this; bool interrupt = (options.getSpecialOption(2)==0); CoinSighandler_t saveSignal=SIG_DFL; if (interrupt) { currentModel = model2; // register signal handler saveSignal = signal(SIGINT,signal_handler); } ClpPresolve pinfo; double timePresolve=0.0; double timeIdiot=0.0; double timeCore=0.0; int savePerturbation=perturbation_; int saveScaling = scalingFlag_; if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) { // network - switch off stuff presolve = ClpSolve::presolveOff; } // For below >0 overrides // 0 means no, -1 means maybe int doIdiot=0; int doCrash=0; int doSprint=0; if (method!=ClpSolve::useDual) { switch (options.getSpecialOption(1)) { case 0: doIdiot=-1; doCrash=-1; doSprint=-1; break; case 1: doIdiot=0; doCrash=1; doSprint=0; break; case 2: doIdiot=1; doCrash=0; doSprint=0; break; case 3: doIdiot=0; doCrash=0; doSprint=1; break; case 4: doIdiot=0; doCrash=0; doSprint=0; break; case 5: doIdiot=0; doCrash=-1; doSprint=-1; break; case 6: doIdiot=-1; doCrash=-1; doSprint=0; break; case 7: doIdiot=-1; doCrash=0; doSprint=-1; break; case 8: doIdiot=-1; doCrash=0; doSprint=0; break; case 9: doIdiot=0; doCrash=0; doSprint=-1; break; default: abort(); } } else { // Dual switch (options.getSpecialOption(0)) { case 0: doIdiot=0; doCrash=0; doSprint=0; break; case 1: doIdiot=0; doCrash=1; doSprint=0; break; case 2: doIdiot=-1; if (options.getExtraInfo(0)>0) doIdiot = options.getExtraInfo(0); doCrash=0; doSprint=0; break; default: abort(); } } // Just do this number of passes in Sprint int maxSprintPass=100; // PreSolve to file - not fully tested bool presolveToFile=false; if (presolve!=ClpSolve::presolveOff) { int numberPasses=5; if (presolve==ClpSolve::presolveNumber) { numberPasses=options.getPresolvePasses(); presolve = ClpSolve::presolveOn; } if (presolveToFile) { printf("***** temp test\n"); pinfo.presolvedModelToFile(*this,"ss.sav",1.0e-8, false,numberPasses); model2=this; } else { model2 = pinfo.presolvedModel(*this,1.0e-8, false,numberPasses,true); } time2 = CoinCpuTime(); timePresolve = time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Presolve"<<timePresolve<<time2-time1 <<CoinMessageEol; timeX=time2; if (model2) { } else { handler_->message(CLP_INFEASIBLE,messages_) <<CoinMessageEol; model2 = this; presolve=ClpSolve::presolveOff; } // We may be better off using original (but if dual leave because of bounds) if (presolve!=ClpSolve::presolveOff&& numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_ &&method!=ClpSolve::useDual&&model2!=this) { delete model2; model2 = this; presolve=ClpSolve::presolveOff; } } if (interrupt) currentModel = model2; // See if worth trying +- one matrix bool plusMinus=false; int numberElements=model2->getNumElements(); if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) { // network - switch off stuff doIdiot=0; doSprint=0; } int numberColumns = model2->numberColumns(); int numberRows = model2->numberRows(); // If not all slack basis - switch off all int number=0; int iRow; for (iRow=0;iRow<numberRows;iRow++) if (model2->getRowStatus(iRow)==basic) number++; if (number<numberRows) { doIdiot=0; doCrash=0; doSprint=0; } if (options.getSpecialOption(3)==0) { if(numberElements>100000) plusMinus=true; if(numberElements>10000&&(doIdiot||doSprint)) plusMinus=true; } if (plusMinus) { saveMatrix = model2->clpMatrix(); ClpPackedMatrix* clpMatrix = dynamic_cast< ClpPackedMatrix*>(saveMatrix); if (clpMatrix) { ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix())); if (newMatrix->getIndices()) { model2->replaceMatrix(newMatrix); } else { handler_->message(CLP_MATRIX_CHANGE,messages_) <<"+- 1" <<CoinMessageEol; saveMatrix=NULL; plusMinus=false; delete newMatrix; } } else { saveMatrix=NULL; plusMinus=false; } } if (model2->factorizationFrequency()==200) { // User did not touch preset model2->setFactorizationFrequency(min(2000,100+model2->numberRows()/200)); } if (method==ClpSolve::usePrimalorSprint) { if (doSprint<0) { if (numberElements<500000) { // Small problem if(numberRows*10>numberColumns||numberColumns<6000 ||(numberRows*20>numberColumns&&!plusMinus)) method=ClpSolve::usePrimal; // switch off sprint } else { // larger problem if(numberRows*8>numberColumns) method=ClpSolve::usePrimal; // switch off sprint // but make lightweight if(numberRows*10>numberColumns||numberColumns<6000 ||(numberRows*20>numberColumns&&!plusMinus)) maxSprintPass=5; } } else if (doSprint==0) { method=ClpSolve::usePrimal; // switch off sprint } } if (method==ClpSolve::useDual) { // pick up number passes int nPasses=0; int numberNotE=0; if ((doIdiot<0&&plusMinus)||doIdiot>0) { // See if candidate for idiot nPasses=0; Idiot info(*model2); // Get average number of elements per column double ratio = ((double) numberElements/(double) numberColumns); // look at rhs int iRow; double largest=0.0; double smallest = 1.0e30; double largestGap=0.0; for (iRow=0;iRow<numberRows;iRow++) { double value1 = model2->rowLower_[iRow]; if (value1&&value1>-1.0e31) { largest = max(largest,fabs(value1)); smallest=min(smallest,fabs(value1)); } double value2 = model2->rowUpper_[iRow]; if (value2&&value2<1.0e31) { largest = max(largest,fabs(value2)); smallest=min(smallest,fabs(value2)); } if (value2>value1) { numberNotE++; if (value2>1.0e31||value1<-1.0e31) largestGap = COIN_DBL_MAX; else largestGap = value2-value1; } } if (doIdiot<0) { if (numberRows>200&&numberColumns>5000&&ratio>=3.0&& largest/smallest<1.1&&!numberNotE) { nPasses = 71; } } if (doIdiot>0) { nPasses=max(nPasses,doIdiot); if (nPasses>70) info.setStartingWeight(1.0e3); } if (nPasses) { info.setReduceIterations(5); doCrash=0; info.crash(nPasses,model2->messageHandler(),model2->messagesPointer()); time2 = CoinCpuTime(); timeIdiot = time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Idiot Crash"<<timeIdiot<<time2-time1 <<CoinMessageEol; timeX=time2; } } if (presolve==ClpSolve::presolveOn) { int numberInfeasibilities = model2->tightenPrimalBounds(); if (numberInfeasibilities) { handler_->message(CLP_INFEASIBLE,messages_) <<CoinMessageEol; model2 = this; presolve=ClpSolve::presolveOff; } } if (options.getSpecialOption(0)==1) model2->crash(1000,1); if (!nPasses) { model2->dual(0); } else if (!numberNotE&&0) { // E so we can do in another way double * pi = model2->dualRowSolution(); int i; int numberColumns = model2->numberColumns(); int numberRows = model2->numberRows(); double * saveObj = new double[numberColumns]; memcpy(saveObj,model2->objective(),numberColumns*sizeof(double)); memcpy(model2->dualColumnSolution(),model2->objective(), numberColumns*sizeof(double)); model2->clpMatrix()->transposeTimes(-1.0,pi,model2->dualColumnSolution()); memcpy(model2->objective(),model2->dualColumnSolution(), numberColumns*sizeof(double)); const double * rowsol = model2->primalRowSolution(); double offset=0.0; for (i=0;i<numberRows;i++) { offset += pi[i]*rowsol[i]; } double value2; model2->getDblParam(ClpObjOffset,value2); printf("Offset %g %g\n",offset,value2); model2->setRowObjective(pi); // zero out pi memset(pi,0,numberRows*sizeof(double)); // Could put some in basis - only partially tested model2->allSlackBasis(); model2->factorization()->maximumPivots(200); //model2->setLogLevel(63); // solve model2->dual(1); memcpy(model2->objective(),saveObj,numberColumns*sizeof(double)); // zero out pi memset(pi,0,numberRows*sizeof(double)); model2->setRowObjective(pi); delete [] saveObj; model2->primal(); } else { // solve model2->dual(1); } time2 = CoinCpuTime(); timeCore = time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Dual"<<timeCore<<time2-time1 <<CoinMessageEol; timeX=time2; } else if (method==ClpSolve::usePrimal) { if (doIdiot) { int nPasses=0; Idiot info(*model2); // Get average number of elements per column double ratio = ((double) numberElements/(double) numberColumns); // look at rhs int iRow; double largest=0.0; double smallest = 1.0e30; double largestGap=0.0; int numberNotE=0; for (iRow=0;iRow<numberRows;iRow++) { double value1 = model2->rowLower_[iRow]; if (value1&&value1>-1.0e31) { largest = max(largest,fabs(value1)); smallest=min(smallest,fabs(value1)); } double value2 = model2->rowUpper_[iRow]; if (value2&&value2<1.0e31) { largest = max(largest,fabs(value2)); smallest=min(smallest,fabs(value2)); } if (value2>value1) { numberNotE++; if (value2>1.0e31||value1<-1.0e31) largestGap = COIN_DBL_MAX; else largestGap = value2-value1; } } if (numberRows>200&&numberColumns>5000&&numberColumns>2*numberRows) { if (plusMinus) { if (largest/smallest>2.0) { nPasses = 10+numberColumns/100000; nPasses = min(nPasses,50); nPasses = max(nPasses,15); if (numberRows>25000&&nPasses>5) { // Might as well go for it nPasses = max(nPasses,71); } else if (numberElements<3*numberColumns) { nPasses=min(nPasses,10); // probably not worh it } if (doIdiot<0) info.setLightweight(1); // say lightweight idiot } else if (largest/smallest>1.01||numberElements<=3*numberColumns) { nPasses = 10+numberColumns/1000; nPasses = min(nPasses,100); nPasses = max(nPasses,30); if (numberRows>25000) { // Might as well go for it nPasses = max(nPasses,71); } if (!largestGap) nPasses *= 2; } else { nPasses = 10+numberColumns/1000; nPasses = min(nPasses,200); nPasses = max(nPasses,100); info.setStartingWeight(1.0e-1); info.setReduceIterations(6); if (!largestGap) nPasses *= 2; //info.setFeasibilityTolerance(1.0e-7); } // If few passes - don't bother if (nPasses<=5) nPasses=0; } else { if (doIdiot<0) info.setLightweight(1); // say lightweight idiot if (largest/smallest>1.01||numberNotE) { if (numberRows>25000||numberColumns>5*numberRows) { nPasses = 50; } else if (numberColumns>4*numberRows) { nPasses = 20; } else { nPasses=5; } } else { if (numberRows>25000||numberColumns>5*numberRows) { nPasses = 50; info.setLightweight(0); // say not lightweight idiot } else if (numberColumns>4*numberRows) { nPasses = 20; } else { nPasses=15; } } if (numberElements<3*numberColumns) { nPasses=(int) ((2.0*(double) nPasses)/ratio); // probably not worh it } else { nPasses = max(nPasses,5); nPasses = (int) (((double) nPasses)*5.0/ratio); // reduce if lots of elements per column } if (numberRows>25000&&nPasses>5) { // Might as well go for it nPasses = max(nPasses,71); } else if (plusMinus) { nPasses *= 2; nPasses=min(nPasses,71); } if (nPasses<=5) nPasses=0; //info.setStartingWeight(1.0e-1); } } if (doIdiot>0) { // pick up number passes nPasses=options.getExtraInfo(1); if (nPasses>70) { info.setStartingWeight(1.0e3); info.setReduceIterations(6); // also be more lenient on infeasibilities info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility()); info.setDropEnoughWeighted(-2.0); } else if (nPasses>=50) { info.setStartingWeight(1.0e3); //info.setReduceIterations(6); } // For experimenting if (nPasses<70&&(nPasses%10)>0&&(nPasses%10)<4) { info.setStartingWeight(1.0e3); info.setLightweight(nPasses%10); // special testing //info.setReduceIterations(6); } } if (nPasses) { doCrash=0; #if 0 double * solution = model2->primalColumnSolution(); int iColumn; double * saveLower = new double[numberColumns]; memcpy(saveLower,model2->columnLower(),numberColumns*sizeof(double)); double * saveUpper = new double[numberColumns]; memcpy(saveUpper,model2->columnUpper(),numberColumns*sizeof(double)); printf("doing tighten before idiot\n"); model2->tightenPrimalBounds(); // Move solution double * columnLower = model2->columnLower(); double * columnUpper = model2->columnUpper(); for (iColumn=0;iColumn<numberColumns;iColumn++) { if (columnLower[iColumn]>0.0) solution[iColumn]=columnLower[iColumn]; else if (columnUpper[iColumn]<0.0) solution[iColumn]=columnUpper[iColumn]; else solution[iColumn]=0.0; } memcpy(columnLower,saveLower,numberColumns*sizeof(double)); memcpy(columnUpper,saveUpper,numberColumns*sizeof(double)); delete [] saveLower; delete [] saveUpper; #else // Allow for crossover info.setStrategy(512|info.getStrategy()); // Allow for scaling info.setStrategy(32|info.getStrategy()); info.crash(nPasses,model2->messageHandler(),model2->messagesPointer()); #endif time2 = CoinCpuTime(); timeIdiot = time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Idiot Crash"<<timeIdiot<<time2-time1 <<CoinMessageEol; timeX=time2; } } // ? if (doCrash) model2->crash(1000,1); model2->primal(1); time2 = CoinCpuTime(); timeCore = time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Primal"<<timeCore<<time2-time1 <<CoinMessageEol; timeX=time2; } else if (method==ClpSolve::usePrimalorSprint) { // Sprint /* This driver implements what I called Sprint when I introduced the idea many years ago. Cplex calls it "sifting" which I think is just as silly. When I thought of this trivial idea it reminded me of an LP code of the 60's called sprint which after every factorization took a subset of the matrix into memory (all 64K words!) and then iterated very fast on that subset. On the problems of those days it did not work very well, but it worked very well on aircrew scheduling problems where there were very large numbers of columns all with the same flavor. */ /* The idea works best if you can get feasible easily. To make it more general we can add in costed slacks */ int originalNumberColumns = model2->numberColumns(); int numberRows = model2->numberRows(); // We will need arrays to choose variables. These are too big but .. double * weight = new double [numberRows+originalNumberColumns]; int * sort = new int [numberRows+originalNumberColumns]; int numberSort=0; // We are going to add slacks to get feasible. // initial list will just be artificials // first we will set all variables as close to zero as possible int iColumn; const double * columnLower = model2->columnLower(); const double * columnUpper = model2->columnUpper(); double * columnSolution = model2->primalColumnSolution(); for (iColumn=0;iColumn<originalNumberColumns;iColumn++) { double value =0.0; if (columnLower[iColumn]>0.0) value = columnLower[iColumn]; else if (columnUpper[iColumn]<0.0) value = columnUpper[iColumn]; columnSolution[iColumn]=value; } // now see what that does to row solution double * rowSolution = model2->primalRowSolution(); memset (rowSolution,0,numberRows*sizeof(double)); model2->times(1.0,columnSolution,rowSolution); int * addStarts = new int [numberRows+1]; int * addRow = new int[numberRows]; double * addElement = new double[numberRows]; const double * lower = model2->rowLower(); const double * upper = model2->rowUpper(); addStarts[0]=0; int numberArtificials=0; double * addCost = new double [numberRows]; const double penalty=1.0e8; int iRow; for (iRow=0;iRow<numberRows;iRow++) { if (lower[iRow]>rowSolution[iRow]) { addRow[numberArtificials]=iRow; addElement[numberArtificials]=1.0; addCost[numberArtificials]=penalty; numberArtificials++; addStarts[numberArtificials]=numberArtificials; } else if (upper[iRow]<rowSolution[iRow]) { addRow[numberArtificials]=iRow; addElement[numberArtificials]=-1.0; addCost[numberArtificials]=penalty; numberArtificials++; addStarts[numberArtificials]=numberArtificials; } } model2->addColumns(numberArtificials,NULL,NULL,addCost, addStarts,addRow,addElement); delete [] addStarts; delete [] addRow; delete [] addElement; delete [] addCost; // look at rhs to see if to perturb double largest=0.0; double smallest = 1.0e30; for (iRow=0;iRow<numberRows;iRow++) { double value; value = fabs(model2->rowLower_[iRow]); if (value&&value<1.0e30) { largest = max(largest,value); smallest=min(smallest,value); } value = fabs(model2->rowUpper_[iRow]); if (value&&value<1.0e30) { largest = max(largest,value); smallest=min(smallest,value); } } double * saveLower = NULL; double * saveUpper = NULL; if (largest<2.01*smallest) { // perturb - so switch off standard model2->setPerturbation(100); saveLower = new double[numberRows]; memcpy(saveLower,model2->rowLower_,numberRows*sizeof(double)); saveUpper = new double[numberRows]; memcpy(saveUpper,model2->rowUpper_,numberRows*sizeof(double)); double * lower = model2->rowLower(); double * upper = model2->rowUpper(); for (iRow=0;iRow<numberRows;iRow++) { double lowerValue=lower[iRow], upperValue=upper[iRow]; double value = CoinDrand48(); if (upperValue>lowerValue+primalTolerance_) { if (lowerValue>-1.0e20&&lowerValue) lowerValue -= value * 1.0e-4*fabs(lowerValue); if (upperValue<1.0e20&&upperValue) upperValue += value * 1.0e-4*fabs(upperValue); } else if (upperValue>0.0) { upperValue -= value * 1.0e-4*fabs(lowerValue); lowerValue -= value * 1.0e-4*fabs(lowerValue); } else if (upperValue<0.0) { upperValue += value * 1.0e-4*fabs(lowerValue); lowerValue += value * 1.0e-4*fabs(lowerValue); } else { } lower[iRow]=lowerValue; upper[iRow]=upperValue; } } int i; // Set up initial list if (numberArtificials) { numberSort=numberArtificials; for (i=0;i<numberSort;i++) sort[i] = i+originalNumberColumns; } else { numberSort = min(numberRows_,numberColumns_); for (i=0;i<numberSort;i++) sort[i] = i; } // redo as will have changed columnLower = model2->columnLower(); columnUpper = model2->columnUpper(); int numberColumns = model2->numberColumns(); double * fullSolution = model2->primalColumnSolution(); // Just do this number of passes in Sprint if (doSprint>0) maxSprintPass=options.getExtraInfo(1); int iPass; double lastObjective=1.0e31; // It will be safe to allow dense model2->setInitialDenseFactorization(true); // Just take this number of columns in small problem int smallNumberColumns = min(3*numberRows,numberColumns); smallNumberColumns = max(smallNumberColumns,3000); //int smallNumberColumns = min(12*numberRows/10,numberColumns); //smallNumberColumns = max(smallNumberColumns,3000); //smallNumberColumns = max(smallNumberColumns,numberRows+1000); // We will be using all rows int * whichRows = new int [numberRows]; for (iRow=0;iRow<numberRows;iRow++) whichRows[iRow]=iRow; double originalOffset; model2->getDblParam(ClpObjOffset,originalOffset); int totalIterations=0; for (iPass=0;iPass<maxSprintPass;iPass++) { //printf("Bug until submodel new version\n"); //CoinSort_2(sort,sort+numberSort,weight); // Create small problem ClpSimplex small(model2,numberRows,whichRows,numberSort,sort); small.setPerturbation(model2->perturbation()); // now see what variables left out do to row solution double * rowSolution = model2->primalRowSolution(); double * sumFixed = new double[numberRows]; memset (sumFixed,0,numberRows*sizeof(double)); int iRow,iColumn; // zero out ones in small problem for (iColumn=0;iColumn<numberSort;iColumn++) { int kColumn = sort[iColumn]; fullSolution[kColumn]=0.0; } // Get objective offset double offset=0.0; const double * objective = model2->objective(); for (iColumn=0;iColumn<numberColumns;iColumn++) offset += fullSolution[iColumn]*objective[iColumn]; small.setDblParam(ClpObjOffset,originalOffset-offset); model2->times(1.0,fullSolution,sumFixed); double * lower = small.rowLower(); double * upper = small.rowUpper(); for (iRow=0;iRow<numberRows;iRow++) { if (lower[iRow]>-1.0e50) lower[iRow] -= sumFixed[iRow]; if (upper[iRow]<1.0e50) upper[iRow] -= sumFixed[iRow]; rowSolution[iRow] -= sumFixed[iRow]; } delete [] sumFixed; // Solve if (interrupt) currentModel = &small; small.primal(); totalIterations += small.numberIterations(); // move solution back const double * solution = small.primalColumnSolution(); for (iColumn=0;iColumn<numberSort;iColumn++) { int kColumn = sort[iColumn]; model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn)); fullSolution[kColumn]=solution[iColumn]; } for (iRow=0;iRow<numberRows;iRow++) model2->setRowStatus(iRow,small.getRowStatus(iRow)); memcpy(model2->primalRowSolution(),small.primalRowSolution(), numberRows*sizeof(double)); // get reduced cost for large problem memcpy(weight,model2->objective(),numberColumns*sizeof(double)); model2->transposeTimes(-1.0,small.dualRowSolution(),weight); int numberNegative=0; double sumNegative = 0.0; // now massage weight so all basic in plus good djs for (iColumn=0;iColumn<numberColumns;iColumn++) { double dj = weight[iColumn]*optimizationDirection_; double value = fullSolution[iColumn]; if (model2->getColumnStatus(iColumn)==ClpSimplex::basic) dj = -1.0e50; else if (dj<0.0&&value<columnUpper[iColumn]) dj = dj; else if (dj>0.0&&value>columnLower[iColumn]) dj = -dj; else if (columnUpper[iColumn]>columnLower[iColumn]) dj = fabs(dj); else dj = 1.0e50; weight[iColumn] = dj; if (dj<-dualTolerance_&&dj>-1.0e50) { numberNegative++; sumNegative -= dj; } sort[iColumn] = iColumn; } handler_->message(CLP_SPRINT,messages_) <<iPass+1<<small.numberIterations()<<small.objectiveValue()<<sumNegative <<numberNegative <<CoinMessageEol; if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5)|| !small.numberIterations()|| iPass==maxSprintPass-1||small.status()==3) { break; // finished } else { lastObjective = small.objectiveValue()*optimizationDirection_; // sort CoinSort_2(weight,weight+numberColumns,sort); numberSort = smallNumberColumns; } } if (interrupt) currentModel = model2; for (i=0;i<numberArtificials;i++) sort[i] = i + originalNumberColumns; model2->deleteColumns(numberArtificials,sort); delete [] weight; delete [] sort; delete [] whichRows; if (saveLower) { // unperturb and clean for (iRow=0;iRow<numberRows;iRow++) { double diffLower = saveLower[iRow]-model2->rowLower_[iRow]; double diffUpper = saveUpper[iRow]-model2->rowUpper_[iRow]; model2->rowLower_[iRow]=saveLower[iRow]; model2->rowUpper_[iRow]=saveUpper[iRow]; if (diffLower) assert (!diffUpper||fabs(diffLower-diffUpper)<1.0e-5); else diffLower = diffUpper; model2->rowActivity_[iRow] += diffLower; } delete [] saveLower; delete [] saveUpper; } model2->primal(0); model2->setPerturbation(savePerturbation); time2 = CoinCpuTime(); timeCore = time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Sprint"<<timeCore<<time2-time1 <<CoinMessageEol; timeX=time2; model2->setNumberIterations(model2->numberIterations()+totalIterations); } else if (method==ClpSolve::useBarrier) { printf("***** experimental pretty crude barrier\n"); ClpInterior barrier(*model2); barrier.primalDual(); time2 = CoinCpuTime(); timeCore = time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Barrier"<<timeCore<<time2-time1 <<CoinMessageEol; timeX=time2; printf("***** crude crossover - I will improve\n"); // move solutions CoinMemcpyN(barrier.primalRowSolution(), model2->numberRows(),model2->primalRowSolution()); CoinMemcpyN(barrier.dualRowSolution(), model2->numberRows(),model2->dualRowSolution()); CoinMemcpyN(barrier.primalColumnSolution(), model2->numberColumns(),model2->primalColumnSolution()); CoinMemcpyN(barrier.dualColumnSolution(), model2->numberColumns(),model2->dualColumnSolution()); // make sure no status left model2->createStatus(); // solve model2->setPerturbation(100); model2->primal(1); model2->setPerturbation(savePerturbation); time2 = CoinCpuTime(); timeCore = time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Crossover"<<timeCore<<time2-time1 <<CoinMessageEol; timeX=time2; } else { assert (method!=ClpSolve::automatic); // later time2=0.0; } if (saveMatrix) { // delete and replace delete model2->clpMatrix(); model2->replaceMatrix(saveMatrix); } numberIterations = model2->numberIterations(); finalStatus=model2->status(); if (presolve==ClpSolve::presolveOn) { int saveLevel = logLevel(); setLogLevel(1); pinfo.postsolve(true); time2 = CoinCpuTime(); timePresolve += time2-timeX; handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Postsolve"<<time2-timeX<<time2-time1 <<CoinMessageEol; timeX=time2; if (!presolveToFile) delete model2; if (interrupt) currentModel = this; checkSolution(); setLogLevel(saveLevel); if (finalStatus!=3&&(finalStatus||status()==-1)) { int savePerturbation = perturbation(); setPerturbation(100); primal(1); setPerturbation(savePerturbation); numberIterations += this->numberIterations(); finalStatus=status(); time2 = CoinCpuTime(); handler_->message(CLP_INTERVAL_TIMING,messages_) <<"Cleanup"<<time2-timeX<<time2-time1 <<CoinMessageEol; timeX=time2; } } setMaximumIterations(saveMaxIterations); std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped"}; assert (finalStatus>=-1&&finalStatus<=3); handler_->message(CLP_TIMING,messages_) <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1; handler_->printing(presolve==ClpSolve::presolveOn) <<timePresolve; handler_->printing(timeIdiot) <<timeIdiot; handler_->message()<<CoinMessageEol; if (interrupt) signal(SIGINT,saveSignal); perturbation_=savePerturbation; scalingFlag_=saveScaling; return finalStatus; }

void ClpSimplex::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 )

Just like the other loadProblem() method except that the matrix is given in a standard column major ordered format (without gaps). Reimplemented from ClpModel. Definition at line 3634 of file ClpSimplex.cpp. References createStatus(), and ClpModel::loadProblem(). { ClpModel::loadProblem(numcols, numrows, start, index, value, collb, colub, obj, rowlb, rowub, rowObjective); createStatus(); }

void ClpSimplex::loadProblem (  const ClpMatrixBase &  matrix, const double *  collb, const double *  colub, const double *  obj, const double *  rowlb, const double *  rowub, const double *  rowObjective = NULL )

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: colub: all columns have upper bound infinity collb: all columns have lower bound 0 rowub: all rows have upper bound infinity rowlb: all rows have lower bound -infinity obj: all variables have 0 objective coefficient Reimplemented from ClpModel. Definition at line 3609 of file ClpSimplex.cpp. References createStatus(), and ClpModel::loadProblem(). Referenced by ClpSimplexPrimalQuadratic::makeQuadratic(). { ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub, rowObjective); createStatus(); }

void ClpSimplex::originalModel (  ClpSimplex *  miniModel  )

This copies back stuff from miniModel and then deletes miniModel. Only to be used with mini constructor Definition at line 5152 of file ClpSimplex.cpp. References ClpNonLinearCost::checkInfeasibilities(), columnActivityWork_, columnLowerWork_, columnScale_, columnUpperWork_, cost_, dj_, ClpModel::getDblParam(), getStatus(), lower_, ClpModel::matrix_, nonLinearCost_, ClpModel::numberColumns_, ClpNonLinearCost::numberInfeasibilities(), ClpModel::numberRows_, objectiveWork_, pivotVariable_, primalColumnPivot_, reducedCostWork_, rowActivityWork_, ClpModel::rowCopy_, rowLowerWork_, rowObjectiveWork_, rowReducedCost_, rowScale_, ClpModel::rowUpper_, rowUpperWork_, savedSolution_, saveStatus_, ClpPrimalColumnPivot::saveWeights(), ClpModel::setDblParam(), solution_, ClpModel::status_, ClpNonLinearCost::sumInfeasibilities(), ClpMatrixBase::times(), and upper_. Referenced by ClpSimplexPrimal::primal(). { int numberSmall = numberColumns_; numberColumns_ = miniModel->numberColumns_; int numberTotal = numberSmall+numberRows_; // copy back int iColumn; int * mapping = (int *) miniModel->rowUpper_; #ifdef CLP_DEBUG for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) printf("mapb %d %d\n",iColumn,mapping[iColumn]); #endif // miniModel actually has full arrays // now see what variables left out do to row solution double * fullSolution = miniModel->solution_; double * sumFixed = new double[numberRows_]; memset (sumFixed,0,numberRows_*sizeof(double)); miniModel->matrix_->times(1.0,fullSolution,sumFixed,rowScale_,miniModel->columnScale_); for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; miniModel->lower_[jColumn]=lower_[iColumn]; miniModel->upper_[jColumn]=upper_[iColumn]; miniModel->cost_[jColumn]=cost_[iColumn]; miniModel->dj_[jColumn]=dj_[iColumn]; miniModel->solution_[jColumn]=solution_[iColumn]; miniModel->status_[jColumn]=status_[iColumn]; #ifdef CLP_DEBUG printf("%d in small -> %d in original\n",iColumn,jColumn); #endif } delete [] lower_; lower_ = miniModel->lower_; delete [] upper_; upper_ = miniModel->upper_; delete [] cost_; cost_ = miniModel->cost_; delete [] dj_; dj_ = miniModel->dj_; delete [] solution_; solution_ = miniModel->solution_; delete [] status_; status_ = miniModel->status_; if (columnScale_) { for (iColumn=0;iColumn<numberSmall;iColumn++) { int jColumn = mapping[iColumn]; miniModel->columnScale_[jColumn]=columnScale_[iColumn]; } delete [] columnScale_; columnScale_ = miniModel->columnScale_; } if (savedSolution_) { if (!miniModel->savedSolution_) { miniModel->savedSolution_ = ClpCopyOfArray(solution_,numberColumns_+numberRows_); } else { for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; miniModel->savedSolution_[jColumn]=savedSolution_[iColumn]; } } delete [] savedSolution_; savedSolution_ = miniModel->savedSolution_; } if (saveStatus_) { if (!miniModel->saveStatus_) { miniModel->saveStatus_ = ClpCopyOfArray(status_,numberColumns_+numberRows_); } else { for (iColumn=0;iColumn<numberTotal;iColumn++) { int jColumn = mapping[iColumn]; miniModel->saveStatus_[jColumn]=saveStatus_[iColumn]; } } delete [] saveStatus_; saveStatus_ = miniModel->saveStatus_; } // Re-define pivotVariable_ int iRow; for (iRow=0;iRow<numberRows_;iRow++) { int iPivot = pivotVariable_[iRow]; #ifdef CLP_DEBUG printf("pb row %d, pivot %d -> %d\n",iRow,iPivot,mapping[iPivot]); #endif pivotVariable_[iRow]=mapping[iPivot]; assert (pivotVariable_[iRow]>=0); } // delete stuff and move back delete matrix_; delete rowCopy_; delete primalColumnPivot_; delete nonLinearCost_; matrix_ = miniModel->matrix_; rowCopy_ = miniModel->rowCopy_; nonLinearCost_ = miniModel->nonLinearCost_; double originalOffset; miniModel->getDblParam(ClpObjOffset,originalOffset); setDblParam(ClpObjOffset,originalOffset); // Redo some stuff reducedCostWork_ = dj_; rowReducedCost_ = dj_+numberColumns_; columnActivityWork_ = solution_; rowActivityWork_ = solution_+numberColumns_; objectiveWork_ = cost_; rowObjectiveWork_ = cost_+numberColumns_; rowLowerWork_ = lower_+numberColumns_; columnLowerWork_ = lower_; rowUpperWork_ = upper_+numberColumns_; columnUpperWork_ = upper_; // Cleanup for (iRow=0;iRow<numberRows_;iRow++) { double value = rowActivityWork_[iRow] + sumFixed[iRow]; rowActivityWork_[iRow] = value; switch(getRowStatus(iRow)) { case basic: break; case atUpperBound: //rowActivityWork_[iRow]=rowUpperWork_[iRow]; break; case ClpSimplex::isFixed: case atLowerBound: //rowActivityWork_[iRow]=rowLowerWork_[iRow]; break; case isFree: break; // superbasic case superBasic: break; } } delete [] sumFixed; nonLinearCost_->checkInfeasibilities(); printf("in original %d infeasibilities summing to %g\n", nonLinearCost_->numberInfeasibilities(),nonLinearCost_->sumInfeasibilities()); // Initialize weights primalColumnPivot_ = new ClpPrimalColumnSteepest(10); primalColumnPivot_->saveWeights(this,2); #ifndef NDEBUG // Check status ClpSimplex * xxxx = this; int nBasic=0; for (iColumn=0;iColumn<xxxx->numberRows_+xxxx->numberColumns_;iColumn++) if (xxxx->getStatus(iColumn)==basic) nBasic++; assert (nBasic==xxxx->numberRows_); for (iRow=0;iRow<xxxx->numberRows_;iRow++) { int iPivot=xxxx->pivotVariable_[iRow]; assert (xxxx->getStatus(iPivot)==basic); } #endif }

int ClpSimplex::perturbation (   )  const [inline]

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 Definition at line 301 of file ClpSimplex.hpp. References perturbation_. Referenced by initialSolve(), and ClpItem::intParameter(). { return perturbation_;};

int ClpSimplex::pivot (   )

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. Definition at line 4153 of file ClpSimplex.cpp. References algorithm_, alpha_, CoinIndexedVector::clear(), columnArray_, CoinIndexedVector::denseVector(), directionIn_, directionOut_, dj_, dualIn_, dualOut_, factorization_, CoinIndexedVector::getIndices(), gutsOfSolution(), housekeeping(), lastGoodIteration_, lower_, lowerIn_, lowerOut_, pivotRow_, pivotVariable_, rowArray_, scalingFlag_, sequenceIn_, sequenceOut_, solution_, statusOfProblem(), ClpMatrixBase::transposeTimes(), unpack(), upper_, upperIn_, upperOut_, valueIn_, and valueOut_. { // scaling not allowed assert (!scalingFlag_); // assume In_ and Out_ are correct and directionOut_ set // (or In_ if flip lowerIn_ = lower_[sequenceIn_]; valueIn_ = solution_[sequenceIn_]; upperIn_ = upper_[sequenceIn_]; dualIn_ = dj_[sequenceIn_]; if (sequenceOut_>=0&&sequenceIn_!=sequenceIn_) { assert (pivotRow_>=0&&pivotRow_<numberRows_); assert (pivotVariable_[pivotRow_]==sequenceOut_); lowerOut_ = lower_[sequenceOut_]; valueOut_ = solution_[sequenceOut_]; upperOut_ = upper_[sequenceOut_]; // for now assume primal is feasible (or in dual) dualOut_ = dj_[sequenceOut_]; assert(fabs(dualOut_)<1.0e-6); } else { assert (pivotRow_<0); } bool roundAgain = true; int returnCode=0; while (roundAgain) { roundAgain=false; unpack(rowArray_[1]); factorization_->updateColumnFT(rowArray_[2],rowArray_[1]); // we are going to subtract movement from current basic double movement; // see where incoming will go to if (sequenceOut_<0||sequenceIn_==sequenceOut_) { // flip so go to bound movement = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_; } else { // get where outgoing needs to get to double outValue = (directionOut_>0) ? upperOut_ : lowerOut_; // solutionOut_ - movement*alpha_ == outValue movement = (outValue-valueOut_)/alpha_; // set directionIn_ correctly directionIn_ = (movement>0) ? 1 :-1; } // update primal solution { int i; int * index = rowArray_[1]->getIndices(); int number = rowArray_[1]->getNumElements(); double * element = rowArray_[1]->denseVector(); for (i=0;i<number;i++) { int ii = index[i]; // get column ii = pivotVariable_[ii]; solution_[ii] -= movement*element[i]; element[i]=0.0; } // see where something went to if (sequenceOut_<0) { if (directionIn_<0) { assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7); solution_[sequenceIn_]=upperIn_; } else { assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7); solution_[sequenceIn_]=lowerIn_; } } else { if (directionOut_<0) { assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7); solution_[sequenceOut_]=upperOut_; } else { assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7); solution_[sequenceOut_]=lowerOut_; } solution_[sequenceIn_]=valueIn_+movement; } } double objectiveChange = dualIn_*movement; // update duals if (pivotRow_>=0) { alpha_ = rowArray_[1]->denseVector()[pivotRow_]; assert (fabs(alpha_)>1.0e-8); double multiplier = dualIn_/alpha_; rowArray_[0]->insert(pivotRow_,multiplier); factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]); // put row of tableau in rowArray[0] and columnArray[0] matrix_->transposeTimes(this,-1.0, rowArray_[0],columnArray_[1],columnArray_[0]); // update column djs int i; int * index = columnArray_[0]->getIndices(); int number = columnArray_[0]->getNumElements(); double * element = columnArray_[0]->denseVector(); for (i=0;i<number;i++) { int ii = index[i]; dj_[ii] += element[ii]; element[ii]=0.0; } columnArray_[0]->setNumElements(0); // and row djs index = rowArray_[0]->getIndices(); number = rowArray_[0]->getNumElements(); element = rowArray_[0]->denseVector(); for (i=0;i<number;i++) { int ii = index[i]; dj_[ii+numberColumns_] += element[ii]; dual_[ii] = dj_[ii+numberColumns_]; element[ii]=0.0; } rowArray_[0]->setNumElements(0); // check incoming assert (fabs(dj_[sequenceIn_])<1.0e-6); } // if stable replace in basis int updateStatus = factorization_->replaceColumn(rowArray_[2], pivotRow_, alpha_); bool takePivot=true; // if no pivots, bad update but reasonable alpha - take and invert if (updateStatus==2&& lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5) updateStatus=4; if (updateStatus==1||updateStatus==4) { // slight error if (factorization_->pivots()>5||updateStatus==4) { returnCode=-1; } } else if (updateStatus==2) { // major error rowArray_[1]->clear(); takePivot=false; if (factorization_->pivots()) { // refactorize here statusOfProblem(); roundAgain=true; } else { returnCode=1; } } else if (updateStatus==3) { // out of memory // increase space if not many iterations if (factorization_->pivots()< 0.5*factorization_->maximumPivots()&& factorization_->pivots()<200) factorization_->areaFactor( factorization_->areaFactor() * 1.1); returnCode =-1; // factorize now } if (takePivot) { int save = algorithm_; // make simple so always primal algorithm_=1; housekeeping(objectiveChange); algorithm_=save; } } if (returnCode == -1) { // refactorize here statusOfProblem(); } else { // just for now - recompute anyway gutsOfSolution(NULL,NULL); } return returnCode; }

int ClpSimplex::primal (  int  ifValuesPass = 0  )

Primal algorithm - see ClpSimplexPrimal.hpp for method Reimplemented in ClpSimplexPrimal. Definition at line 2900 of file ClpSimplex.cpp. References perturbation_, and setInitialDenseFactorization(). Referenced by initialSolve(). { assert (ifValuesPass>=0&&ifValuesPass<2); int returnCode = ((ClpSimplexPrimal *) this)->primal(ifValuesPass); if (problemStatus_==10) { //printf("Cleaning up with dual\n"); int savePerturbation = perturbation_; perturbation_=100; bool denseFactorization = initialDenseFactorization(); // It will be safe to allow dense setInitialDenseFactorization(true); returnCode = ((ClpSimplexDual *) this)->dual(0); setInitialDenseFactorization(denseFactorization); perturbation_=savePerturbation; if (problemStatus_==10) problemStatus_=0; } return returnCode; }

int ClpSimplex::primalPivotResult (   )

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 Definition at line 4323 of file ClpSimplex.cpp. References dj_, dualIn_, lower_, lowerIn_, sequenceIn_, solution_, upper_, upperIn_, and valueIn_. { assert (sequenceIn_>=0); valueIn_=solution_[sequenceIn_]; lowerIn_=lower_[sequenceIn_]; upperIn_=upper_[sequenceIn_]; dualIn_=dj_[sequenceIn_]; int returnCode = ((ClpSimplexPrimal *) this)->pivotResult(); if (returnCode<0&&returnCode>-4) { return 0; } else { printf("Return code of %d from ClpSimplexPrimal::pivotResult\n", returnCode); return -1; } }

int ClpSimplex::quadraticSLP (  int  numberPasses, double  deltaTolerance )

Solves quadratic 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 Definition at line 2924 of file ClpSimplex.cpp. { return ((ClpSimplexPrimalQuadratic *) this)->primalSLP(numberPasses,deltaTolerance); }

int ClpSimplex::restoreModel (  const char *  fileName  )

Restore model from file, returns 0 if success, deletes current model Definition at line 3197 of file ClpSimplex.cpp. References algorithm_, columnArray_, dualBound_, dualRowPivot_, dualTolerance_, factorization_, ClpModel::gutsOfDelete(), infeasibilityCost_, numberDualInfeasibilities_, numberDualInfeasibilitiesWithoutFree_, numberPrimalInfeasibilities_, numberRefinements_, primalColumnPivot_, primalTolerance_, rowArray_, scalingFlag_, specialOptions_, sumDualInfeasibilities_, and sumPrimalInfeasibilities_. { FILE * fp = fopen(fileName,"rb"); if (fp) { // Get rid of current model ClpModel::gutsOfDelete(); gutsOfDelete(0); int i; for (i=0;i<6;i++) { rowArray_[i]=NULL; columnArray_[i]=NULL; } // get an empty factorization so we can set tolerances etc factorization_ = new ClpFactorization(); // Say sparse factorization_->sparseThreshold(1); Clp_scalars scalars; CoinBigIndex numberRead; // get scalars numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp); if (numberRead!=1) return 1; // Fill in scalars optimizationDirection_ = scalars.optimizationDirection; memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double)); objectiveValue_ = scalars.objectiveValue; dualBound_ = scalars.dualBound; dualTolerance_ = scalars.dualTolerance; primalTolerance_ = scalars.primalTolerance; sumDualInfeasibilities_ = scalars.sumDualInfeasibilities; sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities; infeasibilityCost_ = scalars.infeasibilityCost; numberRows_ = scalars.numberRows; numberColumns_ = scalars.numberColumns; memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double)); numberIterations_ = scalars.numberIterations; problemStatus_ = scalars.problemStatus; setMaximumIterations(scalars.maximumIterations); lengthNames_ = scalars.lengthNames; numberDualInfeasibilities_ = scalars.numberDualInfeasibilities; numberDualInfeasibilitiesWithoutFree_ = scalars.numberDualInfeasibilitiesWithoutFree; numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities; numberRefinements_ = scalars.numberRefinements; scalingFlag_ = scalars.scalingFlag; algorithm_ = scalars.algorithm; specialOptions_ = scalars.specialOptions; // strings CoinBigIndex length; for (i=0;i<ClpLastStrParam;i++) { numberRead = fread(&length,sizeof(int),1,fp); if (numberRead!=1) return 1; if (length) { char * array = new char[length+1]; numberRead = fread(array,length,1,fp); if (numberRead!=1) return 1; array[length]='\0'; strParam_[i]=array; delete [] array; } } // arrays - in no particular order if (inDoubleArray(rowActivity_,numberRows_,fp)) return 1; if (inDoubleArray(columnActivity_,numberColumns_,fp)) return 1; if (inDoubleArray(dual_,numberRows_,fp)) return 1; if (inDoubleArray(reducedCost_,numberColumns_,fp)) return 1; if (inDoubleArray(rowLower_,numberRows_,fp)) return 1; if (inDoubleArray(rowUpper_,numberRows_,fp)) return 1; double * objective; if (inDoubleArray(objective,numberColumns_,fp)) return 1; delete objective_; objective_ = new ClpLinearObjective(objective,numberColumns_); delete [] objective; if (inDoubleArray(rowObjective_,numberRows_,fp)) return 1; if (inDoubleArray(columnLower_,numberColumns_,fp)) return 1; if (inDoubleArray(columnUpper_,numberColumns_,fp)) return 1; if (problemStatus_==1) { if (inDoubleArray(ray_,numberRows_,fp)) return 1; } else if (problemStatus_==2) { if (inDoubleArray(ray_,numberColumns_,fp)) return 1; } else { // ray should be null numberRead = fread(&length,sizeof(int),1,fp); if (numberRead!=1) return 1; if (length) return 2; } delete [] status_; status_=NULL; // status region numberRead = fread(&length,sizeof(int),1,fp); if (numberRead!=1) return 1; if (length) { if (length!=numberRows_+numberColumns_) return 1; status_ = new char unsigned[length]; numberRead = fread(status_,sizeof(char),length, fp); if (numberRead!=length) return 1; } if (lengthNames_) { char * array = new char[max(numberRows_,numberColumns_)*(lengthNames_+1)]; char * get = array; numberRead = fread(array,lengthNames_+1,numberRows_,fp); if (numberRead!=numberRows_) return 1; rowNames_ = std::vector<std::string> (); rowNames_.resize(numberRows_); for (i=0;i<numberRows_;i++) { rowNames_.push_back(get); get += lengthNames_+1; } get = array; numberRead = fread(array,lengthNames_+1,numberColumns_,fp); if (numberRead!=numberColumns_) return 1; columnNames_ = std::vector<std::string> (); columnNames_.resize(numberColumns_); for (i=0;i<numberColumns_;i++) { columnNames_.push_back(get); get += lengthNames_+1; } delete [] array; } // Pivot choices assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3); delete dualRowPivot_; switch ((scalars.dualPivotChoice&63)) { default: printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63); case 1: // Dantzig dualRowPivot_ = new ClpDualRowDantzig(); break; case 2: // Steepest - use mode dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6); break; } assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3); delete primalColumnPivot_; switch ((scalars.primalPivotChoice&63)) { default: printf("Need another primalPivot case %d\n", scalars.primalPivotChoice&63); case 1: // Dantzig primalColumnPivot_ = new ClpPrimalColumnDantzig(); break; case 2: // Steepest - use mode primalColumnPivot_ = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6); break; } assert(scalars.matrixStorageChoice==1); delete matrix_; // get arrays numberRead = fread(&length,sizeof(int),1,fp); if (numberRead!=1) return 1; double * elements = new double[length]; int * indices = new int[length]; CoinBigIndex * starts = new CoinBigIndex[numberColumns_+1]; int * lengths = new int[numberColumns_]; numberRead = fread(elements, sizeof(double),length,fp); if (numberRead!=length) return 1; numberRead = fread(indices, sizeof(int),length,fp); if (numberRead!=length) return 1; numberRead = fread(starts, sizeof(int),numberColumns_+1,fp); if (numberRead!=numberColumns_+1) return 1; numberRead = fread(lengths, sizeof(int),numberColumns_,fp); if (numberRead!=numberColumns_) return 1; // assign matrix CoinPackedMatrix * matrix = new CoinPackedMatrix(); matrix->assignMatrix(true, numberRows_, numberColumns_, length, elements, indices, starts, lengths); // and transfer to Clp matrix_ = new ClpPackedMatrix(matrix); // finished fclose(fp); return 0; } else { return -1; } return 0; }

void ClpSimplex::returnModel (  ClpSimplex &  otherModel  )

Return model - updates any scalars Definition at line 4529 of file ClpSimplex.cpp. References alpha_, columnDualInfeasibility_, columnDualSequence_, columnPrimalInfeasibility_, columnPrimalSequence_, directionIn_, directionOut_, dualIn_, dualOut_, largestDualError_, largestPrimalError_, largestSolutionError_, lowerIn_, lowerOut_, numberDualInfeasibilities_, numberDualInfeasibilitiesWithoutFree_, numberPrimalInfeasibilities_, numberTimesOptimal_, pivotRow_, primalToleranceToGetOptimal_, remainingDualInfeasibility_, ClpModel::returnModel(), rowDualInfeasibility_, rowDualSequence_, rowPrimalInfeasibility_, rowPrimalSequence_, sequenceIn_, sequenceOut_, sumDualInfeasibilities_, sumOfRelaxedDualInfeasibilities_, sumOfRelaxedPrimalInfeasibilities_, sumPrimalInfeasibilities_, theta_, upperIn_, upperOut_, valueIn_, and valueOut_. { ClpModel::returnModel(otherModel); otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_; otherModel.columnPrimalSequence_ = columnPrimalSequence_; otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_; otherModel.rowPrimalSequence_ = rowPrimalSequence_; otherModel.columnDualInfeasibility_ = columnDualInfeasibility_; otherModel.columnDualSequence_ = columnDualSequence_; otherModel.rowDualInfeasibility_ = rowDualInfeasibility_; otherModel.rowDualSequence_ = rowDualSequence_; otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_; otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_; otherModel.largestPrimalError_ = largestPrimalError_; otherModel.largestDualError_ = largestDualError_; otherModel.largestSolutionError_ = largestSolutionError_; otherModel.alpha_ = alpha_; otherModel.theta_ = theta_; otherModel.lowerIn_ = lowerIn_; otherModel.valueIn_ = valueIn_; otherModel.upperIn_ = upperIn_; otherModel.dualIn_ = dualIn_; otherModel.sequenceIn_ = sequenceIn_; otherModel.directionIn_ = directionIn_; otherModel.lowerOut_ = lowerOut_; otherModel.valueOut_ = valueOut_; otherModel.upperOut_ = upperOut_; otherModel.dualOut_ = dualOut_; otherModel.sequenceOut_ = sequenceOut_; otherModel.directionOut_ = directionOut_; otherModel.pivotRow_ = pivotRow_; otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_; otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_; otherModel.numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilitiesWithoutFree_; otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_; otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_; otherModel.numberTimesOptimal_ = numberTimesOptimal_; otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_; otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_; }

ClpDataSave ClpSimplex::saveData (   )

Save data. 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 Definition at line 4486 of file ClpSimplex.cpp. References dualBound_, ClpDataSave::dualBound_, factorization_, infeasibilityCost_, ClpDataSave::infeasibilityCost_, perturbation_, ClpDataSave::perturbation_, progress_, and ClpDataSave::sparseThreshold_. Referenced by ClpSimplexDual::dual(), ClpSimplexDual::fastDual(), ClpSimplexPrimal::primal(), and ClpSimplexPrimalQuadratic::primalQuadratic2(). { ClpDataSave saved; saved.dualBound_ = dualBound_; saved.infeasibilityCost_ = infeasibilityCost_; saved.sparseThreshold_ = factorization_->sparseThreshold(); saved.perturbation_ = perturbation_; // Progress indicator delete progress_; progress_ = new ClpSimplexProgress (this); return saved; }

int ClpSimplex::saveModel (  const char *  fileName  )

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. Definition at line 3015 of file ClpSimplex.cpp. References algorithm_, dualBound_, dualRowPivot_, dualTolerance_, ClpMatrixBase::getElements(), ClpMatrixBase::getIndices(), ClpMatrixBase::getNumCols(), ClpMatrixBase::getNumRows(), ClpMatrixBase::getVectorLengths(), ClpMatrixBase::getVectorStarts(), infeasibilityCost_, ClpModel::maximumIterations(), numberDualInfeasibilities_, numberDualInfeasibilitiesWithoutFree_, numberPrimalInfeasibilities_, numberRefinements_, ClpModel::objective(), primalColumnPivot_, primalTolerance_, scalingFlag_, specialOptions_, sumDualInfeasibilities_, sumPrimalInfeasibilities_, ClpMatrixBase::type(), ClpPrimalColumnPivot::type(), and ClpDualRowPivot::type(). { FILE * fp = fopen(fileName,"wb"); if (fp) { Clp_scalars scalars; int i; CoinBigIndex numberWritten; // Fill in scalars scalars.optimizationDirection = optimizationDirection_; memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double)); scalars.objectiveValue = objectiveValue_; scalars.dualBound = dualBound_; scalars.dualTolerance = dualTolerance_; scalars.primalTolerance = primalTolerance_; scalars.sumDualInfeasibilities = sumDualInfeasibilities_; scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_; scalars.infeasibilityCost = infeasibilityCost_; scalars.numberRows = numberRows_; scalars.numberColumns = numberColumns_; memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double)); scalars.numberIterations = numberIterations_; scalars.problemStatus = problemStatus_; scalars.maximumIterations = maximumIterations(); scalars.lengthNames = lengthNames_; scalars.numberDualInfeasibilities = numberDualInfeasibilities_; scalars.numberDualInfeasibilitiesWithoutFree = numberDualInfeasibilitiesWithoutFree_; scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_; scalars.numberRefinements = numberRefinements_; scalars.scalingFlag = scalingFlag_; scalars.algorithm = algorithm_; scalars.specialOptions = specialOptions_; scalars.dualPivotChoice = dualRowPivot_->type(); scalars.primalPivotChoice = primalColumnPivot_->type(); scalars.matrixStorageChoice = matrix_->type(); // put out scalars numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp); if (numberWritten!=1) return 1; // strings CoinBigIndex length; for (i=0;i<ClpLastStrParam;i++) { length = strParam_[i].size(); numberWritten = fwrite(&length,sizeof(int),1,fp); if (numberWritten!=1) return 1; if (length) { numberWritten = fwrite(strParam_[i].c_str(),length,1,fp); if (numberWritten!=1) return 1; } } // arrays - in no particular order if (outDoubleArray(rowActivity_,numberRows_,fp)) return 1; if (outDoubleArray(columnActivity_,numberColumns_,fp)) return 1; if (outDoubleArray(dual_,numberRows_,fp)) return 1; if (outDoubleArray(reducedCost_,numberColumns_,fp)) return 1; if (outDoubleArray(rowLower_,numberRows_,fp)) return 1; if (outDoubleArray(rowUpper_,numberRows_,fp)) return 1; if (outDoubleArray(objective(),numberColumns_,fp)) return 1; if (outDoubleArray(rowObjective_,numberRows_,fp)) return 1; if (outDoubleArray(columnLower_,numberColumns_,fp)) return 1; if (outDoubleArray(columnUpper_,numberColumns_,fp)) return 1; if (ray_) { if (problemStatus_==1) if (outDoubleArray(ray_,numberRows_,fp)) return 1; else if (problemStatus_==2) if (outDoubleArray(ray_,numberColumns_,fp)) return 1; else if (outDoubleArray(NULL,0,fp)) return 1; } else { if (outDoubleArray(NULL,0,fp)) return 1; } if (status_&&(numberRows_+numberColumns_)>0) { length = numberRows_+numberColumns_; numberWritten = fwrite(&length,sizeof(int),1,fp); if (numberWritten!=1) return 1; numberWritten = fwrite(status_,sizeof(char),length, fp); if (numberWritten!=length) return 1; } else { length = 0; numberWritten = fwrite(&length,sizeof(int),1,fp); if (numberWritten!=1) return 1; } if (lengthNames_) { char * array = new char[max(numberRows_,numberColumns_)*(lengthNames_+1)]; char * put = array; assert (numberRows_ == (int) rowNames_.size()); for (i=0;i<numberRows_;i++) { assert((int)rowNames_[i].size()<=lengthNames_); strcpy(put,rowNames_[i].c_str()); put += lengthNames_+1; } numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp); if (numberWritten!=numberRows_) return 1; put=array; assert (numberColumns_ == (int) columnNames_.size()); for (i=0;i<numberColumns_;i++) { assert((int) columnNames_[i].size()<=lengthNames_); strcpy(put,columnNames_[i].c_str()); put += lengthNames_+1; } numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp); if (numberWritten!=numberColumns_) return 1; delete [] array; } // just standard type at present assert (matrix_->type()==1); assert (matrix_->getNumCols() == numberColumns_); assert (matrix_->getNumRows() == numberRows_); // we are going to save with gaps length = matrix_->getVectorStarts()[numberColumns_-1] + matrix_->getVectorLengths()[numberColumns_-1]; numberWritten = fwrite(&length,sizeof(int),1,fp); if (numberWritten!=1) return 1; numberWritten = fwrite(matrix_->getElements(), sizeof(double),length,fp); if (numberWritten!=length) return 1; numberWritten = fwrite(matrix_->getIndices(), sizeof(int),length,fp); if (numberWritten!=length) return 1; numberWritten = fwrite(matrix_->getVectorStarts(), sizeof(int),numberColumns_+1,fp); if (numberWritten!=numberColumns_+1) return 1; numberWritten = fwrite(matrix_->getVectorLengths(), sizeof(int),numberColumns_,fp); if (numberWritten!=numberColumns_) return 1; // finished fclose(fp); return 0; } else { return -1; } }

int ClpSimplex::sequenceIn (   )  const [inline]

Return sequence In or Out Definition at line 608 of file ClpSimplex.hpp. References sequenceIn_. Referenced by housekeeping(), and ClpNetworkBasis::replaceColumn(). {return sequenceIn_;};

void ClpSimplex::setDirectionIn (  int  direction  )  [inline]

Set directionIn or Out Definition at line 623 of file ClpSimplex.hpp. References directionIn_. { directionIn_=direction;};

void ClpSimplex::setInitialDenseFactorization (  bool  onOff  )

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 Definition at line 4929 of file ClpSimplex.cpp. References specialOptions_. Referenced by dual(), and primal(). { if (onOff) specialOptions_ |= 8; else specialOptions_ &= ~8; }

void ClpSimplex::setSequenceIn (  int  sequence  )  [inline]

Set sequenceIn or Out Definition at line 613 of file ClpSimplex.hpp. References sequenceIn_. { sequenceIn_=sequence;};

void ClpSimplex::setValuesPassAction (  float  incomingInfeasibility, float  allowedInfeasibility )

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 Definition at line 4613 of file ClpSimplex.cpp. References incomingInfeasibility_. { incomingInfeasibility_=incomingInfeasibility; allowedInfeasibility_=allowedInfeasibility; assert(incomingInfeasibility_>=0.0); assert(allowedInfeasibility_>=incomingInfeasibility_); }

double* ClpSimplex::solutionRegion (  int  section  )  [inline]

Return row or column sections - not as much needed as it once was. These just map into single arrays Definition at line 572 of file ClpSimplex.hpp. References columnActivityWork_, and rowActivityWork_. Referenced by ClpNonLinearCost::checkInfeasibilities(), and ClpQuadraticInfo::createGradient(). { if (!section) return rowActivityWork_; else return columnActivityWork_;};

int ClpSimplex::startup (  int  ifValuesPass  )

Common bits of coding for dual and primal. Return s0 if okay, 1 if bad matrix, 2 if very bad factorization Definition at line 4369 of file ClpSimplex.cpp. References algorithm_, createRim(), dualTolerance_, factorization_, ClpMatrixBase::getNumElements(), gutsOfSolution(), CoinMessageHandler::message(), nonLinearCost_, numberTimesOptimal_, perturbation_, pivotRow_, primalTolerance_, sequenceIn_, and sequenceOut_. Referenced by ClpSimplexDual::dual(), ClpSimplexPrimal::primal(), and ClpSimplexPrimalQuadratic::primalQuadratic2(). { // sanity check // bad if empty (trap here to avoid using bad matrix_) if (!matrix_||!matrix_->getNumElements()) { handler_->message(CLP_EMPTY_PROBLEM,messages_) <<numberRows_ <<numberColumns_ <<0 <<CoinMessageEol; problemStatus_=4; return 2; } pivotRow_=-1; sequenceIn_=-1; sequenceOut_=-1; secondaryStatus_=0; primalTolerance_=dblParam_[ClpPrimalTolerance]; dualTolerance_=dblParam_[ClpDualTolerance]; if (problemStatus_!=10) numberIterations_=0; // put in standard form (and make row copy) // create modifiable copies of model rim and do optional scaling bool goodMatrix=createRim(7+8+16,true); if (goodMatrix) { // Model looks okay // Do initial factorization // and set certain stuff // We can either set increasing rows so ...IsBasic gives pivot row // or we can just increment iBasic one by one // for now let ...iBasic give pivot row factorization_->increasingRows(2); // row activities have negative sign factorization_->slackValue(-1.0); factorization_->zeroTolerance(1.0e-13); // Switch off dense (unless special option set) int saveThreshold = factorization_->denseThreshold(); factorization_->setDenseThreshold(0); // If values pass then perturb (otherwise may be optimal so leave a bit) if (ifValuesPass) { // do perturbation if asked for if (perturbation_<100) { if (algorithm_>0) { ((ClpSimplexPrimal *) this)->perturb(0); } else if (algorithm_<0) { ((ClpSimplexDual *) this)->perturb(); } } } // for primal we will change bounds using infeasibilityCost_ if (nonLinearCost_==NULL&&algorithm_>0) { // get a valid nonlinear cost function delete nonLinearCost_; nonLinearCost_= new ClpNonLinearCost(this); } // loop round to clean up solution if values pass int numberThrownOut = -1; int totalNumberThrownOut=0; while(numberThrownOut) { int status = internalFactorize(0+10*ifValuesPass); if (status<0) return 1; // some error else numberThrownOut = status; // for this we need clean basis so it is after factorize if (!numberThrownOut) numberThrownOut = gutsOfSolution( NULL,NULL, ifValuesPass!=0); totalNumberThrownOut+= numberThrownOut; } if (totalNumberThrownOut) handler_->message(CLP_SINGULARITIES,messages_) <<totalNumberThrownOut <<CoinMessageEol; // Switch back dense factorization_->setDenseThreshold(saveThreshold); problemStatus_ = -1; // number of times we have declared optimality numberTimesOptimal_=0; return 0; } else { // bad matrix return 2; } }

bool ClpSimplex::statusOfProblem (   )

Factorizes and returns true if optimal. Used by user Definition at line 4511 of file ClpSimplex.cpp. References createRim(), deleteRim(), dualFeasible(), gutsOfSolution(), and primalFeasible(). Referenced by pivot(). { // is factorization okay? assert (internalFactorize(1)==0); // put back original costs and then check // also move to work arrays createRim(4+32); //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double)); //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double)); gutsOfSolution(NULL,NULL); //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double)); //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double)); //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double)); deleteRim(-1); return (primalFeasible()&&dualFeasible()); }

int ClpSimplex::strongBranching (  int  numberVariables, const int *  variables, double *  newLower, double *  newUpper, double **  outputSolution, int *  outputStatus, int *  outputIterations, bool  stopOnFirstInfeasible = true, bool  alwaysFinish = false )

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 Reimplemented in ClpSimplexDual. Definition at line 2939 of file ClpSimplex.cpp. { return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables, newLower, newUpper,outputSolution, outputStatus, outputIterations, stopOnFirstInfeasible, alwaysFinish); }

int ClpSimplex::tightenPrimalBounds (  double  factor = 0.0  )

Tightens primal bounds to make dual faster. Unless fixed, 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) Definition at line 2485 of file ClpSimplex.cpp. References largeValue(), ClpModel::matrix(), CoinMessageHandler::message(), and ClpModel::primalTolerance(). { // Get a row copy in standard format CoinPackedMatrix copy; copy.reverseOrderedCopyOf(*matrix()); // get matrix data pointers const int * column = copy.getIndices(); const CoinBigIndex * rowStart = copy.getVectorStarts(); const int * rowLength = copy.getVectorLengths(); const double * element = copy.getElements(); int numberChanged=1,iPass=0; double large = largeValue(); // treat bounds > this as infinite #ifndef NDEBUG double large2= 1.0e10*large; #endif int numberInfeasible=0; int totalTightened = 0; double tolerance = primalTolerance(); // Save column bounds double * saveLower = new double [numberColumns_]; memcpy(saveLower,columnLower_,numberColumns_*sizeof(double)); double * saveUpper = new double [numberColumns_]; memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double)); int iRow, iColumn; // If wanted - tighten column bounds using solution if (factor) { assert (factor>1.0); double largest=0.0; for (iColumn=0;iColumn<numberColumns_;iColumn++) { if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) { largest = max(largest,fabs(columnActivity_[iColumn])); } } largest *= factor; for (iColumn=0;iColumn<numberColumns_;iColumn++) { if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) { columnUpper_[iColumn] = min(columnUpper_[iColumn],largest); columnLower_[iColumn] = max(columnLower_[iColumn],-largest); } } } #define MAXPASS 10 // Loop round seeing if we can tighten bounds // Would be faster to have a stack of possible rows // and we put altered rows back on stack int numberCheck=-1; while(numberChanged>numberCheck) { numberChanged = 0; // Bounds tightened this pass if (iPass==MAXPASS) break; iPass++; for (iRow = 0; iRow < numberRows_; iRow++) { if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) { // possible row int infiniteUpper = 0; int infiniteLower = 0; double maximumUp = 0.0; double maximumDown = 0.0; double newBound; CoinBigIndex rStart = rowStart[iRow]; CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow]; CoinBigIndex j; // Compute possible lower and upper ranges for (j = rStart; j < rEnd; ++j) { double value=element[j]; iColumn = column[j]; if (value > 0.0) { if (columnUpper_[iColumn] >= large) { ++infiniteUpper; } else { maximumUp += columnUpper_[iColumn] * value; } if (columnLower_[iColumn] <= -large) { ++infiniteLower; } else { maximumDown += columnLower_[iColumn] * value; } } else if (value<0.0) { if (columnUpper_[iColumn] >= large) { ++infiniteLower; } else { maximumDown += columnUpper_[iColumn] * value; } if (columnLower_[iColumn] <= -large) { ++infiniteUpper; } else { maximumUp += columnLower_[iColumn] * value; } } } // Build in a margin of error maximumUp += 1.0e-8*fabs(maximumUp); maximumDown -= 1.0e-8*fabs(maximumDown); double maxUp = maximumUp+infiniteUpper*1.0e31; double maxDown = maximumDown-infiniteLower*1.0e31; if (maxUp <= rowUpper_[iRow] + tolerance && maxDown >= rowLower_[iRow] - tolerance) { // Row is redundant - make totally free // NO - can't do this for postsolve // rowLower_[iRow]=-COIN_DBL_MAX; // rowUpper_[iRow]=COIN_DBL_MAX; //printf("Redundant row in presolveX %d\n",iRow); } else { if (maxUp < rowLower_[iRow] -100.0*tolerance || maxDown > rowUpper_[iRow]+100.0*tolerance) { // problem is infeasible - exit at once numberInfeasible++; break; } double lower = rowLower_[iRow]; double upper = rowUpper_[iRow]; for (j = rStart; j < rEnd; ++j) { double value=element[j]; iColumn = column[j]; double nowLower = columnLower_[iColumn]; double nowUpper = columnUpper_[iColumn]; if (value > 0.0) { // positive value if (lower>-large) { if (!infiniteUpper) { assert(nowUpper < large2); newBound = nowUpper + (lower - maximumUp) / value; // relax if original was large if (fabs(maximumUp)>1.0e8) newBound -= 1.0e-12*fabs(maximumUp); } else if (infiniteUpper==1&&nowUpper>large) { newBound = (lower -maximumUp) / value; // relax if original was large if (fabs(maximumUp)>1.0e8) newBound -= 1.0e-12*fabs(maximumUp); } else { newBound = -COIN_DBL_MAX; } if (newBound > nowLower + 1.0e-12&&newBound>-large) { // Tighten the lower bound columnLower_[iColumn] = newBound; numberChanged++; // check infeasible (relaxed) if (nowUpper - newBound < -100.0*tolerance) { numberInfeasible++; } // adjust double now; if (nowLower<-large) { now=0.0; infiniteLower--; } else { now = nowLower; } maximumDown += (newBound-now) * value; nowLower = newBound; #ifdef DEBUG checkCorrect(this,iRow, element, rowStart, rowLength, column, columnLower_, columnUpper_, infiniteUpper, infiniteLower, maximumUp, maximumDown); #endif } } if (upper <large) { if (!infiniteLower) { assert(nowLower >- large2); newBound = nowLower + (upper - maximumDown) / value; // relax if original was large if (fabs(maximumDown)>1.0e8) newBound += 1.0e-12*fabs(maximumDown); } else if (infiniteLower==1&&nowLower<-large) { newBound = (upper - maximumDown) / value; // relax if original was large if (fabs(maximumDown)>1.0e8) newBound += 1.0e-12*fabs(maximumDown); } else { newBound = COIN_DBL_MAX; } if (newBound < nowUpper - 1.0e-12&&newBound<large) { // Tighten the upper bound columnUpper_[iColumn] = newBound; numberChanged++; // check infeasible (relaxed) if (newBound - nowLower < -100.0*tolerance) { numberInfeasible++; } // adjust double now; if (nowUpper>large) { now=0.0; infiniteUpper--; } else { now = nowUpper; } maximumUp += (newBound-now) * value; nowUpper = newBound; #ifdef DEBUG checkCorrect(this,iRow, element, rowStart, rowLength, column, columnLower_, columnUpper_, infiniteUpper, infiniteLower, maximumUp, maximumDown); #endif } } } else { // negative value if (lower>-large) { if (!infiniteUpper) { assert(nowLower < large2); newBound = nowLower + (lower - maximumUp) / value; // relax if original was large if (fabs(maximumUp)>1.0e8) newBound += 1.0e-12*fabs(maximumUp); } else if (infiniteUpper==1&&nowLower<-large) { newBound = (lower -maximumUp) / value; // relax if original was large if (fabs(maximumUp)>1.0e8) newBound += 1.0e-12*fabs(maximumUp); } else { newBound = COIN_DBL_MAX; } if (newBound < nowUpper - 1.0e-12&&newBound<large) { // Tighten the upper bound columnUpper_[iColumn] = newBound; numberChanged++; // check infeasible (relaxed) if (newBound - nowLower < -100.0*tolerance) { numberInfeasible++; } // adjust double now; if (nowUpper>large) { now=0.0; infiniteLower--; } else { now = nowUpper; } maximumDown += (newBound-now) * value; nowUpper = newBound; #ifdef DEBUG checkCorrect(this,iRow, element, rowStart, rowLength, column, columnLower_, columnUpper_, infiniteUpper, infiniteLower, maximumUp, maximumDown); #endif } } if (upper <large) { if (!infiniteLower) { assert(nowUpper < large2); newBound = nowUpper + (upper - maximumDown) / value; // relax if original was large if (fabs(maximumDown)>1.0e8) newBound -= 1.0e-12*fabs(maximumDown); } else if (infiniteLower==1&&nowUpper>large) { newBound = (upper - maximumDown) / value; // relax if original was large if (fabs(maximumDown)>1.0e8) newBound -= 1.0e-12*fabs(maximumDown); } else { newBound = -COIN_DBL_MAX; } if (newBound > nowLower + 1.0e-12&&newBound>-large) { // Tighten the lower bound columnLower_[iColumn] = newBound; numberChanged++; // check infeasible (relaxed) if (nowUpper - newBound < -100.0*tolerance) { numberInfeasible++; } // adjust double now; if (nowLower<-large) { now=0.0; infiniteUpper--; } else { now = nowLower; } maximumUp += (newBound-now) * value; nowLower = newBound; #ifdef DEBUG checkCorrect(this,iRow, element, rowStart, rowLength, column, columnLower_, columnUpper_, infiniteUpper, infiniteLower, maximumUp, maximumDown); #endif } } } } } } } totalTightened += numberChanged; if (iPass==1) numberCheck=numberChanged>>4; if (numberInfeasible) break; } if (!numberInfeasible) { handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_) <<totalTightened <<CoinMessageEol; // Set bounds slightly loose double useTolerance = 1.0e-3; for (iColumn=0;iColumn<numberColumns_;iColumn++) { if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) { if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e8) { // relax enough so will have correct dj #if 1 columnLower_[iColumn]=max(saveLower[iColumn], columnLower_[iColumn]-100.0*useTolerance); columnUpper_[iColumn]=min(saveUpper[iColumn], columnUpper_[iColumn]+100.0*useTolerance); #else if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) { if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) { columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance; } else { columnLower_[iColumn]=saveLower[iColumn]; columnUpper_[iColumn]=min(saveUpper[iColumn], saveLower[iColumn]+100.0*useTolerance); } } else { if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) { columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance; } else { columnUpper_[iColumn]=saveUpper[iColumn]; columnLower_[iColumn]=max(saveLower[iColumn], saveUpper[iColumn]-100.0*useTolerance); } } #endif } else { if (columnUpper_[iColumn]<saveUpper[iColumn]) { // relax a bit columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*useTolerance, saveUpper[iColumn]); } if (columnLower_[iColumn]>saveLower[iColumn]) { // relax a bit columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*useTolerance, saveLower[iColumn]); } } } } } else { handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_) <<numberInfeasible <<CoinMessageEol; // restore column bounds memcpy(columnLower_,saveLower,numberColumns_*sizeof(double)); memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double)); } delete [] saveLower; delete [] saveUpper; return (numberInfeasible); }

void ClpSimplex::times (  double  scalar, const double *  x, double *  y )  const

Return y + A * x * scalar in y. Precondition: x must be of size numColumns() y must be of size numRows() Definition at line 2369 of file ClpSimplex.cpp. References columnScale_, rowScale_, and ClpMatrixBase::times(). Referenced by cleanStatus(), computePrimals(), ClpSimplexPrimalQuadratic::endQuadratic(), gutsOfSolution(), and ClpSimplexPrimalQuadratic::makeQuadratic(). { if (rowScale_) matrix_->times(scalar,x,y,rowScale_,columnScale_); else matrix_->times(scalar,x,y); }

void ClpSimplex::transposeTimes (  double  scalar, const double *  x, double *  y )  const

Return y + x * scalar * A in y. Precondition: x must be of size numRows() y must be of size numColumns() Definition at line 2378 of file ClpSimplex.cpp. References columnScale_, rowScale_, and ClpMatrixBase::transposeTimes(). Referenced by computeDuals(), and ClpSimplexDual::dual(). { if (rowScale_) matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_); else matrix_->transposeTimes(scalar,x,y); }

void ClpSimplex::unpack (  CoinIndexedVector *  rowArray, int  sequence )  const

Unpacks one column of the matrix into indexed array Slack if sequence>= numberColumns Also applies scaling if needed Definition at line 1912 of file ClpSimplex.cpp. References CoinIndexedVector::clear(), CoinIndexedVector::insert(), and ClpMatrixBase::unpack(). { rowArray->clear(); if (sequence>=numberColumns_) { //slack rowArray->insert(sequence-numberColumns_,-1.0); } else { // column matrix_->unpack(this,rowArray,sequence); } }

void ClpSimplex::unpack (  CoinIndexedVector *  rowArray  )  const

Unpacks one column of the matrix into indexed array Uses sequenceIn_ Also applies scaling if needed Definition at line 1900 of file ClpSimplex.cpp. References CoinIndexedVector::clear(), CoinIndexedVector::insert(), sequenceIn_, and ClpMatrixBase::unpack(). Referenced by ClpSimplexDual::dualRow(), pivot(), ClpSimplexPrimal::pivotResult(), ClpNetworkBasis::replaceColumn(), ClpSimplexDual::statusOfProblemInDual(), ClpSimplexPrimalQuadratic::whileIterating(), and ClpSimplexDual::whileIterating(). { rowArray->clear(); if (sequenceIn_>=numberColumns_) { //slack rowArray->insert(sequenceIn_-numberColumns_,-1.0); } else { // column matrix_->unpack(this,rowArray,sequenceIn_); } }

void ClpSimplex::unpackPacked (  CoinIndexedVector *  rowArray, int  sequence )

Unpacks one column of the matrix into indexed array as packed vector Slack if sequence>= numberColumns Also applies scaling if needed Definition at line 1944 of file ClpSimplex.cpp. References CoinIndexedVector::clear(), CoinIndexedVector::denseVector(), CoinIndexedVector::getIndices(), CoinIndexedVector::setNumElements(), CoinIndexedVector::setPackedMode(), and ClpMatrixBase::unpackPacked(). { rowArray->clear(); if (sequence>=numberColumns_) { //slack int * index = rowArray->getIndices(); double * array = rowArray->denseVector(); array[0]=-1.0; index[0]=sequence-numberColumns_; rowArray->setNumElements(1); rowArray->setPackedMode(true); } else { // column matrix_->unpackPacked(this,rowArray,sequence); } }

void ClpSimplex::unpackPacked (  CoinIndexedVector *  rowArray  )

Unpacks one column of the matrix into indexed array as packed vector Uses sequenceIn_ Also applies scaling if needed Definition at line 1927 of file ClpSimplex.cpp. References CoinIndexedVector::clear(), CoinIndexedVector::denseVector(), CoinIndexedVector::getIndices(), sequenceIn_, CoinIndexedVector::setNumElements(), CoinIndexedVector::setPackedMode(), and ClpMatrixBase::unpackPacked(). Referenced by ClpSimplexPrimal::pivotResult(), and ClpSimplexDual::whileIterating(). { rowArray->clear(); if (sequenceIn_>=numberColumns_) { //slack int * index = rowArray->getIndices(); double * array = rowArray->denseVector(); array[0]=-1.0; index[0]=sequenceIn_-numberColumns_; rowArray->setNumElements(1); rowArray->setPackedMode(true); } else { // column matrix_->unpackPacked(this,rowArray,sequenceIn_); } }

---

Friends And Related Function Documentation

void ClpSimplexUnitTest (  const std::string &  mpsDir, const std::string &  netlibDir )  [friend]

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 Definition at line 371 of file unitTest.cpp. { CoinRelFltEq eq(0.000001); { ClpSimplex solution; // matrix data //deliberate hiccup of 2 between 0 and 1 CoinBigIndex start[5]={0,4,7,8,9}; int length[5]={2,3,1,1,1}; int rows[11]={0,2,-1,-1,0,1,2,0,1,2}; double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1}; CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length); // rim data double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0}; double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10}; double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10}; double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0}; double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0}; // basis 1 int rowBasis1[3]={-1,-1,-1}; int colBasis1[5]={1,1,-1,-1,1}; solution.loadProblem(matrix,colLower,colUpper,objective, rowLower,rowUpper); int i; solution.createStatus(); for (i=0;i<3;i++) { if (rowBasis1[i]<0) { solution.setRowStatus(i,ClpSimplex::atLowerBound); } else { solution.setRowStatus(i,ClpSimplex::basic); } } for (i=0;i<5;i++) { if (colBasis1[i]<0) { solution.setColumnStatus(i,ClpSimplex::atLowerBound); } else { solution.setColumnStatus(i,ClpSimplex::basic); } } solution.setLogLevel(3+4+8+16+32); solution.primal(); for (i=0;i<3;i++) { if (rowBasis1[i]<0) { solution.setRowStatus(i,ClpSimplex::atLowerBound); } else { solution.setRowStatus(i,ClpSimplex::basic); } } for (i=0;i<5;i++) { if (colBasis1[i]<0) { solution.setColumnStatus(i,ClpSimplex::atLowerBound); } else { solution.setColumnStatus(i,ClpSimplex::basic); } } // intricate stuff does not work with scaling solution.scaling(0); assert(!solution.factorize ( )); const double * colsol = solution.primalColumnSolution(); const double * rowsol = solution.primalRowSolution(); solution.getSolution(rowsol,colsol); double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0}; for (i=0;i<5;i++) { assert(eq(colsol[i],colsol1[i])); } // now feed in again without actually doing factorization ClpFactorization factorization2 = *solution.factorization(); ClpSimplex solution2 = solution; solution2.setFactorization(factorization2); solution2.createStatus(); for (i=0;i<3;i++) { if (rowBasis1[i]<0) { solution2.setRowStatus(i,ClpSimplex::atLowerBound); } else { solution2.setRowStatus(i,ClpSimplex::basic); } } for (i=0;i<5;i++) { if (colBasis1[i]<0) { solution2.setColumnStatus(i,ClpSimplex::atLowerBound); } else { solution2.setColumnStatus(i,ClpSimplex::basic); } } // intricate stuff does not work with scaling solution2.scaling(0); solution2.getSolution(rowsol,colsol); colsol = solution2.primalColumnSolution(); rowsol = solution2.primalRowSolution(); for (i=0;i<5;i++) { assert(eq(colsol[i],colsol1[i])); } solution2.setDualBound(0.1); solution2.dual(); objective[2]=-1.0; objective[3]=-0.5; objective[4]=10.0; solution.dual(); for (i=0;i<3;i++) { rowLower[i]=-1.0e50; colUpper[i+2]=0.0; } solution.setLogLevel(3); solution.dual(); double rowObjective[]={1.0,0.5,-10.0}; solution.loadProblem(matrix,colLower,colUpper,objective, rowLower,rowUpper,rowObjective); solution.dual(); } { CoinMpsIO m; std::string fn = mpsDir+"exmip1"; m.readMps(fn.c_str(),"mps"); ClpSimplex solution; solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); solution.dual(); } // Test Message handler { CoinMpsIO m; std::string fn = mpsDir+"exmip1"; //fn = "Test/subGams4"; m.readMps(fn.c_str(),"mps"); ClpSimplex model; model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); // Message handler MyMessageHandler messageHandler(&model); std::cout<<"Testing derived message handler"<<std::endl; model.passInMessageHandler(&messageHandler); model.messagesPointer()->setDetailMessage(1,102); model.setFactorizationFrequency(10); model.primal(); // Write saved solutions int nc = model.getNumCols(); int s; std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints(); int numSavedSolutions = fep.size(); for ( s=0; s<numSavedSolutions; ++s ) { const StdVectorDouble & solnVec = fep[s]; for ( int c=0; c<nc; ++c ) { if (fabs(solnVec[c])>1.0e-8) std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl; } } // Solve again without scaling // and maximize then minimize messageHandler.clearFeasibleExtremePoints(); model.scaling(0); model.setOptimizationDirection(-1); model.primal(); model.setOptimizationDirection(1); model.primal(); fep = messageHandler.getFeasibleExtremePoints(); numSavedSolutions = fep.size(); for ( s=0; s<numSavedSolutions; ++s ) { const StdVectorDouble & solnVec = fep[s]; for ( int c=0; c<nc; ++c ) { if (fabs(solnVec[c])>1.0e-8) std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl; } } } // test steepest edge { CoinMpsIO m; std::string fn = netlibDir+"finnis"; m.readMps(fn.c_str(),"mps"); ClpModel model; model.loadProblem(*m.getMatrixByCol(),m.getColLower(), m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); ClpSimplex solution(model); solution.scaling(1); solution.setDualBound(1.0e8); //solution.factorization()->maximumPivots(1); //solution.setLogLevel(3); solution.setDualTolerance(1.0e-7); // set objective sense, ClpDualRowSteepest steep; solution.setDualRowPivotAlgorithm(steep); solution.setDblParam(ClpObjOffset,m.objectiveOffset()); solution.dual(); } // test normal solution { CoinMpsIO m; std::string fn = netlibDir+"afiro"; m.readMps(fn.c_str(),"mps"); ClpSimplex solution; ClpModel model; // do twice - without and with scaling int iPass; for (iPass=0;iPass<2;iPass++) { // explicit row objective for testing int nr = m.getNumRows(); double * rowObj = new double[nr]; CoinFillN(rowObj,nr,0.0); model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper(),rowObj); delete [] rowObj; solution = ClpSimplex(model); if (iPass) { solution.scaling(); } solution.dual(); solution.dual(); // test optimal assert (solution.status()==0); int numberColumns = solution.numberColumns(); int numberRows = solution.numberRows(); CoinPackedVector colsol(numberColumns,solution.primalColumnSolution()); double * objective = solution.objective(); double objValue = colsol.dotProduct(objective); CoinRelFltEq eq(1.0e-8); assert(eq(objValue,-4.6475314286e+02)); double * lower = solution.columnLower(); double * upper = solution.columnUpper(); double * sol = solution.primalColumnSolution(); double * result = new double[numberColumns]; CoinFillN ( result, numberColumns,0.0); solution.matrix()->transposeTimes(solution.dualRowSolution(), result); int iRow , iColumn; // see if feasible and dual feasible for (iColumn=0;iColumn<numberColumns;iColumn++) { double value = sol[iColumn]; assert(value<upper[iColumn]+1.0e-8); assert(value>lower[iColumn]-1.0e-8); value = objective[iColumn]-result[iColumn]; assert (value>-1.0e-5); if (sol[iColumn]>1.0e-5) assert (value<1.0e-5); } delete [] result; result = new double[numberRows]; CoinFillN ( result, numberRows,0.0); solution.matrix()->times(colsol, result); lower = solution.rowLower(); upper = solution.rowUpper(); sol = solution.primalRowSolution(); for (iRow=0;iRow<numberRows;iRow++) { double value = result[iRow]; assert(eq(value,sol[iRow])); assert(value<upper[iRow]+1.0e-8); assert(value>lower[iRow]-1.0e-8); } delete [] result; // test row objective double * rowObjective = solution.rowObjective(); CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective); CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective); // this sets up all slack basis solution.createStatus(); solution.dual(); CoinFillN(rowObjective,numberRows,0.0); CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective); solution.dual(); } } // test unbounded { CoinMpsIO m; std::string fn = netlibDir+"brandy"; m.readMps(fn.c_str(),"mps"); ClpSimplex solution; // do twice - without and with scaling int iPass; for (iPass=0;iPass<2;iPass++) { solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); if (iPass) solution.scaling(); solution.setOptimizationDirection(-1); // test unbounded and ray #ifdef DUAL solution.setDualBound(100.0); solution.dual(); #else solution.primal(); #endif assert (solution.status()==2); int numberColumns = solution.numberColumns(); int numberRows = solution.numberRows(); double * lower = solution.columnLower(); double * upper = solution.columnUpper(); double * sol = solution.primalColumnSolution(); double * ray = solution.unboundedRay(); double * objective = solution.objective(); double objChange=0.0; int iRow , iColumn; // make sure feasible and columns form ray for (iColumn=0;iColumn<numberColumns;iColumn++) { double value = sol[iColumn]; assert(value<upper[iColumn]+1.0e-8); assert(value>lower[iColumn]-1.0e-8); value = ray[iColumn]; if (value>0.0) assert(upper[iColumn]>1.0e30); else if (value<0.0) assert(lower[iColumn]<-1.0e30); objChange += value*objective[iColumn]; } // make sure increasing objective assert(objChange>0.0); double * result = new double[numberRows]; CoinFillN ( result, numberRows,0.0); solution.matrix()->times(sol, result); lower = solution.rowLower(); upper = solution.rowUpper(); sol = solution.primalRowSolution(); for (iRow=0;iRow<numberRows;iRow++) { double value = result[iRow]; assert(eq(value,sol[iRow])); assert(value<upper[iRow]+2.0e-8); assert(value>lower[iRow]-2.0e-8); } CoinFillN ( result, numberRows,0.0); solution.matrix()->times(ray, result); // there may be small differences (especially if scaled) for (iRow=0;iRow<numberRows;iRow++) { double value = result[iRow]; if (value>1.0e-8) assert(upper[iRow]>1.0e30); else if (value<-1.0e-8) assert(lower[iRow]<-1.0e30); } delete [] result; delete [] ray; } } // test infeasible { CoinMpsIO m; std::string fn = netlibDir+"brandy"; m.readMps(fn.c_str(),"mps"); ClpSimplex solution; // do twice - without and with scaling int iPass; for (iPass=0;iPass<2;iPass++) { solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); if (iPass) solution.scaling(); // test infeasible and ray solution.columnUpper()[0]=0.0; #ifdef DUAL solution.setDualBound(100.0); solution.dual(); #else solution.primal(); #endif assert (solution.status()==1); int numberColumns = solution.numberColumns(); int numberRows = solution.numberRows(); double * lower = solution.rowLower(); double * upper = solution.rowUpper(); double * ray = solution.infeasibilityRay(); assert(ray); // construct proof of infeasibility int iRow , iColumn; double lo=0.0,up=0.0; int nl=0,nu=0; for (iRow=0;iRow<numberRows;iRow++) { if (lower[iRow]>-1.0e20) { lo += ray[iRow]*lower[iRow]; } else { if (ray[iRow]>1.0e-8) nl++; } if (upper[iRow]<1.0e20) { up += ray[iRow]*upper[iRow]; } else { if (ray[iRow]>1.0e-8) nu++; } } if (nl) lo=-1.0e100; if (nu) up=1.0e100; double * result = new double[numberColumns]; double lo2=0.0,up2=0.0; CoinFillN ( result, numberColumns,0.0); solution.matrix()->transposeTimes(ray, result); lower = solution.columnLower(); upper = solution.columnUpper(); nl=nu=0; for (iColumn=0;iColumn<numberColumns;iColumn++) { if (result[iColumn]>1.0e-8) { if (lower[iColumn]>-1.0e20) lo2 += result[iColumn]*lower[iColumn]; else nl++; if (upper[iColumn]<1.0e20) up2 += result[iColumn]*upper[iColumn]; else nu++; } else if (result[iColumn]<-1.0e-8) { if (lower[iColumn]>-1.0e20) up2 += result[iColumn]*lower[iColumn]; else nu++; if (upper[iColumn]<1.0e20) lo2 += result[iColumn]*upper[iColumn]; else nl++; } } if (nl) lo2=-1.0e100; if (nu) up2=1.0e100; // make sure inconsistency assert(lo2>up||up2<lo); delete [] result; delete [] ray; } } // test delete and add { CoinMpsIO m; std::string fn = netlibDir+"brandy"; m.readMps(fn.c_str(),"mps"); ClpSimplex solution; solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); solution.dual(); CoinRelFltEq eq(1.0e-8); assert(eq(solution.objectiveValue(),1.5185098965e+03)); int numberColumns = solution.numberColumns(); int numberRows = solution.numberRows(); double * saveObj = new double [numberColumns]; double * saveLower = new double[numberRows+numberColumns]; double * saveUpper = new double[numberRows+numberColumns]; int * which = new int [numberRows+numberColumns]; int numberElements = m.getMatrixByCol()->getNumElements(); int * starts = new int[numberRows+numberColumns]; int * index = new int[numberElements]; double * element = new double[numberElements]; const int * startM; const int * lengthM; const int * indexM; const double * elementM; int n,nel; // delete non basic columns n=0; nel=0; int iRow , iColumn; double * lower = solution.columnLower(); double * upper = solution.columnUpper(); double * objective = solution.objective(); startM = m.getMatrixByCol()->getVectorStarts(); lengthM = m.getMatrixByCol()->getVectorLengths(); indexM = m.getMatrixByCol()->getIndices(); elementM = m.getMatrixByCol()->getElements(); starts[0]=0; for (iColumn=0;iColumn<numberColumns;iColumn++) { if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) { saveObj[n]=objective[iColumn]; saveLower[n]=lower[iColumn]; saveUpper[n]=upper[iColumn]; int j; for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) { index[nel]=indexM[j]; element[nel++]=elementM[j]; } which[n++]=iColumn; starts[n]=nel; } } solution.deleteColumns(n,which); solution.dual(); // Put back solution.addColumns(n,saveLower,saveUpper,saveObj, starts,index,element); solution.dual(); assert(eq(solution.objectiveValue(),1.5185098965e+03)); // reload with original solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); // delete half rows n=0; nel=0; lower = solution.rowLower(); upper = solution.rowUpper(); startM = m.getMatrixByRow()->getVectorStarts(); lengthM = m.getMatrixByRow()->getVectorLengths(); indexM = m.getMatrixByRow()->getIndices(); elementM = m.getMatrixByRow()->getElements(); starts[0]=0; for (iRow=0;iRow<numberRows;iRow++) { if ((iRow&1)==0) { saveLower[n]=lower[iRow]; saveUpper[n]=upper[iRow]; int j; for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) { index[nel]=indexM[j]; element[nel++]=elementM[j]; } which[n++]=iRow; starts[n]=nel; } } solution.deleteRows(n,which); solution.dual(); // Put back solution.addRows(n,saveLower,saveUpper, starts,index,element); solution.dual(); assert(eq(solution.objectiveValue(),1.5185098965e+03)); // Zero out status array to give some interest memset(solution.statusArray()+numberColumns,0,numberRows); solution.primal(1); assert(eq(solution.objectiveValue(),1.5185098965e+03)); delete [] saveObj; delete [] saveLower; delete [] saveUpper; delete [] which; delete [] starts; delete [] index; delete [] element; } #if 0 // Test barrier { CoinMpsIO m; std::string fn = mpsDir+"exmip1"; m.readMps(fn.c_str(),"mps"); ClpInterior solution; solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); solution.primalDual(); } #endif // test network //#define QUADRATIC #ifndef QUADRATIC if (1) { std::string fn = mpsDir+"input.130"; int numberColumns; int numberRows; FILE * fp = fopen(fn.c_str(),"r"); if (!fp) { // Try in Samples fn = "Samples/input.130"; fp = fopen(fn.c_str(),"r"); } if (!fp) { fprintf(stderr,"Unable to open file input.130 in mpsDir or Samples directory\n"); } else { int problem; char temp[100]; // read and skip fscanf(fp,"%s",temp); assert (!strcmp(temp,"BEGIN")); fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, &numberColumns); // scan down to SUPPLY while (fgets(temp,100,fp)) { if (!strncmp(temp,"SUPPLY",6)) break; } if (strncmp(temp,"SUPPLY",6)) { fprintf(stderr,"Unable to find SUPPLY\n"); exit(2); } // get space for rhs double * lower = new double[numberRows]; double * upper = new double[numberRows]; int i; for (i=0;i<numberRows;i++) { lower[i]=0.0; upper[i]=0.0; } // ***** Remember to convert to C notation while (fgets(temp,100,fp)) { int row; int value; if (!strncmp(temp,"ARCS",4)) break; sscanf(temp,"%d %d",&row,&value); upper[row-1]=-value; lower[row-1]=-value; } if (strncmp(temp,"ARCS",4)) { fprintf(stderr,"Unable to find ARCS\n"); exit(2); } // number of columns may be underestimate int * head = new int[2*numberColumns]; int * tail = new int[2*numberColumns]; int * ub = new int[2*numberColumns]; int * cost = new int[2*numberColumns]; // ***** Remember to convert to C notation numberColumns=0; while (fgets(temp,100,fp)) { int iHead; int iTail; int iUb; int iCost; if (!strncmp(temp,"DEMAND",6)) break; sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb); iHead--; iTail--; head[numberColumns]=iHead; tail[numberColumns]=iTail; ub[numberColumns]=iUb; cost[numberColumns]=iCost; numberColumns++; } if (strncmp(temp,"DEMAND",6)) { fprintf(stderr,"Unable to find DEMAND\n"); exit(2); } // ***** Remember to convert to C notation while (fgets(temp,100,fp)) { int row; int value; if (!strncmp(temp,"END",3)) break; sscanf(temp,"%d %d",&row,&value); upper[row-1]=value; lower[row-1]=value; } if (strncmp(temp,"END",3)) { fprintf(stderr,"Unable to find END\n"); exit(2); } printf("Problem %d has %d rows and %d columns\n",problem, numberRows,numberColumns); fclose(fp); ClpSimplex model; // now build model double * objective =new double[numberColumns]; double * lowerColumn = new double[numberColumns]; double * upperColumn = new double[numberColumns]; double * element = new double [2*numberColumns]; int * start = new int[numberColumns+1]; int * row = new int[2*numberColumns]; start[numberColumns]=2*numberColumns; for (i=0;i<numberColumns;i++) { start[i]=2*i; element[2*i]=-1.0; element[2*i+1]=1.0; row[2*i]=head[i]; row[2*i+1]=tail[i]; lowerColumn[i]=0.0; upperColumn[i]=ub[i]; objective[i]=cost[i]; } // Create Packed Matrix CoinPackedMatrix matrix; int * lengths = NULL; matrix.assignMatrix(true,numberRows,numberColumns, 2*numberColumns,element,row,start,lengths); // load model model.loadProblem(matrix, lowerColumn,upperColumn,objective, lower,upper); model.factorization()->maximumPivots(200+model.numberRows()/100); model.createStatus(); double time1 = CoinCpuTime(); model.dual(); std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl; ClpPlusMinusOneMatrix plusMinus(matrix); assert (plusMinus.getIndices()); // would be zero if not +- one model.loadProblem(plusMinus, lowerColumn,upperColumn,objective, lower,upper); model.factorization()->maximumPivots(200+model.numberRows()/100); model.createStatus(); time1 = CoinCpuTime(); model.dual(); std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl; ClpNetworkMatrix network(numberColumns,head,tail); model.loadProblem(network, lowerColumn,upperColumn,objective, lower,upper); model.factorization()->maximumPivots(200+model.numberRows()/100); model.createStatus(); time1 = CoinCpuTime(); model.dual(); std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl; delete [] lower; delete [] upper; delete [] head; delete [] tail; delete [] ub; delete [] cost; delete [] objective; delete [] lowerColumn; delete [] upperColumn; } } #endif #ifdef QUADRATIC // Test quadratic to solve linear if (1) { CoinMpsIO m; std::string fn = mpsDir+"exmip1"; m.readMps(fn.c_str(),"mps"); ClpSimplex solution; solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); //solution.dual(); // get quadratic part int numberColumns=solution.numberColumns(); int * start=new int [numberColumns+1]; int * column = new int[numberColumns]; double * element = new double[numberColumns]; int i; start[0]=0; int n=0; int kk=numberColumns-1; int kk2=numberColumns-1; for (i=0;i<numberColumns;i++) { if (i>=kk) { column[n]=i; if (i>=kk2) element[n]=1.0e-1; else element[n]=0.0; n++; } start[i+1]=n; } // Load up objective solution.loadQuadraticObjective(numberColumns,start,column,element); delete [] start; delete [] column; delete [] element; //solution.quadraticSLP(50,1.0e-4); double objValue = solution.getObjValue(); CoinRelFltEq eq(1.0e-4); //assert(eq(objValue,820.0)); solution.setLogLevel(63); solution.quadraticPrimal(1); objValue = solution.getObjValue(); //assert(eq(objValue,3.2368421)); //exit(77); } // Test quadratic if (1) { CoinMpsIO m; std::string fn = mpsDir+"share2qp"; //fn = "share2qpb"; m.readMps(fn.c_str(),"mps"); ClpSimplex model; model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); model.dual(); // get quadratic part int * start=NULL; int * column = NULL; double * element = NULL; m.readQuadraticMps(NULL,start,column,element,2); int column2[200]; double element2[200]; int start2[80]; int j; start2[0]=0; int nel=0; bool good=false; for (j=0;j<79;j++) { if (start[j]==start[j+1]) { column2[nel]=j; element2[nel]=0.0; nel++; } else { int i; for (i=start[j];i<start[j+1];i++) { column2[nel]=column[i]; element2[nel++]=element[i]; } } start2[j+1]=nel; } // Load up objective if (good) model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2); else model.loadQuadraticObjective(model.numberColumns(),start,column,element); delete [] start; delete [] column; delete [] element; int numberColumns=model.numberColumns(); #if 0 model.quadraticSLP(50,1.0e-4); #else // Get feasible ClpObjective * saveObjective = model.objectiveAsObject()->clone(); ClpLinearObjective zeroObjective(NULL,numberColumns); model.setObjective(&zeroObjective); model.dual(); model.setObjective(saveObjective); delete saveObjective; #endif model.setLogLevel(63); //exit(77); model.setFactorizationFrequency(10); model.quadraticPrimal(1); double objValue = model.getObjValue(); const double * solution = model.primalColumnSolution(); int i; for (i=0;i<numberColumns;i++) if (solution[i]) printf("%d %g\n",i,solution[i]); CoinRelFltEq eq(1.0e-4); assert(eq(objValue,-400.92)); } if (0) { CoinMpsIO m; std::string fn = "./beale"; //fn = "./jensen"; m.readMps(fn.c_str(),"mps"); ClpSimplex solution; solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), m.getObjCoefficients(), m.getRowLower(),m.getRowUpper()); solution.setDblParam(ClpObjOffset,m.objectiveOffset()); solution.dual(); // get quadratic part int * start=NULL; int * column = NULL; double * element = NULL; m.readQuadraticMps(NULL,start,column,element,2); // Load up objective solution.loadQuadraticObjective(solution.numberColumns(),start,column,element); delete [] start; delete [] column; delete [] element; solution.quadraticPrimal(1); solution.quadraticSLP(50,1.0e-4); double objValue = solution.getObjValue(); CoinRelFltEq eq(1.0e-4); assert(eq(objValue,0.5)); solution.quadraticPrimal(); objValue = solution.getObjValue(); assert(eq(objValue,0.5)); } #endif }

---

Member Data Documentation

int ClpSimplex::forceFactorization_ [protected]

Now for some reliability aids This forces re-factorization early Definition at line 907 of file ClpSimplex.hpp. Referenced by forceFactorization(), gutsOfCopy(), and housekeeping().

float ClpSimplex::incomingInfeasibility_ [protected]

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 Definition at line 947 of file ClpSimplex.hpp. Referenced by gutsOfCopy(), gutsOfSolution(), and setValuesPassAction().

ClpNonLinearCost* ClpSimplex::nonLinearCost_ [protected]

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. Definition at line 922 of file ClpSimplex.hpp. Referenced by ClpSimplex(), createPiecewiseLinearCosts(), gutsOfCopy(), gutsOfDelete(), gutsOfSolution(), nonLinearCost(), operator=(), originalModel(), startup(), and ~ClpSimplex().

int ClpSimplex::perturbation_ [protected]

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 Definition at line 915 of file ClpSimplex.hpp. Referenced by dual(), gutsOfCopy(), initialSolve(), perturbation(), primal(), restoreData(), saveData(), and startup().

unsigned int ClpSimplex::specialOptions_ [protected]

For advanced options 1 - Don't keep changing infeasibility weight 2 - Keep nonLinearCost round solves 4 - Force outgoing variables to exact bound (primal) 8 - Safe to use dense initial factorization Definition at line 929 of file ClpSimplex.hpp. Referenced by ClpSimplex(), createPiecewiseLinearCosts(), gutsOfCopy(), gutsOfDelete(), gutsOfSolution(), operator=(), restoreModel(), saveModel(), and setInitialDenseFactorization().

---

The documentation for this class was generated from the following files:

• ClpSimplex.hpp
• ClpSimplex.cpp
• ClpSolve.cpp

---