ClpPackedMatrix.hpp

Go to the documentation of this file.

/* $Id: ClpPackedMatrix.hpp 2385 2019-01-06 19:43:06Z unxusr $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#ifndef ClpPackedMatrix_H
#define ClpPackedMatrix_H

#include "CoinPragma.hpp"

#include "ClpMatrixBase.hpp"
#include "ClpPrimalColumnSteepest.hpp"

class ClpPackedMatrix2;
class ClpPackedMatrix3;
class CoinDoubleArrayWithLength;
class ClpPackedMatrix : public ClpMatrixBase {

public:
  virtual CoinPackedMatrix *getPackedMatrix() const
  {
    return matrix_;
  }
  virtual bool isColOrdered() const
  {
    return matrix_->isColOrdered();
  }
  virtual CoinBigIndex getNumElements() const
  {
    return matrix_->getNumElements();
  }
  virtual int getNumCols() const
  {
    return matrix_->getNumCols();
  }
  virtual int getNumRows() const
  {
    return matrix_->getNumRows();
  }

  virtual const double *getElements() const
  {
    return matrix_->getElements();
  }
  inline double *getMutableElements() const
  {
    return matrix_->getMutableElements();
  }
  virtual const int *getIndices() const
  {
    return matrix_->getIndices();
  }

  virtual const CoinBigIndex *getVectorStarts() const
  {
    return matrix_->getVectorStarts();
  }
  virtual const int *getVectorLengths() const
  {
    return matrix_->getVectorLengths();
  }
  virtual int getVectorLength(int index) const
  {
    return matrix_->getVectorSize(index);
  }

  virtual void deleteCols(const int numDel, const int *indDel);
  virtual void deleteRows(const int numDel, const int *indDel);
#ifndef CLP_NO_VECTOR
  virtual void appendCols(int number, const CoinPackedVectorBase *const *columns);
  virtual void appendRows(int number, const CoinPackedVectorBase *const *rows);
#endif

  virtual int appendMatrix(int number, int type,
    const CoinBigIndex *starts, const int *index,
    const double *element, int numberOther = -1);
  virtual void replaceVector(const int index,
    const int numReplace, const double *newElements)
  {
    matrix_->replaceVector(index, numReplace, newElements);
  }
  virtual void modifyCoefficient(int row, int column, double newElement,
    bool keepZero = false)
  {
    matrix_->modifyCoefficient(row, column, newElement, keepZero);
  }
  virtual ClpMatrixBase *reverseOrderedCopy() const;
  virtual int countBasis(const int *whichColumn,
    int &numberColumnBasic);
  virtual void fillBasis(ClpSimplex *model,
    const int *whichColumn,
    int &numberColumnBasic,
    int *row, int *start,
    int *rowCount, int *columnCount,
    CoinFactorizationDouble *element);
  virtual int scale(ClpModel *model, ClpSimplex *simplex = NULL) const;
  virtual void scaleRowCopy(ClpModel *model) const;
  void createScaledMatrix(ClpSimplex *model) const;
  virtual ClpMatrixBase *scaledColumnCopy(ClpModel *model) const;
  virtual bool allElementsInRange(ClpModel *model,
    double smallest, double largest,
    int check = 15);
  virtual void rangeOfElements(double &smallestNegative, double &largestNegative,
    double &smallestPositive, double &largestPositive);

  virtual void unpack(const ClpSimplex *model, CoinIndexedVector *rowArray,
    int column) const;
  virtual void unpackPacked(ClpSimplex *model,
    CoinIndexedVector *rowArray,
    int column) const;
  virtual void add(const ClpSimplex *model, CoinIndexedVector *rowArray,
    int column, double multiplier) const;
  virtual void add(const ClpSimplex *model, double *array,
    int column, double multiplier) const;
  virtual void releasePackedMatrix() const {}
  virtual CoinBigIndex *dubiousWeights(const ClpSimplex *model, int *inputWeights) const;
  virtual bool canDoPartialPricing() const;
  virtual void partialPricing(ClpSimplex *model, double start, double end,
    int &bestSequence, int &numberWanted);
  virtual int refresh(ClpSimplex *model);
  // Really scale matrix
  virtual void reallyScale(const double *rowScale, const double *columnScale);
  virtual void setDimensions(int numrows, int numcols);

  virtual void times(double scalar,
    const double *x, double *y) const;
  virtual void times(double scalar,
    const double *x, double *y,
    const double *rowScale,
    const double *columnScale) const;
  virtual void transposeTimes(double scalar,
    const double *x, double *y) const;
  virtual void transposeTimes(double scalar,
    const double *x, double *y,
    const double *rowScale,
    const double *columnScale,
    double *spare = NULL) const;
  void transposeTimesSubset(int number,
    const int *which,
    const double *pi, double *y,
    const double *rowScale,
    const double *columnScale,
    double *spare = NULL) const;
  virtual void transposeTimes(const ClpSimplex *model, double scalar,
    const CoinIndexedVector *x,
    CoinIndexedVector *y,
    CoinIndexedVector *z) const;
  void transposeTimesByColumn(const ClpSimplex *model, double scalar,
    const CoinIndexedVector *x,
    CoinIndexedVector *y,
    CoinIndexedVector *z) const;
  virtual void transposeTimesByRow(const ClpSimplex *model, double scalar,
    const CoinIndexedVector *x,
    CoinIndexedVector *y,
    CoinIndexedVector *z) const;
  virtual void subsetTransposeTimes(const ClpSimplex *model,
    const CoinIndexedVector *x,
    const CoinIndexedVector *y,
    CoinIndexedVector *z) const;
  virtual bool canCombine(const ClpSimplex *model,
    const CoinIndexedVector *pi) const;
  virtual int transposeTimes2(const ClpSimplex *model,
    const CoinIndexedVector *pi1, CoinIndexedVector *dj1,
    const CoinIndexedVector *pi2,
    CoinIndexedVector *spare,
    double *infeas, double *reducedCost,
    double referenceIn, double devex,
    // Array for exact devex to say what is in reference framework
    unsigned int *reference,
    double *weights, double scaleFactor);
  virtual void subsetTimes2(const ClpSimplex *model,
    CoinIndexedVector *dj1,
    const CoinIndexedVector *pi2, CoinIndexedVector *dj2,
    double referenceIn, double devex,
    // Array for exact devex to say what is in reference framework
    unsigned int *reference,
    double *weights, double scaleFactor);
  void useEffectiveRhs(ClpSimplex *model);
#if COIN_LONG_WORK
  // For long double versions
  virtual void times(CoinWorkDouble scalar,
    const CoinWorkDouble *x, CoinWorkDouble *y) const;
  virtual void transposeTimes(CoinWorkDouble scalar,
    const CoinWorkDouble *x, CoinWorkDouble *y) const;
#endif


  inline CoinPackedMatrix *matrix() const
  {
    return matrix_;
  }
  inline void setMatrixNull()
  {
    matrix_ = NULL;
  }
  inline void makeSpecialColumnCopy()
  {
    flags_ |= 16;
  }
  void releaseSpecialColumnCopy();
  inline bool zeros() const
  {
    return ((flags_ & 1) != 0);
  }
  inline bool wantsSpecialColumnCopy() const
  {
    return ((flags_ & 16) != 0);
  }
  inline int flags() const
  {
    return flags_;
  }
  inline void checkGaps()
  {
    flags_ = (matrix_->hasGaps()) ? (flags_ | 2) : (flags_ & (~2));
  }
  inline int numberActiveColumns() const
  {
    return numberActiveColumns_;
  }
  inline void setNumberActiveColumns(int value)
  {
    numberActiveColumns_ = value;
  }

  ClpPackedMatrix();
  virtual ~ClpPackedMatrix();

  ClpPackedMatrix(const ClpPackedMatrix &);
  ClpPackedMatrix(const CoinPackedMatrix &);
  ClpPackedMatrix(const ClpPackedMatrix &wholeModel,
    int numberRows, const int *whichRows,
    int numberColumns, const int *whichColumns);
  ClpPackedMatrix(const CoinPackedMatrix &wholeModel,
    int numberRows, const int *whichRows,
    int numberColumns, const int *whichColumns);

  ClpPackedMatrix(CoinPackedMatrix *matrix);

  ClpPackedMatrix &operator=(const ClpPackedMatrix &);
  virtual ClpMatrixBase *clone() const;
  virtual void copy(const ClpPackedMatrix *from);
  virtual ClpMatrixBase *subsetClone(
    int numberRows, const int *whichRows,
    int numberColumns, const int *whichColumns) const;
  void specialRowCopy(ClpSimplex *model, const ClpMatrixBase *rowCopy);
  void specialColumnCopy(ClpSimplex *model);
  virtual void correctSequence(const ClpSimplex *model, int &sequenceIn, int &sequenceOut);
private:
  int gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
    int *COIN_RESTRICT index,
    double *COIN_RESTRICT array,
    const double tolerance) const;
  int gutsOfTransposeTimesScaled(const double *COIN_RESTRICT pi,
    const double *COIN_RESTRICT columnScale,
    int *COIN_RESTRICT index,
    double *COIN_RESTRICT array,
    const double tolerance) const;
  int gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
    int *COIN_RESTRICT index,
    double *COIN_RESTRICT array,
    const unsigned char *status,
    const double tolerance) const;
  int gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi,
    int *COIN_RESTRICT index,
    double *COIN_RESTRICT array,
    const unsigned char *status,
    int *COIN_RESTRICT spareIndex,
    double *COIN_RESTRICT spareArray,
    const double *COIN_RESTRICT reducedCost,
    double &upperTheta,
    double acceptablePivot,
    double dualTolerance,
    int &numberRemaining,
    const double zeroTolerance) const;
  int gutsOfTransposeTimesScaled(const double *COIN_RESTRICT pi,
    const double *COIN_RESTRICT columnScale,
    int *COIN_RESTRICT index,
    double *COIN_RESTRICT array,
    const unsigned char *status,
    const double tolerance) const;
  int gutsOfTransposeTimesByRowGEK(const CoinIndexedVector *COIN_RESTRICT piVector,
    int *COIN_RESTRICT index,
    double *COIN_RESTRICT output,
    int numberColumns,
    const double tolerance,
    const double scalar) const;
  int gutsOfTransposeTimesByRowGE3(const CoinIndexedVector *COIN_RESTRICT piVector,
    int *COIN_RESTRICT index,
    double *COIN_RESTRICT output,
    double *COIN_RESTRICT array2,
    const double tolerance,
    const double scalar) const;
  int gutsOfTransposeTimesByRowGE3a(const CoinIndexedVector *COIN_RESTRICT piVector,
    int *COIN_RESTRICT index,
    double *COIN_RESTRICT output,
    int *COIN_RESTRICT lookup,
    char *COIN_RESTRICT marked,
    const double tolerance,
    const double scalar) const;
  void gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector *piVector, CoinIndexedVector *output,
    CoinIndexedVector *spareVector, const double tolerance, const double scalar) const;
  void gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector *piVector, CoinIndexedVector *output,
    const double tolerance, const double scalar) const;
  void clearCopies();

protected:
  void checkFlags(int type) const;
  CoinPackedMatrix *matrix_;
  int numberActiveColumns_;
  mutable int flags_;
  ClpPackedMatrix2 *rowCopy_;
  ClpPackedMatrix3 *columnCopy_;
};
#ifdef THREAD
#include <pthread.h>
typedef struct {
  double acceptablePivot;
  const ClpSimplex *model;
  double *spare;
  int *spareIndex;
  double *arrayTemp;
  int *indexTemp;
  int *numberInPtr;
  //double * bestPossiblePtr;
  double *upperThetaPtr;
  int *posFreePtr;
  double *freePivotPtr;
  int *numberOutPtr;
  const unsigned short *count;
  const double *pi;
  const CoinBigIndex *rowStart;
  const double *element;
  const unsigned short *column;
  int offset;
  int numberInRowArray;
  int numberLook;
} dualColumn0Struct;
#endif
class ClpPackedMatrix2 {

public:
  void transposeTimes(const ClpSimplex *model,
    const CoinPackedMatrix *rowCopy,
    const CoinIndexedVector *x,
    CoinIndexedVector *spareArray,
    CoinIndexedVector *z) const;
  inline bool usefulInfo() const
  {
    return rowStart_ != NULL;
  }

  ClpPackedMatrix2();
  ClpPackedMatrix2(ClpSimplex *model, const CoinPackedMatrix *rowCopy);
  virtual ~ClpPackedMatrix2();

  ClpPackedMatrix2(const ClpPackedMatrix2 &);
  ClpPackedMatrix2 &operator=(const ClpPackedMatrix2 &);

protected:
  int numberBlocks_;
  int numberRows_;
  int *offset_;
  mutable unsigned short *count_;
  mutable CoinBigIndex *rowStart_;
  unsigned short *column_;
  double *work_;
#ifdef THREAD
  pthread_t *threadId_;
  dualColumn0Struct *info_;
#endif

};
typedef struct {
  CoinBigIndex startElements_; // point to data
  CoinBigIndex startRows_; // point to data later
  int startIndices_; // point to column_
  int numberInBlock_;
  int numberScan_; // i.e. miss out basic and fixed
  /* order is -
     free or superbasic
     at lower
     at upper
     fixed or basic */
  int firstAtLower_;
  int firstAtUpper_;
  int firstBasic_; // or fixed
  int numberElements_; // number elements per column
  int numberOnes_; // later
} blockStruct;
class ClpPackedMatrix3 {

public:
  void transposeTimes(const ClpSimplex *model,
    const double *pi,
    CoinIndexedVector *output) const;
  void transposeTimes(const ClpSimplex *model,
    const double *pi,
    CoinIndexedVector *output,
    CoinIndexedVector *candidate,
    const CoinIndexedVector *rowArray) const;
  void transposeTimes2(const ClpSimplex *model,
    const double *pi, CoinIndexedVector *dj1,
    const double *piWeight,
    double *COIN_RESTRICT infeas,
    double *COIN_RESTRICT reducedCost,
    double referenceIn, double devex,
    // Array for exact devex to say what is in reference framework
    unsigned int *reference,
    double *weights, double scaleFactor);

