1 /* $Id: ClpPackedMatrix.hpp 2385 2019-01-06 19:43:06Z unxusr $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others.  All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5 
6 #ifndef ClpPackedMatrix_H
7 #define ClpPackedMatrix_H
8 
9 #include "CoinPragma.hpp"
10 
11 #include "ClpMatrixBase.hpp"
12 #include "ClpPrimalColumnSteepest.hpp"
13 
14 /** This implements CoinPackedMatrix as derived from ClpMatrixBase.
15 
16     It adds a few methods that know about model as well as matrix
17 
18     For details see CoinPackedMatrix */
19 
20 class ClpPackedMatrix2;
21 class ClpPackedMatrix3;
22 class CoinDoubleArrayWithLength;
23 class ClpPackedMatrix : public ClpMatrixBase {
24 
25 public:
26   /**@name Useful methods */
27   //@{
28   /// Return a complete CoinPackedMatrix
getPackedMatrix() const29   virtual CoinPackedMatrix *getPackedMatrix() const
30   {
31     return matrix_;
32   }
33   /** Whether the packed matrix is column major ordered or not. */
isColOrdered() const34   virtual bool isColOrdered() const
35   {
36     return matrix_->isColOrdered();
37   }
38   /** Number of entries in the packed matrix. */
getNumElements() const39   virtual CoinBigIndex getNumElements() const
40   {
41     return matrix_->getNumElements();
42   }
43   /** Number of columns. */
getNumCols() const44   virtual int getNumCols() const
45   {
46     return matrix_->getNumCols();
47   }
48   /** Number of rows. */
getNumRows() const49   virtual int getNumRows() const
50   {
51     return matrix_->getNumRows();
52   }
53 
54   /** A vector containing the elements in the packed matrix. Note that there
55          might be gaps in this list, entries that do not belong to any
56          major-dimension vector. To get the actual elements one should look at
57          this vector together with vectorStarts and vectorLengths. */
getElements() const58   virtual const double *getElements() const
59   {
60     return matrix_->getElements();
61   }
62   /// Mutable elements
getMutableElements() const63   inline double *getMutableElements() const
64   {
65     return matrix_->getMutableElements();
66   }
67   /** A vector containing the minor indices of the elements in the packed
68           matrix. Note that there might be gaps in this list, entries that do not
69           belong to any major-dimension vector. To get the actual elements one
70           should look at this vector together with vectorStarts and
71           vectorLengths. */
getIndices() const72   virtual const int *getIndices() const
73   {
74     return matrix_->getIndices();
75   }
76 
getVectorStarts() const77   virtual const CoinBigIndex *getVectorStarts() const
78   {
79     return matrix_->getVectorStarts();
80   }
81   /** The lengths of the major-dimension vectors. */
getVectorLengths() const82   virtual const int *getVectorLengths() const
83   {
84     return matrix_->getVectorLengths();
85   }
86   /** The length of a single major-dimension vector. */
getVectorLength(int index) const87   virtual int getVectorLength(int index) const
88   {
89     return matrix_->getVectorSize(index);
90   }
91 
92   /** Delete the columns whose indices are listed in <code>indDel</code>. */
93   virtual void deleteCols(const int numDel, const int *indDel);
94   /** Delete the rows whose indices are listed in <code>indDel</code>. */
95   virtual void deleteRows(const int numDel, const int *indDel);
96 #ifndef CLP_NO_VECTOR
97   /// Append Columns
98   virtual void appendCols(int number, const CoinPackedVectorBase *const *columns);
99   /// Append Rows
100   virtual void appendRows(int number, const CoinPackedVectorBase *const *rows);
101 #endif
102   /** Append a set of rows/columns to the end of the matrix. Returns number of errors
103          i.e. if any of the new rows/columns contain an index that's larger than the
104          number of columns-1/rows-1 (if numberOther>0) or duplicates
105          If 0 then rows, 1 if columns */
106   virtual int appendMatrix(int number, int type,
107     const CoinBigIndex *starts, const int *index,
108     const double *element, int numberOther = -1);
109   /** Replace the elements of a vector.  The indices remain the same.
110          This is only needed if scaling and a row copy is used.
111          At most the number specified will be replaced.
112          The index is between 0 and major dimension of matrix */
replaceVector(const int index,const int numReplace,const double * newElements)113   virtual void replaceVector(const int index,
114     const int numReplace, const double *newElements)
115   {
116     matrix_->replaceVector(index, numReplace, newElements);
117   }
118   /** Modify one element of packed matrix.  An element may be added.
119          This works for either ordering If the new element is zero it will be
120          deleted unless keepZero true */
modifyCoefficient(int row,int column,double newElement,bool keepZero=false)121   virtual void modifyCoefficient(int row, int column, double newElement,
122     bool keepZero = false)
123   {
124     matrix_->modifyCoefficient(row, column, newElement, keepZero);
125   }
126   /** Returns a new matrix in reverse order without gaps */
127   virtual ClpMatrixBase *reverseOrderedCopy() const;
128   /// Returns number of elements in column part of basis
129   virtual int countBasis(const int *whichColumn,
130     int &numberColumnBasic);
131   /// Fills in column part of basis
132   virtual void fillBasis(ClpSimplex *model,
133     const int *whichColumn,
134     int &numberColumnBasic,
135     int *row, int *start,
136     int *rowCount, int *columnCount,
137     CoinFactorizationDouble *element);
138   /** Creates scales for column copy (rowCopy in model may be modified)
139          returns non-zero if no scaling done */
140   virtual int scale(ClpModel *model, ClpSimplex *simplex = NULL) const;
141   /** Scales rowCopy if column copy scaled
142          Only called if scales already exist */
143   virtual void scaleRowCopy(ClpModel *model) const;
144   /// Creates scaled column copy if scales exist
145   void createScaledMatrix(ClpSimplex *model) const;
146   /** Realy really scales column copy
147          Only called if scales already exist.
148          Up to user ro delete */
149   virtual ClpMatrixBase *scaledColumnCopy(ClpModel *model) const;
150   /** Checks if all elements are in valid range.  Can just
151          return true if you are not paranoid.  For Clp I will
152          probably expect no zeros.  Code can modify matrix to get rid of
153          small elements.
154          check bits (can be turned off to save time) :
155          1 - check if matrix has gaps
156          2 - check if zero elements
157          4 - check and compress duplicates
158          8 - report on large and small
159      */
160   virtual bool allElementsInRange(ClpModel *model,
161     double smallest, double largest,
162     int check = 15);
163   /** Returns largest and smallest elements of both signs.
164          Largest refers to largest absolute value.
165      */
166   virtual void rangeOfElements(double &smallestNegative, double &largestNegative,
167     double &smallestPositive, double &largestPositive);
168 
169   /** Unpacks a column into an CoinIndexedvector
170       */
171   virtual void unpack(const ClpSimplex *model, CoinIndexedVector *rowArray,
172     int column) const;
173   /** Unpacks a column into an CoinIndexedvector
174       ** in packed foramt
175          Note that model is NOT const.  Bounds and objective could
176          be modified if doing column generation (just for this variable) */
177   virtual void unpackPacked(ClpSimplex *model,
178     CoinIndexedVector *rowArray,
179     int column) const;
180   /** Adds multiple of a column into an CoinIndexedvector
181          You can use quickAdd to add to vector */
182   virtual void add(const ClpSimplex *model, CoinIndexedVector *rowArray,
183     int column, double multiplier) const;
184   /** Adds multiple of a column into an array */
185   virtual void add(const ClpSimplex *model, double *array,
186     int column, double multiplier) const;
187   /// Allow any parts of a created CoinPackedMatrix to be deleted
releasePackedMatrix() const188   virtual void releasePackedMatrix() const {}
189   /** Given positive integer weights for each row fills in sum of weights
190          for each column (and slack).
191          Returns weights vector
192      */
193   virtual CoinBigIndex *dubiousWeights(const ClpSimplex *model, int *inputWeights) const;
194   /// Says whether it can do partial pricing
195   virtual bool canDoPartialPricing() const;
196   /// Partial pricing
197   virtual void partialPricing(ClpSimplex *model, double start, double end,
198     int &bestSequence, int &numberWanted);
199   /// makes sure active columns correct
200   virtual int refresh(ClpSimplex *model);
201   // Really scale matrix
202   virtual void reallyScale(const double *rowScale, const double *columnScale);
203   /** Set the dimensions of the matrix. In effect, append new empty
204          columns/rows to the matrix. A negative number for either dimension
205          means that that dimension doesn't change. Otherwise the new dimensions
206          MUST be at least as large as the current ones otherwise an exception
207          is thrown. */
208   virtual void setDimensions(int numrows, int numcols);
209   //@}
210 
211   /**@name Matrix times vector methods */
212   //@{
213   /** Return <code>y + A * scalar *x</code> in <code>y</code>.
214          @pre <code>x</code> must be of size <code>numColumns()</code>
215          @pre <code>y</code> must be of size <code>numRows()</code> */
216   virtual void times(double scalar,
217     const double *x, double *y) const;
218   /// And for scaling
219   virtual void times(double scalar,
220     const double *x, double *y,
221     const double *rowScale,
222     const double *columnScale) const;
223   /** Return <code>y + x * scalar * A</code> in <code>y</code>.
224          @pre <code>x</code> must be of size <code>numRows()</code>
225          @pre <code>y</code> must be of size <code>numColumns()</code> */
226   virtual void transposeTimes(double scalar,
227     const double *x, double *y) const;
228   /// And for scaling
229   virtual void transposeTimes(double scalar,
230     const double *x, double *y,
231     const double *rowScale,
232     const double *columnScale,
233     double *spare = NULL) const;
234   /** Return <code>y - pi * A</code> in <code>y</code>.
235          @pre <code>pi</code> must be of size <code>numRows()</code>
236          @pre <code>y</code> must be of size <code>numColumns()</code>
237      This just does subset (but puts in correct place in y) */
238   void transposeTimesSubset(int number,
239     const int *which,
240     const double *pi, double *y,
241     const double *rowScale,
242     const double *columnScale,
243     double *spare = NULL) const;
244   /** Return <code>x * scalar * A in <code>z</code>.
245      Can use y as temporary array (will be empty at end)
246      Note - If x packed mode - then z packed mode
247      Squashes small elements and knows about ClpSimplex */
248   virtual void transposeTimes(const ClpSimplex *model, double scalar,
249     const CoinIndexedVector *x,
250     CoinIndexedVector *y,
251     CoinIndexedVector *z) const;
252   /** Return <code>x * scalar * A in <code>z</code>.
253      Note - If x packed mode - then z packed mode
254      This does by column and knows no gaps
255      Squashes small elements and knows about ClpSimplex */
256   void transposeTimesByColumn(const ClpSimplex *model, double scalar,
257     const CoinIndexedVector *x,
258     CoinIndexedVector *y,
259     CoinIndexedVector *z) const;
260   /** Return <code>x * scalar * A in <code>z</code>.
261      Can use y as temporary array (will be empty at end)
262      Note - If x packed mode - then z packed mode
263      Squashes small elements and knows about ClpSimplex.
264      This version uses row copy*/
265   virtual void transposeTimesByRow(const ClpSimplex *model, double scalar,
266     const CoinIndexedVector *x,
267     CoinIndexedVector *y,
268     CoinIndexedVector *z) const;
269   /** Return <code>x *A</code> in <code>z</code> but
270      just for indices in y.
271      Note - z always packed mode */
272   virtual void subsetTransposeTimes(const ClpSimplex *model,
273     const CoinIndexedVector *x,
274     const CoinIndexedVector *y,
275     CoinIndexedVector *z) const;
276   /** Returns true if can combine transposeTimes and subsetTransposeTimes
277          and if it would be faster */
278   virtual bool canCombine(const ClpSimplex *model,
279     const CoinIndexedVector *pi) const;
280   /** Updates two arrays for steepest and does devex weights
281 	 Returns nonzero if updates reduced cost and infeas -
282 	 new infeas in dj1 */
283   virtual int transposeTimes2(const ClpSimplex *model,
284     const CoinIndexedVector *pi1, CoinIndexedVector *dj1,
285     const CoinIndexedVector *pi2,
286     CoinIndexedVector *spare,
287     double *infeas, double *reducedCost,
288     double referenceIn, double devex,
289     // Array for exact devex to say what is in reference framework
290     unsigned int *reference,
291     double *weights, double scaleFactor);
292   /// Updates second array for steepest and does devex weights
293   virtual void subsetTimes2(const ClpSimplex *model,
294     CoinIndexedVector *dj1,
295     const CoinIndexedVector *pi2, CoinIndexedVector *dj2,
296     double referenceIn, double devex,
297     // Array for exact devex to say what is in reference framework
298     unsigned int *reference,
299     double *weights, double scaleFactor);
300   /// Sets up an effective RHS
301   void useEffectiveRhs(ClpSimplex *model);
302 #if COIN_LONG_WORK
303   // For long double versions
304   virtual void times(CoinWorkDouble scalar,
305     const CoinWorkDouble *x, CoinWorkDouble *y) const;
306   virtual void transposeTimes(CoinWorkDouble scalar,
307     const CoinWorkDouble *x, CoinWorkDouble *y) const;
308 #endif
309   //@}
310 
311   /**@name Other */
312   //@{
313   /// Returns CoinPackedMatrix (non const)
matrix() const314   inline CoinPackedMatrix *matrix() const
315   {
316     return matrix_;
317   }
318   /** Just sets matrix_ to NULL so it can be used elsewhere.
319          used in GUB
320      */
setMatrixNull()321   inline void setMatrixNull()
322   {
323     matrix_ = NULL;
324   }
325   /// Say we want special column copy
makeSpecialColumnCopy()326   inline void makeSpecialColumnCopy()
327   {
328     flags_ |= 16;
329   }
330   /// Say we don't want special column copy
331   void releaseSpecialColumnCopy();
332   /// Are there zeros?
zeros() const333   inline bool zeros() const
334   {
335     return ((flags_ & 1) != 0);
336   }
337   /// Do we want special column copy
wantsSpecialColumnCopy() const338   inline bool wantsSpecialColumnCopy() const
339   {
340     return ((flags_ & 16) != 0);
341   }
342   /// Flags
flags() const343   inline int flags() const
344   {
345     return flags_;
346   }
347   /// Sets flags_ correctly
checkGaps()348   inline void checkGaps()
349   {
350     flags_ = (matrix_->hasGaps()) ? (flags_ | 2) : (flags_ & (~2));
351   }
352   /// number of active columns (normally same as number of columns)
numberActiveColumns() const353   inline int numberActiveColumns() const
354   {
355     return numberActiveColumns_;
356   }
357   /// Set number of active columns (normally same as number of columns)
setNumberActiveColumns(int value)358   inline void setNumberActiveColumns(int value)
359   {
360     numberActiveColumns_ = value;
361   }
362   //@}
363 
364   /**@name Constructors, destructor */
365   //@{
366   /** Default constructor. */
367   ClpPackedMatrix();
368   /** Destructor */
369   virtual ~ClpPackedMatrix();
370   //@}
371 
372   /**@name Copy method */
373   //@{
374   /** The copy constructor. */
375   ClpPackedMatrix(const ClpPackedMatrix &);
376   /** The copy constructor from an CoinPackedMatrix. */
377   ClpPackedMatrix(const CoinPackedMatrix &);
378   /** Subset constructor (without gaps).  Duplicates are allowed
379          and order is as given */
380   ClpPackedMatrix(const ClpPackedMatrix &wholeModel,
381     int numberRows, const int *whichRows,
382     int numberColumns, const int *whichColumns);
383   ClpPackedMatrix(const CoinPackedMatrix &wholeModel,
384     int numberRows, const int *whichRows,
385     int numberColumns, const int *whichColumns);
386 
387   /** This takes over ownership (for space reasons) */
388   ClpPackedMatrix(CoinPackedMatrix *matrix);
389 
390   ClpPackedMatrix &operator=(const ClpPackedMatrix &);
391   /// Clone
392   virtual ClpMatrixBase *clone() const;
393   /// Copy contents - resizing if necessary - otherwise re-use memory
394   virtual void copy(const ClpPackedMatrix *from);
395   /** Subset clone (without gaps).  Duplicates are allowed
396          and order is as given */
397   virtual ClpMatrixBase *subsetClone(
398     int numberRows, const int *whichRows,
399     int numberColumns, const int *whichColumns) const;
400   /// make special row copy
401   void specialRowCopy(ClpSimplex *model, const ClpMatrixBase *rowCopy);
402   /// make special column copy
403   void specialColumnCopy(ClpSimplex *model);
404   /// Correct sequence in and out to give true value
405   virtual void correctSequence(const ClpSimplex *model, int &sequenceIn, int &sequenceOut);
406   //@}
407 private:
408   /// Meat of transposeTimes by column when not scaled
409   int gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
410     int *COIN_RESTRICT index,
411     double *COIN_RESTRICT array,
412     const double tolerance) const;
413   /// Meat of transposeTimes by column when scaled
414   int gutsOfTransposeTimesScaled(const double *COIN_RESTRICT pi,
415     const double *COIN_RESTRICT columnScale,
416     int *COIN_RESTRICT index,
417     double *COIN_RESTRICT array,
418     const double tolerance) const;
419   /// Meat of transposeTimes by column when not scaled and skipping
420   int gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
421     int *COIN_RESTRICT index,
422     double *COIN_RESTRICT array,
423     const unsigned char *status,
424     const double tolerance) const;
425   /** Meat of transposeTimes by column when not scaled and skipping
426          and doing part of dualColumn */
427   int gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
428     int *COIN_RESTRICT index,
429     double *COIN_RESTRICT array,
430     const unsigned char *status,
431     int *COIN_RESTRICT spareIndex,
432     double *COIN_RESTRICT spareArray,
433     const double *COIN_RESTRICT reducedCost,
434     double &upperTheta,
435     double acceptablePivot,
436     double dualTolerance,
437     int &numberRemaining,
438     const double zeroTolerance) const;
439   /// Meat of transposeTimes by column when scaled and skipping
440   int gutsOfTransposeTimesScaled(const double *COIN_RESTRICT pi,
441     const double *COIN_RESTRICT columnScale,
442     int *COIN_RESTRICT index,
443     double *COIN_RESTRICT array,
444     const unsigned char *status,
445     const double tolerance) const;
446   /// Meat of transposeTimes by row n > K if packed - returns number nonzero
447   int gutsOfTransposeTimesByRowGEK(const CoinIndexedVector *COIN_RESTRICT piVector,
448     int *COIN_RESTRICT index,
449     double *COIN_RESTRICT output,
450     int numberColumns,
451     const double tolerance,
452     const double scalar) const;
453   /// Meat of transposeTimes by row n > 2 if packed - returns number nonzero
454   int gutsOfTransposeTimesByRowGE3(const CoinIndexedVector *COIN_RESTRICT piVector,
455     int *COIN_RESTRICT index,
456     double *COIN_RESTRICT output,
457     double *COIN_RESTRICT array2,
458     const double tolerance,
459     const double scalar) const;
460   /// Meat of transposeTimes by row n > 2 if packed - returns number nonzero
461   int gutsOfTransposeTimesByRowGE3a(const CoinIndexedVector *COIN_RESTRICT piVector,
462     int *COIN_RESTRICT index,
463     double *COIN_RESTRICT output,
464     int *COIN_RESTRICT lookup,
465     char *COIN_RESTRICT marked,
466     const double tolerance,
467     const double scalar) const;
468   /// Meat of transposeTimes by row n == 2 if packed
469   void gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector *piVector, CoinIndexedVector *output,
470     CoinIndexedVector *spareVector, const double tolerance, const double scalar) const;
471   /// Meat of transposeTimes by row n == 1 if packed
472   void gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector *piVector, CoinIndexedVector *output,
473     const double tolerance, const double scalar) const;
474   /// Gets rid of special copies
475   void clearCopies();
476 
477 protected:
478   /// Check validity
479   void checkFlags(int type) const;
480   /**@name Data members
481         The data members are protected to allow access for derived classes. */
482   //@{
483   /// Data
484   CoinPackedMatrix *matrix_;
485   /// number of active columns (normally same as number of columns)
486   int numberActiveColumns_;
487   /** Flags -
488          1 - has zero elements
489          2 - has gaps
490          4 - has special row copy
491          8 - has special column copy
492          16 - wants special column copy
493      */
494   mutable int flags_;
495   /// Special row copy
496   ClpPackedMatrix2 *rowCopy_;
497   /// Special column copy
498   ClpPackedMatrix3 *columnCopy_;
499   //@}
500 };
501 #ifdef THREAD
502 #include <pthread.h>
503 typedef struct {
504   double acceptablePivot;
505   const ClpSimplex *model;
506   double *spare;
507   int *spareIndex;
508   double *arrayTemp;
509   int *indexTemp;
510   int *numberInPtr;
511   //double * bestPossiblePtr;
512   double *upperThetaPtr;
513   int *posFreePtr;
514   double *freePivotPtr;
515   int *numberOutPtr;
516   const unsigned short *count;
517   const double *pi;
518   const CoinBigIndex *rowStart;
519   const double *element;
520   const unsigned short *column;
521   int offset;
522   int numberInRowArray;
523   int numberLook;
524 } dualColumn0Struct;
525 #endif
526 class ClpPackedMatrix2 {
527 
528 public:
529   /**@name Useful methods */
530   //@{
531   /** Return <code>x * -1 * A in <code>z</code>.
532      Note - x packed and z will be packed mode
533      Squashes small elements and knows about ClpSimplex */
534   void transposeTimes(const ClpSimplex *model,
535     const CoinPackedMatrix *rowCopy,
536     const CoinIndexedVector *x,
537     CoinIndexedVector *spareArray,
538     CoinIndexedVector *z) const;
539   /// Returns true if copy has useful information
usefulInfo() const540   inline bool usefulInfo() const
541   {
542     return rowStart_ != NULL;
543   }
544   //@}
545 
546   /**@name Constructors, destructor */
547   //@{
548   /** Default constructor. */
549   ClpPackedMatrix2();
550   /** Constructor from copy. */
551   ClpPackedMatrix2(ClpSimplex *model, const CoinPackedMatrix *rowCopy);
552   /** Destructor */
553   virtual ~ClpPackedMatrix2();
554   //@}
555 
556   /**@name Copy method */
557   //@{
558   /** The copy constructor. */
559   ClpPackedMatrix2(const ClpPackedMatrix2 &);
560   ClpPackedMatrix2 &operator=(const ClpPackedMatrix2 &);
561   //@}
562 
563 protected:
564   /**@name Data members
565         The data members are protected to allow access for derived classes. */
566   //@{
567   /// Number of blocks
568   int numberBlocks_;
569   /// Number of rows
570   int numberRows_;
571   /// Column offset for each block (plus one at end)
572   int *offset_;
573   /// Counts of elements in each part of row
574   mutable unsigned short *count_;
575   /// Row starts
576   mutable CoinBigIndex *rowStart_;
577   /// columns within block
578   unsigned short *column_;
579   /// work arrays
580   double *work_;
581 #ifdef THREAD
582   pthread_t *threadId_;
583   dualColumn0Struct *info_;
584 #endif
585   //@}
586 };
587 typedef struct {
588   CoinBigIndex startElements_; // point to data
589   CoinBigIndex startRows_; // point to data later
590   int startIndices_; // point to column_
591   int numberInBlock_;
592   int numberScan_; // i.e. miss out basic and fixed
593   /* order is -
594      free or superbasic
595      at lower
596      at upper
597      fixed or basic */
598   int firstAtLower_;
599   int firstAtUpper_;
600   int firstBasic_; // or fixed
601   int numberElements_; // number elements per column
602   int numberOnes_; // later
603 } blockStruct;
604 class ClpPackedMatrix3 {
605 
606 public:
607   /**@name Useful methods */
608   //@{
609   /** Return <code>x * -1 * A in <code>z</code>.
610      Note - x packed and z will be packed mode
611      Squashes small elements and knows about ClpSimplex */
612   void transposeTimes(const ClpSimplex *model,
613     const double *pi,
614     CoinIndexedVector *output) const;
615   /// This version does dualColumn0
616   /// Updates two arrays for steepest
617   void transposeTimes(const ClpSimplex *model,
618     const double *pi,
619     CoinIndexedVector *output,
620     CoinIndexedVector *candidate,
621     const CoinIndexedVector *rowArray) const;
622   void transposeTimes2(const ClpSimplex *model,
623     const double *pi, CoinIndexedVector *dj1,
624     const double *piWeight,
625     double *COIN_RESTRICT infeas,
626     double *COIN_RESTRICT reducedCost,
627     double referenceIn, double devex,
628     // Array for exact devex to say what is in reference framework
629     unsigned int *reference,
630     double *weights, double scaleFactor);
631   //@}
632 
633   /**@name Constructors, destructor */
634   //@{
635   /** Default constructor. */
636   ClpPackedMatrix3();
637   /** Constructor from copy. */
638   ClpPackedMatrix3(ClpSimplex *model, const CoinPackedMatrix *columnCopy);
639   /** Destructor */
640   virtual ~ClpPackedMatrix3();
641   //@}
642 
643   /**@name Copy method */
644   //@{
645   /** The copy constructor. */
646   ClpPackedMatrix3(const ClpPackedMatrix3 &);
647   ClpPackedMatrix3 &operator=(const ClpPackedMatrix3 &);
648   //@}
649   /**@name Sort methods */
650   //@{
651   /** Sort blocks */
652   void sortBlocks(const ClpSimplex *model);
653   /// Swap one variable
654   void swapOne(const ClpSimplex *model, const ClpPackedMatrix *matrix,
655     int iColumn);
656   /// Part of above
657   void swapOne(int iBlock, int kA, int kB);
658   /** Debug - check blocks */
659   void checkBlocks(const ClpSimplex *model);
660   /**
661 	type - 1 redo infeasible, 2 choose sequenceIn, 3 both
662 	returns sequenceIn (or -1) for type 2
663      */
664   int redoInfeasibilities(const ClpSimplex *model,
665     ClpPrimalColumnSteepest *pivotChoose,
666     int type);
667   /// Get temporary array (aligned)
668   //@}
669 
670 protected:
671   /**@name Data members
672         The data members are protected to allow access for derived classes. */
673   //@{
674   /// Number of blocks
675   int numberBlocks_;
676   /// Number of columns
677   int numberColumns_;
678   /// Number of columns including gaps
679   int numberColumnsWithGaps_;
680 #if ABOCA_LITE
681   /// Number of chunks
682   int numberChunks_;
683 #endif
684   /// Number of elements (including gaps)
685   CoinBigIndex numberElements_;
686   /// Maximum size of any block
687   int maxBlockSize_;
688   /// Column indices and reverse lookup (within block)
689   int *column_;
690   /// Starts for odd/long vectors??
691   CoinBigIndex *start_;
692   /// Rows
693   int *row_;
694   /// Elements
695   double *element_;
696   /// Temporary work area (aligned)
697   CoinDoubleArrayWithLength *temporary_;
698 #if ABOCA_LITE
699   /// Chunk ends (could have more than cpus)
700   int endChunk_[2 * ABOCA_LITE + 1];
701 #endif
702   /// Blocks (ordinary start at 0 and go to first block)
703   blockStruct *block_;
704   /// If active
705   int ifActive_;
706   //@}
707 };
708 #elif INCLUDE_MATRIX3_PRICING
709 int iColumn = *column;
710 column++;
711 if (fabs(value) > zeroTolerance) {
712   double thisWeight = weights[iColumn];
713   double pivot = value * scaleFactor;
714   double pivotSquared = pivot * pivot;
715   thisWeight += pivotSquared * devex + pivot * modification;
716   if (thisWeight < DEVEX_TRY_NORM) {
717     if (referenceIn < 0.0) {
718       // steepest
719       thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
720     } else {
721       // exact
722       thisWeight = referenceIn * pivotSquared;
723       if (reference(iColumn))
724         thisWeight += 1.0;
725       thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
726     }
727   }
728   // out basic or fixed
729   weights[iColumn] = thisWeight;
730   value = reducedCost[iColumn] - value;
731   reducedCost[iColumn] = value;
732   unsigned char thisStatus = status[iColumn] & 7;
733   assert(thisStatus != 0 && thisStatus != 4);
734   if (thisStatus == 3) {
735     //} else if ((thisStatus&1)!=0) {
736     // basic or fixed
737     //value=0.0;
738   } else {
739     assert(thisStatus == 2);
740     value = -value;
741   }
742   if (value < dualTolerance) {
743     value *= value;
744     if (value > bestRatio * weights[iColumn]) {
745       bestSequence = iColumn;
746       bestRatio = value / weights[iColumn];
747 #if NO_CHANGE_MULTIPLIER != 1
748       bestRatio2 = bestRatio * NO_CHANGE_MULTIPLIER;
749 #endif
750     }
751   }
752 } else {
753   // interesting - was faster without this?!
754   value = reducedCost[iColumn];
755   unsigned char thisStatus = status[iColumn] & 7;
756   assert(thisStatus != 0 && thisStatus != 4);
757   if (thisStatus == 3) {
758   } else if ((thisStatus & 1) != 0) {
759     // basic or fixed
760     value = 0.0;
761   } else {
762     value = -value;
763   }
764   if (value < dualTolerance) {
765     value *= value;
766     if (value > bestRatio2 * weights[iColumn]) {
767       bestSequence = iColumn;
768       bestRatio2 = value / weights[iColumn];
769 #if NO_CHANGE_MULTIPLIER != 1
770       bestRatio = bestRatio2 * INVERSE_MULTIPLIER;
771 #endif
772     }
773   }
774 }
775 #endif
776 
777 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
778 */
779