source: branches/devel-1/include/PresolveMatrix.hpp @ 39

Last change on this file since 39 was 39, checked in by forrest, 17 years ago

Messing about with presolve

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