source: trunk/Clp/src/AbcMatrix.hpp @ 1958

Last change on this file since 1958 was 1910, checked in by stefan, 7 years ago
  • add configure option --enable-aboca={1,2,3,4,yes,no}
  • compile Aboca source only if --enable-aboca set (instead of compiling empty source files)
  • fix svn properties
  • Property svn:keywords set to Id
File size: 21.2 KB
Line 
1/* $Id: AbcMatrix.hpp 1910 2013-01-27 02:00:13Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#ifndef AbcMatrix_H
7#define AbcMatrix_H
8
9#include "CoinPragma.hpp"
10
11#include "ClpMatrixBase.hpp"
12#include "AbcSimplex.hpp"
13#include "CoinAbcHelperFunctions.hpp"
14/** This implements a scaled version of CoinPackedMatrix
15    It may have THREE! copies
16    1) scaled CoinPackedMatrix without gaps
17    2) row copy non-basic,basic, fixed
18    3) vector copy
19*/
20class AbcMatrix2;
21class AbcMatrix3;
22class AbcMatrix {
23 
24public:
25  /**@name Useful methods */
26  //@{
27  /// Return a complete CoinPackedMatrix
28  inline CoinPackedMatrix * getPackedMatrix() const {
29    return matrix_;
30  }
31  /** Whether the packed matrix is column major ordered or not. */
32  inline bool isColOrdered() const {
33    return true;
34  }
35  /** Number of entries in the packed matrix. */
36  inline CoinBigIndex getNumElements() const {
37    return matrix_->getNumElements();
38  }
39  /** Number of columns. */
40  inline int getNumCols() const {
41    assert(matrix_->getNumCols()==model_->numberColumns());return matrix_->getNumCols();
42  }
43  /** Number of rows. */
44  inline int getNumRows() const {
45    assert(matrix_->getNumRows()==model_->numberRows());return matrix_->getNumRows();
46  }
47  /// Sets model
48  void setModel(AbcSimplex * model);
49  /// A vector containing the elements in the packed matrix.
50  inline const double * getElements() const {
51    return matrix_->getElements();
52  }
53  /// Mutable elements
54  inline double * getMutableElements() const {
55    return matrix_->getMutableElements();
56  }
57  /// A vector containing the minor indices of the elements in the packed matrix.
58  inline const int * getIndices() const {
59    return matrix_->getIndices();
60  }
61  /// A vector containing the minor indices of the elements in the packed matrix.
62  inline int * getMutableIndices() const {
63    return matrix_->getMutableIndices();
64  }
65  /// Starts
66  inline const CoinBigIndex * getVectorStarts() const {
67    return matrix_->getVectorStarts();
68  }
69  inline  CoinBigIndex * getMutableVectorStarts() const {
70    return matrix_->getMutableVectorStarts();
71  }
72  /** The lengths of the major-dimension vectors. */
73  inline const int * getVectorLengths() const {
74    return matrix_->getVectorLengths();
75  }
76  /** The lengths of the major-dimension vectors. */
77  inline int * getMutableVectorLengths() const {
78    return matrix_->getMutableVectorLengths();
79  }
80  /// Row starts
81  CoinBigIndex * rowStart() const;
82  /// Row ends
83  CoinBigIndex * rowEnd() const;
84  /// Row elements
85  double * rowElements() const;
86  /// Row columns
87  CoinSimplexInt * rowColumns() const;
88  /** Returns a new matrix in reverse order without gaps */
89  CoinPackedMatrix * reverseOrderedCopy() const;
90  /// Returns number of elements in column part of basis
91  CoinBigIndex countBasis(const int * whichColumn,
92                          int & numberColumnBasic);
93  /// Fills in column part of basis
94  void fillBasis(const int * whichColumn,
95                 int & numberColumnBasic,
96                 int * row, int * start,
97                 int * rowCount, int * columnCount,
98                 CoinSimplexDouble * element);
99  /// Fills in column part of basis
100  void fillBasis(const int * whichColumn,
101                 int & numberColumnBasic,
102                 int * row, int * start,
103                 int * rowCount, int * columnCount,
104                 long double * element);
105  /** Scales and creates row copy
106  */
107  void scale(int numberRowsAlreadyScaled);
108  /// Creates row copy
109  void createRowCopy();
110  /// Take out of useful
111  void takeOutOfUseful(int sequence,CoinIndexedVector & spare);
112  /// Put into useful
113  void putIntofUseful(int sequence,CoinIndexedVector & spare);
114  /// Put in and out for useful
115  void inOutUseful(int sequenceIn,int sequenceOut);
116  /// Make all useful
117  void makeAllUseful(CoinIndexedVector & spare);
118  /// Sort into useful
119  void sortUseful(CoinIndexedVector & spare);
120  /// Move largest in column to beginning (not used as doesn't help factorization)
121  void moveLargestToStart();
122 
123  /** Unpacks a column into an CoinIndexedVector
124   */
125  void unpack(CoinIndexedVector & rowArray,
126              int column) const ;
127  /** Adds multiple of a column (or slack) into an CoinIndexedvector
128      You can use quickAdd to add to vector */
129  void add(CoinIndexedVector & rowArray, int column, double multiplier) const ;
130  //@}
131 
132  /**@name Matrix times vector methods */
133  //@{
134  /** Return <code>y + A * scalar *x</code> in <code>y</code>.
135      @pre <code>x</code> must be of size <code>numColumns()</code>
136      @pre <code>y</code> must be of size <code>numRows()</code> */
137  void timesModifyExcludingSlacks(double scalar,
138             const double * x, double * y) const;
139  /** Return <code>y + A * scalar(+-1) *x</code> in <code>y</code>.
140      @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
141      @pre <code>y</code> must be of size <code>numRows()</code> */
142  void timesModifyIncludingSlacks(double scalar,
143             const double * x, double * y) const;
144  /** Return <code>A * scalar(+-1) *x</code> in <code>y</code>.
145      @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
146      @pre <code>y</code> must be of size <code>numRows()</code> */
147  void timesIncludingSlacks(double scalar,
148             const double * x, double * y) const;
149  /** Return A * scalar(+-1) *x + y</code> in <code>y</code>.
150      @pre <code>x</code> must be of size <code>numRows()</code>
151      @pre <code>y</code> must be of size <code>numRows()+numColumns()</code> */
152  void transposeTimesNonBasic(double scalar,
153                            const double * x, double * y) const;
154  /** Return y - A * x</code> in <code>y</code>.
155      @pre <code>x</code> must be of size <code>numRows()</code>
156      @pre <code>y</code> must be of size <code>numRows()+numColumns()</code> */
157  void transposeTimesAll(const double * x, double * y) const;
158  /** Return y + A * scalar(+-1) *x</code> in <code>y</code>.
159      @pre <code>x</code> must be of size <code>numRows()</code>
160      @pre <code>y</code> must be of size <code>numRows()</code> */
161  void transposeTimesBasic(double scalar,
162                            const double * x, double * y) const;
163  /** Return <code>x * scalar * A/code> in <code>z</code>.
164      Note - x unpacked mode - z packed mode including slacks
165      All these return atLo/atUp first then free/superbasic
166      number of first set returned
167      pivotVariable is extended to have that order
168      reversePivotVariable used to update that list
169      free/superbasic only stored in normal format
170      can use spare array to get this effect
171      may put djs alongside atLo/atUp
172      Squashes small elements and knows about AbcSimplex */
173  int transposeTimesNonBasic(double scalar,
174                      const CoinIndexedVector & x,
175                      CoinIndexedVector & z) const;
176  /// gets sorted tableau row and a possible value of theta
177  double dualColumn1(const CoinIndexedVector & update,
178                     CoinPartitionedVector & tableauRow,
179                     CoinPartitionedVector & candidateList) const;
180  /// gets sorted tableau row and a possible value of theta
181  double dualColumn1Row(int iBlock, double upperThetaSlack, int & freeSequence,
182                        const CoinIndexedVector & update,
183                        CoinPartitionedVector & tableauRow,
184                        CoinPartitionedVector & candidateList) const;
185  /// gets sorted tableau row and a possible value of theta
186  double dualColumn1RowFew(int iBlock, double upperThetaSlack, int & freeSequence,
187                         const CoinIndexedVector & update,
188                         CoinPartitionedVector & tableauRow,
189                         CoinPartitionedVector & candidateList) const;
190  /// gets sorted tableau row and a possible value of theta
191  double dualColumn1Row2(double upperThetaSlack, int & freeSequence,
192                         const CoinIndexedVector & update,
193                         CoinPartitionedVector & tableauRow,
194                         CoinPartitionedVector & candidateList) const;
195  /// gets sorted tableau row and a possible value of theta
196  double dualColumn1Row1(double upperThetaSlack, int & freeSequence,
197                         const CoinIndexedVector & update,
198                         CoinPartitionedVector & tableauRow,
199                         CoinPartitionedVector & candidateList) const;
200  /** gets sorted tableau row and a possible value of theta
201      On input first,last give what to scan
202      On output is number in tableauRow and candidateList */
203  void dualColumn1Part(int iBlock,int & sequenceIn, double & upperTheta,
204                         const CoinIndexedVector & update,
205                         CoinPartitionedVector & tableauRow,
206                         CoinPartitionedVector & candidateList) const;
207  /// rebalance for parallel
208  void rebalance() const;
209  /// Get sequenceIn when Dantzig
210  int pivotColumnDantzig(const CoinIndexedVector & updates,
211                         CoinPartitionedVector & spare) const; 
212  /// Get sequenceIn when Dantzig (One block)
213  int pivotColumnDantzig(int iBlock,bool doByRow,const CoinIndexedVector & updates,
214                         CoinPartitionedVector & spare,
215                         double & bestValue) const; 
216  /// gets tableau row - returns number of slacks in block
217  int primalColumnRow(int iBlock,bool doByRow,const CoinIndexedVector & update,
218                          CoinPartitionedVector & tableauRow) const;
219  /// gets tableau row and dj row - returns number of slacks in block
220  int primalColumnRowAndDjs(int iBlock,const CoinIndexedVector & updateTableau,
221                          const CoinIndexedVector & updateDjs,
222                            CoinPartitionedVector & tableauRow) const;
223  /** Chooses best weighted dj
224   */
225  int chooseBestDj(int iBlock,const CoinIndexedVector & infeasibilities, 
226                   const double * weights) const;
227  /** does steepest edge double or triple update
228      If scaleFactor!=0 then use with tableau row to update djs
229      otherwise use updateForDjs
230      Returns best sequence
231   */
232  int primalColumnDouble(int iBlock,CoinPartitionedVector & updateForTableauRow,
233                          CoinPartitionedVector & updateForDjs,
234                          const CoinIndexedVector & updateForWeights,
235                          CoinPartitionedVector & spareColumn1,
236                          double * infeasibilities, 
237                          double referenceIn, double devex,
238                          // Array for exact devex to say what is in reference framework
239                          unsigned int * reference,
240                          double * weights, double scaleFactor) const;
241  /** does steepest edge double or triple update
242      If scaleFactor!=0 then use with tableau row to update djs
243      otherwise use updateForDjs
244      Returns best sequence
245   */
246  int primalColumnSparseDouble(int iBlock,CoinPartitionedVector & updateForTableauRow,
247                          CoinPartitionedVector & updateForDjs,
248                          const CoinIndexedVector & updateForWeights,
249                          CoinPartitionedVector & spareColumn1,
250                          double * infeasibilities, 
251                          double referenceIn, double devex,
252                          // Array for exact devex to say what is in reference framework
253                          unsigned int * reference,
254                          double * weights, double scaleFactor) const;
255  /** does steepest edge double or triple update
256      If scaleFactor!=0 then use with tableau row to update djs
257      otherwise use updateForDjs
258      Returns best sequence
259   */
260  int primalColumnDouble(CoinPartitionedVector & updateForTableauRow,
261                         CoinPartitionedVector & updateForDjs,
262                         const CoinIndexedVector & updateForWeights,
263                         CoinPartitionedVector & spareColumn1,
264                         CoinIndexedVector & infeasible, 
265                         double referenceIn, double devex,
266                         // Array for exact devex to say what is in reference framework
267                         unsigned int * reference,
268                         double * weights, double scaleFactor) const;
269  /// gets subset updates
270  void primalColumnSubset(int iBlock,const CoinIndexedVector & update,
271                          const CoinPartitionedVector & tableauRow,
272                          CoinPartitionedVector & weights) const;
273  /// Partial pricing
274  void partialPricing(double startFraction, double endFraction,
275                      int & bestSequence, int & numberWanted);
276  /** Return <code>x *A</code> in <code>z</code> but
277      just for indices Already in z.
278      Note - z always packed mode */
279  void subsetTransposeTimes(const CoinIndexedVector & x,
280                            CoinIndexedVector & z) const;
281  /// Return <code>-x *A</code> in <code>z</code>
282  void transposeTimes(const CoinIndexedVector & x,
283                            CoinIndexedVector & z) const;
284  //@}
285 
286  /**@name Other */
287  //@{
288  /// Returns CoinPackedMatrix (non const)
289  inline CoinPackedMatrix * matrix() const {
290    return matrix_;
291  }
292  /** Partial pricing tuning parameter - minimum number of "objects" to scan.
293      e.g. number of Gub sets but could be number of variables */
294  inline int minimumObjectsScan() const {
295    return minimumObjectsScan_;
296  }
297  inline void setMinimumObjectsScan(int value) {
298    minimumObjectsScan_ = value;
299  }
300  /// Partial pricing tuning parameter - minimum number of negative reduced costs to get
301  inline int minimumGoodReducedCosts() const {
302    return minimumGoodReducedCosts_;
303  }
304  inline void setMinimumGoodReducedCosts(int value) {
305    minimumGoodReducedCosts_ = value;
306  }
307  /// Current start of search space in matrix (as fraction)
308  inline double startFraction() const {
309    return startFraction_;
310  }
311  inline void setStartFraction(double value) {
312    startFraction_ = value;
313  }
314  /// Current end of search space in matrix (as fraction)
315  inline double endFraction() const {
316    return endFraction_;
317  }
318  inline void setEndFraction(double value) {
319    endFraction_ = value;
320  }
321  /// Current best reduced cost
322  inline double savedBestDj() const {
323    return savedBestDj_;
324  }
325  inline void setSavedBestDj(double value) {
326    savedBestDj_ = value;
327  }
328  /// Initial number of negative reduced costs wanted
329  inline int originalWanted() const {
330    return originalWanted_;
331  }
332  inline void setOriginalWanted(int value) {
333    originalWanted_ = value;
334  }
335  /// Current number of negative reduced costs which we still need
336  inline int currentWanted() const {
337    return currentWanted_;
338  }
339  inline void setCurrentWanted(int value) {
340    currentWanted_ = value;
341  }
342  /// Current best sequence
343  inline int savedBestSequence() const {
344    return savedBestSequence_;
345  }
346  inline void setSavedBestSequence(int value) {
347    savedBestSequence_ = value;
348  }
349  /// Start of each column block
350  inline int * startColumnBlock() const
351  {return startColumnBlock_;}
352  /// Start of each block (in stored)
353  inline const int * blockStart() const
354  { return blockStart_;}
355  inline bool gotRowCopy() const
356  { return rowStart_!=0;}
357  /// Start of each block (in stored)
358  inline int blockStart(int block) const
359  { return blockStart_[block];}
360  /// Number of actual column blocks
361  inline int numberColumnBlocks() const
362  { return numberColumnBlocks_;}
363  /// Number of actual row blocks
364  inline int numberRowBlocks() const
365  { return numberRowBlocks_;}
366  //@}
367 
368 
369  /**@name Constructors, destructor */
370  //@{
371  /** Default constructor. */
372  AbcMatrix();
373  /** Destructor */
374  ~AbcMatrix();
375  //@}
376 
377  /**@name Copy method */
378  //@{
379  /** The copy constructor. */
380  AbcMatrix(const AbcMatrix&);
381  /** The copy constructor from an CoinPackedMatrix. */
382  AbcMatrix(const CoinPackedMatrix&);
383  /** Subset constructor (without gaps).  Duplicates are allowed
384      and order is as given */
385  AbcMatrix (const AbcMatrix & wholeModel,
386             int numberRows, const int * whichRows,
387             int numberColumns, const int * whichColumns);
388  AbcMatrix (const CoinPackedMatrix & wholeModel,
389             int numberRows, const int * whichRows,
390             int numberColumns, const int * whichColumns);
391 
392  AbcMatrix& operator=(const AbcMatrix&);
393  /// Copy contents - resizing if necessary - otherwise re-use memory
394  void copy(const AbcMatrix * from);
395  //@}
396private:
397 
398protected:
399  /**@name Data members
400     The data members are protected to allow access for derived classes. */
401  //@{
402  /// Data
403  CoinPackedMatrix * matrix_;
404  /// Model
405  mutable AbcSimplex * model_;
406#if ABC_PARALLEL==0
407#define NUMBER_ROW_BLOCKS 1
408#define NUMBER_COLUMN_BLOCKS 1
409#elif ABC_PARALLEL==1
410#define NUMBER_ROW_BLOCKS 4
411#define NUMBER_COLUMN_BLOCKS 4
412#else
413#define NUMBER_ROW_BLOCKS 8
414#define NUMBER_COLUMN_BLOCKS 8
415#endif
416  /** Start of each row (per block) - last lot are useless
417      first all row starts for block 0, then for block2
418      so NUMBER_ROW_BLOCKS+2 times number rows */
419  CoinBigIndex * rowStart_;
420  /// Values by row
421  double * element_;
422  /// Columns
423  int * column_;
424  /// Start of each column block
425  mutable int startColumnBlock_[NUMBER_COLUMN_BLOCKS+1];
426  /// Start of each block (in stored)
427  int blockStart_[NUMBER_ROW_BLOCKS+1];
428  /// Number of actual column blocks
429  mutable int numberColumnBlocks_;
430  /// Number of actual row blocks
431  int numberRowBlocks_;
432  //#define COUNT_COPY
433#ifdef COUNT_COPY
434#define MAX_COUNT 13
435  /// Start in elements etc
436  CoinBigIndex countStart_[MAX_COUNT+1];
437  /// First column
438  int countFirst_[MAX_COUNT+1];
439  // later int countEndUseful_[MAX_COUNT+1];
440  int * countRealColumn_;
441  // later int * countInverseRealColumn_;
442  CoinBigIndex * countStartLarge_;
443  int * countRow_;
444  double * countElement_;
445  int smallestCount_;
446  int largestCount_;
447#endif
448  /// Special row copy
449  //AbcMatrix2 * rowCopy_;
450  /// Special column copy
451  //AbcMatrix3 * columnCopy_;
452  /// Current start of search space in matrix (as fraction)
453  double startFraction_;
454  /// Current end of search space in matrix (as fraction)
455  double endFraction_;
456  /// Best reduced cost so far
457  double savedBestDj_;
458  /// Initial number of negative reduced costs wanted
459  int originalWanted_;
460  /// Current number of negative reduced costs which we still need
461  int currentWanted_;
462  /// Saved best sequence in pricing
463  int savedBestSequence_;
464  /// Partial pricing tuning parameter - minimum number of "objects" to scan
465  int minimumObjectsScan_;
466  /// Partial pricing tuning parameter - minimum number of negative reduced costs to get
467  int minimumGoodReducedCosts_;
468  //@}
469};
470#ifdef THREAD
471#include <pthread.h>
472typedef struct {
473  double acceptablePivot;
474  const AbcSimplex * model;
475  double * spare;
476  int * spareIndex;
477  double * arrayTemp;
478  int * indexTemp;
479  int * numberInPtr;
480  double * bestPossiblePtr;
481  double * upperThetaPtr;
482  int * posFreePtr;
483  double * freePivotPtr;
484  int * numberOutPtr;
485  const unsigned short * count;
486  const double * pi;
487  const CoinBigIndex * rowStart;
488  const double * element;
489  const unsigned short * column;
490  int offset;
491  int numberInRowArray;
492  int numberLook;
493} dualColumn0Struct;
494#endif
495class AbcMatrix2 {
496 
497public:
498  /**@name Useful methods */
499  //@{
500  /** Return <code>x * -1 * A in <code>z</code>.
501      Note - x packed and z will be packed mode
502      Squashes small elements and knows about AbcSimplex */
503  void transposeTimes(const AbcSimplex * model,
504                      const CoinPackedMatrix * rowCopy,
505                      const CoinIndexedVector & x,
506                      CoinIndexedVector & spareArray,
507                      CoinIndexedVector & z) const;
508  /// Returns true if copy has useful information
509  inline bool usefulInfo() const {
510    return rowStart_ != NULL;
511  }
512  //@}
513 
514 
515  /**@name Constructors, destructor */
516  //@{
517  /** Default constructor. */
518  AbcMatrix2();
519  /** Constructor from copy. */
520  AbcMatrix2(AbcSimplex * model, const CoinPackedMatrix * rowCopy);
521  /** Destructor */
522  ~AbcMatrix2();
523  //@}
524 
525  /**@name Copy method */
526  //@{
527  /** The copy constructor. */
528  AbcMatrix2(const AbcMatrix2&);
529  AbcMatrix2& operator=(const AbcMatrix2&);
530  //@}
531 
532 
533protected:
534  /**@name Data members
535     The data members are protected to allow access for derived classes. */
536  //@{
537  /// Number of blocks
538  int numberBlocks_;
539  /// Number of rows
540  int numberRows_;
541  /// Column offset for each block (plus one at end)
542  int * offset_;
543  /// Counts of elements in each part of row
544  mutable unsigned short * count_;
545  /// Row starts
546  mutable CoinBigIndex * rowStart_;
547  /// columns within block
548  unsigned short * column_;
549  /// work arrays
550  double * work_;
551#ifdef THREAD
552  pthread_t * threadId_;
553  dualColumn0Struct * info_;
554#endif
555  //@}
556};
557typedef struct {
558  CoinBigIndex startElements_; // point to data
559  int startIndices_; // point to column_
560  int numberInBlock_;
561  int numberPrice_; // at beginning
562  int numberElements_; // number elements per column
563} blockStruct3;
564class AbcMatrix3 {
565 
566public:
567  /**@name Useful methods */
568  //@{
569  /** Return <code>x * -1 * A in <code>z</code>.
570      Note - x packed and z will be packed mode
571      Squashes small elements and knows about AbcSimplex */
572  void transposeTimes(const AbcSimplex * model,
573                      const double * pi,
574                      CoinIndexedVector & output) const;
575  /// Updates two arrays for steepest
576  void transposeTimes2(const AbcSimplex * model,
577                       const double * pi, CoinIndexedVector & dj1,
578                       const double * piWeight,
579                       double referenceIn, double devex,
580                       // Array for exact devex to say what is in reference framework
581                       unsigned int * reference,
582                       double * weights, double scaleFactor);
583  //@}
584 
585 
586  /**@name Constructors, destructor */
587  //@{
588  /** Default constructor. */
589  AbcMatrix3();
590  /** Constructor from copy. */
591  AbcMatrix3(AbcSimplex * model, const CoinPackedMatrix * columnCopy);
592  /** Destructor */
593  ~AbcMatrix3();
594  //@}
595 
596  /**@name Copy method */
597  //@{
598  /** The copy constructor. */
599  AbcMatrix3(const AbcMatrix3&);
600  AbcMatrix3& operator=(const AbcMatrix3&);
601  //@}
602  /**@name Sort methods */
603  //@{
604  /** Sort blocks */
605  void sortBlocks(const AbcSimplex * model);
606  /// Swap one variable
607  void swapOne(const AbcSimplex * model, const AbcMatrix * matrix,
608               int iColumn);
609  //@}
610 
611 
612protected:
613  /**@name Data members
614     The data members are protected to allow access for derived classes. */
615  //@{
616  /// Number of blocks
617  int numberBlocks_;
618  /// Number of columns
619  int numberColumns_;
620  /// Column indices and reverse lookup (within block)
621  int * column_;
622  /// Starts for odd/long vectors
623  CoinBigIndex * start_;
624  /// Rows
625  int * row_;
626  /// Elements
627  double * element_;
628  /// Blocks (ordinary start at 0 and go to first block)
629  blockStruct * block_;
630  //@}
631};
632
633#endif
Note: See TracBrowser for help on using the repository browser.