source: trunk/Clp/src/ClpMatrixBase.hpp @ 1367

Last change on this file since 1367 was 1367, checked in by forrest, 11 years ago

try and improve stability - also long double interior point

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.4 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#ifndef ClpMatrixBase_H
4#define ClpMatrixBase_H
5
6#include "CoinPragma.hpp"
7#include "CoinFinite.hpp"
8
9#include "CoinPackedMatrix.hpp"
10class CoinIndexedVector;
11class ClpSimplex;
12class ClpModel;
13
14/** Abstract base class for Clp Matrices
15   
16Since this class is abstract, no object of this type can be created.
17
18If a derived class provides all methods then all Clp algorithms
19should work.  Some can be very inefficient e.g. getElements etc is
20only used for tightening bounds for dual and the copies are
21deleted.  Many methods can just be dummy i.e. abort(); if not
22all features are being used.  So if column generation was being done
23then it makes no sense to do steepest edge so there would be
24no point providing subsetTransposeTimes.
25*/
26
27class ClpMatrixBase  {
28 
29public:
30  /**@name Virtual methods that the derived classes must provide */
31  //@{
32  /// Return a complete CoinPackedMatrix
33  virtual CoinPackedMatrix * getPackedMatrix() const = 0;
34  /** Whether the packed matrix is column major ordered or not. */
35  virtual bool isColOrdered() const = 0;
36  /** Number of entries in the packed matrix. */
37  virtual CoinBigIndex getNumElements() const = 0;
38  /** Number of columns. */
39  virtual int getNumCols() const = 0;
40  /** Number of rows. */
41  virtual int getNumRows() const = 0;
42 
43  /** A vector containing the elements in the packed matrix. Note that there
44      might be gaps in this list, entries that do not belong to any
45      major-dimension vector. To get the actual elements one should look at
46      this vector together with vectorStarts and vectorLengths. */
47  virtual const double * getElements() const = 0;
48  /** A vector containing the minor indices of the elements in the packed
49      matrix. Note that there might be gaps in this list, entries that do not
50      belong to any major-dimension vector. To get the actual elements one
51      should look at this vector together with vectorStarts and
52      vectorLengths. */
53  virtual const int * getIndices() const = 0;
54 
55  virtual const CoinBigIndex * getVectorStarts() const = 0;
56  /** The lengths of the major-dimension vectors. */
57  virtual const int * getVectorLengths() const = 0 ;
58  /** The length of a single major-dimension vector. */
59  virtual int getVectorLength(int index) const ;
60  /** Delete the columns whose indices are listed in <code>indDel</code>. */
61  virtual void deleteCols(const int numDel, const int * indDel) = 0;
62  /** Delete the rows whose indices are listed in <code>indDel</code>. */
63  virtual void deleteRows(const int numDel, const int * indDel) = 0;
64#ifndef CLP_NO_VECTOR
65  /// Append Columns
66  virtual void appendCols(int number, const CoinPackedVectorBase * const * columns);
67  /// Append Rows
68  virtual void appendRows(int number, const CoinPackedVectorBase * const * rows);
69#endif
70  /** Modify one element of packed matrix.  An element may be added.
71      This works for either ordering If the new element is zero it will be
72      deleted unless keepZero true */
73  virtual void modifyCoefficient(int row, int column, double newElement,
74                                 bool keepZero=false);
75  /** Append a set of rows/columns to the end of the matrix. Returns number of errors
76      i.e. if any of the new rows/columns contain an index that's larger than the
77      number of columns-1/rows-1 (if numberOther>0) or duplicates
78      If 0 then rows, 1 if columns */
79  virtual int appendMatrix(int number, int type,
80                           const CoinBigIndex * starts, const int * index,
81                           const double * element, int numberOther=-1);
82 
83  /** Returns a new matrix in reverse order without gaps
84      Is allowed to return NULL if doesn't want to have row copy */
85  virtual ClpMatrixBase * reverseOrderedCopy() const {return NULL;}
86 
87  /// Returns number of elements in column part of basis
88  virtual CoinBigIndex countBasis(ClpSimplex * model,
89                                 const int * whichColumn, 
90                                 int numberRowBasic,
91                                  int & numberColumnBasic)=0;
92  /// Fills in column part of basis
93  virtual void fillBasis(ClpSimplex * model,
94                                 const int * whichColumn, 
95                                 int & numberColumnBasic,
96                                 int * row, int * start,
97                                 int * rowCount, int * columnCount,
98                                 CoinFactorizationDouble * element)=0;
99  /** Creates scales for column copy (rowCopy in model may be modified)
100      default does not allow scaling
101      returns non-zero if no scaling done */
102  virtual int scale(ClpModel * model, const ClpSimplex * baseModel=NULL) const 
103  { return 1;}
104  /** Scales rowCopy if column copy scaled
105      Only called if scales already exist */
106  virtual void scaleRowCopy(ClpModel * model) const 
107  { }
108  /// Returns true if can create row copy
109  virtual bool canGetRowCopy() const
110  { return true;}
111  /** Realy really scales column copy
112      Only called if scales already exist.
113      Up to user to delete */
114  inline virtual ClpMatrixBase * scaledColumnCopy(ClpModel * model) const 
115  { return this->clone();}
116 
117  /** Checks if all elements are in valid range.  Can just
118      return true if you are not paranoid.  For Clp I will
119      probably expect no zeros.  Code can modify matrix to get rid of
120      small elements.
121      check bits (can be turned off to save time) :
122      1 - check if matrix has gaps
123      2 - check if zero elements
124      4 - check and compress duplicates
125      8 - report on large and small
126  */
127  virtual bool allElementsInRange(ClpModel * model,
128                                  double smallest, double largest,
129                                  int check=15)
130  { return true;}
131  /** Set the dimensions of the matrix. In effect, append new empty
132      columns/rows to the matrix. A negative number for either dimension
133      means that that dimension doesn't change. Otherwise the new dimensions
134      MUST be at least as large as the current ones otherwise an exception
135      is thrown. */
136  virtual void setDimensions(int numrows, int numcols);
137  /** Returns largest and smallest elements of both signs.
138      Largest refers to largest absolute value.
139      If returns zeros then can't tell anything */
140  virtual void rangeOfElements(double & smallestNegative, double & largestNegative,
141                               double & smallestPositive, double & largestPositive);
142 
143  /** Unpacks a column into an CoinIndexedvector
144   */
145  virtual void unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
146                      int column) const =0;
147  /** Unpacks a column into an CoinIndexedvector
148   ** in packed format
149   Note that model is NOT const.  Bounds and objective could
150   be modified if doing column generation (just for this variable) */
151  virtual void unpackPacked(ClpSimplex * model,
152                            CoinIndexedVector * rowArray,
153                            int column) const =0;
154  /** Purely for column generation and similar ideas.  Allows
155      matrix and any bounds or costs to be updated (sensibly).
156      Returns non-zero if any changes.
157  */
158  virtual int refresh(ClpSimplex * model)
159  { return 0;}
160 
161  // Really scale matrix
162  virtual void reallyScale(const double * rowScale, const double * columnScale);
163  /** Given positive integer weights for each row fills in sum of weights
164      for each column (and slack).
165      Returns weights vector
166      Default returns vector of ones
167  */
168  virtual CoinBigIndex * dubiousWeights(const ClpSimplex * model,int * inputWeights) const;
169  /** Adds multiple of a column into an CoinIndexedvector
170      You can use quickAdd to add to vector */
171  virtual void add(const ClpSimplex * model,CoinIndexedVector * rowArray,
172                   int column, double multiplier) const =0;
173  /** Adds multiple of a column into an array */
174  virtual void add(const ClpSimplex * model,double * array,
175                   int column, double multiplier) const =0;
176  /// Allow any parts of a created CoinPackedMatrix to be deleted
177  virtual void releasePackedMatrix() const =0;
178  /// Says whether it can do partial pricing
179  virtual bool canDoPartialPricing() const;
180  /// Returns number of hidden rows e.g. gub
181  virtual int hiddenRows() const;
182  /// Partial pricing
183  virtual void partialPricing(ClpSimplex * model, double start, double end,
184                              int & bestSequence, int & numberWanted);
185  /** expands an updated column to allow for extra rows which the main
186      solver does not know about and returns number added.
187     
188      This will normally be a no-op - it is in for GUB but may get extended to
189      general non-overlapping and embedded networks.
190     
191      mode 0 - extend
192      mode 1 - delete etc
193  */
194  virtual int extendUpdated(ClpSimplex * model,CoinIndexedVector * update,int mode);
195  /**
196     utility primal function for dealing with dynamic constraints
197     mode=0  - Set up before "update" and "times" for primal solution using extended rows
198     mode=1  - Cleanup primal solution after "times" using extended rows.
199     mode=2  - Check (or report on) primal infeasibilities
200  */
201  virtual void primalExpanded(ClpSimplex * model,int mode);
202  /**
203      utility dual function for dealing with dynamic constraints
204      mode=0  - Set up before "updateTranspose" and "transposeTimes" for duals using extended
205                updates array (and may use other if dual values pass)
206      mode=1  - Update dual solution after "transposeTimes" using extended rows.
207      mode=2  - Compute all djs and compute key dual infeasibilities
208      mode=3  - Report on key dual infeasibilities
209      mode=4  - Modify before updateTranspose in partial pricing
210  */
211  virtual void dualExpanded(ClpSimplex * model,CoinIndexedVector * array,
212                            double * other,int mode);
213  /**
214      general utility function for dealing with dynamic constraints
215      mode=0  - Create list of non-key basics in pivotVariable_ using
216                number as numberBasic in and out
217      mode=1  - Set all key variables as basic
218      mode=2  - return number extra rows needed, number gives maximum number basic
219      mode=3  - before replaceColumn
220      mode=4  - return 1 if can do primal, 2 if dual, 3 if both
221      mode=5  - save any status stuff (when in good state)
222      mode=6  - restore status stuff
223      mode=7  - flag given variable (normally sequenceIn)
224      mode=8  - unflag all variables
225      mode=9  - synchronize costs and bounds
226      mode=10  - return 1 if there may be changing bounds on variable (column generation)
227      mode=11  - make sure set is clean (used when a variable rejected - but not flagged)
228      mode=12  - after factorize but before permute stuff
229      mode=13  - at end of simplex to delete stuff
230
231  */
232  virtual int generalExpanded(ClpSimplex * model,int mode,int & number);
233  /**
234     update information for a pivot (and effective rhs)
235  */
236  virtual int updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue);
237  /** Creates a variable.  This is called after partial pricing and may modify matrix.
238      May update bestSequence.
239  */
240  virtual void createVariable(ClpSimplex * model, int & bestSequence);
241  /** Just for debug if odd type matrix.
242      Returns number of primal infeasibilities. */
243  virtual int checkFeasible(ClpSimplex * model,double & sum) const ;
244  /// Returns reduced cost of a variable
245  double reducedCost(ClpSimplex * model,int sequence) const;
246  /// Correct sequence in and out to give true value (if both -1 maybe do whole matrix)
247  virtual void correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) ;
248  //@}
249 
250  //---------------------------------------------------------------------------
251  /**@name Matrix times vector methods
252     They can be faster if scalar is +- 1
253     Also for simplex I am not using basic/non-basic split */
254  //@{
255  /** Return <code>y + A * x * scalar</code> in <code>y</code>.
256      @pre <code>x</code> must be of size <code>numColumns()</code>
257      @pre <code>y</code> must be of size <code>numRows()</code> */
258  virtual void times(double scalar,
259                     const double * x, double * y) const=0;
260  /** And for scaling - default aborts for when scaling not supported
261      (unless pointers NULL when as normal)
262  */
263  virtual void times(double scalar,
264                     const double * x, double * y,
265                     const double * rowScale, 
266                     const double * columnScale) const;
267  /** Return <code>y + x * scalar * A</code> in <code>y</code>.
268      @pre <code>x</code> must be of size <code>numRows()</code>
269      @pre <code>y</code> must be of size <code>numColumns()</code> */
270  virtual void transposeTimes(double scalar,
271                              const double * x, double * y) const = 0;
272  /** And for scaling - default aborts for when scaling not supported
273      (unless pointers NULL when as normal)
274  */
275  virtual void transposeTimes(double scalar,
276                              const double * x, double * y,
277                              const double * rowScale, 
278                              const double * columnScale,
279                              double * spare=NULL) const;
280#if COIN_LONG_WORK
281  // For long double versions (aborts if not supported)
282  virtual void times(CoinWorkDouble scalar,
283                     const CoinWorkDouble * x, CoinWorkDouble * y) const ;
284  virtual void transposeTimes(CoinWorkDouble scalar,
285                              const CoinWorkDouble * x, CoinWorkDouble * y) const ;
286#endif
287  /** Return <code>x * scalar *A + y</code> in <code>z</code>.
288      Can use y as temporary array (will be empty at end)
289      Note - If x packed mode - then z packed mode
290      Squashes small elements and knows about ClpSimplex */
291  virtual void transposeTimes(const ClpSimplex * model, double scalar,
292                              const CoinIndexedVector * x,
293                              CoinIndexedVector * y,
294                              CoinIndexedVector * z) const = 0;
295  /** Return <code>x *A</code> in <code>z</code> but
296      just for indices in y.
297      This is only needed for primal steepest edge.
298      Note - z always packed mode */
299  virtual void subsetTransposeTimes(const ClpSimplex * model,
300                                    const CoinIndexedVector * x,
301                                    const CoinIndexedVector * y,
302                                    CoinIndexedVector * z) const = 0;
303  /** Returns true if can combine transposeTimes and subsetTransposeTimes
304      and if it would be faster */
305  virtual bool canCombine(const ClpSimplex * model,
306                          const CoinIndexedVector * pi) const {return false;}
307  /// Updates two arrays for steepest and does devex weights (need not be coded)
308  virtual void transposeTimes2(const ClpSimplex * model,
309                               const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
310                               const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
311                               CoinIndexedVector * spare,
312                               double referenceIn, double devex,
313                               // Array for exact devex to say what is in reference framework
314                               unsigned int * reference,
315                               double * weights, double scaleFactor);
316  /// Updates second array for steepest and does devex weights (need not be coded)
317  virtual void subsetTimes2(const ClpSimplex * model,
318                            CoinIndexedVector * dj1,
319                               const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
320                               double referenceIn, double devex,
321                               // Array for exact devex to say what is in reference framework
322                               unsigned int * reference,
323                               double * weights, double scaleFactor);
324  /** Return <code>x *A</code> in <code>z</code> but
325      just for number indices in y.
326      Default cheats with fake CoinIndexedVector and
327      then calls subsetTransposeTimes */
328  virtual void listTransposeTimes(const ClpSimplex * model,
329                                  double * x,
330                                  int * y,
331                                  int number,
332                                  double * z) const;
333  //@}
334  //@{
335  ///@name Other
336  /// Clone
337  virtual ClpMatrixBase * clone() const = 0;
338  /** Subset clone (without gaps).  Duplicates are allowed
339      and order is as given.
340      Derived classes need not provide this as it may not always make
341      sense */
342  virtual ClpMatrixBase * subsetClone (
343                                       int numberRows, const int * whichRows,
344                                       int numberColumns, const int * whichColumns) const;
345  /// Gets rid of any mutable by products
346  virtual void backToBasics() {}
347  /** Returns type.
348      The types which code may need to know about are:
349      1  - ClpPackedMatrix
350      11 - ClpNetworkMatrix
351      12 - ClpPlusMinusOneMatrix
352  */
353  inline int type() const
354  { return type_;}
355  /// Sets type
356  void setType(int type) {type_=type;}
357  /// Sets up an effective RHS
358  void useEffectiveRhs(ClpSimplex * model);
359  /** Returns effective RHS offset if it is being used.  This is used for long problems
360      or big gub or anywhere where going through full columns is
361      expensive.  This may re-compute */
362  virtual double * rhsOffset(ClpSimplex * model,bool forceRefresh=false,
363                                bool check=false);
364  /// If rhsOffset used this is iteration last refreshed
365  inline int lastRefresh() const
366  { return lastRefresh_;}
367  /// If rhsOffset used this is refresh frequency (0==off)
368  inline int refreshFrequency() const
369  { return refreshFrequency_;}
370  inline void setRefreshFrequency(int value)
371  { refreshFrequency_=value;}
372  /// whether to skip dual checks most of time
373  inline bool skipDualCheck() const
374  { return skipDualCheck_;}
375  inline void setSkipDualCheck(bool yes)
376  { skipDualCheck_=yes;}
377  /** Partial pricing tuning parameter - minimum number of "objects" to scan.
378      e.g. number of Gub sets but could be number of variables */
379  inline int minimumObjectsScan() const
380  { return minimumObjectsScan_;}
381  inline void setMinimumObjectsScan(int value)
382  { minimumObjectsScan_=value;}
383  /// Partial pricing tuning parameter - minimum number of negative reduced costs to get
384  inline int minimumGoodReducedCosts() const
385  { return minimumGoodReducedCosts_;}
386  inline void setMinimumGoodReducedCosts(int value)
387  { minimumGoodReducedCosts_=value;}
388  /// Current start of search space in matrix (as fraction)
389  inline double startFraction() const
390  { return startFraction_;}
391  inline void setStartFraction(double value) 
392  { startFraction_ = value;}
393  /// Current end of search space in matrix (as fraction)
394  inline double endFraction() const
395  { return endFraction_;}
396  inline void setEndFraction(double value) 
397  { endFraction_ = value;}
398  /// Current best reduced cost
399  inline double savedBestDj() const
400  { return savedBestDj_;}
401  inline void setSavedBestDj(double value) 
402  { savedBestDj_ = value;}
403  /// Initial number of negative reduced costs wanted
404  inline int originalWanted() const
405  { return originalWanted_;}
406  inline void setOriginalWanted(int value) 
407  { originalWanted_ = value;}
408  /// Current number of negative reduced costs which we still need
409  inline int currentWanted() const
410  { return currentWanted_;}
411  inline void setCurrentWanted(int value) 
412  { currentWanted_ = value;}
413  /// Current best sequence
414  inline int savedBestSequence() const
415  { return savedBestSequence_;}
416  inline void setSavedBestSequence(int value) 
417  { savedBestSequence_ = value;}
418  //@}
419 
420 
421protected:
422 
423  /**@name Constructors, destructor<br>
424     <strong>NOTE</strong>: All constructors are protected. There's no need
425     to expose them, after all, this is an abstract class. */
426  //@{
427  /** Default constructor. */
428  ClpMatrixBase();
429  /** Destructor (has to be public) */
430public:
431  virtual ~ClpMatrixBase();
432protected:
433  // Copy
434  ClpMatrixBase(const ClpMatrixBase&);
435  // Assignment
436  ClpMatrixBase& operator=(const ClpMatrixBase&);
437  //@}
438 
439 
440protected:
441  /**@name Data members
442     The data members are protected to allow access for derived classes. */
443  //@{
444  /** Effective RHS offset if it is being used.  This is used for long problems
445      or big gub or anywhere where going through full columns is
446      expensive */
447  double * rhsOffset_;
448  /// Current start of search space in matrix (as fraction)
449  double startFraction_;
450  /// Current end of search space in matrix (as fraction)
451  double endFraction_;
452  /// Best reduced cost so far
453  double savedBestDj_;
454  /// Initial number of negative reduced costs wanted
455  int originalWanted_;
456  /// Current number of negative reduced costs which we still need
457  int currentWanted_;
458  /// Saved best sequence in pricing
459  int savedBestSequence_;
460  /// type (may be useful)
461  int type_;
462  /// If rhsOffset used this is iteration last refreshed
463  int lastRefresh_;
464  /// If rhsOffset used this is refresh frequency (0==off)
465  int refreshFrequency_;
466  /// Partial pricing tuning parameter - minimum number of "objects" to scan
467  int minimumObjectsScan_;
468  /// Partial pricing tuning parameter - minimum number of negative reduced costs to get
469  int minimumGoodReducedCosts_;
470  /// True sequence in (i.e. from larger problem)
471  int trueSequenceIn_;
472  /// True sequence out (i.e. from larger problem)
473  int trueSequenceOut_;
474  /// whether to skip dual checks most of time
475  bool skipDualCheck_;
476  //@}
477};
478// bias for free variables
479#define FREE_BIAS 1.0e1
480// Acceptance criteria for free variables
481#define FREE_ACCEPT 1.0e2
482
483#endif
Note: See TracBrowser for help on using the repository browser.