Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (3 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpPESimplex.cpp

    r2271 r2385  
    1515#include "ClpMessage.hpp"
    1616
    17 
    1817/** SHARED METHODS FOR USEFUL ALGEBRAIC OPERATIONS */
    1918
    20 /** inner product between a coin vector and a pointer */
    21  double PEdot(CoinIndexedVector& v1, const double* v2){
    22    double sum = 0;
    23    int size = v1.getNumElements();
    24    int* indices = v1.getIndices();
    25    
    26    for (int i = 0; i < size; i++)
    27      sum += v1[indices[i]] * v2[indices[i]];
    28    return sum;
    29  }
     19/** inner product between a coin vector and a pointer */
     20double PEdot(CoinIndexedVector &v1, const double *v2)
     21{
     22  double sum = 0;
     23  int size = v1.getNumElements();
     24  int *indices = v1.getIndices();
     25
     26  for (int i = 0; i < size; i++)
     27    sum += v1[indices[i]] * v2[indices[i]];
     28  return sum;
     29}
    3030
    3131/** inner product between two coin vectors
    32     call the function with the sparser vector first for efficiency */
    33  double PEdot(CoinIndexedVector& v1, CoinIndexedVector& v2){
    34    double sum = 0;
    35    int size = v1.getNumElements();
    36    int* indices = v1.getIndices();
    37 
    38    for (int i = 0; i < size; i++)
    39      sum += v1[indices[i]] * v2[indices[i]];
    40    return sum;
    41  }
    42 
    43   /** compute the product y + x^T*[A I] for the indices "which" of [A I] and store it into y */
    44   void PEtransposeTimesSubsetAll( ClpSimplex* model, int number, const int * which,
    45                                   const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
    46                                   const double * COIN_RESTRICT rowScale,
    47                                   const double * COIN_RESTRICT columnScale)
    48   {
    49     // get the packed matrix
    50     CoinPackedMatrix* clpMatrix = model->matrix();
    51 
    52    // get the matrix data pointers
    53    const int *  row = clpMatrix->getIndices();
    54    const CoinBigIndex *  columnStart = clpMatrix->getVectorStarts();
    55    const int * columnLength = clpMatrix->getVectorLengths();
    56    const double *  elementByColumn = clpMatrix->getElements();
    57 
    58    // there is scaling iff rowScale is not null
    59    if (rowScale) {
    60      for (int jColumn=0;jColumn<number;jColumn++) {
    61        int iColumn = which[jColumn];
    62        CoinBigIndex j;
    63        CoinBigIndex start=columnStart[iColumn];
    64        CoinBigIndex next=start+columnLength[iColumn];
    65        double value=0.0;
    66 
    67        // if looking at a slack variable, there is no scaling (rowScale=1/columnScale)
    68        if (iColumn > model->getNumCols()){
    69          int jRow = iColumn - model->getNumCols();
    70          y[iColumn] = x[jRow]* (-1);
    71        }
    72        // don't forget the scaling for the decision variables
    73        else{
    74          for (j=start;j<next;j++) {
    75            int jRow=row[j];
    76            value += x[jRow]*elementByColumn[j]*rowScale[jRow];
    77          }
    78          y[iColumn] += value*columnScale[iColumn];
    79        }
    80      }
    81    } else {
    82      for (int jColumn=0;jColumn<number;jColumn++) {
    83        int iColumn = which[jColumn];
    84        CoinBigIndex j;
    85        CoinBigIndex start=columnStart[iColumn];
    86        CoinBigIndex next=start+columnLength[iColumn];
    87        double value=0.0;
    88        if (iColumn > model->getNumCols()){
    89          int jRow = iColumn - model->getNumCols();
    90          value = x[jRow]* (-1);
    91        }
    92        else{
    93          for (j=start;j<next;j++) {
    94            int jRow=row[j];
    95            value += x[jRow]*elementByColumn[j];
    96          }
    97        }
    98        y[iColumn] += value;
    99      }
    100    }
    101   }
    102 
     32    call the function with the sparser vector first for efficiency */
     33double PEdot(CoinIndexedVector &v1, CoinIndexedVector &v2)
     34{
     35  double sum = 0;
     36  int size = v1.getNumElements();
     37  int *indices = v1.getIndices();
     38
     39  for (int i = 0; i < size; i++)
     40    sum += v1[indices[i]] * v2[indices[i]];
     41  return sum;
     42}
     43
     44/** compute the product y + x^T*[A I] for the indices "which" of [A I] and store it into y */
     45void PEtransposeTimesSubsetAll(ClpSimplex *model, int number, const int *which,
     46  const double *COIN_RESTRICT x, double *COIN_RESTRICT y,
     47  const double *COIN_RESTRICT rowScale,
     48  const double *COIN_RESTRICT columnScale)
     49{
     50  // get the packed matrix
     51  CoinPackedMatrix *clpMatrix = model->matrix();
     52
     53  // get the matrix data pointers
     54  const int *row = clpMatrix->getIndices();
     55  const CoinBigIndex *columnStart = clpMatrix->getVectorStarts();
     56  const int *columnLength = clpMatrix->getVectorLengths();
     57  const double *elementByColumn = clpMatrix->getElements();
     58
     59  // there is scaling iff rowScale is not null
     60  if (rowScale) {
     61    for (int jColumn = 0; jColumn < number; jColumn++) {
     62      int iColumn = which[jColumn];
     63      CoinBigIndex j;
     64      CoinBigIndex start = columnStart[iColumn];
     65      CoinBigIndex next = start + columnLength[iColumn];
     66      double value = 0.0;
     67
     68      // if looking at a slack variable, there is no scaling (rowScale=1/columnScale)
     69      if (iColumn > model->getNumCols()) {
     70        int jRow = iColumn - model->getNumCols();
     71        y[iColumn] = x[jRow] * (-1);
     72      }
     73      // don't forget the scaling for the decision variables
     74      else {
     75        for (j = start; j < next; j++) {
     76          int jRow = row[j];
     77          value += x[jRow] * elementByColumn[j] * rowScale[jRow];
     78        }
     79        y[iColumn] += value * columnScale[iColumn];
     80      }
     81    }
     82  } else {
     83    for (int jColumn = 0; jColumn < number; jColumn++) {
     84      int iColumn = which[jColumn];
     85      CoinBigIndex j;
     86      CoinBigIndex start = columnStart[iColumn];
     87      CoinBigIndex next = start + columnLength[iColumn];
     88      double value = 0.0;
     89      if (iColumn > model->getNumCols()) {
     90        int jRow = iColumn - model->getNumCols();
     91        value = x[jRow] * (-1);
     92      } else {
     93        for (j = start; j < next; j++) {
     94          int jRow = row[j];
     95          value += x[jRow] * elementByColumn[j];
     96        }
     97      }
     98      y[iColumn] += value;
     99    }
     100  }
     101}
    103102
    104103/** BASE CLASS FOR THE IMPROVED SIMPLEX
    105104*/
    106105
    107  /** Constructor */
    108 ClpPESimplex::ClpPESimplex(ClpSimplex *model): coPrimalDegenerates_(0), primalDegenerates_(NULL), isPrimalDegenerate_(NULL),
    109 coDualDegenerates_(0), dualDegenerates_(NULL), isDualDegenerate_(NULL),
    110 coCompatibleCols_(0),  isCompatibleCol_(NULL), coCompatibleRows_(0),isCompatibleRow_(NULL),
    111 model_(model), epsDegeneracy_(1.0e-07), epsCompatibility_(1.0e-07),
    112 tempRandom_(NULL),
    113 coPrimalDegeneratesAvg_(0), coDualDegeneratesAvg_(0), coCompatibleColsAvg_(0), coCompatibleRowsAvg_(0),
    114 coUpdateDegenerates_(0), coIdentifyCompatibles_(0),
    115 coDegeneratePivots_(0), coCompatiblePivots_(0),
    116 coDegenerateCompatiblePivots_(0),coDegeneratePivotsConsecutive_(0),
    117 coPriorityPivots_(0),doStatistics_(0),
    118 lastObjectiveValue_(COIN_DBL_MAX), isLastPivotCompatible_(false),
    119 timeCompatibility_(0.0), timeMultRandom_(0.0), timeLinearSystem_(0.0), timeTmp_(0.0) {
     106/** Constructor */
     107ClpPESimplex::ClpPESimplex(ClpSimplex *model)
     108  : coPrimalDegenerates_(0)
     109  , primalDegenerates_(NULL)
     110  , isPrimalDegenerate_(NULL)
     111  , coDualDegenerates_(0)
     112  , dualDegenerates_(NULL)
     113  , isDualDegenerate_(NULL)
     114  , coCompatibleCols_(0)
     115  , isCompatibleCol_(NULL)
     116  , coCompatibleRows_(0)
     117  , isCompatibleRow_(NULL)
     118  , model_(model)
     119  , epsDegeneracy_(1.0e-07)
     120  , epsCompatibility_(1.0e-07)
     121  , tempRandom_(NULL)
     122  , coPrimalDegeneratesAvg_(0)
     123  , coDualDegeneratesAvg_(0)
     124  , coCompatibleColsAvg_(0)
     125  , coCompatibleRowsAvg_(0)
     126  , coUpdateDegenerates_(0)
     127  , coIdentifyCompatibles_(0)
     128  , coDegeneratePivots_(0)
     129  , coCompatiblePivots_(0)
     130  , coDegenerateCompatiblePivots_(0)
     131  , coDegeneratePivotsConsecutive_(0)
     132  , coPriorityPivots_(0)
     133  , doStatistics_(0)
     134  , lastObjectiveValue_(COIN_DBL_MAX)
     135  , isLastPivotCompatible_(false)
     136  , timeCompatibility_(0.0)
     137  , timeMultRandom_(0.0)
     138  , timeLinearSystem_(0.0)
     139  , timeTmp_(0.0)
     140{
    120141
    121142  // the improved simplex object should only be implemented after loading the model
     
    124145  // size of the original model
    125146  numberColumns_ = model_->numberColumns();
    126   numberRows_    = model_->numberRows();
     147  numberRows_ = model_->numberRows();
    127148
    128149  // allocate memory to the arrays
    129   primalDegenerates_ = reinterpret_cast<int *> ( malloc(numberRows_*sizeof(int)));
    130   isPrimalDegenerate_ = reinterpret_cast<bool *> ( malloc( (numberRows_+numberColumns_)*sizeof(bool) ));
    131 
    132   dualDegenerates_ = reinterpret_cast<int *> ( malloc(numberColumns_*sizeof(int)));
    133   isDualDegenerate_ = reinterpret_cast<bool *> ( malloc( (numberRows_+numberColumns_)*sizeof(bool) ));
    134 
    135   compatibilityCol_ = reinterpret_cast<double *> ( malloc( (numberRows_+numberColumns_)*sizeof(double)) );
    136   isCompatibleCol_ = reinterpret_cast<bool *> ( malloc( (numberRows_+numberColumns_)*sizeof(bool) ));
    137   std::fill(isCompatibleCol_, isCompatibleCol_+numberRows_+numberColumns_, false);
    138 
    139 
    140   compatibilityRow_ = reinterpret_cast<double *> ( malloc( numberRows_*sizeof(double) ));
    141   isCompatibleRow_ = reinterpret_cast<bool *> ( malloc( numberRows_*sizeof(bool) ));
    142   std::fill(isCompatibleRow_, isCompatibleRow_+numberRows_, false);
    143 
    144   tempRandom_ = reinterpret_cast<double *> ( malloc(CoinMax(numberColumns_,numberRows_)*sizeof(double)));
     150  primalDegenerates_ = reinterpret_cast< int * >(malloc(numberRows_ * sizeof(int)));
     151  isPrimalDegenerate_ = reinterpret_cast< bool * >(malloc((numberRows_ + numberColumns_) * sizeof(bool)));
     152
     153  dualDegenerates_ = reinterpret_cast< int * >(malloc(numberColumns_ * sizeof(int)));
     154  isDualDegenerate_ = reinterpret_cast< bool * >(malloc((numberRows_ + numberColumns_) * sizeof(bool)));
     155
     156  compatibilityCol_ = reinterpret_cast< double * >(malloc((numberRows_ + numberColumns_) * sizeof(double)));
     157  isCompatibleCol_ = reinterpret_cast< bool * >(malloc((numberRows_ + numberColumns_) * sizeof(bool)));
     158  std::fill(isCompatibleCol_, isCompatibleCol_ + numberRows_ + numberColumns_, false);
     159
     160  compatibilityRow_ = reinterpret_cast< double * >(malloc(numberRows_ * sizeof(double)));
     161  isCompatibleRow_ = reinterpret_cast< bool * >(malloc(numberRows_ * sizeof(bool)));
     162  std::fill(isCompatibleRow_, isCompatibleRow_ + numberRows_, false);
     163
     164  tempRandom_ = reinterpret_cast< double * >(malloc(CoinMax(numberColumns_, numberRows_) * sizeof(double)));
    145165  // fill
    146166  CoinThreadRandom generator = *model_->randomNumberGenerator();
    147   for (int i=0;i<CoinMax(numberColumns_,numberRows_);i++) {
     167  for (int i = 0; i < CoinMax(numberColumns_, numberRows_); i++) {
    148168    double random;
    149169    do
    150      random = static_cast<int>(generator.randomDouble()*1.0e6)-5.0e5;
    151     while (random==0.0);
     170      random = static_cast< int >(generator.randomDouble() * 1.0e6) - 5.0e5;
     171    while (random == 0.0);
    152172    //random = static_cast<int>(model_->randomNumberGenerator()->randomDouble()*1.0e7)-5.0e6;
    153173    tempRandom_[i] = random;
    154174  }
    155   if (model_->logLevel()>2)
    156     doStatistics_=model_->logLevel();
     175  if (model_->logLevel() > 2)
     176    doStatistics_ = model_->logLevel();
    157177}
    158178
    159179/** Destructor */
    160 ClpPESimplex::~ClpPESimplex() {
     180ClpPESimplex::~ClpPESimplex()
     181{
    161182
    162183  // free memory
    163   if (primalDegenerates_)  free(primalDegenerates_);
    164   if (isPrimalDegenerate_) free(isPrimalDegenerate_);
    165   if (dualDegenerates_)  free(dualDegenerates_);
    166   if (isDualDegenerate_) free(isDualDegenerate_);
    167   if (isCompatibleCol_) free(isCompatibleCol_);
    168   if (compatibilityCol_) free(compatibilityCol_);
    169   if (isCompatibleRow_) free(isCompatibleRow_);
    170   if (compatibilityRow_) free(compatibilityRow_);
    171   if (tempRandom_)      free(tempRandom_);
    172 
    173 
    174   if (doStatistics_&&model_&&model_->numberIterations()) {
     184  if (primalDegenerates_)
     185    free(primalDegenerates_);
     186  if (isPrimalDegenerate_)
     187    free(isPrimalDegenerate_);
     188  if (dualDegenerates_)
     189    free(dualDegenerates_);
     190  if (isDualDegenerate_)
     191    free(isDualDegenerate_);
     192  if (isCompatibleCol_)
     193    free(isCompatibleCol_);
     194  if (compatibilityCol_)
     195    free(compatibilityCol_);
     196  if (isCompatibleRow_)
     197    free(isCompatibleRow_);
     198  if (compatibilityRow_)
     199    free(compatibilityRow_);
     200  if (tempRandom_)
     201    free(tempRandom_);
     202
     203  if (doStatistics_ && model_ && model_->numberIterations()) {
    175204    char generalPrint[200];
    176205    sprintf(generalPrint, "Degenerate pivots   : %d, compatibility time %.2f",
    177             coDegeneratePivots(),
    178             timeCompatibility());
     206      coDegeneratePivots(),
     207      timeCompatibility());
    179208    model_->messageHandler()->message(CLP_GENERAL,
    180                                       *model_->messagesPointer())
     209      *model_->messagesPointer())
    181210      << generalPrint << CoinMessageEol;
    182211
    183212#if 1
    184     int numberPivots=model_->numberIterations();
    185     if ( coDualDegeneratesAvg() ){
     213    int numberPivots = model_->numberIterations();
     214    if (coDualDegeneratesAvg()) {
    186215      sprintf(generalPrint, "coDegenAvg/rows %g coCompatAvg/rows %g",
    187               coDualDegeneratesAvg()/numberRows_,
    188               coCompatibleRowsAvg()/numberRows_);
     216        coDualDegeneratesAvg() / numberRows_,
     217        coCompatibleRowsAvg() / numberRows_);
    189218      model_->messageHandler()->message(CLP_GENERAL,
    190                                         *model_->messagesPointer())
    191         << generalPrint << CoinMessageEol;
     219        *model_->messagesPointer())
     220        << generalPrint << CoinMessageEol;
    192221    } else if (coPrimalDegeneratesAvg()) {
    193       sprintf(generalPrint, "coDegenAvg/columns %g coCompatAvg/columns %g", 
    194               coPrimalDegeneratesAvg()/numberColumns_,
    195               coCompatibleColsAvg()/numberColumns_);
     222      sprintf(generalPrint, "coDegenAvg/columns %g coCompatAvg/columns %g",
     223        coPrimalDegeneratesAvg() / numberColumns_,
     224        coCompatibleColsAvg() / numberColumns_);
    196225      model_->messageHandler()->message(CLP_GENERAL,
    197                                         *model_->messagesPointer())
    198         << generalPrint << CoinMessageEol;
    199     }   
    200     if (numberPivots-coCompatiblePivots()) {
    201       sprintf(generalPrint,"(coDegeneratePivots()-coDegenerateCompatiblePivots())/( (numberPivots-coCompatiblePivots()) %g",(static_cast<double> (coDegeneratePivots()-coDegenerateCompatiblePivots()))/(static_cast<double> (numberPivots-coCompatiblePivots())));
     226        *model_->messagesPointer())
     227        << generalPrint << CoinMessageEol;
     228    }
     229    if (numberPivots - coCompatiblePivots()) {
     230      sprintf(generalPrint, "(coDegeneratePivots()-coDegenerateCompatiblePivots())/( (numberPivots-coCompatiblePivots()) %g", (static_cast< double >(coDegeneratePivots() - coDegenerateCompatiblePivots())) / (static_cast< double >(numberPivots - coCompatiblePivots())));
    202231      model_->messageHandler()->message(CLP_GENERAL,
    203                                         *model_->messagesPointer())
    204         << generalPrint << CoinMessageEol;
     232        *model_->messagesPointer())
     233        << generalPrint << CoinMessageEol;
    205234    }
    206235    if (coCompatiblePivots()) {
    207       sprintf(generalPrint,  "coDegenerateCompatiblePivots()/coCompatiblePivots() %g",static_cast<double> ( coDegenerateCompatiblePivots())/static_cast<double> (coCompatiblePivots()));
     236      sprintf(generalPrint, "coDegenerateCompatiblePivots()/coCompatiblePivots() %g", static_cast< double >(coDegenerateCompatiblePivots()) / static_cast< double >(coCompatiblePivots()));
    208237      model_->messageHandler()->message(CLP_GENERAL,
    209                                         *model_->messagesPointer())
    210         << generalPrint << CoinMessageEol;
    211     }
    212     sprintf(generalPrint,  "coDegeneratePivots()/ numberPivots %g",static_cast<double> ( coDegeneratePivots())/static_cast<double> ( numberPivots));
     238        *model_->messagesPointer())
     239        << generalPrint << CoinMessageEol;
     240    }
     241    sprintf(generalPrint, "coDegeneratePivots()/ numberPivots %g", static_cast< double >(coDegeneratePivots()) / static_cast< double >(numberPivots));
    213242    model_->messageHandler()->message(CLP_GENERAL,
    214                                       *model_->messagesPointer())
     243      *model_->messagesPointer())
    215244      << generalPrint << CoinMessageEol;
    216     sprintf(generalPrint, "coCompatiblePivots() %d coPriorityPivots() %d", 
    217             coCompatiblePivots(),coPriorityPivots());
     245    sprintf(generalPrint, "coCompatiblePivots() %d coPriorityPivots() %d",
     246      coCompatiblePivots(), coPriorityPivots());
    218247    model_->messageHandler()->message(CLP_GENERAL,
    219                                       *model_->messagesPointer())
     248      *model_->messagesPointer())
    220249      << generalPrint << CoinMessageEol;
    221250    //sprintf(generalPrint, numberPivots);
     
    227256
    228257/** Updates the set of variables that are not at their bounds */
    229 void ClpPESimplex::updatePrimalDegenerates() {
    230 
     258void ClpPESimplex::updatePrimalDegenerates()
     259{
    231260
    232261  coPrimalDegenerates_ = 0;
    233   epsDegeneracy_ = 1.0e-04;//std::min(1.0e-2, 1.0e-7*fabs(model_->objectiveValue()));
    234   std::fill(isPrimalDegenerate_, isPrimalDegenerate_+numberRows_+numberColumns_, false);
     262  epsDegeneracy_ = 1.0e-04; //std::min(1.0e-2, 1.0e-7*fabs(model_->objectiveValue()));
     263  std::fill(isPrimalDegenerate_, isPrimalDegenerate_ + numberRows_ + numberColumns_, false);
    235264  int *pivotVariable = model_->pivotVariable();
    236265
     
    238267  // set to one of their bounds (or zero if the variable is free)
    239268  // An epsDegeneracy_ tolerance is used to detect variables at a bound
    240   for (int i = 0 ; i < numberRows_ ; i++){
     269  for (int i = 0; i < numberRows_; i++) {
    241270    int iVar = pivotVariable[i]; // index of basic variable at row i
    242    
     271
    243272    // std::cout << "solution[" << iVar << "] = " << model_->solution(iVar) << " ; lb = " << model_->lower(iVar) << " ; ub = "<< model_->upper(iVar) << "\n" ;
    244273
     
    247276    double dLb = model_->lower(iVar);
    248277    // if (fabs(dVal) <= epsDegeneracy_) {
    249     if ((dLb > -COIN_DBL_MAX && fabs(dVal-dLb) <= std::max(1.0, fabs(dLb)) * epsDegeneracy_)
    250         || (dUb < COIN_DBL_MAX && fabs(dVal-dUb) <= std::max(1.0, fabs(dUb)) * epsDegeneracy_)) {
     278    if ((dLb > -COIN_DBL_MAX && fabs(dVal - dLb) <= std::max(1.0, fabs(dLb)) * epsDegeneracy_)
     279      || (dUb < COIN_DBL_MAX && fabs(dVal - dUb) <= std::max(1.0, fabs(dUb)) * epsDegeneracy_)) {
    251280      primalDegenerates_[coPrimalDegenerates_++] = i;
    252281      isPrimalDegenerate_[iVar] = true;
     
    257286
    258287/** Updates the set of dual degenerate variables */
    259 void ClpPESimplex::updateDualDegenerates() {
     288void ClpPESimplex::updateDualDegenerates()
     289{
    260290
    261291  coDualDegenerates_ = 0;
    262   std::fill(isDualDegenerate_, isDualDegenerate_+numberRows_+numberColumns_, false);
     292  std::fill(isDualDegenerate_, isDualDegenerate_ + numberRows_ + numberColumns_, false);
    263293
    264294  // The dual degenerate variables are the nonbasic variables with a zero reduced costs
    265295  // An epsDegeneracy_ tolerance is used to detect zero reduced costs
    266   epsDegeneracy_ = 1.0e-04;//std::min(1.0e-2, 1.0e-7*fabs(model_->objectiveValue()));
     296  epsDegeneracy_ = 1.0e-04; //std::min(1.0e-2, 1.0e-7*fabs(model_->objectiveValue()));
    267297  double maxDegen = 0.0;
    268   for (int i = 0 ; i < numberColumns_+numberRows_ ; i++){
    269 
    270     if ( model_->getStatus(i) != ClpSimplex::basic && fabs(model_->reducedCost(i)) <= epsDegeneracy_  ) {
     298  for (int i = 0; i < numberColumns_ + numberRows_; i++) {
     299
     300    if (model_->getStatus(i) != ClpSimplex::basic && fabs(model_->reducedCost(i)) <= epsDegeneracy_) {
    271301      dualDegenerates_[coDualDegenerates_++] = i;
    272302      isDualDegenerate_[i] = true;
     
    280310  The input argument is a temporary array that is needed for the Clp's BTRAN
    281311  if which is NULL, every variable is tested */
    282 void ClpPESimplex::identifyCompatibleCols(int number, const int * which,
    283                                           CoinIndexedVector * spareRow2,
    284                                          CoinIndexedVector* wPrimal) {
     312void ClpPESimplex::identifyCompatibleCols(int number, const int *which,
     313  CoinIndexedVector *spareRow2,
     314  CoinIndexedVector *wPrimal)
     315{
    285316
    286317  // the update of variables within their bounds must have been called at least once
     
    289320  // initialize every variable to incompatible
    290321  coCompatibleCols_ = 0;
    291   std::fill(isCompatibleCol_, isCompatibleCol_+numberRows_+numberColumns_, false);
    292   std::fill(compatibilityCol_, compatibilityCol_+numberRows_+numberColumns_, -1.0);
    293 
    294   // treat the two obvious cases: 
    295   // - no primal degenerate variable => every column is compatible 
     322  std::fill(isCompatibleCol_, isCompatibleCol_ + numberRows_ + numberColumns_, false);
     323  std::fill(compatibilityCol_, compatibilityCol_ + numberRows_ + numberColumns_, -1.0);
     324
     325  // treat the two obvious cases:
     326  // - no primal degenerate variable => every column is compatible
    296327  // - only primal degenerate variables => no compatible column
    297   if (coPrimalDegenerates_ == 0)  {
     328  if (coPrimalDegenerates_ == 0) {
    298329    if (!which) {
    299       std::fill(isCompatibleCol_, isCompatibleCol_+numberRows_+numberColumns_, true);
    300       coCompatibleCols_ = numberRows_+numberColumns_;
    301     }
    302     else {
     330      std::fill(isCompatibleCol_, isCompatibleCol_ + numberRows_ + numberColumns_, true);
     331      coCompatibleCols_ = numberRows_ + numberColumns_;
     332    } else {
    303333      for (int j = 0; j < number; j++) {
    304334        isCompatibleCol_[which[j]] = true;
     
    307337    }
    308338    return;
    309   }
    310   else if (coPrimalDegenerates_ == numberRows_) {
     339  } else if (coPrimalDegenerates_ == numberRows_) {
    311340    return;
    312341  }
    313 
    314342
    315343  // fill the elements of w corresponding to the degenerate variables
     
    318346  //wPrimal_->clear();
    319347  wPrimal->checkClear();
    320   assert (coPrimalDegenerates_<=CoinMax(numberColumns_,numberRows_));
    321   for (int i = 0; i < coPrimalDegenerates_; i++)  {
     348  assert(coPrimalDegenerates_ <= CoinMax(numberColumns_, numberRows_));
     349  for (int i = 0; i < coPrimalDegenerates_; i++) {
    322350#if 0
    323351    double random;
     
    343371  // the zero elements (with a tolerance epsCompatibility_) are the compatibles
    344372  coCompatibleCols_ = 0;
    345   number = which ? number : numberColumns_+ numberRows_;
     373  number = which ? number : numberColumns_ + numberRows_;
    346374  int jColumn;
    347   assert (!wPrimal->packedMode());
     375  assert(!wPrimal->packedMode());
    348376  double *values = wPrimal->denseVector();
    349   const double * rowScale = model_->rowScale();
    350   CoinPackedMatrix* clpMatrix = model_->matrix();
    351   const int * row = clpMatrix->getIndices();
    352   const CoinBigIndex * columnStart = clpMatrix->getVectorStarts();
    353   const int * columnLength = clpMatrix->getVectorLengths();
    354   const double * elementByColumn = clpMatrix->getElements();
    355   for (int j = 0; j < number; j++)  {
     377  const double *rowScale = model_->rowScale();
     378  CoinPackedMatrix *clpMatrix = model_->matrix();
     379  const int *row = clpMatrix->getIndices();
     380  const CoinBigIndex *columnStart = clpMatrix->getVectorStarts();
     381  const int *columnLength = clpMatrix->getVectorLengths();
     382  const double *elementByColumn = clpMatrix->getElements();
     383  for (int j = 0; j < number; j++) {
    356384    if (which)
    357385      jColumn = which[j];
     
    359387      jColumn = j;
    360388
    361     if ( model_->getStatus(jColumn) == ClpSimplex::basic
    362        /* && !isPrimalDegenerate_[jColumn]*/) {
     389    if (model_->getStatus(jColumn) == ClpSimplex::basic
     390      /* && !isPrimalDegenerate_[jColumn]*/) {
    363391      isCompatibleCol_[jColumn] = false;
    364392      /*coCompatibleCols_++;*/
    365     }
    366     else {
    367       double dotProduct=0.0;
    368       if (jColumn<numberColumns_) {
    369         if (!rowScale) {
    370           for (CoinBigIndex i = columnStart[jColumn]; i < columnStart[jColumn] + columnLength[jColumn]; i++) {
    371               int iRow = row[i];
    372               dotProduct += values[iRow]*elementByColumn[i];
     393    } else {
     394      double dotProduct = 0.0;
     395      if (jColumn < numberColumns_) {
     396        if (!rowScale) {
     397          for (CoinBigIndex i = columnStart[jColumn]; i < columnStart[jColumn] + columnLength[jColumn]; i++) {
     398            int iRow = row[i];
     399            dotProduct += values[iRow] * elementByColumn[i];
    373400          }
    374         }
    375         else {
     401        } else {
    376402          // apply scaling
    377403          double scale = model_->columnScale()[jColumn];
    378404          for (CoinBigIndex i = columnStart[jColumn]; i < columnStart[jColumn] + columnLength[jColumn]; i++) {
    379             int iRow = row[i];
    380             dotProduct += values[iRow]*elementByColumn[i] * rowScale[iRow];
     405            int iRow = row[i];
     406            dotProduct += values[iRow] * elementByColumn[i] * rowScale[iRow];
    381407          }
    382           dotProduct *= scale;
    383         }
     408          dotProduct *= scale;
     409        }
    384410      } else {
    385         // slack
    386         dotProduct = values[jColumn-numberColumns_];
     411        // slack
     412        dotProduct = values[jColumn - numberColumns_];
    387413      }
    388414      dotProduct = fabs(dotProduct);
     
    392418      double dotProduct2 = fabs( PEdot(*tempColumn_, *wPrimal) );
    393419      assert (fabs(dotProduct-dotProduct2)<1.0e-6);
    394 #endif 
     420#endif
    395421      compatibilityCol_[jColumn] = dotProduct;
    396422
    397       if ( dotProduct < epsCompatibility_) {
     423      if (dotProduct < epsCompatibility_) {
    398424        isCompatibleCol_[jColumn] = true;
    399425        coCompatibleCols_++;
     
    405431
    406432/** Identify the dual compatible rows */
    407 void ClpPESimplex::identifyCompatibleRows(CoinIndexedVector * spare,
    408                                           CoinIndexedVector * wDual) {
     433void ClpPESimplex::identifyCompatibleRows(CoinIndexedVector *spare,
     434  CoinIndexedVector *wDual)
     435{
    409436
    410437  // the update of dual degenerate variables must have been called at least once
    411438  assert(dualDegenerates_);
    412439
    413   // only one case is obvious here: 
    414   // - no dual degenerate variable => every row is compatible 
    415   if (coDualDegenerates_ == 0)  {
    416     std::fill(isCompatibleRow_, isCompatibleRow_+numberRows_, false);
     440  // only one case is obvious here:
     441  // - no dual degenerate variable => every row is compatible
     442  if (coDualDegenerates_ == 0) {
     443    std::fill(isCompatibleRow_, isCompatibleRow_ + numberRows_, false);
    417444    coCompatibleRows_ = numberRows_;
    418445    return;
    419446  }
    420   assert (coDualDegenerates_<=CoinMax(numberColumns_,numberRows_));
     447  assert(coDualDegenerates_ <= CoinMax(numberColumns_, numberRows_));
    421448#if 0
    422449  // fill the elements of tempRandom with as many random elements as dual degenerate variables
     
    432459  wDual->checkClear();
    433460  //wDual_->clear();
    434   double timeTmp=0.0;
     461  double timeTmp = 0.0;
    435462  if (doStatistics_)
    436463    timeTmp = CoinCpuTime();
    437464  double *values = wDual->denseVector();
    438   const double * rowScale = model_->rowScale();
    439   CoinPackedMatrix* clpMatrix = model_->matrix();
    440   const int * row = clpMatrix->getIndices();
    441   const CoinBigIndex * columnStart = clpMatrix->getVectorStarts();
    442   const int * columnLength = clpMatrix->getVectorLengths();
    443   const double * elementByColumn = clpMatrix->getElements();
    444   for (int j = 0; j < coDualDegenerates_; j++)  {
     465  const double *rowScale = model_->rowScale();
     466  CoinPackedMatrix *clpMatrix = model_->matrix();
     467  const int *row = clpMatrix->getIndices();
     468  const CoinBigIndex *columnStart = clpMatrix->getVectorStarts();
     469  const int *columnLength = clpMatrix->getVectorLengths();
     470  const double *elementByColumn = clpMatrix->getElements();
     471  for (int j = 0; j < coDualDegenerates_; j++) {
    445472    // I try and save time to avoid calling unpack
    446473    int sequence = dualDegenerates_[j];
    447     if (sequence >= numberColumns_ ) {
    448         //slack
    449         values[sequence-numberColumns_] -= tempRandom_[j];
    450     }
    451     else {
    452           // column
     474    if (sequence >= numberColumns_) {
     475      //slack
     476      values[sequence - numberColumns_] -= tempRandom_[j];
     477    } else {
     478      // column
    453479      CoinBigIndex i;
    454480      if (!rowScale) {
    455481        for (i = columnStart[sequence]; i < columnStart[sequence] + columnLength[sequence]; i++) {
    456               int iRow = row[i];
    457               values[iRow] += tempRandom_[j]*elementByColumn[i];
    458           }
    459       }
    460       else {
    461           // apply scaling
    462           double scale = model_->columnScale()[sequence];
    463           for (i = columnStart[sequence]; i < columnStart[sequence] + columnLength[sequence]; i++) {
    464                int iRow = row[i];
    465               values[iRow] += tempRandom_[j]*elementByColumn[i]*scale * rowScale[iRow];
    466           }
     482          int iRow = row[i];
     483          values[iRow] += tempRandom_[j] * elementByColumn[i];
     484        }
     485      } else {
     486        // apply scaling
     487        double scale = model_->columnScale()[sequence];
     488        for (i = columnStart[sequence]; i < columnStart[sequence] + columnLength[sequence]; i++) {
     489          int iRow = row[i];
     490          values[iRow] += tempRandom_[j] * elementByColumn[i] * scale * rowScale[iRow];
     491        }
    467492      }
    468493    }
     
    480505  }
    481506  int n = 0;
    482   int * indices = wDual->getIndices();
    483   for (int i=0;i<numberRows_;i++) {
     507  int *indices = wDual->getIndices();
     508  for (int i = 0; i < numberRows_; i++) {
    484509    if (values[i])
    485       indices[n++]=i;
     510      indices[n++] = i;
    486511  }
    487512  wDual->setNumElements(n);
     
    504529#endif
    505530
    506   // the zero elements of wDual (with a tolerance epsCompatibility_) are the 
     531  // the zero elements of wDual (with a tolerance epsCompatibility_) are the
    507532  // dual-compatible variables
    508   assert (!wDual->packedMode());
     533  assert(!wDual->packedMode());
    509534  int size = wDual->getNumElements();
    510   std::fill(isCompatibleRow_, isCompatibleRow_+numberRows_, true);
     535  std::fill(isCompatibleRow_, isCompatibleRow_ + numberRows_, true);
    511536  coCompatibleRows_ = numberRows_;
    512   for (int i = 0; i < size; i++)  {
     537  for (int i = 0; i < size; i++) {
    513538    int iRow = indices[i];
    514539    double value = values[iRow];
    515     if ( fabs(value) >= epsCompatibility_*100.0) {
     540    if (fabs(value) >= epsCompatibility_ * 100.0) {
    516541      isCompatibleRow_[iRow] = false;
    517542      coCompatibleRows_--;
     
    525550
    526551*/
    527 void ClpPESimplex::updatePrimalDegeneratesAvg(int coPivots){
    528         int totalPivots = model_->numberIterations()+1;
    529         double fracPivots = static_cast<double> ( coPivots)/totalPivots;
    530         coPrimalDegeneratesAvg_ = floor( (1.0-fracPivots)*
    531                                          (coPrimalDegeneratesAvg_ + fracPivots*coPrimalDegenerates_));
     552void ClpPESimplex::updatePrimalDegeneratesAvg(int coPivots)
     553{
     554  int totalPivots = model_->numberIterations() + 1;
     555  double fracPivots = static_cast< double >(coPivots) / totalPivots;
     556  coPrimalDegeneratesAvg_ = floor((1.0 - fracPivots) * (coPrimalDegeneratesAvg_ + fracPivots * coPrimalDegenerates_));
    532557}
    533558
     
    536561
    537562*/
    538 void ClpPESimplex::updateDualDegeneratesAvg(int coPivots){
    539         int totalPivots = model_->numberIterations()+1;
    540         double fracPivots = static_cast<double> ( coPivots)/totalPivots;
    541         coDualDegeneratesAvg_ = floor((1.0-fracPivots)*coDualDegeneratesAvg_ + fracPivots*coDualDegenerates_);
     563void ClpPESimplex::updateDualDegeneratesAvg(int coPivots)
     564{
     565  int totalPivots = model_->numberIterations() + 1;
     566  double fracPivots = static_cast< double >(coPivots) / totalPivots;
     567  coDualDegeneratesAvg_ = floor((1.0 - fracPivots) * coDualDegeneratesAvg_ + fracPivots * coDualDegenerates_);
    542568}
    543569
     
    545571        The input is the number of pivots since last update.
    546572*/
    547 void ClpPESimplex::updateCompatibleColsAvg(int coPivots) {
    548         int totalPivots = model_->numberIterations()+1;
    549         double fracPivots = static_cast<double> ( coPivots)/totalPivots;
    550         coCompatibleColsAvg_ =  floor((1.0-fracPivots)* coCompatibleColsAvg_+fracPivots*(coCompatibleCols_/*-(numberRows_-coPrimalDegenerates_)*/));
     573void ClpPESimplex::updateCompatibleColsAvg(int coPivots)
     574{
     575  int totalPivots = model_->numberIterations() + 1;
     576  double fracPivots = static_cast< double >(coPivots) / totalPivots;
     577  coCompatibleColsAvg_ = floor((1.0 - fracPivots) * coCompatibleColsAvg_ + fracPivots * (coCompatibleCols_ /*-(numberRows_-coPrimalDegenerates_)*/));
    551578}
    552579
     
    554581        The input is the number of pivots since last update.
    555582*/
    556 void ClpPESimplex::updateCompatibleRowsAvg(int coPivots) {
    557         int totalPivots = model_->numberIterations()+1;
    558         double fracPivots = static_cast<double> ( coPivots)/totalPivots;
    559         coCompatibleRowsAvg_ =  floor((1.0-fracPivots)* coCompatibleRowsAvg_+fracPivots*coCompatibleRows_);
     583void ClpPESimplex::updateCompatibleRowsAvg(int coPivots)
     584{
     585  int totalPivots = model_->numberIterations() + 1;
     586  double fracPivots = static_cast< double >(coPivots) / totalPivots;
     587  coCompatibleRowsAvg_ = floor((1.0 - fracPivots) * coCompatibleRowsAvg_ + fracPivots * coCompatibleRows_);
    560588}
    561589
     
    563591
    564592#if PE_DEBUG >= 1
    565 /** Print the set of variables within their bounds */
    566 void ClpPESimplex::printPrimalDegenerates() {
     593/** Print the set of variables within their bounds */
     594void ClpPESimplex::printPrimalDegenerates()
     595{
    567596  assert(primalDegenerates_);
    568597
    569   if (! coPrimalDegenerates_) {
     598  if (!coPrimalDegenerates_) {
    570599    std::cout << "There is no primal degenerate variable!" << std::endl;
    571600    return;
     
    575604
    576605  std::cout << "List of primal degenerate variables:";
    577   for (int i = 0; i < coPrimalDegenerates_; i++)  {
     606  for (int i = 0; i < coPrimalDegenerates_; i++) {
    578607    std::cout << " " << pivotVariable[primalDegenerates_[i]] << " ;";
    579608  }
     
    581610}
    582611
    583 /** Print the set of primal compatible variables */
    584 void ClpPESimplex::printCompatibleCols() {
    585     assert(isCompatibleCol_);
    586 
    587   if (! coCompatibleCols_) {
     612/** Print the set of primal compatible variables */
     613void ClpPESimplex::printCompatibleCols()
     614{
     615  assert(isCompatibleCol_);
     616
     617  if (!coCompatibleCols_) {
    588618    std::cout << "There is no compatible column!" << std::endl;
    589619    return;
     
    591621
    592622  std::cout << "List of compatible columns:";
    593   for (int i = 0; i < numberColumns_+numberRows_; i++) {
     623  for (int i = 0; i < numberColumns_ + numberRows_; i++) {
    594624    if (isCompatibleCol_[i])
    595625      std::cout << " " << i << " ;";
     
    599629
    600630/** Check that a nonbasic variable is indeed contemptible */
    601 bool ClpPESimplex::checkCompatibilityCol(int sequence,
    602                                          CoinIndexedVector* spareRow2) {
     631bool ClpPESimplex::checkCompatibilityCol(int sequence,
     632  CoinIndexedVector *spareRow2)
     633{
    603634
    604635  bool isCompatible = true;
     
    615646#endif
    616647
    617   for (int j = 0; j < coPrimalDegenerates_; j++)  {
    618     if ( fabs((*tempColumn_)[primalDegenerates_[j]]) > epsDegeneracy_ ) {
    619         return false;
    620       }
     648  for (int j = 0; j < coPrimalDegenerates_; j++) {
     649    if (fabs((*tempColumn_)[primalDegenerates_[j]]) > epsDegeneracy_) {
     650      return false;
     651    }
    621652  }
    622653  delete tempColumn_;
     
    625656#endif
    626657/** Check that a basic row is indeed compatible */
    627 bool ClpPESimplex::checkCompatibilityRow(int pivotRow) {
     658bool ClpPESimplex::checkCompatibilityRow(int pivotRow)
     659{
    628660
    629661  bool isCompatible = true;
     
    639671  model_->clpMatrix()->transposeTimes(model_, -1.0, model_->rowArray(0), model_->rowArray(1), model_->columnArray(0));
    640672
    641   CoinIndexedVector* columnArray = model_->columnArray(0);
    642   CoinIndexedVector* rowArray = model_->rowArray(0);
     673  CoinIndexedVector *columnArray = model_->columnArray(0);
     674  CoinIndexedVector *rowArray = model_->rowArray(0);
    643675  int nzCol = columnArray->getNumElements();
    644   int* indCol =columnArray->getIndices();
    645   double* valCol = columnArray->denseVector();
     676  int *indCol = columnArray->getIndices();
     677  double *valCol = columnArray->denseVector();
    646678  int nzRow = rowArray->getNumElements();
    647   int* indRow = rowArray->getIndices();
    648   double* valRow = rowArray->denseVector();
    649 
    650   if ( columnArray->packedMode() )  {
    651     for (int j = 0; j < nzCol; j++)  {
    652         int iCol = indCol[j];
    653         if ( isDualDegenerate_[iCol] && fabs(valCol[j]) > epsDegeneracy_ ) {
    654           std::cout << "Dual degenerate column: " << valCol[j] <<std::endl;
    655         }
    656     }
    657   }
    658   else {
    659     for (int j = 0; j < nzCol; j++)  {
    660         int iCol = indCol[j];
    661         if ( isDualDegenerate_[iCol] && fabs(valCol[iCol]) > epsDegeneracy_ ) {
    662           std::cout << "Dual degenerate column: " << valCol[iCol] <<std::endl;
    663         }
    664     }
    665   }
    666     if ( rowArray->packedMode() )  {
    667         for (int j = 0; j < nzRow; j++)  {
    668             int iRow = indRow[j];
    669             if ( isDualDegenerate_[iRow+numberColumns_] && fabs(valRow[j]) > epsDegeneracy_ ) {
    670                 std::cout << "Dual degenerate row: " << valRow[j] <<std::endl;
    671             }
    672         }
    673     }
    674     else {
    675         for (int j = 0; j < nzRow; j++)  {
    676             int iRow = indRow[j];
    677             if ( isDualDegenerate_[iRow+numberColumns_] && fabs(valRow[iRow]) > epsDegeneracy_ ) {
    678                 std::cout << "Dual degenerate row: " << valRow[iRow] <<std::endl;
    679             }
    680         }
    681     }
     679  int *indRow = rowArray->getIndices();
     680  double *valRow = rowArray->denseVector();
     681
     682  if (columnArray->packedMode()) {
     683    for (int j = 0; j < nzCol; j++) {
     684      int iCol = indCol[j];
     685      if (isDualDegenerate_[iCol] && fabs(valCol[j]) > epsDegeneracy_) {
     686        std::cout << "Dual degenerate column: " << valCol[j] << std::endl;
     687      }
     688    }
     689  } else {
     690    for (int j = 0; j < nzCol; j++) {
     691      int iCol = indCol[j];
     692      if (isDualDegenerate_[iCol] && fabs(valCol[iCol]) > epsDegeneracy_) {
     693        std::cout << "Dual degenerate column: " << valCol[iCol] << std::endl;
     694      }
     695    }
     696  }
     697  if (rowArray->packedMode()) {
     698    for (int j = 0; j < nzRow; j++) {
     699      int iRow = indRow[j];
     700      if (isDualDegenerate_[iRow + numberColumns_] && fabs(valRow[j]) > epsDegeneracy_) {
     701        std::cout << "Dual degenerate row: " << valRow[j] << std::endl;
     702      }
     703    }
     704  } else {
     705    for (int j = 0; j < nzRow; j++) {
     706      int iRow = indRow[j];
     707      if (isDualDegenerate_[iRow + numberColumns_] && fabs(valRow[iRow]) > epsDegeneracy_) {
     708        std::cout << "Dual degenerate row: " << valRow[iRow] << std::endl;
     709      }
     710    }
     711  }
    682712
    683713  return isCompatible;
    684714}
    685715// checks size
    686 bool
    687 ClpPESimplex::checkSize()
    688 {
    689   return (numberRows_==model_->numberRows()&&
    690           numberColumns_==model_->numberColumns());
     716bool ClpPESimplex::checkSize()
     717{
     718  return (numberRows_ == model_->numberRows() && numberColumns_ == model_->numberColumns());
    691719}
    692720/* Update the dual compatible rows */
    693 void
    694 ClpPESimplex::updateCompatibleRows(int iColumn)
    695 {
    696   if (iColumn<numberColumns_) {
     721void ClpPESimplex::updateCompatibleRows(int iColumn)
     722{
     723  if (iColumn < numberColumns_) {
    697724    // get the packed matrix
    698     CoinPackedMatrix* clpMatrix = model_->matrix();
     725    CoinPackedMatrix *clpMatrix = model_->matrix();
    699726
    700727    // get the matrix data pointers
    701     const int *  row = clpMatrix->getIndices();
    702     const CoinBigIndex *  columnStart = clpMatrix->getVectorStarts();
    703     const int * columnLength = clpMatrix->getVectorLengths();
    704     CoinBigIndex start=columnStart[iColumn];
    705     CoinBigIndex next=start+columnLength[iColumn];
    706     for (CoinBigIndex j=start;j<next;j++) {
    707       int jRow=row[j];
     728    const int *row = clpMatrix->getIndices();
     729    const CoinBigIndex *columnStart = clpMatrix->getVectorStarts();
     730    const int *columnLength = clpMatrix->getVectorLengths();
     731    CoinBigIndex start = columnStart[iColumn];
     732    CoinBigIndex next = start + columnLength[iColumn];
     733    for (CoinBigIndex j = start; j < next; j++) {
     734      int jRow = row[j];
    708735      if (isCompatibleRow_[jRow]) {
    709         isCompatibleRow_[jRow]=false;
    710         coCompatibleRows_--;
     736        isCompatibleRow_[jRow] = false;
     737        coCompatibleRows_--;
    711738      }
    712739    }
    713740  } else {
    714     int jRow=iColumn-numberColumns_;
     741    int jRow = iColumn - numberColumns_;
    715742    if (isCompatibleRow_[jRow]) {
    716       isCompatibleRow_[jRow]=false;
     743      isCompatibleRow_[jRow] = false;
    717744      coCompatibleRows_--;
    718745    }
     
    721748
    722749/** Print the total recorded time */
    723 void ClpPESimplex::printTimer(std::ostream &out) {
    724   out << "Cpu in compatibility: " << timeCompatibility_ << " s"<< std::endl;
    725 }
    726 
    727 
    728 
    729 
     750void ClpPESimplex::printTimer(std::ostream &out)
     751{
     752  out << "Cpu in compatibility: " << timeCompatibility_ << " s" << std::endl;
     753}
     754
     755/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     756*/
Note: See TracChangeset for help on using the changeset viewer.