source: trunk/Clp/src/ClpPESimplex.cpp @ 2449

Last change on this file since 2449 was 2385, checked in by unxusr, 11 months ago

formatting

File size: 25.4 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3/*
4   Authors
5
6   Jeremy Omer
7
8   Last update: april 10, 2015
9
10 */
11
12#include <iostream>
13#include "ClpPESimplex.hpp"
14#include "ClpFactorization.hpp"
15#include "ClpMessage.hpp"
16
17/** SHARED METHODS FOR USEFUL ALGEBRAIC OPERATIONS */
18
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}
30
31/** inner product between two coin vectors
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}
102
103/** BASE CLASS FOR THE IMPROVED SIMPLEX
104*/
105
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{
141
142  // the improved simplex object should only be implemented after loading the model
143  assert(model_->numberColumns() > 0);
144
145  // size of the original model
146  numberColumns_ = model_->numberColumns();
147  numberRows_ = model_->numberRows();
148
149  // allocate memory to the arrays
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)));
165  // fill
166  CoinThreadRandom generator = *model_->randomNumberGenerator();
167  for (int i = 0; i < CoinMax(numberColumns_, numberRows_); i++) {
168    double random;
169    do
170      random = static_cast< int >(generator.randomDouble() * 1.0e6) - 5.0e5;
171    while (random == 0.0);
172    //random = static_cast<int>(model_->randomNumberGenerator()->randomDouble()*1.0e7)-5.0e6;
173    tempRandom_[i] = random;
174  }
175  if (model_->logLevel() > 2)
176    doStatistics_ = model_->logLevel();
177}
178
179/** Destructor */
180ClpPESimplex::~ClpPESimplex()
181{
182
183  // free memory
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()) {
204    char generalPrint[200];
205    sprintf(generalPrint, "Degenerate pivots   : %d, compatibility time %.2f",
206      coDegeneratePivots(),
207      timeCompatibility());
208    model_->messageHandler()->message(CLP_GENERAL,
209      *model_->messagesPointer())
210      << generalPrint << CoinMessageEol;
211
212#if 1
213    int numberPivots = model_->numberIterations();
214    if (coDualDegeneratesAvg()) {
215      sprintf(generalPrint, "coDegenAvg/rows %g coCompatAvg/rows %g",
216        coDualDegeneratesAvg() / numberRows_,
217        coCompatibleRowsAvg() / numberRows_);
218      model_->messageHandler()->message(CLP_GENERAL,
219        *model_->messagesPointer())
220        << generalPrint << CoinMessageEol;
221    } else if (coPrimalDegeneratesAvg()) {
222      sprintf(generalPrint, "coDegenAvg/columns %g coCompatAvg/columns %g",
223        coPrimalDegeneratesAvg() / numberColumns_,
224        coCompatibleColsAvg() / numberColumns_);
225      model_->messageHandler()->message(CLP_GENERAL,
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())));
231      model_->messageHandler()->message(CLP_GENERAL,
232        *model_->messagesPointer())
233        << generalPrint << CoinMessageEol;
234    }
235    if (coCompatiblePivots()) {
236      sprintf(generalPrint, "coDegenerateCompatiblePivots()/coCompatiblePivots() %g", static_cast< double >(coDegenerateCompatiblePivots()) / static_cast< double >(coCompatiblePivots()));
237      model_->messageHandler()->message(CLP_GENERAL,
238        *model_->messagesPointer())
239        << generalPrint << CoinMessageEol;
240    }
241    sprintf(generalPrint, "coDegeneratePivots()/ numberPivots %g", static_cast< double >(coDegeneratePivots()) / static_cast< double >(numberPivots));
242    model_->messageHandler()->message(CLP_GENERAL,
243      *model_->messagesPointer())
244      << generalPrint << CoinMessageEol;
245    sprintf(generalPrint, "coCompatiblePivots() %d coPriorityPivots() %d",
246      coCompatiblePivots(), coPriorityPivots());
247    model_->messageHandler()->message(CLP_GENERAL,
248      *model_->messagesPointer())
249      << generalPrint << CoinMessageEol;
250    //sprintf(generalPrint, numberPivots);
251    // sprintf(generalPrint, timeMultRandom());
252    //sprintf(generalPrint, timeLinearSystem());
253#endif
254  }
255}
256
257/** Updates the set of variables that are not at their bounds */
258void ClpPESimplex::updatePrimalDegenerates()
259{
260
261  coPrimalDegenerates_ = 0;
262  epsDegeneracy_ = 1.0e-04; //std::min(1.0e-2, 1.0e-7*fabs(model_->objectiveValue()));
263  std::fill(isPrimalDegenerate_, isPrimalDegenerate_ + numberRows_ + numberColumns_, false);
264  int *pivotVariable = model_->pivotVariable();
265
266  // Only go over the basic variables, since non basic varaibles are
267  // set to one of their bounds (or zero if the variable is free)
268  // An epsDegeneracy_ tolerance is used to detect variables at a bound
269  for (int i = 0; i < numberRows_; i++) {
270    int iVar = pivotVariable[i]; // index of basic variable at row i
271
272    // std::cout << "solution[" << iVar << "] = " << model_->solution(iVar) << " ; lb = " << model_->lower(iVar) << " ; ub = "<< model_->upper(iVar) << "\n" ;
273
274    double dVal = model_->solution(iVar);
275    double dUb = model_->upper(iVar);
276    double dLb = model_->lower(iVar);
277    // if (fabs(dVal) <= 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_)) {
280      primalDegenerates_[coPrimalDegenerates_++] = i;
281      isPrimalDegenerate_[iVar] = true;
282    }
283  }
284  coUpdateDegenerates_++;
285}
286
287/** Updates the set of dual degenerate variables */
288void ClpPESimplex::updateDualDegenerates()
289{
290
291  coDualDegenerates_ = 0;
292  std::fill(isDualDegenerate_, isDualDegenerate_ + numberRows_ + numberColumns_, false);
293
294  // The dual degenerate variables are the nonbasic variables with a zero reduced costs
295  // An epsDegeneracy_ tolerance is used to detect zero reduced costs
296  epsDegeneracy_ = 1.0e-04; //std::min(1.0e-2, 1.0e-7*fabs(model_->objectiveValue()));
297  double maxDegen = 0.0;
298  for (int i = 0; i < numberColumns_ + numberRows_; i++) {
299
300    if (model_->getStatus(i) != ClpSimplex::basic && fabs(model_->reducedCost(i)) <= epsDegeneracy_) {
301      dualDegenerates_[coDualDegenerates_++] = i;
302      isDualDegenerate_[i] = true;
303      maxDegen = std::max(maxDegen, fabs(model_->reducedCost(i)));
304    }
305  }
306  coUpdateDegenerates_++;
307}
308
309/** Identify the primal compatible columns in a subset of columns
310  The input argument is a temporary array that is needed for the Clp's BTRAN
311  if which is NULL, every variable is tested */
312void ClpPESimplex::identifyCompatibleCols(int number, const int *which,
313  CoinIndexedVector *spareRow2,
314  CoinIndexedVector *wPrimal)
315{
316
317  // the update of variables within their bounds must have been called at least once
318  assert(primalDegenerates_);
319
320  // initialize every variable to incompatible
321  coCompatibleCols_ = 0;
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
327  // - only primal degenerate variables => no compatible column
328  if (coPrimalDegenerates_ == 0) {
329    if (!which) {
330      std::fill(isCompatibleCol_, isCompatibleCol_ + numberRows_ + numberColumns_, true);
331      coCompatibleCols_ = numberRows_ + numberColumns_;
332    } else {
333      for (int j = 0; j < number; j++) {
334        isCompatibleCol_[which[j]] = true;
335      }
336      coCompatibleCols_ = number;
337    }
338    return;
339  } else if (coPrimalDegenerates_ == numberRows_) {
340    return;
341  }
342
343  // fill the elements of w corresponding to the degenerate variables
344  // with random values and leave the other to zero
345  // no longer using wPrimal_
346  //wPrimal_->clear();
347  wPrimal->checkClear();
348  assert(coPrimalDegenerates_ <= CoinMax(numberColumns_, numberRows_));
349  for (int i = 0; i < coPrimalDegenerates_; i++) {
350#if 0
351    double random;
352      do
353        random = static_cast<double> ((rand() %(static_cast<int> (1.0e6)+1))) -5.0e5;
354      while (random == 0.0);
355    wPrimal->quickInsert(primalDegenerates_[i], 0.5+random);
356#else
357    wPrimal->quickInsert(primalDegenerates_[i], tempRandom_[i]);
358#endif
359  }
360
361  // compute wTran * Binv and store it into wTran
362#ifdef PE_TEST
363  model_->factorization()->doStatistics(false);
364#endif
365  model_->factorization()->updateColumnTranspose(spareRow2, wPrimal);
366#ifdef PE_TEST
367  model_->factorization()->doStatistics(true);
368#endif
369
370  // compute wTran * AN, where AN is the matrix of nonbasic variables
371  // the zero elements (with a tolerance epsCompatibility_) are the compatibles
372  coCompatibleCols_ = 0;
373  number = which ? number : numberColumns_ + numberRows_;
374  int jColumn;
375  assert(!wPrimal->packedMode());
376  double *values = wPrimal->denseVector();
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++) {
384    if (which)
385      jColumn = which[j];
386    else
387      jColumn = j;
388
389    if (model_->getStatus(jColumn) == ClpSimplex::basic
390      /* && !isPrimalDegenerate_[jColumn]*/) {
391      isCompatibleCol_[jColumn] = false;
392      /*coCompatibleCols_++;*/
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];
400          }
401        } else {
402          // apply scaling
403          double scale = model_->columnScale()[jColumn];
404          for (CoinBigIndex i = columnStart[jColumn]; i < columnStart[jColumn] + columnLength[jColumn]; i++) {
405            int iRow = row[i];
406            dotProduct += values[iRow] * elementByColumn[i] * rowScale[iRow];
407          }
408          dotProduct *= scale;
409        }
410      } else {
411        // slack
412        dotProduct = values[jColumn - numberColumns_];
413      }
414      dotProduct = fabs(dotProduct);
415#if 0
416      model_->unpack(tempColumn_, jColumn);
417      // perform the inner product <tempColumn_,wPrimal_>
418      double dotProduct2 = fabs( PEdot(*tempColumn_, *wPrimal) );
419      assert (fabs(dotProduct-dotProduct2)<1.0e-6);
420#endif
421      compatibilityCol_[jColumn] = dotProduct;
422
423      if (dotProduct < epsCompatibility_) {
424        isCompatibleCol_[jColumn] = true;
425        coCompatibleCols_++;
426      }
427    }
428  }
429  wPrimal->clear();
430}
431
432/** Identify the dual compatible rows */
433void ClpPESimplex::identifyCompatibleRows(CoinIndexedVector *spare,
434  CoinIndexedVector *wDual)
435{
436
437  // the update of dual degenerate variables must have been called at least once
438  assert(dualDegenerates_);
439
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);
444    coCompatibleRows_ = numberRows_;
445    return;
446  }
447  assert(coDualDegenerates_ <= CoinMax(numberColumns_, numberRows_));
448#if 0
449  // fill the elements of tempRandom with as many random elements as dual degenerate variables
450  for (int j = 0; j < coDualDegenerates_; j++) {
451    //do
452    tempRandom_[j] = static_cast<double> ((rand() %(static_cast<int> (1.0e7)+1))) -5.0e6;
453    //while (tempRandom_[j] == 0.0);
454  }
455#endif
456
457  // compute the product A_D * tempRandom and store it into wDual
458  // No longer using wDual_
459  wDual->checkClear();
460  //wDual_->clear();
461  double timeTmp = 0.0;
462  if (doStatistics_)
463    timeTmp = CoinCpuTime();
464  double *values = wDual->denseVector();
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++) {
472    // I try and save time to avoid calling unpack
473    int sequence = dualDegenerates_[j];
474    if (sequence >= numberColumns_) {
475      //slack
476      values[sequence - numberColumns_] -= tempRandom_[j];
477    } else {
478      // column
479      CoinBigIndex i;
480      if (!rowScale) {
481        for (i = columnStart[sequence]; i < columnStart[sequence] + columnLength[sequence]; i++) {
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        }
492      }
493    }
494    //       matrix_->unpack(this, rowArray, sequenceIn_);
495    //  }
496    // model_->unpack(tempColumn_, dualDegenerates_[j]);
497
498    // int nz = tempColumn_->getNumElements();
499    // int *ind = tempColumn_->getIndices();
500    // double *val = tempColumn_->denseVector();
501
502    // for (int k = 0; k < nz ; k++) {
503    //   values[ind[k]] = values[ind[k]] + tempRandom_[j]*val[ind[k]];
504    // }
505  }
506  int n = 0;
507  int *indices = wDual->getIndices();
508  for (int i = 0; i < numberRows_; i++) {
509    if (values[i])
510      indices[n++] = i;
511  }
512  wDual->setNumElements(n);
513  wDual->setPackedMode(false);
514
515  // compute Binv*wDual and store it into wDual
516#ifdef PE_TEST
517  if (doStatistics_) {
518    timeMultRandom_ += CoinCpuTime() - timeTmp;
519    timeTmp = CoinCpuTime();
520  }
521  model_->factorization()->doStatistics(false);
522#endif
523  model_->factorization()->updateColumn(spare, wDual);
524#ifdef PE_TEST
525  model_->factorization()->doStatistics(true);
526  if (doStatistics_) {
527    timeLinearSystem_ += CoinCpuTime() - timeTmp;
528  }
529#endif
530
531  // the zero elements of wDual (with a tolerance epsCompatibility_) are the
532  // dual-compatible variables
533  assert(!wDual->packedMode());
534  int size = wDual->getNumElements();
535  std::fill(isCompatibleRow_, isCompatibleRow_ + numberRows_, true);
536  coCompatibleRows_ = numberRows_;
537  for (int i = 0; i < size; i++) {
538    int iRow = indices[i];
539    double value = values[iRow];
540    if (fabs(value) >= epsCompatibility_ * 100.0) {
541      isCompatibleRow_[iRow] = false;
542      coCompatibleRows_--;
543    }
544  }
545  wDual->clear();
546}
547
548/* Update the average number of primal degenerate columns 
549        The input is the number of pivots since last update.
550
551*/
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_));
557}
558
559/* Update the average number of dual degenerate columns
560        The input is the number of pivots since last update.
561
562*/
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_);
568}
569
570/* Update the average number of compatible columns.
571        The input is the number of pivots since last update.
572*/
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_)*/));
578}
579
580/* Update the average number of compatible rows.
581        The input is the number of pivots since last update.
582*/
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_);
588}
589
590/** DEBUG METHODS */
591
592#if PE_DEBUG >= 1
593/** Print the set of variables within their bounds */
594void ClpPESimplex::printPrimalDegenerates()
595{
596  assert(primalDegenerates_);
597
598  if (!coPrimalDegenerates_) {
599    std::cout << "There is no primal degenerate variable!" << std::endl;
600    return;
601  }
602
603  int *pivotVariable = model_->pivotVariable();
604
605  std::cout << "List of primal degenerate variables:";
606  for (int i = 0; i < coPrimalDegenerates_; i++) {
607    std::cout << " " << pivotVariable[primalDegenerates_[i]] << " ;";
608  }
609  std::cout << std::endl;
610}
611
612/** Print the set of primal compatible variables */
613void ClpPESimplex::printCompatibleCols()
614{
615  assert(isCompatibleCol_);
616
617  if (!coCompatibleCols_) {
618    std::cout << "There is no compatible column!" << std::endl;
619    return;
620  }
621
622  std::cout << "List of compatible columns:";
623  for (int i = 0; i < numberColumns_ + numberRows_; i++) {
624    if (isCompatibleCol_[i])
625      std::cout << " " << i << " ;";
626  }
627  std::cout << std::endl;
628}
629
630/** Check that a nonbasic variable is indeed contemptible */
631bool ClpPESimplex::checkCompatibilityCol(int sequence,
632  CoinIndexedVector *spareRow2)
633{
634
635  bool isCompatible = true;
636  tempColumn_ = new CoinIndexedVector();
637  tempColumn_->reserve(numberRows_);
638
639  model_->unpack(tempColumn_, sequence);
640#ifdef PE_TEST
641  model_->factorization()->doStatistics(false);
642#endif
643  model_->factorization()->updateColumn(spareRow2, tempColumn_);
644#ifdef PE_TEST
645  model_->factorization()->doStatistics(true);
646#endif
647
648  for (int j = 0; j < coPrimalDegenerates_; j++) {
649    if (fabs((*tempColumn_)[primalDegenerates_[j]]) > epsDegeneracy_) {
650      return false;
651    }
652  }
653  delete tempColumn_;
654  return isCompatible;
655}
656#endif
657/** Check that a basic row is indeed compatible */
658bool ClpPESimplex::checkCompatibilityRow(int pivotRow)
659{
660
661  bool isCompatible = true;
662  double direction = 1;
663  model_->rowArray(0)->createPacked(1, &pivotRow, &direction);
664#ifdef PE_TEST
665  model_->factorization()->doStatistics(false);
666#endif
667  model_->factorization()->updateColumnTranspose(model_->rowArray(1), model_->rowArray(0));
668#ifdef PE_TEST
669  model_->factorization()->doStatistics(true);
670#endif
671  model_->clpMatrix()->transposeTimes(model_, -1.0, model_->rowArray(0), model_->rowArray(1), model_->columnArray(0));
672
673  CoinIndexedVector *columnArray = model_->columnArray(0);
674  CoinIndexedVector *rowArray = model_->rowArray(0);
675  int nzCol = columnArray->getNumElements();
676  int *indCol = columnArray->getIndices();
677  double *valCol = columnArray->denseVector();
678  int nzRow = rowArray->getNumElements();
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  }
712
713  return isCompatible;
714}
715// checks size
716bool ClpPESimplex::checkSize()
717{
718  return (numberRows_ == model_->numberRows() && numberColumns_ == model_->numberColumns());
719}
720/* Update the dual compatible rows */
721void ClpPESimplex::updateCompatibleRows(int iColumn)
722{
723  if (iColumn < numberColumns_) {
724    // get the packed matrix
725    CoinPackedMatrix *clpMatrix = model_->matrix();
726
727    // get the matrix data pointers
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];
735      if (isCompatibleRow_[jRow]) {
736        isCompatibleRow_[jRow] = false;
737        coCompatibleRows_--;
738      }
739    }
740  } else {
741    int jRow = iColumn - numberColumns_;
742    if (isCompatibleRow_[jRow]) {
743      isCompatibleRow_[jRow] = false;
744      coCompatibleRows_--;
745    }
746  }
747}
748
749/** Print the total recorded time */
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 TracBrowser for help on using the repository browser.