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