  ClpPackedMatrix3();
  ClpPackedMatrix3(ClpSimplex *model, const CoinPackedMatrix *columnCopy);
  virtual ~ClpPackedMatrix3();

  ClpPackedMatrix3(const ClpPackedMatrix3 &);
  ClpPackedMatrix3 &operator=(const ClpPackedMatrix3 &);

  void sortBlocks(const ClpSimplex *model);
  void swapOne(const ClpSimplex *model, const ClpPackedMatrix *matrix,
    int iColumn);
  void swapOne(int iBlock, int kA, int kB);
  void checkBlocks(const ClpSimplex *model);
  int redoInfeasibilities(const ClpSimplex *model,
    ClpPrimalColumnSteepest *pivotChoose,
    int type);


protected:
  int numberBlocks_;
  int numberColumns_;
  int numberColumnsWithGaps_;
#if ABOCA_LITE
  int numberChunks_;
#endif
  CoinBigIndex numberElements_;
  int maxBlockSize_;
  int *column_;
  CoinBigIndex *start_;
  int *row_;
  double *element_;
  CoinDoubleArrayWithLength *temporary_;
#if ABOCA_LITE
  int endChunk_[2 * ABOCA_LITE + 1];
#endif
  blockStruct *block_;
  int ifActive_;
};
#elif INCLUDE_MATRIX3_PRICING
int iColumn = *column;
column++;
if (fabs(value) > zeroTolerance) {
  double thisWeight = weights[iColumn];
  double pivot = value * scaleFactor;
  double pivotSquared = pivot * pivot;
  thisWeight += pivotSquared * devex + pivot * modification;
  if (thisWeight < DEVEX_TRY_NORM) {
    if (referenceIn < 0.0) {
      // steepest
      thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
    } else {
      // exact
      thisWeight = referenceIn * pivotSquared;
      if (reference(iColumn))
        thisWeight += 1.0;
      thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
    }
  }
  // out basic or fixed
  weights[iColumn] = thisWeight;
  value = reducedCost[iColumn] - value;
  reducedCost[iColumn] = value;
  unsigned char thisStatus = status[iColumn] & 7;
  assert(thisStatus != 0 && thisStatus != 4);
  if (thisStatus == 3) {
    //} else if ((thisStatus&1)!=0) {
    // basic or fixed
    //value=0.0;
  } else {
    assert(thisStatus == 2);
    value = -value;
  }
  if (value < dualTolerance) {
    value *= value;
    if (value > bestRatio * weights[iColumn]) {
      bestSequence = iColumn;
      bestRatio = value / weights[iColumn];
#if NO_CHANGE_MULTIPLIER != 1
      bestRatio2 = bestRatio * NO_CHANGE_MULTIPLIER;
#endif
    }
  }
} else {
  // interesting - was faster without this?!
  value = reducedCost[iColumn];
  unsigned char thisStatus = status[iColumn] & 7;
  assert(thisStatus != 0 && thisStatus != 4);
  if (thisStatus == 3) {
  } else if ((thisStatus & 1) != 0) {
    // basic or fixed
    value = 0.0;
  } else {
    value = -value;
  }
  if (value < dualTolerance) {
    value *= value;
    if (value > bestRatio2 * weights[iColumn]) {
      bestSequence = iColumn;
      bestRatio2 = value / weights[iColumn];
#if NO_CHANGE_MULTIPLIER != 1
      bestRatio = bestRatio2 * INVERSE_MULTIPLIER;
#endif
    }
  }
}
#endif

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/