1 /* $Id: CoinPresolveMatrix.hpp 2083 2019-01-06 19:38:09Z unxusr $ */
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 #ifndef CoinPresolveMatrix_H
7 #define CoinPresolveMatrix_H
8
9 #include "CoinPragma.hpp"
10 #include "CoinPackedMatrix.hpp"
11 #include "CoinMessage.hpp"
12 #include "CoinTime.hpp"
13
14 #include <cmath>
15 #include <cassert>
16 #include <cfloat>
17 #include <cassert>
18 #include <cstdlib>
19
20 //# define COIN_PRESOLVE_TUNING 2
21 #if PRESOLVE_DEBUG > 0
22 #include "CoinFinite.hpp"
23 #endif
24
25 /*! \file
26
27 Declarations for CoinPresolveMatrix and CoinPostsolveMatrix and their
28 common base class CoinPrePostsolveMatrix. Also declarations for
29 CoinPresolveAction and a number of non-member utility functions.
30 */
31
32 #if defined(_MSC_VER)
33 // Avoid MS Compiler problem in recognizing type to delete
34 // by casting to type.
35 // Is this still necessary? -- lh, 111202 --
36 #define deleteAction(array, type) delete[]((type)array)
37 #else
38 #define deleteAction(array, type) delete[] array
39 #endif
40
41 /*
42 Define PRESOLVE_DEBUG and PRESOLVE_CONSISTENCY on the configure command
43 line or in a Makefile! See comments in CoinPresolvePsdebug.hpp.
44 */
45 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
46
47 #define PRESOLVE_STMT(s) s
48
49 #define PRESOLVEASSERT(x) \
50 ((x) ? 1 : ((std::cerr << "FAILED ASSERTION at line " << __LINE__ << ": " #x "\n"), abort(), 0))
51
DIE(const char * s)52 inline void DIE(const char *s)
53 {
54 std::cout << s;
55 abort();
56 }
57
58 /*! \brief Indicate column or row present at start of postsolve
59
60 This code is used during postsolve in [cr]done to indicate columns and rows
61 that are present in the presolved system (i.e., present at the start of
62 postsolve processing).
63
64 \todo
65 There are a bunch of these code definitions, scattered through presolve
66 files. They should be collected in one place.
67 */
68 #define PRESENT_IN_REDUCED '\377'
69
70 #else
71
72 #define PRESOLVEASSERT(x) \
73 { \
74 }
75 #define PRESOLVE_STMT(s) \
76 { \
77 }
78
DIE(const char *)79 inline void DIE(const char *) {}
80
81 #endif
82
83 /*
84 Unclear why these are separate from standard debug.
85 */
86 #ifndef PRESOLVE_DETAIL
87 #define PRESOLVE_DETAIL_PRINT(s) \
88 { \
89 }
90 #else
91 #define PRESOLVE_DETAIL_PRINT(s) s
92 #endif
93
94 /*! \brief Zero tolerance
95
96 OSL had a fixed zero tolerance; we still use that here.
97 */
98 const double ZTOLDP = 1e-12;
99 /*! \brief Alternate zero tolerance
100
101 Use a different one if we are doing doubletons, etc.
102 */
103 const double ZTOLDP2 = 1e-10;
104
105 /// The usual finite infinity
106 #define PRESOLVE_INF COIN_DBL_MAX
107 /// And a small infinity
108 #define PRESOLVE_SMALL_INF 1.0e20
109 /// Check for infinity using finite infinity
110 #define PRESOLVEFINITE(n) (-PRESOLVE_INF < (n) && (n) < PRESOLVE_INF)
111
112 class CoinPostsolveMatrix;
113
114 /*! \class CoinPresolveAction
115 \brief Abstract base class of all presolve routines.
116
117 The details will make more sense after a quick overview of the grand plan:
118 A presolve object is handed a problem object, which it is expected to
119 modify in some useful way. Assuming that it succeeds, the presolve object
120 should create a postsolve object, <i>i.e.</i>, an object that contains
121 instructions for backing out the presolve transform to recover the original
122 problem. These postsolve objects are accumulated in a linked list, with each
123 successive presolve action adding its postsolve action to the head of the
124 list. The end result of all this is a presolved problem object, and a list
125 of postsolve objects. The presolved problem object is then handed to a
126 solver for optimization, and the problem object augmented with the
127 results. The list of postsolve objects is then traversed. Each of them
128 (un)modifies the problem object, with the end result being the original
129 problem, augmented with solution information.
130
131 The problem object representation is CoinPrePostsolveMatrix and subclasses.
132 Check there for details. The \c CoinPresolveAction class and subclasses
133 represent the presolve and postsolve objects.
134
135 In spite of the name, the only information held in a \c CoinPresolveAction
136 object is the information needed to postsolve (<i>i.e.</i>, the information
137 needed to back out the presolve transformation). This information is not
138 expected to change, so the fields are all \c const.
139
140 A subclass of \c CoinPresolveAction, implementing a specific pre/postsolve
141 action, is expected to declare a static function that attempts to perform a
142 presolve transformation. This function will be handed a CoinPresolveMatrix
143 to transform, and a pointer to the head of the list of postsolve objects.
144 If the transform is successful, the function will create a new
145 \c CoinPresolveAction object, link it at the head of the list of postsolve
146 objects, and return a pointer to the postsolve object it has just created.
147 Otherwise, it should return 0. It is expected that these static functions
148 will be the only things that can create new \c CoinPresolveAction objects;
149 this is expressed by making each subclass' constructor(s) private.
150
151 Every subclass must also define a \c postsolve method.
152 This function will be handed a CoinPostsolveMatrix to transform.
153
154 It is the client's responsibility to implement presolve and postsolve driver
155 routines. See OsiPresolve for examples.
156
157 \note Since the only fields in a \c CoinPresolveAction are \c const, anything
158 one can do with a variable declared \c CoinPresolveAction* can also be
159 done with a variable declared \c const \c CoinPresolveAction* It is
160 expected that all derived subclasses of \c CoinPresolveAction also have
161 this property.
162 */
163 class CoinPresolveAction {
164 public:
165 /*! \brief Stub routine to throw exceptions.
166
167 Exceptions are inefficient, particularly with g++. Even with xlC, the
168 use of exceptions adds a long prologue to a routine. Therefore, rather
169 than use throw directly in the routine, I use it in a stub routine.
170 */
throwCoinError(const char * error,const char * ps_routine)171 static void throwCoinError(const char *error, const char *ps_routine)
172 {
173 throw CoinError(error, ps_routine, "CoinPresolve");
174 }
175
176 /*! \brief The next presolve transformation
177
178 Set at object construction.
179 */
180 const CoinPresolveAction *next;
181
182 /*! \brief Construct a postsolve object and add it to the transformation list.
183
184 This is an `add to head' operation. This object will point to the
185 one passed as the parameter.
186 */
CoinPresolveAction(const CoinPresolveAction * next)187 CoinPresolveAction(const CoinPresolveAction *next)
188 : next(next)
189 {
190 }
191 /// modify next (when building rather than passing)
setNext(const CoinPresolveAction * nextAction)192 inline void setNext(const CoinPresolveAction *nextAction)
193 {
194 next = nextAction;
195 }
196
197 /*! \brief A name for debug printing.
198
199 It is expected that the name is not stored in the transform itself.
200 */
201 virtual const char *name() const = 0;
202
203 /*! \brief Apply the postsolve transformation for this particular
204 presolve action.
205 */
206 virtual void postsolve(CoinPostsolveMatrix *prob) const = 0;
207
208 /*! \brief Virtual destructor. */
~CoinPresolveAction()209 virtual ~CoinPresolveAction() {}
210 };
211
212 /*
213 These are needed for OSI-aware constructors associated with
214 CoinPrePostsolveMatrix, CoinPresolveMatrix, and CoinPostsolveMatrix.
215 */
216 class ClpSimplex;
217 class OsiSolverInterface;
218
219 /*
220 CoinWarmStartBasis is required for methods in CoinPrePostsolveMatrix
221 that accept/return a CoinWarmStartBasis object.
222 */
223 class CoinWarmStartBasis;
224
225 /*! \class CoinPrePostsolveMatrix
226 \brief Collects all the information about the problem that is needed
227 in both presolve and postsolve.
228
229 In a bit more detail, a column-major representation of the constraint
230 matrix and upper and lower bounds on variables and constraints, plus row
231 and column solutions, reduced costs, and status. There's also a set of
232 arrays holding the original row and column numbers.
233
234 As presolve and postsolve transform the matrix, it will occasionally be
235 necessary to expand the number of entries in a column. There are two
236 aspects:
237 <ul>
238 <li> During postsolve, the constraint system is expected to grow as
239 the smaller presolved system is transformed back to the original
240 system.
241 <li> During both pre- and postsolve, transforms can increase the number
242 of coefficients in a row or column. (See the
243 variable substitution, doubleton, and tripleton transforms.)
244 </ul>
245
246 The first is addressed by the members #ncols0_, #nrows0_, and #nelems0_.
247 These should be set (via constructor parameters) to values large enough
248 for the largest size taken on by the constraint system. Typically, this
249 will be the size of the original constraint system.
250
251 The second is addressed by a generous allocation of extra (empty) space
252 for the arrays used to hold coefficients and row indices. When columns
253 must be expanded, they are moved into the empty space. When it is used up,
254 the arrays are compacted. When compaction fails to produce sufficient
255 space, presolve/postsolve will fail.
256
257 CoinPrePostsolveMatrix isn't really intended to be used `bare' --- the
258 expectation is that it'll be used through CoinPresolveMatrix or
259 CoinPostsolveMatrix. Some of the functions needed to load a problem are
260 defined in the derived classes.
261
262 When CoinPresolve is applied when reoptimising, we need to be prepared to
263 accept a basis and modify it in step with the presolve actions (otherwise
264 we throw away all the advantages of warm start for reoptimization). But
265 other solution components (#acts_, #rowduals_, #sol_, and #rcosts_) are
266 needed only for postsolve, where they're used in places to determine the
267 proper action(s) when restoring rows or columns. If presolve is provided
268 with a solution, it will modify it in step with the presolve actions.
269 Moving the solution components from CoinPrePostsolveMatrix to
270 CoinPostsolveMatrix would break a lot of code. It's not clear that it's
271 worth it, and it would preclude upgrades to the presolve side that might
272 make use of any of these. -- lh, 080501 --
273
274 The constructors that take an OSI or ClpSimplex as a parameter really should
275 not be here, but for historical reasons they will likely remain for the
276 forseeable future. -- lh, 111202 --
277 */
278
279 class CoinPrePostsolveMatrix {
280 public:
281 /*! \name Constructors & Destructors */
282
283 //@{
284 /*! \brief `Native' constructor
285
286 This constructor creates an empty object which must then be loaded. On
287 the other hand, it doesn't assume that the client is an
288 OsiSolverInterface.
289 */
290 CoinPrePostsolveMatrix(int ncols_alloc, int nrows_alloc,
291 CoinBigIndex nelems_alloc);
292
293 /*! \brief Generic OSI constructor
294
295 See OSI code for the definition.
296 */
297 CoinPrePostsolveMatrix(const OsiSolverInterface *si,
298 int ncols_,
299 int nrows_,
300 CoinBigIndex nelems_);
301
302 /*! ClpOsi constructor
303
304 See Clp code for the definition.
305 */
306 CoinPrePostsolveMatrix(const ClpSimplex *si,
307 int ncols_,
308 int nrows_,
309 CoinBigIndex nelems_,
310 double bulkRatio);
311
312 /// Destructor
313 ~CoinPrePostsolveMatrix();
314 //@}
315
316 /*! \brief Enum for status of various sorts
317
318 Matches CoinWarmStartBasis::Status and adds superBasic. Most code that
319 converts between CoinPrePostsolveMatrix::Status and
320 CoinWarmStartBasis::Status will break if this correspondence is broken.
321
322 superBasic is an unresolved problem: there's no analogue in
323 CoinWarmStartBasis::Status.
324 */
325 enum Status {
326 isFree = 0x00,
327 basic = 0x01,
328 atUpperBound = 0x02,
329 atLowerBound = 0x03,
330 superBasic = 0x04
331 };
332
333 /*! \name Functions to work with variable status
334
335 Functions to work with the CoinPrePostsolveMatrix::Status enum and
336 related vectors.
337
338 \todo
339 Why are we futzing around with three bit status? A holdover from the
340 packed arrays of CoinWarmStartBasis? Big swaths of the presolve code
341 manipulates colstat_ and rowstat_ as unsigned char arrays using simple
342 assignment to set values.
343 */
344 //@{
345
346 /// Set row status (<i>i.e.</i>, status of artificial for this row)
setRowStatus(int sequence,Status status)347 inline void setRowStatus(int sequence, Status status)
348 {
349 unsigned char &st_byte = rowstat_[sequence];
350 st_byte = static_cast< unsigned char >(st_byte & (~7));
351 st_byte = static_cast< unsigned char >(st_byte | status);
352 }
353 /// Get row status
getRowStatus(int sequence) const354 inline Status getRowStatus(int sequence) const
355 {
356 return static_cast< Status >(rowstat_[sequence] & 7);
357 }
358 /// Check if artificial for this row is basic
rowIsBasic(int sequence) const359 inline bool rowIsBasic(int sequence) const
360 {
361 return (static_cast< Status >(rowstat_[sequence] & 7) == basic);
362 }
363 /// Set column status (<i>i.e.</i>, status of primal variable)
setColumnStatus(int sequence,Status status)364 inline void setColumnStatus(int sequence, Status status)
365 {
366 unsigned char &st_byte = colstat_[sequence];
367 st_byte = static_cast< unsigned char >(st_byte & (~7));
368 st_byte = static_cast< unsigned char >(st_byte | status);
369
370 #ifdef PRESOLVE_DEBUG
371 switch (status) {
372 case isFree: {
373 if (clo_[sequence] > -PRESOLVE_INF || cup_[sequence] < PRESOLVE_INF) {
374 std::cout << "Bad status: Var " << sequence
375 << " isFree, lb = " << clo_[sequence]
376 << ", ub = " << cup_[sequence] << std::endl;
377 }
378 break;
379 }
380 case basic: {
381 break;
382 }
383 case atUpperBound: {
384 if (cup_[sequence] >= PRESOLVE_INF) {
385 std::cout << "Bad status: Var " << sequence
386 << " atUpperBound, lb = " << clo_[sequence]
387 << ", ub = " << cup_[sequence] << std::endl;
388 }
389 break;
390 }
391 case atLowerBound: {
392 if (clo_[sequence] <= -PRESOLVE_INF) {
393 std::cout << "Bad status: Var " << sequence
394 << " atLowerBound, lb = " << clo_[sequence]
395 << ", ub = " << cup_[sequence] << std::endl;
396 }
397 break;
398 }
399 case superBasic: {
400 if (clo_[sequence] <= -PRESOLVE_INF && cup_[sequence] >= PRESOLVE_INF) {
401 std::cout << "Bad status: Var " << sequence
402 << " superBasic, lb = " << clo_[sequence]
403 << ", ub = " << cup_[sequence] << std::endl;
404 }
405 break;
406 }
407 default: {
408 assert(false);
409 break;
410 }
411 }
412 #endif
413 }
414 /// Get column (structural variable) status
getColumnStatus(int sequence) const415 inline Status getColumnStatus(int sequence) const
416 {
417 return static_cast< Status >(colstat_[sequence] & 7);
418 }
419 /// Check if column (structural variable) is basic
columnIsBasic(int sequence) const420 inline bool columnIsBasic(int sequence) const
421 {
422 return (static_cast< Status >(colstat_[sequence] & 7) == basic);
423 }
424 /*! \brief Set status of row (artificial variable) to the correct nonbasic
425 status given bounds and current value
426 */
427 void setRowStatusUsingValue(int iRow);
428 /*! \brief Set status of column (structural variable) to the correct
429 nonbasic status given bounds and current value
430 */
431 void setColumnStatusUsingValue(int iColumn);
432 /*! \brief Set column (structural variable) status vector */
433 void setStructuralStatus(const char *strucStatus, int lenParam);
434 /*! \brief Set row (artificial variable) status vector */
435 void setArtificialStatus(const char *artifStatus, int lenParam);
436 /*! \brief Set the status of all variables from a basis */
437 void setStatus(const CoinWarmStartBasis *basis);
438 /*! \brief Get status in the form of a CoinWarmStartBasis */
439 CoinWarmStartBasis *getStatus();
440 /*! \brief Return a print string for status of a column (structural
441 variable)
442 */
443 const char *columnStatusString(int j) const;
444 /*! \brief Return a print string for status of a row (artificial
445 variable)
446 */
447 const char *rowStatusString(int i) const;
448 //@}
449
450 /*! \name Functions to load problem and solution information
451
452 These functions can be used to load portions of the problem definition
453 and solution. See also the CoinPresolveMatrix and CoinPostsolveMatrix
454 classes.
455 */
456 //@{
457 /// Set the objective function offset for the original system.
458 void setObjOffset(double offset);
459 /*! \brief Set the objective sense (max/min)
460
461 Coded as 1.0 for min, -1.0 for max.
462 Yes, there's a method, and a matching attribute. No, you really
463 don't want to set this to maximise.
464 */
465 void setObjSense(double objSense);
466 /// Set the primal feasibility tolerance
467 void setPrimalTolerance(double primTol);
468 /// Set the dual feasibility tolerance
469 void setDualTolerance(double dualTol);
470 /// Set column lower bounds
471 void setColLower(const double *colLower, int lenParam);
472 /// Set column upper bounds
473 void setColUpper(const double *colUpper, int lenParam);
474 /// Set column solution
475 void setColSolution(const double *colSol, int lenParam);
476 /// Set objective coefficients
477 void setCost(const double *cost, int lenParam);
478 /// Set reduced costs
479 void setReducedCost(const double *redCost, int lenParam);
480 /// Set row lower bounds
481 void setRowLower(const double *rowLower, int lenParam);
482 /// Set row upper bounds
483 void setRowUpper(const double *rowUpper, int lenParam);
484 /// Set row solution
485 void setRowPrice(const double *rowSol, int lenParam);
486 /// Set row activity
487 void setRowActivity(const double *rowAct, int lenParam);
488 //@}
489
490 /*! \name Functions to retrieve problem and solution information */
491 //@{
492 /// Get current number of columns
getNumCols() const493 inline int getNumCols() const
494 {
495 return (ncols_);
496 }
497 /// Get current number of rows
getNumRows() const498 inline int getNumRows() const
499 {
500 return (nrows_);
501 }
502 /// Get current number of non-zero coefficients
getNumElems() const503 inline CoinBigIndex getNumElems() const
504 {
505 return (nelems_);
506 }
507 /// Get column start vector for column-major packed matrix
getColStarts() const508 inline const CoinBigIndex *getColStarts() const
509 {
510 return (mcstrt_);
511 }
512 /// Get column length vector for column-major packed matrix
getColLengths() const513 inline const int *getColLengths() const
514 {
515 return (hincol_);
516 }
517 /// Get vector of row indices for column-major packed matrix
getRowIndicesByCol() const518 inline const int *getRowIndicesByCol() const
519 {
520 return (hrow_);
521 }
522 /// Get vector of elements for column-major packed matrix
getElementsByCol() const523 inline const double *getElementsByCol() const
524 {
525 return (colels_);
526 }
527 /// Get column lower bounds
getColLower() const528 inline const double *getColLower() const
529 {
530 return (clo_);
531 }
532 /// Get column upper bounds
getColUpper() const533 inline const double *getColUpper() const
534 {
535 return (cup_);
536 }
537 /// Get objective coefficients
getCost() const538 inline const double *getCost() const
539 {
540 return (cost_);
541 }
542 /// Get row lower bounds
getRowLower() const543 inline const double *getRowLower() const
544 {
545 return (rlo_);
546 }
547 /// Get row upper bounds
getRowUpper() const548 inline const double *getRowUpper() const
549 {
550 return (rup_);
551 }
552 /// Get column solution (primal variable values)
getColSolution() const553 inline const double *getColSolution() const
554 {
555 return (sol_);
556 }
557 /// Get row activity (constraint lhs values)
getRowActivity() const558 inline const double *getRowActivity() const
559 {
560 return (acts_);
561 }
562 /// Get row solution (dual variables)
getRowPrice() const563 inline const double *getRowPrice() const
564 {
565 return (rowduals_);
566 }
567 /// Get reduced costs
getReducedCost() const568 inline const double *getReducedCost() const
569 {
570 return (rcosts_);
571 }
572 /// Count empty columns
countEmptyCols()573 inline int countEmptyCols()
574 {
575 int empty = 0;
576 for (int i = 0; i < ncols_; i++)
577 if (hincol_[i] == 0)
578 empty++;
579 return (empty);
580 }
581 //@}
582
583 /*! \name Message handling */
584 //@{
585 /// Return message handler
messageHandler() const586 inline CoinMessageHandler *messageHandler() const
587 {
588 return handler_;
589 }
590 /*! \brief Set message handler
591
592 The client retains responsibility for the handler --- it will not be
593 destroyed with the \c CoinPrePostsolveMatrix object.
594 */
setMessageHandler(CoinMessageHandler * handler)595 inline void setMessageHandler(CoinMessageHandler *handler)
596 {
597 if (defaultHandler_ == true) {
598 delete handler_;
599 defaultHandler_ = false;
600 }
601 handler_ = handler;
602 }
603 /// Return messages
messages() const604 inline CoinMessages messages() const
605 {
606 return messages_;
607 }
608 //@}
609
610 /*! \name Current and Allocated Size
611
612 During pre- and postsolve, the matrix will change in size. During presolve
613 it will shrink; during postsolve it will grow. Hence there are two sets of
614 size variables, one for the current size and one for the allocated size.
615 (See the general comments for the CoinPrePostsolveMatrix class for more
616 information.)
617 */
618 //@{
619
620 /// current number of columns
621 int ncols_;
622 /// current number of rows
623 int nrows_;
624 /// current number of coefficients
625 CoinBigIndex nelems_;
626
627 /// Allocated number of columns
628 int ncols0_;
629 /// Allocated number of rows
630 int nrows0_;
631 /// Allocated number of coefficients
632 CoinBigIndex nelems0_;
633 /*! \brief Allocated size of bulk storage for row indices and coefficients
634
635 This is the space allocated for hrow_ and colels_. This must be large
636 enough to allow columns to be copied into empty space when they need to
637 be expanded. For efficiency (to minimize the number of times the
638 representation must be compressed) it's recommended that this be at least
639 2*nelems0_.
640 */
641 CoinBigIndex bulk0_;
642 /// Ratio of bulk0_ to nelems0_; default is 2.
643 double bulkRatio_;
644 //@}
645
646 /*! \name Problem representation
647
648 The matrix is the common column-major format: A pair of vectors with
649 positional correspondence to hold coefficients and row indices, and a
650 second pair of vectors giving the starting position and length of each
651 column in the first pair.
652 */
653 //@{
654 /// Vector of column start positions in #hrow_, #colels_
655 CoinBigIndex *mcstrt_;
656 /// Vector of column lengths
657 int *hincol_;
658 /// Row indices (positional correspondence with #colels_)
659 int *hrow_;
660 /// Coefficients (positional correspondence with #hrow_)
661 double *colels_;
662
663 /// Objective coefficients
664 double *cost_;
665 /// Original objective offset
666 double originalOffset_;
667
668 /// Column (primal variable) lower bounds
669 double *clo_;
670 /// Column (primal variable) upper bounds
671 double *cup_;
672
673 /// Row (constraint) lower bounds
674 double *rlo_;
675 /// Row (constraint) upper bounds
676 double *rup_;
677
678 /*! \brief Original column numbers
679
680 Over the current range of column numbers in the presolved problem,
681 the entry for column j will contain the index of the corresponding
682 column in the original problem.
683 */
684 int *originalColumn_;
685 /*! \brief Original row numbers
686
687 Over the current range of row numbers in the presolved problem, the
688 entry for row i will contain the index of the corresponding row in
689 the original problem.
690 */
691 int *originalRow_;
692
693 /// Primal feasibility tolerance
694 double ztolzb_;
695 /// Dual feasibility tolerance
696 double ztoldj_;
697
698 /*! \brief Maximization/minimization
699
700 Yes, there's a variable here. No, you really don't want to set this to
701 maximise. See the main notes for CoinPresolveMatrix.
702 */
703 double maxmin_;
704 //@}
705
706 /*! \name Problem solution information
707
708 The presolve phase will work without any solution information
709 (appropriate for initial optimisation) or with solution information
710 (appropriate for reoptimisation). When solution information is supplied,
711 presolve will maintain it to the best of its ability. #colstat_ is
712 checked to determine the presence/absence of status information. #sol_ is
713 checked for primal solution information, and #rowduals_ for dual solution
714 information.
715
716 The postsolve phase requires the complete solution information from the
717 presolved problem (status, primal and dual solutions). It will be
718 transformed into a correct solution for the original problem.
719 */
720 //@{
721 /*! \brief Vector of primal variable values
722
723 If #sol_ exists, it is assumed that primal solution information should be
724 updated and that #acts_ also exists.
725 */
726 double *sol_;
727 /*! \brief Vector of dual variable values
728
729 If #rowduals_ exists, it is assumed that dual solution information should
730 be updated and that #rcosts_ also exists.
731 */
732 double *rowduals_;
733 /*! \brief Vector of constraint left-hand-side values (row activity)
734
735 Produced by evaluating constraints according to #sol_. Updated iff
736 #sol_ exists.
737 */
738 double *acts_;
739 /*! \brief Vector of reduced costs
740
741 Produced by evaluating dual constraints according to #rowduals_. Updated
742 iff #rowduals_ exists.
743 */
744 double *rcosts_;
745
746 /*! \brief Status of primal variables
747
748 Coded with CoinPrePostSolveMatrix::Status, one code per char. colstat_ and
749 #rowstat_ <b>MUST</b> be allocated as a single vector. This is to maintain
750 compatibility with ClpPresolve and OsiPresolve, which do it this way.
751 */
752 unsigned char *colstat_;
753
754 /*! \brief Status of constraints
755
756 More accurately, the status of the logical variable associated with the
757 constraint. Coded with CoinPrePostSolveMatrix::Status, one code per char.
758 Note that this must be allocated as a single vector with #colstat_.
759 */
760 unsigned char *rowstat_;
761 //@}
762
763 /*! \name Message handling
764
765 Uses the standard COIN approach: a default handler is installed, and the
766 CoinPrePostsolveMatrix object takes responsibility for it. If the client
767 replaces the handler with one of their own, it becomes their
768 responsibility.
769 */
770 //@{
771 /// Message handler
772 CoinMessageHandler *handler_;
773 /// Indicates if the current #handler_ is default (true) or not (false).
774 bool defaultHandler_;
775 /// Standard COIN messages
776 CoinMessage messages_;
777 //@}
778 };
779
780 /*! \relates CoinPrePostsolveMatrix
781 \brief Generate a print string for a status code.
782 */
783 const char *statusName(CoinPrePostsolveMatrix::Status status);
784
785 /*! \class presolvehlink
786 \brief Links to aid in packed matrix modification
787
788 Currently, the matrices held by the CoinPrePostsolveMatrix and
789 CoinPresolveMatrix objects are represented in the same way as a
790 CoinPackedMatrix. In the course of presolve and postsolve transforms, it
791 will happen that a major-dimension vector needs to increase in size. In
792 order to check whether there is enough room to add another coefficient in
793 place, it helps to know the next vector (in memory order) in the bulk
794 storage area. To do that, a linked list of major-dimension vectors is
795 maintained; the "pre" and "suc" fields give the previous and next vector,
796 in memory order (that is, the vector whose mcstrt_ or mrstrt_ entry is
797 next smaller or larger).
798
799 Consider a column-major matrix with ncols columns. By definition,
800 presolvehlink[ncols].pre points to the column in the last occupied
801 position of the bulk storage arrays. There is no easy way to find the
802 column which occupies the first position (there is no presolvehlink[-1] to
803 consult). If the column that initially occupies the first position is
804 moved for expansion, there is no way to reclaim the space until the bulk
805 storage is compacted. The same holds for the last and first rows of a
806 row-major matrix, of course.
807 */
808
809 class presolvehlink {
810 public:
811 int pre, suc;
812 };
813
814 #define NO_LINK -66666666
815
816 /*! \relates presolvehlink
817 \brief unlink vector i
818
819 Remove vector i from the ordering.
820 */
PRESOLVE_REMOVE_LINK(presolvehlink * link,int i)821 inline void PRESOLVE_REMOVE_LINK(presolvehlink *link, int i)
822 {
823 int ipre = link[i].pre;
824 int isuc = link[i].suc;
825 if (ipre >= 0) {
826 link[ipre].suc = isuc;
827 }
828 if (isuc >= 0) {
829 link[isuc].pre = ipre;
830 }
831 link[i].pre = NO_LINK, link[i].suc = NO_LINK;
832 }
833
834 /*! \relates presolvehlink
835 \brief insert vector i after vector j
836
837 Insert vector i between j and j.suc.
838 */
PRESOLVE_INSERT_LINK(presolvehlink * link,int i,int j)839 inline void PRESOLVE_INSERT_LINK(presolvehlink *link, int i, int j)
840 {
841 int isuc = link[j].suc;
842 link[j].suc = i;
843 link[i].pre = j;
844 if (isuc >= 0) {
845 link[isuc].pre = i;
846 }
847 link[i].suc = isuc;
848 }
849
850 /*! \relates presolvehlink
851 \brief relink vector j in place of vector i
852
853 Replace vector i in the ordering with vector j. This is equivalent to
854 <pre>
855 int pre = link[i].pre;
856 PRESOLVE_REMOVE_LINK(link,i);
857 PRESOLVE_INSERT_LINK(link,j,pre);
858 </pre>
859 But, this routine will work even if i happens to be first in the order.
860 */
PRESOLVE_MOVE_LINK(presolvehlink * link,int i,int j)861 inline void PRESOLVE_MOVE_LINK(presolvehlink *link, int i, int j)
862 {
863 int ipre = link[i].pre;
864 int isuc = link[i].suc;
865 if (ipre >= 0) {
866 link[ipre].suc = j;
867 }
868 if (isuc >= 0) {
869 link[isuc].pre = j;
870 }
871 link[i].pre = NO_LINK, link[i].suc = NO_LINK;
872 }
873
874 /*! \class CoinPresolveMatrix
875 \brief Augments CoinPrePostsolveMatrix with information about the problem
876 that is only needed during presolve.
877
878 For problem manipulation, this class adds a row-major matrix
879 representation, linked lists that allow for easy manipulation of the matrix
880 when applying presolve transforms, and vectors to track row and column
881 processing status (changed, needs further processing, change prohibited)
882
883 For problem representation, this class adds information about variable type
884 (integer or continuous), an objective offset, and a feasibility tolerance.
885
886 <b>NOTE</b> that the #anyInteger_ and #anyProhibited_ flags are independent
887 of the vectors used to track this information for individual variables
888 (#integerType_ and #rowChanged_ and #colChanged_, respectively).
889
890 <b>NOTE</b> also that at the end of presolve the column-major and row-major
891 matrix representations are loosely packed (<i>i.e.</i>, there may be gaps
892 between columns in the bulk storage arrays).
893
894 <b>NOTE</b> that while you might think that CoinPresolve is prepared to
895 handle minimisation or maximisation, it's unlikely that this still works.
896 This is a good thing: better to convert objective coefficients and duals
897 once, before starting presolve, rather than doing it over and over in
898 each transform that considers dual variables.
899
900 The constructors that take an OSI or ClpSimplex as a parameter really should
901 not be here, but for historical reasons they will likely remain for the
902 forseeable future. -- lh, 111202 --
903 */
904
905 class CoinPresolveMatrix : public CoinPrePostsolveMatrix {
906 public:
907 /*! \brief `Native' constructor
908
909 This constructor creates an empty object which must then be loaded.
910 On the other hand, it doesn't assume that the client is an
911 OsiSolverInterface.
912 */
913 CoinPresolveMatrix(int ncols_alloc, int nrows_alloc,
914 CoinBigIndex nelems_alloc);
915
916 /*! \brief Clp OSI constructor
917
918 See Clp code for the definition.
919 */
920 CoinPresolveMatrix(int ncols0,
921 double maxmin,
922 // end prepost members
923
924 ClpSimplex *si,
925
926 // rowrep
927 int nrows,
928 CoinBigIndex nelems,
929 bool doStatus,
930 double nonLinearVariable,
931 double bulkRatio);
932
933 /*! \brief Update the model held by a Clp OSI */
934 void update_model(ClpSimplex *si,
935 int nrows0,
936 int ncols0,
937 CoinBigIndex nelems0);
938 /*! \brief Generic OSI constructor
939
940 See OSI code for the definition.
941 */
942 CoinPresolveMatrix(int ncols0,
943 double maxmin,
944 // end prepost members
945 OsiSolverInterface *si,
946 // rowrep
947 int nrows,
948 CoinBigIndex nelems,
949 bool doStatus,
950 double nonLinearVariable,
951 const char *prohibited,
952 const char *rowProhibited = NULL);
953
954 /*! \brief Update the model held by a generic OSI */
955 void update_model(OsiSolverInterface *si,
956 int nrows0,
957 int ncols0,
958 CoinBigIndex nelems0);
959
960 /// Destructor
961 ~CoinPresolveMatrix();
962
963 /*! \brief Initialize a CoinPostsolveMatrix object, destroying the
964 CoinPresolveMatrix object.
965
966 See CoinPostsolveMatrix::assignPresolveToPostsolve.
967 */
968 friend void assignPresolveToPostsolve(CoinPresolveMatrix *&preObj);
969
970 /*! \name Functions to load the problem representation
971 */
972 //@{
973 /*! \brief Load the cofficient matrix.
974
975 Load the coefficient matrix before loading the other vectors (bounds,
976 objective, variable type) required to define the problem.
977 */
978 void setMatrix(const CoinPackedMatrix *mtx);
979
980 /// Count number of empty rows
countEmptyRows()981 inline int countEmptyRows()
982 {
983 int empty = 0;
984 for (int i = 0; i < nrows_; i++)
985 if (hinrow_[i] == 0)
986 empty++;
987 return (empty);
988 }
989
990 /*! \brief Set variable type information for a single variable
991
992 Set \p variableType to 0 for continous, 1 for integer.
993 Does not manipulate the #anyInteger_ flag.
994 */
setVariableType(int i,int variableType)995 inline void setVariableType(int i, int variableType)
996 {
997 if (integerType_ == 0)
998 integerType_ = new unsigned char[ncols0_];
999 integerType_[i] = static_cast< unsigned char >(variableType);
1000 }
1001
1002 /*! \brief Set variable type information for all variables
1003
1004 Set \p variableType[i] to 0 for continuous, 1 for integer.
1005 Does not manipulate the #anyInteger_ flag.
1006 */
1007 void setVariableType(const unsigned char *variableType, int lenParam);
1008
1009 /*! \brief Set the type of all variables
1010
1011 allIntegers should be true to set the type to integer, false to set the
1012 type to continuous.
1013 */
1014 void setVariableType(bool allIntegers, int lenParam);
1015
1016 /// Set a flag for presence (true) or absence (false) of integer variables
setAnyInteger(bool anyInteger=true)1017 inline void setAnyInteger(bool anyInteger = true)
1018 {
1019 anyInteger_ = anyInteger;
1020 }
1021 //@}
1022
1023 /*! \name Functions to retrieve problem information
1024 */
1025 //@{
1026
1027 /// Get row start vector for row-major packed matrix
getRowStarts() const1028 inline const CoinBigIndex *getRowStarts() const
1029 {
1030 return (mrstrt_);
1031 }
1032 /// Get vector of column indices for row-major packed matrix
getColIndicesByRow() const1033 inline const int *getColIndicesByRow() const
1034 {
1035 return (hcol_);
1036 }
1037 /// Get vector of elements for row-major packed matrix
getElementsByRow() const1038 inline const double *getElementsByRow() const
1039 {
1040 return (rowels_);
1041 }
1042
1043 /*! \brief Check for integrality of the specified variable.
1044
1045 Consults the #integerType_ vector if present; fallback is the
1046 #anyInteger_ flag.
1047 */
isInteger(int i) const1048 inline bool isInteger(int i) const
1049 {
1050 if (integerType_ == 0) {
1051 return (anyInteger_);
1052 } else if (integerType_[i] == 1) {
1053 return (true);
1054 } else {
1055 return (false);
1056 }
1057 }
1058
1059 /*! \brief Check if there are any integer variables
1060
1061 Consults the #anyInteger_ flag
1062 */
anyInteger() const1063 inline bool anyInteger() const
1064 {
1065 return (anyInteger_);
1066 }
1067 /// Picks up any special options
presolveOptions() const1068 inline int presolveOptions() const
1069 {
1070 return presolveOptions_;
1071 }
1072 /// Sets any special options (see #presolveOptions_)
setPresolveOptions(int value)1073 inline void setPresolveOptions(int value)
1074 {
1075 presolveOptions_ = value;
1076 }
1077 //@}
1078
1079 /*! \name Matrix storage management links
1080
1081 Linked lists, modelled after the linked lists used in OSL
1082 factorization. They are used for management of the bulk coefficient
1083 and minor index storage areas.
1084 */
1085 //@{
1086 /// Linked list for the column-major representation.
1087 presolvehlink *clink_;
1088 /// Linked list for the row-major representation.
1089 presolvehlink *rlink_;
1090 //@}
1091
1092 /// Objective function offset introduced during presolve
1093 double dobias_;
1094
1095 /// Adjust objective function constant offset
change_bias(double change_amount)1096 inline void change_bias(double change_amount)
1097 {
1098 dobias_ += change_amount;
1099 #if PRESOLVE_DEBUG > 2
1100 assert(fabs(change_amount) < 1.0e50);
1101 if (change_amount)
1102 PRESOLVE_STMT(printf("changing bias by %g to %g\n",
1103 change_amount, dobias_));
1104 #endif
1105 }
1106
1107 /*! \name Row-major representation
1108
1109 Common row-major format: A pair of vectors with positional
1110 correspondence to hold coefficients and column indices, and a second pair
1111 of vectors giving the starting position and length of each row in
1112 the first pair.
1113 */
1114 //@{
1115 /// Vector of row start positions in #hcol, #rowels_
1116 CoinBigIndex *mrstrt_;
1117 /// Vector of row lengths
1118 int *hinrow_;
1119 /// Coefficients (positional correspondence with #hcol_)
1120 double *rowels_;
1121 /// Column indices (positional correspondence with #rowels_)
1122 int *hcol_;
1123 //@}
1124
1125 /// Tracks integrality of columns (1 for integer, 0 for continuous)
1126 unsigned char *integerType_;
1127 /*! \brief Flag to say if any variables are integer
1128
1129 Note that this flag is <i>not</i> manipulated by the various
1130 \c setVariableType routines.
1131 */
1132 bool anyInteger_;
1133 /// Print statistics for tuning
1134 bool tuning_;
1135 /// Say we want statistics - also set time
1136 void statistics();
1137 /// Start time of presolve
1138 double startTime_;
1139
1140 /// Bounds can be moved by this to retain feasibility
1141 double feasibilityTolerance_;
1142 /// Return feasibility tolerance
feasibilityTolerance()1143 inline double feasibilityTolerance()
1144 {
1145 return (feasibilityTolerance_);
1146 }
1147 /// Set feasibility tolerance
setFeasibilityTolerance(double val)1148 inline void setFeasibilityTolerance(double val)
1149 {
1150 feasibilityTolerance_ = val;
1151 }
1152
1153 /*! \brief Output status: 0 = feasible, 1 = infeasible, 2 = unbounded
1154
1155 Actually implemented as single bit flags: 1^0 = infeasible, 1^1 =
1156 unbounded.
1157 */
1158 int status_;
1159 /// Returns problem status (0 = feasible, 1 = infeasible, 2 = unbounded)
status()1160 inline int status()
1161 {
1162 return (status_);
1163 }
1164 /// Set problem status
setStatus(int status)1165 inline void setStatus(int status)
1166 {
1167 status_ = (status & 0x3);
1168 }
1169
1170 /*! \brief Presolve pass number
1171
1172 Should be incremented externally by the method controlling application of
1173 presolve transforms.
1174 Used to control the execution of testRedundant (evoked by the
1175 implied_free transform).
1176 */
1177 int pass_;
1178 /// Set pass number
setPass(int pass=0)1179 inline void setPass(int pass = 0)
1180 {
1181 pass_ = pass;
1182 }
1183
1184 /*! \brief Maximum substitution level
1185
1186 Used to control the execution of subst from implied_free
1187 */
1188 int maxSubstLevel_;
1189 /// Set Maximum substitution level (normally 3)
setMaximumSubstitutionLevel(int level)1190 inline void setMaximumSubstitutionLevel(int level)
1191 {
1192 maxSubstLevel_ = level;
1193 }
1194
1195 /*! \name Row and column processing status
1196
1197 Information used to determine if rows or columns can be changed and
1198 if they require further processing due to changes.
1199
1200 There are four major lists: the [row,col]ToDo list, and the
1201 [row,col]NextToDo list. In general, a transform processes entries from
1202 the ToDo list and adds entries to the NextToDo list.
1203
1204 There are two vectors, [row,col]Changed, which track the status of
1205 individual rows and columns.
1206 */
1207 //@{
1208 /*! \brief Column change status information
1209
1210 Coded using the following bits:
1211 <ul>
1212 <li> 0x01: Column has changed
1213 <li> 0x02: preprocessing prohibited
1214 <li> 0x04: Column has been used
1215 <li> 0x08: Column originally had infinite ub
1216 </ul>
1217 */
1218 unsigned char *colChanged_;
1219 /// Input list of columns to process
1220 int *colsToDo_;
1221 /// Length of #colsToDo_
1222 int numberColsToDo_;
1223 /// Output list of columns to process next
1224 int *nextColsToDo_;
1225 /// Length of #nextColsToDo_
1226 int numberNextColsToDo_;
1227
1228 /*! \brief Row change status information
1229
1230 Coded using the following bits:
1231 <ul>
1232 <li> 0x01: Row has changed
1233 <li> 0x02: preprocessing prohibited
1234 <li> 0x04: Row has been used
1235 </ul>
1236 */
1237 unsigned char *rowChanged_;
1238 /// Input list of rows to process
1239 int *rowsToDo_;
1240 /// Length of #rowsToDo_
1241 int numberRowsToDo_;
1242 /// Output list of rows to process next
1243 int *nextRowsToDo_;
1244 /// Length of #nextRowsToDo_
1245 int numberNextRowsToDo_;
1246 /*! \brief Fine control over presolve actions
1247
1248 Set/clear the following bits to allow or suppress actions:
1249 - 0x01 allow duplicate column tests for integer variables
1250 - 0x02 not used
1251 - 0x04 set to inhibit x+y+z=1 mods
1252 - 0x08 not used
1253 - 0x10 set to allow stuff which won't unroll easily (overlapping
1254 duplicate rows; opportunistic fixing of variables from bound
1255 propagation).
1256 - 0x04000 allow presolve transforms to arbitrarily ignore infeasibility
1257 and set arbitrary feasible bounds.
1258 - 0x10000 instructs implied_free_action to be `more lightweight'; will
1259 return without doing anything after 15 presolve passes.
1260 - 0x(2,4,6)0000 instructs implied_free_action to remove small created elements
1261 - 0x80000000 set by presolve to say dupcol_action compressed columns
1262 */
1263 int presolveOptions_;
1264 /*! Flag to say if any rows or columns are marked as prohibited
1265
1266 Note that this flag is <i>not</i> manipulated by any of the
1267 various \c set*Prohibited routines.
1268 */
1269 bool anyProhibited_;
1270 //@}
1271
1272 /*! \name Scratch work arrays
1273
1274 Preallocated work arrays are useful to avoid having to allocate and free
1275 work arrays in individual presolve methods.
1276
1277 All are allocated from #setMatrix by #initializeStuff, freed from
1278 #~CoinPresolveMatrix. You can use #deleteStuff followed by
1279 #initializeStuff to remove and recreate them.
1280 */
1281 //@{
1282 /// Preallocated scratch work array, 3*nrows_
1283 int *usefulRowInt_;
1284 /// Preallocated scratch work array, 2*nrows_
1285 double *usefulRowDouble_;
1286 /// Preallocated scratch work array, 2*ncols_
1287 int *usefulColumnInt_;
1288 /// Preallocated scratch work array, ncols_
1289 double *usefulColumnDouble_;
1290 /// Array of random numbers (max row,column)
1291 double *randomNumber_;
1292
1293 /// Work array for count of infinite contributions to row lhs upper bound
1294 int *infiniteUp_;
1295 /// Work array for sum of finite contributions to row lhs upper bound
1296 double *sumUp_;
1297 /// Work array for count of infinite contributions to row lhs lower bound
1298 int *infiniteDown_;
1299 /// Work array for sum of finite contributions to row lhs lower bound
1300 double *sumDown_;
1301 //@}
1302
1303 /*! \brief Recompute row lhs bounds
1304
1305 Calculate finite contributions to row lhs upper and lower bounds
1306 and count infinite contributions. Returns the number of rows found
1307 to be infeasible.
1308
1309 If \p whichRow < 0, bounds are recomputed for all rows.
1310
1311 As of 110611, this seems to be a work in progress in the sense that it's
1312 barely used by the existing presolve code.
1313 */
1314 int recomputeSums(int whichRow);
1315
1316 /// Allocate scratch arrays
1317 void initializeStuff();
1318 /// Free scratch arrays
1319 void deleteStuff();
1320
1321 /*! \name Functions to manipulate row and column processing status */
1322 //@{
1323
1324 /*! \brief Initialise the column ToDo lists
1325
1326 Places all columns in the #colsToDo_ list except for columns marked
1327 as prohibited (<i>viz.</i> #colChanged_).
1328 */
1329 void initColsToDo();
1330
1331 /*! \brief Step column ToDo lists
1332
1333 Moves columns on the #nextColsToDo_ list to the #colsToDo_ list, emptying
1334 #nextColsToDo_. Returns the number of columns transferred.
1335 */
1336 int stepColsToDo();
1337
1338 /// Return the number of columns on the #colsToDo_ list
numberColsToDo()1339 inline int numberColsToDo()
1340 {
1341 return (numberColsToDo_);
1342 }
1343
1344 /// Has column been changed?
colChanged(int i) const1345 inline bool colChanged(int i) const
1346 {
1347 return (colChanged_[i] & 1) != 0;
1348 }
1349 /// Mark column as not changed
unsetColChanged(int i)1350 inline void unsetColChanged(int i)
1351 {
1352 colChanged_[i] = static_cast< unsigned char >(colChanged_[i] & (~1));
1353 }
1354 /// Mark column as changed.
setColChanged(int i)1355 inline void setColChanged(int i)
1356 {
1357 colChanged_[i] = static_cast< unsigned char >(colChanged_[i] | (1));
1358 }
1359 /// Mark column as changed and add to list of columns to process next
addCol(int i)1360 inline void addCol(int i)
1361 {
1362 if ((colChanged_[i] & 1) == 0) {
1363 colChanged_[i] = static_cast< unsigned char >(colChanged_[i] | (1));
1364 nextColsToDo_[numberNextColsToDo_++] = i;
1365 }
1366 }
1367 /// Test if column is eligible for preprocessing
colProhibited(int i) const1368 inline bool colProhibited(int i) const
1369 {
1370 return (colChanged_[i] & 2) != 0;
1371 }
1372 /*! \brief Test if column is eligible for preprocessing
1373
1374 The difference between this method and #colProhibited() is that this
1375 method first tests #anyProhibited_ before examining the specific entry
1376 for the specified column.
1377 */
colProhibited2(int i) const1378 inline bool colProhibited2(int i) const
1379 {
1380 if (!anyProhibited_)
1381 return false;
1382 else
1383 return (colChanged_[i] & 2) != 0;
1384 }
1385 /// Mark column as ineligible for preprocessing
setColProhibited(int i)1386 inline void setColProhibited(int i)
1387 {
1388 colChanged_[i] = static_cast< unsigned char >(colChanged_[i] | (2));
1389 }
1390 /*! \brief Test if column is marked as used
1391
1392 This is for doing faster lookups to see where two columns have entries
1393 in common.
1394 */
colUsed(int i) const1395 inline bool colUsed(int i) const
1396 {
1397 return (colChanged_[i] & 4) != 0;
1398 }
1399 /// Mark column as used
setColUsed(int i)1400 inline void setColUsed(int i)
1401 {
1402 colChanged_[i] = static_cast< unsigned char >(colChanged_[i] | (4));
1403 }
1404 /// Mark column as unused
unsetColUsed(int i)1405 inline void unsetColUsed(int i)
1406 {
1407 colChanged_[i] = static_cast< unsigned char >(colChanged_[i] & (~4));
1408 }
1409 /// Has column infinite ub (originally)
colInfinite(int i) const1410 inline bool colInfinite(int i) const
1411 {
1412 return (colChanged_[i] & 8) != 0;
1413 }
1414 /// Mark column as not infinite ub (originally)
unsetColInfinite(int i)1415 inline void unsetColInfinite(int i)
1416 {
1417 colChanged_[i] = static_cast< unsigned char >(colChanged_[i] & (~8));
1418 }
1419 /// Mark column as infinite ub (originally)
setColInfinite(int i)1420 inline void setColInfinite(int i)
1421 {
1422 colChanged_[i] = static_cast< unsigned char >(colChanged_[i] | (8));
1423 }
1424
1425 /*! \brief Initialise the row ToDo lists
1426
1427 Places all rows in the #rowsToDo_ list except for rows marked
1428 as prohibited (<i>viz.</i> #rowChanged_).
1429 */
1430 void initRowsToDo();
1431
1432 /*! \brief Step row ToDo lists
1433
1434 Moves rows on the #nextRowsToDo_ list to the #rowsToDo_ list, emptying
1435 #nextRowsToDo_. Returns the number of rows transferred.
1436 */
1437 int stepRowsToDo();
1438
1439 /// Return the number of rows on the #rowsToDo_ list
numberRowsToDo()1440 inline int numberRowsToDo()
1441 {
1442 return (numberRowsToDo_);
1443 }
1444
1445 /// Has row been changed?
rowChanged(int i) const1446 inline bool rowChanged(int i) const
1447 {
1448 return (rowChanged_[i] & 1) != 0;
1449 }
1450 /// Mark row as not changed
unsetRowChanged(int i)1451 inline void unsetRowChanged(int i)
1452 {
1453 rowChanged_[i] = static_cast< unsigned char >(rowChanged_[i] & (~1));
1454 }
1455 /// Mark row as changed
setRowChanged(int i)1456 inline void setRowChanged(int i)
1457 {
1458 rowChanged_[i] = static_cast< unsigned char >(rowChanged_[i] | (1));
1459 }
1460 /// Mark row as changed and add to list of rows to process next
addRow(int i)1461 inline void addRow(int i)
1462 {
1463 if ((rowChanged_[i] & 1) == 0) {
1464 rowChanged_[i] = static_cast< unsigned char >(rowChanged_[i] | (1));
1465 nextRowsToDo_[numberNextRowsToDo_++] = i;
1466 }
1467 }
1468 /// Test if row is eligible for preprocessing
rowProhibited(int i) const1469 inline bool rowProhibited(int i) const
1470 {
1471 return (rowChanged_[i] & 2) != 0;
1472 }
1473 /*! \brief Test if row is eligible for preprocessing
1474
1475 The difference between this method and #rowProhibited() is that this
1476 method first tests #anyProhibited_ before examining the specific entry
1477 for the specified row.
1478 */
rowProhibited2(int i) const1479 inline bool rowProhibited2(int i) const
1480 {
1481 if (!anyProhibited_)
1482 return false;
1483 else
1484 return (rowChanged_[i] & 2) != 0;
1485 }
1486 /// Mark row as ineligible for preprocessing
setRowProhibited(int i)1487 inline void setRowProhibited(int i)
1488 {
1489 rowChanged_[i] = static_cast< unsigned char >(rowChanged_[i] | (2));
1490 }
1491 /*! \brief Test if row is marked as used
1492
1493 This is for doing faster lookups to see where two rows have entries
1494 in common. It can be used anywhere as long as it ends up zeroed out.
1495 */
rowUsed(int i) const1496 inline bool rowUsed(int i) const
1497 {
1498 return (rowChanged_[i] & 4) != 0;
1499 }
1500 /// Mark row as used
setRowUsed(int i)1501 inline void setRowUsed(int i)
1502 {
1503 rowChanged_[i] = static_cast< unsigned char >(rowChanged_[i] | (4));
1504 }
1505 /// Mark row as unused
unsetRowUsed(int i)1506 inline void unsetRowUsed(int i)
1507 {
1508 rowChanged_[i] = static_cast< unsigned char >(rowChanged_[i] & (~4));
1509 }
1510
1511 /// Check if there are any prohibited rows or columns
anyProhibited() const1512 inline bool anyProhibited() const
1513 {
1514 return anyProhibited_;
1515 }
1516 /// Set a flag for presence of prohibited rows or columns
setAnyProhibited(bool val=true)1517 inline void setAnyProhibited(bool val = true)
1518 {
1519 anyProhibited_ = val;
1520 }
1521 //@}
1522 };
1523
1524 /*! \class CoinPostsolveMatrix
1525 \brief Augments CoinPrePostsolveMatrix with information about the problem
1526 that is only needed during postsolve.
1527
1528 The notable point is that the matrix representation is threaded. The
1529 representation is column-major and starts with the standard two pairs of
1530 arrays: one pair to hold the row indices and coefficients, the second pair
1531 to hold the column starting positions and lengths. But the row indices and
1532 coefficients for a column do not necessarily occupy a contiguous block in
1533 their respective arrays. Instead, a link array gives the position of the
1534 next (row index,coefficient) pair. If the row index and value of a
1535 coefficient a<p,j> occupy position kp in their arrays, then the position of
1536 the next coefficient a<q,j> is found as kq = link[kp].
1537
1538 This threaded representation allows for efficient expansion of columns as
1539 rows are reintroduced during postsolve transformations. The basic packed
1540 structures are allocated to the expected size of the postsolved matrix,
1541 and as new coefficients are added, their location is simply added to the
1542 thread for the column.
1543
1544 There is no provision to convert the threaded representation to a packed
1545 representation. In the context of postsolve, it's not required. (You did
1546 keep a copy of the original matrix, eh?)
1547
1548 The constructors that take an OSI or ClpSimplex as a parameter really should
1549 not be here, but for historical reasons they will likely remain for the
1550 forseeable future. -- lh, 111202 --
1551 */
1552 class CoinPostsolveMatrix : public CoinPrePostsolveMatrix {
1553 public:
1554 /*! \brief `Native' constructor
1555
1556 This constructor creates an empty object which must then be loaded.
1557 On the other hand, it doesn't assume that the client is an
1558 OsiSolverInterface.
1559 */
1560 CoinPostsolveMatrix(int ncols_alloc, int nrows_alloc,
1561 CoinBigIndex nelems_alloc);
1562
1563 /*! \brief Clp OSI constructor
1564
1565 See Clp code for the definition.
1566 */
1567 CoinPostsolveMatrix(ClpSimplex *si,
1568
1569 int ncols0,
1570 int nrows0,
1571 CoinBigIndex nelems0,
1572
1573 double maxmin_,
1574 // end prepost members
1575
1576 double *sol,
1577 double *acts,
1578
1579 unsigned char *colstat,
1580 unsigned char *rowstat);
1581
1582 /*! \brief Generic OSI constructor
1583
1584 See OSI code for the definition.
1585 */
1586 CoinPostsolveMatrix(OsiSolverInterface *si,
1587
1588 int ncols0,
1589 int nrows0,
1590 CoinBigIndex nelems0,
1591
1592 double maxmin_,
1593 // end prepost members
1594
1595 double *sol,
1596 double *acts,
1597
1598 unsigned char *colstat,
1599 unsigned char *rowstat);
1600
1601 /*! \brief Load an empty CoinPostsolveMatrix from a CoinPresolveMatrix
1602
1603 This routine transfers the contents of the CoinPrePostsolveMatrix
1604 object from the CoinPresolveMatrix object to the CoinPostsolveMatrix
1605 object and completes initialisation of the CoinPostsolveMatrix object.
1606 The empty shell of the CoinPresolveMatrix object is destroyed.
1607
1608 The routine expects an empty CoinPostsolveMatrix object. If handed a loaded
1609 object, a lot of memory will leak.
1610 */
1611 void assignPresolveToPostsolve(CoinPresolveMatrix *&preObj);
1612
1613 /// Destructor
1614 ~CoinPostsolveMatrix();
1615
1616 /*! \name Column thread structures
1617
1618 As mentioned in the class documentation, the entries for a given column
1619 do not necessarily occupy a contiguous block of space. The #link_ array
1620 is used to maintain the threading. There is one thread for each column,
1621 and a single thread for all free entries in #hrow_ and #colels_.
1622
1623 The allocated size of #link_ must be at least as large as the allocated
1624 size of #hrow_ and #colels_.
1625 */
1626 //@{
1627
1628 /*! \brief First entry in free entries thread */
1629 CoinBigIndex free_list_;
1630 /// Allocated size of #link_
1631 CoinBigIndex maxlink_;
1632 /*! \brief Thread array
1633
1634 Within a thread, link_[k] points to the next entry in the thread.
1635 */
1636 CoinBigIndex *link_;
1637
1638 //@}
1639
1640 /*! \name Debugging aids
1641
1642 These arrays are allocated only when CoinPresolve is compiled with
1643 PRESOLVE_DEBUG defined. They hold codes which track the reason that
1644 a column or row is added to the problem during postsolve.
1645 */
1646 //@{
1647 char *cdone_;
1648 char *rdone_;
1649 //@}
1650
1651 /// debug
1652 void check_nbasic();
1653 };
1654
1655 /*! \defgroup MtxManip Presolve Matrix Manipulation Functions
1656
1657 Functions to work with the loosely packed and threaded packed matrix
1658 structures used during presolve and postsolve.
1659 */
1660 //@{
1661
1662 /*! \relates CoinPrePostsolveMatrix
1663 \brief Initialise linked list for major vector order in bulk storage
1664 */
1665
1666 void presolve_make_memlists(/*CoinBigIndex *starts,*/ int *lengths,
1667 presolvehlink *link, int n);
1668
1669 /*! \relates CoinPrePostsolveMatrix
1670 \brief Make sure a major-dimension vector k has room for one more
1671 coefficient.
1672
1673 You can use this directly, or use the inline wrappers presolve_expand_col
1674 and presolve_expand_row
1675 */
1676 bool presolve_expand_major(CoinBigIndex *majstrts, double *majels,
1677 int *minndxs, int *majlens,
1678 presolvehlink *majlinks, int nmaj, int k);
1679
1680 /*! \relates CoinPrePostsolveMatrix
1681 \brief Make sure a column (colx) in a column-major matrix has room for
1682 one more coefficient
1683 */
1684
presolve_expand_col(CoinBigIndex * mcstrt,double * colels,int * hrow,int * hincol,presolvehlink * clink,int ncols,int colx)1685 inline bool presolve_expand_col(CoinBigIndex *mcstrt, double *colels,
1686 int *hrow, int *hincol,
1687 presolvehlink *clink, int ncols, int colx)
1688 {
1689 return presolve_expand_major(mcstrt, colels,
1690 hrow, hincol, clink, ncols, colx);
1691 }
1692
1693 /*! \relates CoinPrePostsolveMatrix
1694 \brief Make sure a row (rowx) in a row-major matrix has room for one
1695 more coefficient
1696 */
1697
presolve_expand_row(CoinBigIndex * mrstrt,double * rowels,int * hcol,int * hinrow,presolvehlink * rlink,int nrows,int rowx)1698 inline bool presolve_expand_row(CoinBigIndex *mrstrt, double *rowels,
1699 int *hcol, int *hinrow,
1700 presolvehlink *rlink, int nrows, int rowx)
1701 {
1702 return presolve_expand_major(mrstrt, rowels,
1703 hcol, hinrow, rlink, nrows, rowx);
1704 }
1705
1706 /*! \relates CoinPrePostsolveMatrix
1707 \brief Find position of a minor index in a major vector.
1708
1709 The routine returns the position \c k in \p minndxs for the specified
1710 minor index \p tgt. It will abort if the entry does not exist. Can be
1711 used directly or via the inline wrappers presolve_find_row and
1712 presolve_find_col.
1713 */
presolve_find_minor(int tgt,CoinBigIndex ks,CoinBigIndex ke,const int * minndxs)1714 inline CoinBigIndex presolve_find_minor(int tgt,
1715 CoinBigIndex ks, CoinBigIndex ke,
1716 const int *minndxs)
1717 {
1718 CoinBigIndex k;
1719 for (k = ks; k < ke; k++)
1720 #ifndef NDEBUG
1721 {
1722 if (minndxs[k] == tgt)
1723 return (k);
1724 }
1725 DIE("FIND_MINOR");
1726
1727 abort();
1728 return -1;
1729 #else
1730 {
1731 if (minndxs[k] == tgt)
1732 break;
1733 }
1734 return (k);
1735 #endif
1736 }
1737
1738 /*! \relates CoinPrePostsolveMatrix
1739 \brief Find position of a row in a column in a column-major matrix.
1740
1741 The routine returns the position \c k in \p hrow for the specified \p row.
1742 It will abort if the entry does not exist.
1743 */
presolve_find_row(int row,CoinBigIndex kcs,CoinBigIndex kce,const int * hrow)1744 inline CoinBigIndex presolve_find_row(int row, CoinBigIndex kcs,
1745 CoinBigIndex kce, const int *hrow)
1746 {
1747 return presolve_find_minor(row, kcs, kce, hrow);
1748 }
1749
1750 /*! \relates CoinPostsolveMatrix
1751 \brief Find position of a column in a row in a row-major matrix.
1752
1753 The routine returns the position \c k in \p hcol for the specified \p col.
1754 It will abort if the entry does not exist.
1755 */
presolve_find_col(int col,CoinBigIndex krs,CoinBigIndex kre,const int * hcol)1756 inline CoinBigIndex presolve_find_col(int col, CoinBigIndex krs,
1757 CoinBigIndex kre, const int *hcol)
1758 {
1759 return presolve_find_minor(col, krs, kre, hcol);
1760 }
1761
1762 /*! \relates CoinPrePostsolveMatrix
1763 \brief Find position of a minor index in a major vector.
1764
1765 The routine returns the position \c k in \p minndxs for the specified
1766 minor index \p tgt. A return value of \p ke means the entry does not
1767 exist. Can be used directly or via the inline wrappers
1768 presolve_find_row1 and presolve_find_col1.
1769 */
1770 CoinBigIndex presolve_find_minor1(int tgt, CoinBigIndex ks, CoinBigIndex ke,
1771 const int *minndxs);
1772
1773 /*! \relates CoinPrePostsolveMatrix
1774 \brief Find position of a row in a column in a column-major matrix.
1775
1776 The routine returns the position \c k in \p hrow for the specified \p row.
1777 A return value of \p kce means the entry does not exist.
1778 */
presolve_find_row1(int row,CoinBigIndex kcs,CoinBigIndex kce,const int * hrow)1779 inline CoinBigIndex presolve_find_row1(int row, CoinBigIndex kcs,
1780 CoinBigIndex kce, const int *hrow)
1781 {
1782 return presolve_find_minor1(row, kcs, kce, hrow);
1783 }
1784
1785 /*! \relates CoinPrePostsolveMatrix
1786 \brief Find position of a column in a row in a row-major matrix.
1787
1788 The routine returns the position \c k in \p hcol for the specified \p col.
1789 A return value of \p kre means the entry does not exist.
1790 */
presolve_find_col1(int col,CoinBigIndex krs,CoinBigIndex kre,const int * hcol)1791 inline CoinBigIndex presolve_find_col1(int col, CoinBigIndex krs,
1792 CoinBigIndex kre, const int *hcol)
1793 {
1794 return presolve_find_minor1(col, krs, kre, hcol);
1795 }
1796
1797 /*! \relates CoinPostsolveMatrix
1798 \brief Find position of a minor index in a major vector in a threaded
1799 matrix.
1800
1801 The routine returns the position \c k in \p minndxs for the specified
1802 minor index \p tgt. It will abort if the entry does not exist. Can be
1803 used directly or via the inline wrapper presolve_find_row2.
1804 */
1805 CoinBigIndex presolve_find_minor2(int tgt, CoinBigIndex ks, int majlen,
1806 const int *minndxs,
1807 const CoinBigIndex *majlinks);
1808
1809 /*! \relates CoinPostsolveMatrix
1810 \brief Find position of a row in a column in a column-major threaded
1811 matrix.
1812
1813 The routine returns the position \c k in \p hrow for the specified \p row.
1814 It will abort if the entry does not exist.
1815 */
presolve_find_row2(int row,CoinBigIndex kcs,int collen,const int * hrow,const CoinBigIndex * clinks)1816 inline CoinBigIndex presolve_find_row2(int row, CoinBigIndex kcs, int collen,
1817 const int *hrow,
1818 const CoinBigIndex *clinks)
1819 {
1820 return presolve_find_minor2(row, kcs, collen, hrow, clinks);
1821 }
1822
1823 /*! \relates CoinPostsolveMatrix
1824 \brief Find position of a minor index in a major vector in a threaded
1825 matrix.
1826
1827 The routine returns the position \c k in \p minndxs for the specified
1828 minor index \p tgt. It will return -1 if the entry does not exist.
1829 Can be used directly or via the inline wrappers presolve_find_row3.
1830 */
1831 CoinBigIndex presolve_find_minor3(int tgt, CoinBigIndex ks, int majlen,
1832 const int *minndxs,
1833 const CoinBigIndex *majlinks);
1834
1835 /*! \relates CoinPostsolveMatrix
1836 \brief Find position of a row in a column in a column-major threaded
1837 matrix.
1838
1839 The routine returns the position \c k in \p hrow for the specified \p row.
1840 It will return -1 if the entry does not exist.
1841 */
presolve_find_row3(int row,CoinBigIndex kcs,int collen,const int * hrow,const CoinBigIndex * clinks)1842 inline CoinBigIndex presolve_find_row3(int row, CoinBigIndex kcs, int collen,
1843 const int *hrow,
1844 const CoinBigIndex *clinks)
1845 {
1846 return presolve_find_minor3(row, kcs, collen, hrow, clinks);
1847 }
1848
1849 /*! \relates CoinPrePostsolveMatrix
1850 \brief Delete the entry for a minor index from a major vector.
1851
1852 Deletes the entry for \p minndx from the major vector \p majndx.
1853 Specifically, the relevant entries are removed from the minor index
1854 (\p minndxs) and coefficient (\p els) arrays and the vector length (\p
1855 majlens) is decremented. Loose packing is maintained by swapping the last
1856 entry in the row into the position occupied by the deleted entry.
1857 */
presolve_delete_from_major(int majndx,int minndx,const CoinBigIndex * majstrts,int * majlens,int * minndxs,double * els)1858 inline void presolve_delete_from_major(int majndx, int minndx,
1859 const CoinBigIndex *majstrts,
1860 int *majlens, int *minndxs, double *els)
1861 {
1862 const CoinBigIndex ks = majstrts[majndx];
1863 const CoinBigIndex ke = ks + majlens[majndx];
1864
1865 const CoinBigIndex kmi = presolve_find_minor(minndx, ks, ke, minndxs);
1866
1867 minndxs[kmi] = minndxs[ke - 1];
1868 els[kmi] = els[ke - 1];
1869 majlens[majndx]--;
1870
1871 return;
1872 }
1873
1874 /*! \relates CoinPrePostsolveMatrix
1875 \brief Delete marked entries
1876
1877 Removes the entries specified in \p marked, compressing the major vector
1878 to maintain loose packing. \p marked is cleared in the process.
1879 */
presolve_delete_many_from_major(int majndx,char * marked,const CoinBigIndex * majstrts,int * majlens,int * minndxs,double * els)1880 inline void presolve_delete_many_from_major(int majndx, char *marked,
1881 const CoinBigIndex *majstrts,
1882 int *majlens, int *minndxs, double *els)
1883 {
1884 const CoinBigIndex ks = majstrts[majndx];
1885 const CoinBigIndex ke = ks + majlens[majndx];
1886 CoinBigIndex put = ks;
1887 for (CoinBigIndex k = ks; k < ke; k++) {
1888 int iMinor = minndxs[k];
1889 if (!marked[iMinor]) {
1890 minndxs[put] = iMinor;
1891 els[put++] = els[k];
1892 } else {
1893 marked[iMinor] = 0;
1894 }
1895 }
1896 majlens[majndx] = static_cast< int >(put - ks);
1897 return;
1898 }
1899
1900 /*! \relates CoinPrePostsolveMatrix
1901 \brief Delete the entry for row \p row from column \p col in a
1902 column-major matrix
1903
1904 Deletes the entry for \p row from the major vector for \p col.
1905 Specifically, the relevant entries are removed from the row index (\p
1906 hrow) and coefficient (\p colels) arrays and the vector length (\p
1907 hincol) is decremented. Loose packing is maintained by swapping the last
1908 entry in the row into the position occupied by the deleted entry.
1909 */
presolve_delete_from_col(int row,int col,const CoinBigIndex * mcstrt,int * hincol,int * hrow,double * colels)1910 inline void presolve_delete_from_col(int row, int col,
1911 const CoinBigIndex *mcstrt,
1912 int *hincol, int *hrow, double *colels)
1913 {
1914 presolve_delete_from_major(col, row, mcstrt, hincol, hrow, colels);
1915 }
1916
1917 /*! \relates CoinPrePostsolveMatrix
1918 \brief Delete the entry for column \p col from row \p row in a
1919 row-major matrix
1920
1921 Deletes the entry for \p col from the major vector for \p row.
1922 Specifically, the relevant entries are removed from the column index (\p
1923 hcol) and coefficient (\p rowels) arrays and the vector length (\p
1924 hinrow) is decremented. Loose packing is maintained by swapping the last
1925 entry in the column into the position occupied by the deleted entry.
1926 */
presolve_delete_from_row(int row,int col,const CoinBigIndex * mrstrt,int * hinrow,int * hcol,double * rowels)1927 inline void presolve_delete_from_row(int row, int col,
1928 const CoinBigIndex *mrstrt,
1929 int *hinrow, int *hcol, double *rowels)
1930 {
1931 presolve_delete_from_major(row, col, mrstrt, hinrow, hcol, rowels);
1932 }
1933
1934 /*! \relates CoinPostsolveMatrix
1935 \brief Delete the entry for a minor index from a major vector in a
1936 threaded matrix.
1937
1938 Deletes the entry for \p minndx from the major vector \p majndx.
1939 Specifically, the relevant entries are removed from the minor index (\p
1940 minndxs) and coefficient (\p els) arrays and the vector length (\p
1941 majlens) is decremented. The thread for the major vector is relinked
1942 around the deleted entry and the space is returned to the free list.
1943 */
1944 void presolve_delete_from_major2(int majndx, int minndx,
1945 CoinBigIndex *majstrts, int *majlens,
1946 int *minndxs, CoinBigIndex *majlinks,
1947 CoinBigIndex *free_listp);
1948
1949 /*! \relates CoinPostsolveMatrix
1950 \brief Delete the entry for row \p row from column \p col in a
1951 column-major threaded matrix
1952
1953 Deletes the entry for \p row from the major vector for \p col.
1954 Specifically, the relevant entries are removed from the row index (\p
1955 hrow) and coefficient (\p colels) arrays and the vector length (\p
1956 hincol) is decremented. The thread for the major vector is relinked
1957 around the deleted entry and the space is returned to the free list.
1958 */
presolve_delete_from_col2(int row,int col,CoinBigIndex * mcstrt,int * hincol,int * hrow,CoinBigIndex * clinks,CoinBigIndex * free_listp)1959 inline void presolve_delete_from_col2(int row, int col, CoinBigIndex *mcstrt,
1960 int *hincol, int *hrow,
1961 CoinBigIndex *clinks, CoinBigIndex *free_listp)
1962 {
1963 presolve_delete_from_major2(col, row, mcstrt, hincol, hrow, clinks, free_listp);
1964 }
1965
1966 //@}
1967
1968 /*! \defgroup PresolveUtilities Presolve Utility Functions
1969
1970 Utilities used by multiple presolve transform objects.
1971 */
1972 //@{
1973
1974 /*! \brief Duplicate a major-dimension vector; optionally omit the entry
1975 with minor index \p tgt.
1976
1977 Designed to copy a major-dimension vector from the paired coefficient
1978 (\p elems) and minor index (\p indices) arrays used in the standard
1979 packed matrix representation. Copies \p length entries starting at
1980 \p offset.
1981
1982 If \p tgt is specified, the entry with minor index == \p tgt is
1983 omitted from the copy.
1984 */
1985 double *presolve_dupmajor(const double *elems, const int *indices,
1986 int length, CoinBigIndex offset, int tgt = -1);
1987
1988 /// Initialize a vector with random numbers
1989 void coin_init_random_vec(double *work, int n);
1990
1991 //@}
1992
1993 #endif
1994
1995 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1996 */
1997