source: trunk/include/PresolveMatrix.hpp @ 173

Last change on this file since 173 was 173, checked in by forrest, 16 years ago

Improvements to presolve

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#ifndef PresolveMatrix_H
4#define PresolveMatrix_H
5
6#include "CoinPragma.hpp"
7#include "ClpSimplex.hpp"
8
9#include <cmath>
10#include <cfloat>
11
12#if defined(_MSC_VER)
13// Avoid MS Compiler problem in recognizing type to delete
14// by casting to type.
15#define deleteAction(array,type) delete [] ((type) array)
16#else
17#define deleteAction(array,type) delete [] array
18#endif
19
20
21// OSL had a fixed zero tolerance; we still use that here.
22const double ZTOLDP      = 1e-12;
23
24double *presolve_duparray(const double *d, int n, int n2);
25double *presolve_duparray(const double *d, int n);
26int *presolve_duparray(const int *d, int n, int n2);
27int *presolve_duparray(const int *d, int n);
28char *presolve_duparray(const char *d, int n, int n2);
29char *presolve_duparray(const char *d, int n);
30// This one saves in one go to save [] memory
31double * presolve_duparray(const double * element, const int * index,
32                           int length, int offset);
33
34void presolve_delete_from_row(int row, int col /* thing to delete */,
35                     const CoinBigIndex *mrstrt,
36                     int *hinrow, int *hcol, double *dels);
37
38int presolve_find_row(int row, CoinBigIndex kcs, CoinBigIndex kce, const int *hrow);
39int presolve_find_row1(int row, CoinBigIndex kcs, CoinBigIndex kce, const int *hrow);
40
41//#define       DEBUG_PRESOLVE  1
42
43#ifdef  DEBUG_PRESOLVE
44inline void     DIE(const char *s)      { std::cout<<s; abort(); }
45#else
46  inline void   DIE(const char *s)      {}
47#endif
48
49#if     DEBUG_PRESOLVE
50#define PRESOLVE_STMT(s)        s
51#define PRESOLVEASSERT(x)       ((x) ? 1 : ((std::cerr<< "FAILED ASSERTION at line "<< __LINE__ << ":  " #x "\n"), abort(), 0))
52#else
53#define PRESOLVEASSERT(x) assert(x)
54#define PRESOLVE_STMT(s)
55#endif
56
57struct dropped_zero {
58  int row;
59  int col;
60};
61
62inline int ALIGN(int n, int m)  { return (((n + m - 1) / m) * m); }
63inline int ALIGN_DOUBLE(int n)  { return ALIGN(n,sizeof(double)); }
64
65class PostsolveMatrix;
66
67// Note 77
68// "Members and bases are constructed in ordre of declation
69//  in the class and destroyed in the reverse order."  C++PL 3d Ed. p. 307
70//
71// That's why I put integer members (such as ncols) before the array members;
72// I like to use those integer values during initialization.
73// NOT ANYMORE
74
75//
76// This is the abstract base class of all presolve routines.
77// The object just stores the information needed to postsolve;
78// this information is not expected to be changed, so the fields are
79// all const.
80// It is expected that subclasses will declare static functions that
81// attempt to perform the presolve; if it succeeds, the function creates a
82// new presolve object and returns it, otherwise it returns 0.
83// It is expected that these static functions will be the only things
84// that can create new presolve_action objects;
85// this is expressed by making the constructor(s) private.
86// Every subclass must define a postsolve routine.  This gets two records,
87// one that contains information that is also used in presolving (prob)
88// and one with information that is only used in postsolving (prob2).
89// (suggestions for better names welcome).
90// Finally, there is a single postsolve driver (the friend function)
91// that just calls the postsolve member function of all the postsolve objects
92// in order.
93
94// note that since the only field in a presolve_action is const,
95// anything one can do with a variable declared:
96//      PresolveAction*
97// can also be done with a variable declared:
98//      const PresolveAction*
99//
100// It is expected that all derived subclasses of PresolveAction also
101// have this property.
102class PresolveAction {
103 public:
104  // Exceptions are inefficient, particularly with g++.
105  // Even with xlC, the use of exceptions adds a long prologue to a routine.
106  // Therefore, rather than use throw directly in the routine,
107  // I use it in a stub routine.
108  static void throwCoinError(const char *error, const char *ps_routine);
109
110  // The only thing the base class does is point to the next
111  // presolve transform in the list.
112  const PresolveAction *next;
113 
114  PresolveAction(const PresolveAction *next) : next(next) {}
115
116  // A name for debug printing.
117  // It is expected that the name is not stored in the transform itself.
118  virtual const char *name() const = 0;
119
120  // postsolve this particular presolve action
121  virtual void postsolve(PostsolveMatrix *prob) const = 0;
122
123  virtual ~PresolveAction() {}
124};
125
126// This collects all the information about the problem that is needed
127// both in presolve and postsolve.
128class PrePostsolveMatrix {
129 public:
130  /// enums for status of various sorts (matches CoinWarmStartBasis)
131  enum Status {
132    isFree = 0x00,
133    basic = 0x01,
134    atUpperBound = 0x02,
135    atLowerBound = 0x03,
136    superBasic = 0x04
137  };
138  double *sol_;
139  double *rowduals_;
140  double *acts_;
141
142  double *rcosts_;
143  unsigned char *colstat_;
144  unsigned char *rowstat_;
145
146
147  // Original objective offset
148  double originalOffset_;
149  // Message handler
150   CoinMessageHandler *  handler_; 
151   /// Messages
152   CoinMessages messages_; 
153
154   inline CoinMessageHandler * messageHandler() const 
155  { return handler_; }
156   /// Return messages
157   inline CoinMessages messages() const 
158  { return messages_; }
159  // colrep
160  int ncols_;
161  const int ncols0_;
162
163  CoinBigIndex nelems_;
164
165  CoinBigIndex *mcstrt_;
166  int *hincol_;
167  int *hrow_;
168  double *colels_;
169
170  double *cost_;
171
172  double *clo_;
173  double *cup_;
174  double *rlo_;
175  double *rup_;
176
177  // Original column numbers
178  int * originalColumn_;
179
180  const double ztolzb_;
181  const double ztoldj_;
182
183  double maxmin_;
184
185
186  PrePostsolveMatrix(const ClpSimplex& si,
187                        int ncols_,
188                        int nrows_,
189                        CoinBigIndex nelems_);
190
191  ~PrePostsolveMatrix();
192  // Status stuff
193 
194  inline void setRowStatus(int sequence, Status status)
195  {
196    unsigned char & st_byte = rowstat_[sequence];
197    st_byte &= ~7;
198    st_byte |= status;
199  };
200  inline Status getRowStatus(int sequence) const
201  {return static_cast<Status> (rowstat_[sequence]&7);};
202  inline bool rowIsBasic(int sequence) const
203  {return (static_cast<Status> (rowstat_[sequence]&7)==basic);};
204  inline void setColumnStatus(int sequence, Status status)
205  {
206    unsigned char & st_byte = colstat_[sequence];
207    st_byte &= ~7;
208    st_byte |= status;
209  };
210  inline Status getColumnStatus(int sequence) const
211  {return static_cast<Status> (colstat_[sequence]&7);};
212  inline bool columnIsBasic(int sequence) const
213  {return (static_cast<Status> (colstat_[sequence]&7)==basic);};
214  /// Sets status (non -basic ) using value
215  void setRowStatusUsingValue(int iRow);
216  void setColumnStatusUsingValue(int iColumn);
217
218};
219
220
221
222
223
224/*
225 * Currently, the matrix is represented the same way an CoinPackedMatrix is.
226 * Occasionally columns increase in size.
227 * In order to check whether there is enough room for the column
228 * where it sits, I wanted to know what the next column (in memory order)
229 * in the representation was.
230 * To do that, I maintain a linked list of columns; the "pre" and "suc"
231 * fields give the previous and next columns, in memory order (that is,
232 * the column whose mcstrt entry is next smaller or larger).
233 * The same thing for the row representation.
234 *
235 * This is all likely to change, but I'm leaving it as it is for now.
236 */
237//  static const int    NO_LINK = -66666666;
238#define NO_LINK -66666666
239// Plus infinity
240#ifndef COIN_DBL_MAX
241#define COIN_DBL_MAX DBL_MAX
242#endif
243#define PRESOLVE_INF COIN_DBL_MAX
244
245class presolvehlink {
246public:
247  int pre, suc;
248};
249 
250static inline void PRESOLVE_REMOVE_LINK(presolvehlink *link, int i)
251{ 
252  int ipre = link[i].pre;
253  int isuc = link[i].suc;
254  if (ipre >= 0) {
255    link[ipre].suc = isuc;
256  }
257  if (isuc >= 0) {
258    link[isuc].pre = ipre;
259  }
260  link[i].pre = NO_LINK, link[i].suc = NO_LINK;
261}
262
263// inserts i after pos
264static inline void PRESOLVE_INSERT_LINK(presolvehlink *link, int i, int pos)
265{
266  int isuc = link[pos].suc;
267  link[pos].suc = i;
268  link[i].pre = pos;
269  if (isuc >= 0) {
270    link[isuc].pre = i;
271  }
272  link[i].suc = isuc;
273}
274
275// rename i to j
276// that is, position j should be unused, and i will take its place
277// should be equivalent to:
278//   int pre = link[i].pre;
279//   PRESOLVE_REMOVE_LINK(link, i);
280//   PRESOLVE_INSERT_LINK(link, j, pre);
281// if pre is not NO_LINK (otherwise -- ??)
282static inline void PRESOLVE_MOVE_LINK(presolvehlink *link, int i, int j)
283{ 
284  int ipre = link[i].pre;
285  int isuc = link[i].suc;
286  if (ipre >= 0) {
287    link[ipre].suc = j;
288  }
289  if (isuc >= 0) {
290    link[isuc].pre = j;
291  }
292  link[i].pre = NO_LINK, link[i].suc = NO_LINK;
293}
294
295
296
297// this really should never happen.
298// it will if there isn't enough space to postsolve the matrix.
299// see the note below.
300static inline void check_free_list(int free_list)
301{
302  if (free_list < 0) {
303    printf("RAN OUT OF LINKS!!\n");
304    abort();
305  }
306}
307
308
309  // This collects all the information about the problem that is only
310  // needed during presolve.
311class PresolveMatrix : public PrePostsolveMatrix {
312 public:
313
314  // Crude linked lists, modelled after the linked lists used in OSL factorization.
315  presolvehlink *clink_;
316  presolvehlink *rlink_;
317
318  double dobias_;
319
320  // rowrep
321  int nrows_;   // note 77
322  CoinBigIndex *mrstrt_;
323  int *hinrow_;
324  double *rowels_;
325  int *hcol_;
326
327  char *integerType_;
328  // bounds can be moved by this to stay feasible
329  double feasibilityTolerance_;
330  // Output status 0=feasible, 1 infeasible, 2 unbounded
331  int status_;
332  // Should use templates ?
333  // Rows
334  // Bits to say if row changed
335  unsigned int * rowChanged_;
336  // Input list
337  int * rowsToDo_;
338  int numberRowsToDo_;
339  // Output list
340  int * nextRowsToDo_;
341  int numberNextRowsToDo_;
342
343  inline bool rowChanged(int i) const {
344    return (rowChanged_[i>>5]>>(i&31))!=0;
345  }
346  inline void setRowChanged(int i) {
347    unsigned int & value = rowChanged_[i>>5];
348    int bit = i&31;
349    value |= (1<<bit);
350  }
351  inline void addRow(int i) {
352    unsigned int & value = rowChanged_[i>>5];
353    int bit = i&31;
354    if ((value&(1<<bit))==0) {
355      value |= (1<<bit);
356      nextRowsToDo_[numberNextRowsToDo_++] = i;
357    }
358  }
359  inline void unsetRowChanged(int i) {
360    unsigned int & value = rowChanged_[i>>5];
361    int bit = i&31;
362    value &= ~(1<<bit);
363  }
364
365  // Columns
366  // Bits to say if column changed
367  unsigned int * colChanged_;
368  // Input list
369  int * colsToDo_;
370  int numberColsToDo_;
371  // Output list
372  int * nextColsToDo_;
373  int numberNextColsToDo_;
374
375  inline bool colChanged(int i) const {
376    return (colChanged_[i>>5]>>(i&31))!=0;
377  }
378  inline void setColChanged(int i) {
379    unsigned int & value = colChanged_[i>>5];
380    int bit = i&31;
381    value |= (1<<bit);
382  }
383  inline void addCol(int i) {
384    unsigned int & value = colChanged_[i>>5];
385    int bit = i&31;
386    if ((value&(1<<bit))==0) {
387      value |= (1<<bit);
388      nextColsToDo_[numberNextColsToDo_++] = i;
389    }
390  }
391  inline void unsetColChanged(int i) {
392    unsigned int & value = colChanged_[i>>5];
393    int bit = i&31;
394    value &= ~(1<<bit);
395  }
396  // Bits to say if row can not be touched
397  unsigned int * rowProhibited_;
398  inline bool anyProhibited() const
399  { return rowProhibited_!=NULL;};
400
401  inline bool rowProhibited(int i) const {
402    return (rowProhibited_[i>>5]>>(i&31))!=0;
403  }
404  // This one for lazy testing
405  inline bool rowProhibited2(int i) const {
406    if (!rowProhibited_)
407      return false;
408    else
409      return (rowProhibited_[i>>5]>>(i&31))!=0;
410  }
411  inline void setRowProhibited(int i) {
412    unsigned int & value = rowProhibited_[i>>5];
413    int bit = i&31;
414    value |= (1<<bit);
415  }
416  // Columns
417  // Bits to say if column can not be touched
418  unsigned int * colProhibited_;
419  inline bool colProhibited(int i) const {
420    return (colProhibited_[i>>5]>>(i&31))!=0;
421  }
422  // This one for lazy testing
423  inline bool colProhibited2(int i) const {
424    if (!colProhibited_)
425      return false;
426    else
427      return (colProhibited_[i>>5]>>(i&31))!=0;
428  }
429  inline void setColProhibited(int i) {
430    unsigned int & value = colProhibited_[i>>5];
431    int bit = i&31;
432    value |= (1<<bit);
433  }
434  PresolveMatrix(int ncols0,
435                    double maxmin,
436                    // end prepost members
437
438                    ClpSimplex &si,
439
440                    // rowrep
441                    int nrows,
442                    CoinBigIndex nelems,
443                 bool doStatus,
444                 double nonLinearVariable);
445
446  ~PresolveMatrix();
447
448  void consistent(bool testvals = true);
449
450  inline void change_bias(double change_amount);
451
452  void update_model(ClpSimplex& si,
453                            int nrows0,
454                            int ncols0,
455                            CoinBigIndex nelems0);
456};
457
458
459  // This collects all the information about the problem that is needed
460  // only in postsolve.
461class PostsolveMatrix : public PrePostsolveMatrix {
462 public:
463
464
465  CoinBigIndex free_list_;
466  int maxlink_;
467  int *link_;
468
469  // debug
470  char *cdone_;
471  char *rdone_;
472  int nrows_;
473
474  // needed for presolve_empty
475  int nrows0_;
476
477  PostsolveMatrix(ClpSimplex& si,
478
479                   int ncols0,
480                   int nrows0,
481                   CoinBigIndex nelems0,
482                     
483                   double maxmin_,
484                   // end prepost members
485
486                   double *sol,
487                   double *acts,
488
489                   unsigned char *colstat,
490                   unsigned char *rowstat);
491
492  ~PostsolveMatrix();
493
494  // debugging
495  void check_nbasic();
496
497  ////virtual void postsolve(const PresolveAction *paction);
498};
499
500void PresolveMatrix::change_bias(double change_amount)
501{
502  dobias_ += change_amount;
503#if DEBUG_PRESOLVE
504  assert(fabs(change_amount)<1.0e50);
505#endif
506  if (change_amount)
507    PRESOLVE_STMT(printf("changing bias by %g to %g\n",
508                    change_amount, dobias_)); 
509}
510
511// useful functions
512inline void swap(int &x, int &y) { int temp = x; x = y; y = temp; }
513inline void swap(double &x, double &y) { double temp = x; x = y; y = temp; }
514inline void swap(long &x, long &y) { long temp = x; x = y; y = temp; }
515
516inline void swap(int *&x, int *&y) { int *temp = x; x = y; y = temp; }
517inline void swap(double *&x, double *&y) { double *temp = x; x = y; y = temp; }
518// This returns a non const array filled with input from scalar
519// or actual array
520template <class T> inline T*
521copyOfArray( const T * array, const int size, T value)
522{
523  T * arrayNew = new T[size];
524  if (array) {
525    memcpy(arrayNew,array,size*sizeof(T));
526  } else {
527    int i;
528    for (i=0;i<size;i++) 
529      arrayNew[i] = value;
530  }
531  return arrayNew;
532}
533
534// This returns a non const array filled with actual array (or NULL)
535template <class T> inline T*
536copyOfArray( const T * array, const int size)
537{
538  if (array) {
539    T * arrayNew = new T[size];
540    memcpy(arrayNew,array,size*sizeof(T));
541    return arrayNew;
542  } else {
543    return NULL;
544  }
545}
546#define PRESOLVEFINITE(n)       (-PRESOLVE_INF < (n) && (n) < PRESOLVE_INF)
547int presolve_find_row2(int irow, CoinBigIndex ks, int nc, const int *hrow, const int *link);
548void presolve_make_memlists(CoinBigIndex *starts, int *lengths,
549                       presolvehlink *link,
550                       int n);
551int presolve_find_row3(int irow, CoinBigIndex ks, int nc, const int *hrow, const int *link);
552void presolve_delete_from_row(int row, int col /* thing to delete */,
553                     const CoinBigIndex *mrstrt,
554                              int *hinrow, int *hcol, double *dels);
555void presolve_delete_from_row2(int row, int col /* thing to delete */,
556                      CoinBigIndex *mrstrt,
557                               int *hinrow, int *hcol, double *dels, int *link, CoinBigIndex *free_listp);
558#endif
Note: See TracBrowser for help on using the repository browser.