source: trunk/Clp/src/ClpPackedMatrix.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: 21.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#ifndef ClpPackedMatrix_H
4#define ClpPackedMatrix_H
5
6#include "CoinPragma.hpp"
7
8#include "ClpMatrixBase.hpp"
9
10/** This implements CoinPackedMatrix as derived from ClpMatrixBase.
11
12    It adds a few methods that know about model as well as matrix
13
14    For details see CoinPackedMatrix */
15
16class ClpPackedMatrix2;
17class ClpPackedMatrix3;
18class ClpPackedMatrix : public ClpMatrixBase {
19 
20public:
21  /**@name Useful methods */
22   //@{
23   /// Return a complete CoinPackedMatrix
24  virtual CoinPackedMatrix * getPackedMatrix() const { return matrix_;}
25  /** Whether the packed matrix is column major ordered or not. */
26  virtual bool isColOrdered() const { return matrix_->isColOrdered(); }
27   /** Number of entries in the packed matrix. */
28  virtual  CoinBigIndex getNumElements() const 
29  { return matrix_->getNumElements(); }
30  /** Number of columns. */
31  virtual int getNumCols() const { return matrix_->getNumCols(); }
32  /** Number of rows. */
33  virtual int getNumRows() const { return matrix_->getNumRows(); }
34 
35  /** A vector containing the elements in the packed matrix. Note that there
36      might be gaps in this list, entries that do not belong to any
37      major-dimension vector. To get the actual elements one should look at
38      this vector together with vectorStarts and vectorLengths. */
39  virtual const double * getElements() const 
40  { return matrix_->getElements();}
41  /// Mutable elements
42  inline double * getMutableElements() const
43  { return matrix_->getMutableElements();}
44   /** A vector containing the minor indices of the elements in the packed
45        matrix. Note that there might be gaps in this list, entries that do not
46        belong to any major-dimension vector. To get the actual elements one
47        should look at this vector together with vectorStarts and
48        vectorLengths. */
49   virtual const int * getIndices() const 
50  { return matrix_->getIndices();}
51
52   virtual const CoinBigIndex * getVectorStarts() const 
53  { return matrix_->getVectorStarts();}
54   /** The lengths of the major-dimension vectors. */
55   virtual const int * getVectorLengths() const 
56  { return matrix_->getVectorLengths();} 
57  /** The length of a single major-dimension vector. */
58  virtual int getVectorLength(int index) const 
59  { return matrix_->getVectorSize(index);}
60
61    /** Delete the columns whose indices are listed in <code>indDel</code>. */
62  virtual void deleteCols(const int numDel, const int * indDel);
63    /** Delete the rows whose indices are listed in <code>indDel</code>. */
64  virtual void deleteRows(const int numDel, const int * indDel);
65#ifndef CLP_NO_VECTOR
66  /// Append Columns
67  virtual void appendCols(int number, const CoinPackedVectorBase * const * columns);
68  /// Append Rows
69  virtual void appendRows(int number, const CoinPackedVectorBase * const * rows);
70#endif
71  /** Append a set of rows/columns to the end of the matrix. Returns number of errors
72      i.e. if any of the new rows/columns contain an index that's larger than the
73      number of columns-1/rows-1 (if numberOther>0) or duplicates
74      If 0 then rows, 1 if columns */
75  virtual int appendMatrix(int number, int type,
76                           const CoinBigIndex * starts, const int * index,
77                           const double * element, int numberOther=-1);
78  /** Replace the elements of a vector.  The indices remain the same.
79      This is only needed if scaling and a row copy is used.
80      At most the number specified will be replaced.
81      The index is between 0 and major dimension of matrix */
82  virtual void replaceVector(const int index,
83                       const int numReplace, const double * newElements)
84      {matrix_->replaceVector(index,numReplace,newElements);}
85  /** Modify one element of packed matrix.  An element may be added.
86      This works for either ordering If the new element is zero it will be
87      deleted unless keepZero true */
88  virtual void modifyCoefficient(int row, int column, double newElement,
89                           bool keepZero=false)
90       {matrix_->modifyCoefficient(row, column, newElement, keepZero);}
91  /** Returns a new matrix in reverse order without gaps */
92  virtual ClpMatrixBase * reverseOrderedCopy() const;
93  /// Returns number of elements in column part of basis
94  virtual CoinBigIndex countBasis(ClpSimplex * model,
95                                 const int * whichColumn, 
96                                 int numberRowBasic,
97                                  int & numberColumnBasic);
98  /// Fills in column part of basis
99  virtual void fillBasis(ClpSimplex * model,
100                                 const int * whichColumn, 
101                                 int & numberColumnBasic,
102                                 int * row, int * start,
103                                 int * rowCount, int * columnCount,
104                                 CoinFactorizationDouble * element);
105  /** Creates scales for column copy (rowCopy in model may be modified)
106      returns non-zero if no scaling done */
107  virtual int scale(ClpModel * model,const ClpSimplex * baseModel=NULL) const ;
108  /** Scales rowCopy if column copy scaled
109      Only called if scales already exist */
110  virtual void scaleRowCopy(ClpModel * model) const ;
111  /// Creates scaled column copy if scales exist
112  void createScaledMatrix(ClpSimplex * model) const;
113  /** Realy really scales column copy
114      Only called if scales already exist.
115      Up to user ro delete */
116  virtual ClpMatrixBase * scaledColumnCopy(ClpModel * model) const ;
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  /** Returns largest and smallest elements of both signs.
131      Largest refers to largest absolute value.
132  */
133  virtual void rangeOfElements(double & smallestNegative, double & largestNegative,
134                       double & smallestPositive, double & largestPositive);
135
136  /** Unpacks a column into an CoinIndexedvector
137   */
138  virtual void unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
139                   int column) const ;
140  /** Unpacks a column into an CoinIndexedvector
141   ** in packed foramt
142      Note that model is NOT const.  Bounds and objective could
143      be modified if doing column generation (just for this variable) */
144  virtual void unpackPacked(ClpSimplex * model,
145                            CoinIndexedVector * rowArray,
146                            int column) const;
147  /** Adds multiple of a column into an CoinIndexedvector
148      You can use quickAdd to add to vector */
149  virtual void add(const ClpSimplex * model,CoinIndexedVector * rowArray,
150                   int column, double multiplier) const ;
151  /** Adds multiple of a column into an array */
152  virtual void add(const ClpSimplex * model,double * array,
153                   int column, double multiplier) const;
154   /// Allow any parts of a created CoinPackedMatrix to be deleted
155   virtual void releasePackedMatrix() const { }
156  /** Given positive integer weights for each row fills in sum of weights
157      for each column (and slack).
158      Returns weights vector
159  */
160  virtual CoinBigIndex * dubiousWeights(const ClpSimplex * model,int * inputWeights) const;
161  /// Says whether it can do partial pricing
162  virtual bool canDoPartialPricing() const;
163  /// Partial pricing
164  virtual void partialPricing(ClpSimplex * model, double start, double end,
165                      int & bestSequence, int & numberWanted);
166  /// makes sure active columns correct
167  virtual int refresh(ClpSimplex * model);
168  // Really scale matrix
169  virtual void reallyScale(const double * rowScale, const double * columnScale);
170  /** Set the dimensions of the matrix. In effect, append new empty
171      columns/rows to the matrix. A negative number for either dimension
172      means that that dimension doesn't change. Otherwise the new dimensions
173      MUST be at least as large as the current ones otherwise an exception
174      is thrown. */
175  virtual void setDimensions(int numrows, int numcols);
176   //@}
177
178  /**@name Matrix times vector methods */
179  //@{
180    /** Return <code>y + A * scalar *x</code> in <code>y</code>.
181        @pre <code>x</code> must be of size <code>numColumns()</code>
182        @pre <code>y</code> must be of size <code>numRows()</code> */
183  virtual void times(double scalar,
184                       const double * x, double * y) const;
185  /// And for scaling
186  virtual void times(double scalar,
187                     const double * x, double * y,
188                     const double * rowScale, 
189                     const double * columnScale) const;
190    /** Return <code>y + x * scalar * A</code> in <code>y</code>.
191        @pre <code>x</code> must be of size <code>numRows()</code>
192        @pre <code>y</code> must be of size <code>numColumns()</code> */
193    virtual void transposeTimes(double scalar,
194                                const double * x, double * y) const;
195  /// And for scaling
196    virtual void transposeTimes(double scalar,
197                                const double * x, double * y,
198                                const double * rowScale, 
199                                const double * columnScale,
200                                double * spare=NULL) const;
201    /** Return <code>y - pi * A</code> in <code>y</code>.
202        @pre <code>pi</code> must be of size <code>numRows()</code>
203        @pre <code>y</code> must be of size <code>numColumns()</code>
204        This just does subset (but puts in correct place in y) */
205  void transposeTimesSubset( int number,
206                             const int * which,
207                             const double * pi, double * y,
208                             const double * rowScale, 
209                             const double * columnScale,
210                             double * spare=NULL) const;
211    /** Return <code>x * scalar * A + y</code> in <code>z</code>.
212        Can use y as temporary array (will be empty at end)
213        Note - If x packed mode - then z packed mode
214        Squashes small elements and knows about ClpSimplex */
215  virtual void transposeTimes(const ClpSimplex * model, double scalar,
216                              const CoinIndexedVector * x,
217                              CoinIndexedVector * y,
218                              CoinIndexedVector * z) const;
219    /** Return <code>x * scalar * A + y</code> in <code>z</code>.
220        Note - If x packed mode - then z packed mode
221        This does by column and knows no gaps
222        Squashes small elements and knows about ClpSimplex */
223  void transposeTimesByColumn(const ClpSimplex * model, double scalar,
224                              const CoinIndexedVector * x,
225                              CoinIndexedVector * y,
226                              CoinIndexedVector * z) const;
227    /** Return <code>x * scalar * A + y</code> in <code>z</code>.
228        Can use y as temporary array (will be empty at end)
229        Note - If x packed mode - then z packed mode
230        Squashes small elements and knows about ClpSimplex.
231    This version uses row copy*/
232  virtual void transposeTimesByRow(const ClpSimplex * model, double scalar,
233                              const CoinIndexedVector * x,
234                              CoinIndexedVector * y,
235                              CoinIndexedVector * z) const;
236    /** Return <code>x *A</code> in <code>z</code> but
237        just for indices in y.
238        Note - z always packed mode */
239  virtual void subsetTransposeTimes(const ClpSimplex * model,
240                                    const CoinIndexedVector * x,
241                                    const CoinIndexedVector * y,
242                                    CoinIndexedVector * z) const;
243  /** Returns true if can combine transposeTimes and subsetTransposeTimes
244      and if it would be faster */
245  virtual bool canCombine(const ClpSimplex * model,
246                          const CoinIndexedVector * pi) const;
247  /// Updates two arrays for steepest
248  virtual void transposeTimes2(const ClpSimplex * model,
249                               const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
250                               const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
251                               CoinIndexedVector * spare,
252                               double referenceIn, double devex,
253                               // Array for exact devex to say what is in reference framework
254                               unsigned int * reference,
255                               double * weights, double scaleFactor);
256  /// Updates second array for steepest and does devex weights
257  virtual void subsetTimes2(const ClpSimplex * model,
258                                CoinIndexedVector * dj1,
259                               const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
260                               double referenceIn, double devex,
261                               // Array for exact devex to say what is in reference framework
262                               unsigned int * reference,
263                               double * weights, double scaleFactor);
264  /// Sets up an effective RHS
265  void useEffectiveRhs(ClpSimplex * model);
266#if COIN_LONG_WORK
267  // For long double versions
268  virtual void times(CoinWorkDouble scalar,
269                     const CoinWorkDouble * x, CoinWorkDouble * y) const ;
270  virtual void transposeTimes(CoinWorkDouble scalar,
271                              const CoinWorkDouble * x, CoinWorkDouble * y) const ;
272#endif
273//@}
274
275  /**@name Other */
276   //@{
277  /// Returns CoinPackedMatrix (non const)
278  inline CoinPackedMatrix * matrix() const { return matrix_;}
279  /** Just sets matrix_ to NULL so it can be used elsewhere.
280      used in GUB
281  */
282  inline void setMatrixNull()
283  { matrix_=NULL;}
284  /// Say we want special column copy
285  inline void makeSpecialColumnCopy()
286  { flags_ |= 16;}
287  /// Say we don't want special column copy
288  void releaseSpecialColumnCopy();
289  /// Are there zeros?
290  inline bool zeros() const
291  { return ((flags_&1)!=0);}
292  /// Do we want special column copy
293  inline bool wantsSpecialColumnCopy() const
294  { return ((flags_&16)!=0);}
295  /// Flags
296  inline int flags() const
297  { return flags_;}
298  /// Sets flags_ correctly
299  inline void checkGaps()
300  { flags_ = (matrix_->hasGaps()) ? (flags_|2) : (flags_&(~2));} 
301   //@}
302
303
304  /**@name Constructors, destructor */
305   //@{
306   /** Default constructor. */
307   ClpPackedMatrix();
308   /** Destructor */
309   virtual ~ClpPackedMatrix();
310   //@}
311
312   /**@name Copy method */
313   //@{
314   /** The copy constructor. */
315   ClpPackedMatrix(const ClpPackedMatrix&);
316   /** The copy constructor from an CoinPackedMatrix. */
317   ClpPackedMatrix(const CoinPackedMatrix&);
318  /** Subset constructor (without gaps).  Duplicates are allowed
319      and order is as given */
320  ClpPackedMatrix (const ClpPackedMatrix & wholeModel,
321                    int numberRows, const int * whichRows,
322                    int numberColumns, const int * whichColumns);
323  ClpPackedMatrix (const CoinPackedMatrix & wholeModel,
324                    int numberRows, const int * whichRows,
325                    int numberColumns, const int * whichColumns);
326
327  /** This takes over ownership (for space reasons) */
328   ClpPackedMatrix(CoinPackedMatrix * matrix);
329
330   ClpPackedMatrix& operator=(const ClpPackedMatrix&);
331  /// Clone
332  virtual ClpMatrixBase * clone() const ;
333  /// Copy contents - resizing if necessary - otherwise re-use memory
334  virtual void copy(const ClpPackedMatrix * from);
335  /** Subset clone (without gaps).  Duplicates are allowed
336      and order is as given */
337  virtual ClpMatrixBase * subsetClone (
338                    int numberRows, const int * whichRows,
339                    int numberColumns, const int * whichColumns) const ;
340  /// make special row copy
341  void specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy); 
342  /// make special column copy
343  void specialColumnCopy(ClpSimplex * model); 
344  /// Correct sequence in and out to give true value
345  virtual void correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) ;
346   //@}
347private:
348  /// Meat of transposeTimes by column when not scaled
349  int gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
350                                    int * COIN_RESTRICT index, 
351                                    double * COIN_RESTRICT array,
352                                    const double tolerance) const;
353  /// Meat of transposeTimes by column when scaled
354  int gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
355                                 const double * COIN_RESTRICT columnScale,
356                                 int * COIN_RESTRICT index, 
357                                 double * COIN_RESTRICT array,
358                                 const double tolerance) const;
359  /// Meat of transposeTimes by row n > K if packed - returns number nonzero
360  int gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector, 
361                                   int * COIN_RESTRICT index, 
362                                   double * COIN_RESTRICT output,
363                                   int numberColumns,
364                                   const double tolerance, 
365                                   const double scalar) const;
366  /// Meat of transposeTimes by row n > 2 if packed - returns number nonzero
367  int gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector, 
368                                   int * COIN_RESTRICT index, 
369                                   double * COIN_RESTRICT output,
370                                   int * COIN_RESTRICT lookup,
371                                   char * COIN_RESTRICT marked,
372                                   const double tolerance, 
373                                   const double scalar) const;
374  /// Meat of transposeTimes by row n == 2 if packed
375  void gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
376                                   CoinIndexedVector * spareVector, const double tolerance, const double scalar) const;
377  /// Meat of transposeTimes by row n == 1 if packed
378  void gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
379                                   const double tolerance, const double scalar) const;
380  /// Gets rid of special copies
381  void clearCopies();
382   
383   
384protected:
385  /// Check validity
386  void checkFlags(int type) const;
387   /**@name Data members
388      The data members are protected to allow access for derived classes. */
389   //@{
390  /// Data
391  CoinPackedMatrix * matrix_;
392  /// number of active columns (normally same as number of columns)
393  int numberActiveColumns_;
394  /** Flags -
395      1 - has zero elements
396      2 - has gaps
397      4 - has special row copy
398      8 - has special column copy
399      16 - wants special column copy
400  */
401  int flags_;
402  /// Special row copy
403  ClpPackedMatrix2 * rowCopy_;
404  /// Special column copy
405  ClpPackedMatrix3 * columnCopy_;
406   //@}
407};
408#ifdef THREAD
409#include <pthread.h>
410typedef struct {
411  double acceptablePivot;
412  const ClpSimplex * model;
413  double * spare;
414  int * spareIndex;
415  double * arrayTemp;
416  int * indexTemp;
417  int * numberInPtr;
418  double * bestPossiblePtr;
419  double * upperThetaPtr;
420  int * posFreePtr;
421  double * freePivotPtr;
422  int * numberOutPtr;
423  const unsigned short * count;
424  const double * pi;
425  const CoinBigIndex * rowStart;
426  const double * element;
427  const unsigned short * column;
428  int offset;
429  int numberInRowArray;
430  int numberLook;
431} dualColumn0Struct;
432#endif
433class ClpPackedMatrix2 {
434 
435public:
436  /**@name Useful methods */
437  //@{
438    /** Return <code>x * -1 * A in <code>z</code>.
439        Note - x packed and z will be packed mode
440        Squashes small elements and knows about ClpSimplex */
441  void transposeTimes(const ClpSimplex * model,
442                      const CoinPackedMatrix * rowCopy,
443                      const CoinIndexedVector * x,
444                      CoinIndexedVector * spareArray,
445                      CoinIndexedVector * z) const;
446  /// Returns true if copy has useful information
447  inline bool usefulInfo() const
448  { return rowStart_!=NULL;}
449  //@}
450
451
452  /**@name Constructors, destructor */
453   //@{
454   /** Default constructor. */
455   ClpPackedMatrix2();
456   /** Constructor from copy. */
457   ClpPackedMatrix2(ClpSimplex * model,const CoinPackedMatrix * rowCopy);
458   /** Destructor */
459   virtual ~ClpPackedMatrix2();
460   //@}
461
462   /**@name Copy method */
463   //@{
464   /** The copy constructor. */
465   ClpPackedMatrix2(const ClpPackedMatrix2&);
466   ClpPackedMatrix2& operator=(const ClpPackedMatrix2&);
467   //@}
468   
469   
470protected:
471   /**@name Data members
472      The data members are protected to allow access for derived classes. */
473   //@{
474  /// Number of blocks
475  int numberBlocks_;
476  /// Number of rows
477  int numberRows_;
478  /// Column offset for each block (plus one at end)
479  int * offset_;
480  /// Counts of elements in each part of row
481  mutable unsigned short * count_;
482  /// Row starts
483  mutable CoinBigIndex * rowStart_;
484  /// columns within block
485  unsigned short * column_;
486  /// work arrays
487  double * work_;
488#ifdef THREAD
489  pthread_t * threadId_;
490  dualColumn0Struct * info_;
491#endif
492   //@}
493};
494typedef struct {
495  CoinBigIndex startElements_; // point to data
496  int startIndices_; // point to column_
497  int numberInBlock_; 
498  int numberPrice_; // at beginning
499  int numberElements_; // number elements per column
500} blockStruct;
501class ClpPackedMatrix3 {
502 
503public:
504  /**@name Useful methods */
505  //@{
506    /** Return <code>x * -1 * A in <code>z</code>.
507        Note - x packed and z will be packed mode
508        Squashes small elements and knows about ClpSimplex */
509  void transposeTimes(const ClpSimplex * model,
510                      const double * pi,
511                      CoinIndexedVector * output) const;
512  /// Updates two arrays for steepest
513  void transposeTimes2(const ClpSimplex * model,
514                       const double * pi, CoinIndexedVector * dj1,
515                       const double * piWeight,
516                       double referenceIn, double devex,
517                       // Array for exact devex to say what is in reference framework
518                       unsigned int * reference,
519                       double * weights, double scaleFactor);
520  //@}
521
522
523  /**@name Constructors, destructor */
524   //@{
525   /** Default constructor. */
526   ClpPackedMatrix3();
527   /** Constructor from copy. */
528   ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy);
529   /** Destructor */
530   virtual ~ClpPackedMatrix3();
531   //@}
532
533   /**@name Copy method */
534   //@{
535   /** The copy constructor. */
536   ClpPackedMatrix3(const ClpPackedMatrix3&);
537   ClpPackedMatrix3& operator=(const ClpPackedMatrix3&);
538   //@}
539   /**@name Sort methods */
540   //@{
541   /** Sort blocks */
542  void sortBlocks(const ClpSimplex * model);
543  /// Swap one variable
544  void swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
545                int iColumn);
546   //@}
547   
548   
549protected:
550   /**@name Data members
551      The data members are protected to allow access for derived classes. */
552   //@{
553  /// Number of blocks
554  int numberBlocks_;
555  /// Number of columns
556  int numberColumns_;
557  /// Column indices and reverse lookup (within block)
558  int * column_;
559  /// Starts for odd/long vectors
560  CoinBigIndex * start_;
561  /// Rows
562  int * row_;
563  /// Elements
564  double * element_;
565  /// Blocks (ordinary start at 0 and go to first block)
566  blockStruct * block_;
567   //@}
568};
569
570#endif
Note: See TracBrowser for help on using the repository browser.