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