Changeset 137 for trunk


Ignore:
Timestamp:
Mar 17, 2003 10:37:23 AM (17 years ago)
Author:
forrest
Message:

faster Clp and fix memory leaks

Location:
trunk
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpDualRowSteepest.cpp

    r97 r137  
    140140  double * upper = model_->upperRegion();
    141141  // do last pivot row one here
     142  //#define COLUMN_BIAS 4.0
     143  //#define FIXED_BIAS 10.0
    142144  if (lastPivotRow>=0) {
     145#ifdef COLUMN_BIAS
     146    int numberColumns = model_->numberColumns();
     147#endif
    143148    int iPivot=pivotVariable[lastPivotRow];
    144149    double value = solution[iPivot];
     
    146151    double upper = model_->upper(iPivot);
    147152    if (value>upper+tolerance) {
     153      value -= upper;
     154      value *= value;
     155#ifdef COLUMN_BIAS
     156      if (iPivot<numberColumns)
     157        value *= COLUMN_BIAS; // bias towards columns
     158#endif
    148159      // store square in list
    149160      if (infeas[lastPivotRow])
    150         infeas[lastPivotRow] = (value-upper)*(value-upper); // already there
     161        infeas[lastPivotRow] = value; // already there
    151162      else
    152         infeasible_->quickAdd(lastPivotRow,(value-upper)*(value-upper));
     163        infeasible_->quickAdd(lastPivotRow,value);
    153164    } else if (value<lower-tolerance) {
     165      value -= lower;
     166      value *= value;
     167#ifdef COLUMN_BIAS
     168      if (iPivot<numberColumns)
     169        value *= COLUMN_BIAS; // bias towards columns
     170#endif
    154171      // store square in list
    155172      if (infeas[lastPivotRow])
    156         infeas[lastPivotRow] = (value-lower)*(value-lower); // already there
     173        infeas[lastPivotRow] = value; // already there
    157174      else
    158         infeasible_->add(lastPivotRow,(value-lower)*(value-lower));
     175        infeasible_->add(lastPivotRow,value);
    159176    } else {
    160177      // feasible - was it infeasible - if so set tiny
     
    325342  int pivotRow = model_->pivotRow();
    326343  double * solution = model_->solutionRegion();
     344#ifdef COLUMN_BIAS
     345  int numberColumns = model_->numberColumns();
     346#endif
    327347  for (i=0;i<number;i++) {
    328348    int iRow=which[i];
     
    345365      value = model_->valueIncomingDual();
    346366    }
    347     if (value>upper+tolerance) {
     367    if (value<lower-tolerance) {
     368      value -= lower;
     369      value *= value;
     370#ifdef COLUMN_BIAS
     371      if (iPivot<numberColumns)
     372        value *= COLUMN_BIAS; // bias towards columns
     373#endif
     374#ifdef FIXED_BIAS
     375      if (lower==upper)
     376        value *= FIXED_BIAS; // bias towards taking out fixed variables
     377#endif
    348378      // store square in list
    349379      if (infeas[iRow])
    350         infeas[iRow] = (value-upper)*(value-upper); // already there
     380        infeas[iRow] = value; // already there
    351381      else
    352         infeasible_->quickAdd(iRow,(value-upper)*(value-upper));
    353     } else if (value<lower-tolerance) {
     382        infeasible_->quickAdd(iRow,value);
     383    } else if (value>upper+tolerance) {
     384      value -= upper;
     385      value *= value;
     386#ifdef COLUMN_BIAS
     387      if (iPivot<numberColumns)
     388        value *= COLUMN_BIAS; // bias towards columns
     389#endif
     390#ifdef FIXED_BIAS
     391      if (lower==upper)
     392        value *= FIXED_BIAS; // bias towards taking out fixed variables
     393#endif
    354394      // store square in list
    355395      if (infeas[iRow])
    356         infeas[iRow] = (value-lower)*(value-lower); // already there
     396        infeas[iRow] = value; // already there
    357397      else
    358         infeasible_->quickAdd(iRow,(value-lower)*(value-lower));
     398        infeasible_->quickAdd(iRow,value);
    359399    } else {
    360400      // feasible - was it infeasible - if so set tiny
     
    483523    const int * pivotVariable = model_->pivotVariable();
    484524    double tolerance=model_->currentPrimalTolerance();
    485     for (iRow=0;iRow<model_->numberRows();iRow++) {
     525    for (iRow=0;iRow<numberRows;iRow++) {
    486526      int iPivot=pivotVariable[iRow];
    487527      double value = model_->solution(iPivot);
    488528      double lower = model_->lower(iPivot);
    489529      double upper = model_->upper(iPivot);
    490       if (value>upper+tolerance) {
     530      if (value<lower-tolerance) {
     531        value -= lower;
     532        value *= value;
     533#ifdef COLUMN_BIAS
     534        if (iPivot<numberColumns)
     535          value *= COLUMN_BIAS; // bias towards columns
     536#endif
     537#ifdef FIXED_BIAS
     538        if (lower==upper)
     539          value *= FIXED_BIAS; // bias towards taking out fixed variables
     540#endif
     541        // store square in list
     542        infeasible_->quickAdd(iRow,value);
     543      } else if (value>upper+tolerance) {
    491544        value -= upper;
     545        value *= value;
     546#ifdef COLUMN_BIAS
     547        if (iPivot<numberColumns)
     548          value *= COLUMN_BIAS; // bias towards columns
     549#endif
     550#ifdef FIXED_BIAS
     551        if (lower==upper)
     552          value *= FIXED_BIAS; // bias towards taking out fixed variables
     553#endif
    492554        // store square in list
    493         infeasible_->quickAdd(iRow,value*value);
    494       } else if (value<lower-tolerance) {
    495         value -= lower;
    496         // store square in list
    497         infeasible_->quickAdd(iRow,value*value);
     555        infeasible_->quickAdd(iRow,value);
    498556      }
    499557    }
  • trunk/ClpMessage.cpp

    r102 r137  
    7272  {CLP_EMPTY_PROBLEM,3002,0,"Not solving empty problem - %d rows, %d columns and %d elements"},
    7373  {CLP_CRASH,28,1,"Crash put %d variables in basis, %d dual infeasibilities"},
     74  {CLP_END_VALUES_PASS,29,1,"End of values pass after %d iterations"},
    7475  {CLP_DUMMY_END,999999,0,""}
    7576};
  • trunk/ClpNonLinearCost.cpp

    r123 r137  
    342342    lowerValue = lower_[iRange];
    343343    upperValue = lower_[iRange+1];
    344     switch(model_->getStatus(iSequence)) {
     344    ClpSimplex::Status status = model_->getStatus(iSequence);
     345    if (upperValue==lowerValue) {
     346      if (status != ClpSimplex::basic)
     347        model_->setStatus(iSequence,ClpSimplex::isFixed);
     348    }
     349    switch(status) {
    345350     
    346351    case ClpSimplex::basic:
     
    415420      }
    416421      break;
     422    case ClpSimplex::isFixed:
     423      break;
    417424    }
    418425  }
     
    617624  lower = lower_[iRange];
    618625  upper = lower_[iRange+1];
    619   switch(model_->getStatus(iPivot)) {
     626  ClpSimplex::Status status = model_->getStatus(iPivot);
     627  if (upper==lower) {
     628    if (status != ClpSimplex::basic) {
     629      model_->setStatus(iPivot,ClpSimplex::isFixed);
     630      status = ClpSimplex::basic; // so will skip
     631    }
     632  }
     633  switch(status) {
    620634     
    621635  case ClpSimplex::basic:
     
    625639  case ClpSimplex::atUpperBound:
    626640  case ClpSimplex::atLowerBound:
     641  case ClpSimplex::isFixed:
    627642    // set correctly
    628643    if (fabs(value-lower)<=primalTolerance*1.001){
     
    636651    break;
    637652  }
    638   if (upper-lower<1.0e-8)
    639     model_->setFixed(iPivot);
    640   else
    641     model_->clearFixed(iPivot);
    642653  double difference = cost-cost_[iRange];
    643654  cost = cost_[iRange];
  • trunk/ClpPackedMatrix.cpp

    r119 r137  
    222222  ClpPackedMatrix* rowCopy =
    223223    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
    224   if (numberInRowArray>0.3*numberRows||!rowCopy) {
     224  if (numberInRowArray>0.333*numberRows||!rowCopy) {
    225225    // do by column
    226226    int iColumn;
     
    230230    const int * columnLength = matrix_->getVectorLengths();
    231231    const double * elementByColumn = matrix_->getElements();
    232     double * markVector = y->denseVector(); // probably empty
    233232    const double * rowScale = model->rowScale();
    234233    int numberColumns = model->numberColumns();
    235     if (!rowScale) {
    236       for (iColumn=0;iColumn<numberColumns;iColumn++) {
    237         double value = markVector[iColumn];
    238         markVector[iColumn]=0.0;
    239         double value2 = 0.0;
    240         CoinBigIndex j;
    241         for (j=columnStart[iColumn];
    242              j<columnStart[iColumn]+columnLength[iColumn];j++) {
    243           int iRow = row[j];
    244           value2 += pi[iRow]*elementByColumn[j];
    245         }
    246         value += value2*scalar;
    247         if (fabs(value)>zeroTolerance) {
    248           index[numberNonZero++]=iColumn;
    249           array[iColumn]=value;
     234    if (!y->getNumElements()) {
     235      if (!rowScale) {
     236        if (scalar==-1.0) {
     237          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     238            double value = 0.0;
     239            CoinBigIndex j;
     240            for (j=columnStart[iColumn];
     241                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     242              int iRow = row[j];
     243              value += pi[iRow]*elementByColumn[j];
     244            }
     245            if (fabs(value)>zeroTolerance) {
     246              index[numberNonZero++]=iColumn;
     247              array[iColumn]=-value;
     248            }
     249          }
     250        } else if (scalar==1.0) {
     251          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     252            double value = 0.0;
     253            CoinBigIndex j;
     254            for (j=columnStart[iColumn];
     255                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     256              int iRow = row[j];
     257              value += pi[iRow]*elementByColumn[j];
     258            }
     259            if (fabs(value)>zeroTolerance) {
     260              index[numberNonZero++]=iColumn;
     261              array[iColumn]=value;
     262            }
     263          }
     264        } else {
     265          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     266            double value = 0.0;
     267            CoinBigIndex j;
     268            for (j=columnStart[iColumn];
     269                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     270              int iRow = row[j];
     271              value += pi[iRow]*elementByColumn[j];
     272            }
     273            value *= scalar;
     274            if (fabs(value)>zeroTolerance) {
     275              index[numberNonZero++]=iColumn;
     276              array[iColumn]=value;
     277            }
     278          }
     279        }
     280      } else {
     281        // scaled
     282        if (scalar==-1.0) {
     283          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     284            double value = 0.0;
     285            CoinBigIndex j;
     286            const double * columnScale = model->columnScale();
     287            for (j=columnStart[iColumn];
     288                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     289              int iRow = row[j];
     290              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     291            }
     292            value *= columnScale[iColumn];
     293            if (fabs(value)>zeroTolerance) {
     294              index[numberNonZero++]=iColumn;
     295              array[iColumn]=-value;
     296            }
     297          }
     298        } else if (scalar==1.0) {
     299          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     300            double value = 0.0;
     301            CoinBigIndex j;
     302            const double * columnScale = model->columnScale();
     303            for (j=columnStart[iColumn];
     304                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     305              int iRow = row[j];
     306              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     307            }
     308            value *= columnScale[iColumn];
     309            if (fabs(value)>zeroTolerance) {
     310              index[numberNonZero++]=iColumn;
     311              array[iColumn]=value;
     312            }
     313          }
     314        } else {
     315          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     316            double value = 0.0;
     317            CoinBigIndex j;
     318            const double * columnScale = model->columnScale();
     319            for (j=columnStart[iColumn];
     320                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     321              int iRow = row[j];
     322              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     323            }
     324            value *= scalar*columnScale[iColumn];
     325            if (fabs(value)>zeroTolerance) {
     326              index[numberNonZero++]=iColumn;
     327              array[iColumn]=value;
     328            }
     329          }
    250330        }
    251331      }
    252332    } else {
    253       // scaled
    254       for (iColumn=0;iColumn<numberColumns;iColumn++) {
    255         double value = markVector[iColumn];
    256         markVector[iColumn]=0.0;
    257         CoinBigIndex j;
    258         const double * columnScale = model->columnScale();
    259         double value2 = 0.0;
    260         for (j=columnStart[iColumn];
    261              j<columnStart[iColumn]+columnLength[iColumn];j++) {
    262           int iRow = row[j];
    263           value2 += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    264         }
    265         value += value2*scalar*columnScale[iColumn];
    266         if (fabs(value)>zeroTolerance) {
    267           index[numberNonZero++]=iColumn;
    268           array[iColumn]=value;
     333      double * markVector = y->denseVector(); // not empty
     334      if (!rowScale) {
     335        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     336          double value = markVector[iColumn];
     337          markVector[iColumn]=0.0;
     338          double value2 = 0.0;
     339          CoinBigIndex j;
     340          for (j=columnStart[iColumn];
     341               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     342            int iRow = row[j];
     343            value2 += pi[iRow]*elementByColumn[j];
     344          }
     345          value += value2*scalar;
     346          if (fabs(value)>zeroTolerance) {
     347            index[numberNonZero++]=iColumn;
     348            array[iColumn]=value;
     349          }
     350        }
     351      } else {
     352        // scaled
     353        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     354          double value = markVector[iColumn];
     355          markVector[iColumn]=0.0;
     356          CoinBigIndex j;
     357          const double * columnScale = model->columnScale();
     358          double value2 = 0.0;
     359          for (j=columnStart[iColumn];
     360               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     361            int iRow = row[j];
     362            value2 += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     363          }
     364          value += value2*scalar*columnScale[iColumn];
     365          if (fabs(value)>zeroTolerance) {
     366            index[numberNonZero++]=iColumn;
     367            array[iColumn]=value;
     368          }
    269369        }
    270370      }
  • trunk/ClpPrimalColumnDantzig.cpp

    r73 r137  
    140140  for (iSequence=0;iSequence<number;iSequence++) {
    141141    double value = reducedCost[iSequence];
    142     if (!model_->fixed(iSequence)) {
    143       ClpSimplex::Status status = model_->getStatus(iSequence);
     142    ClpSimplex::Status status = model_->getStatus(iSequence);
     143   
     144    switch(status) {
    144145     
    145       switch(status) {
    146 
    147       case ClpSimplex::basic:
    148         break;
    149       case ClpSimplex::isFree:
    150       case ClpSimplex::superBasic:
    151         if (fabs(value)>bestFreeDj) {
    152           bestFreeDj = fabs(value);
    153           bestFreeSequence = iSequence;
    154         }
    155         break;
    156       case ClpSimplex::atUpperBound:
    157         if (value>bestDj) {
    158           bestDj = value;
    159           bestSequence = iSequence;
    160         }
    161         break;
    162       case ClpSimplex::atLowerBound:
    163         if (value<-bestDj) {
    164           bestDj = -value;
    165           bestSequence = iSequence;
    166         }
     146    case ClpSimplex::basic:
     147    case ClpSimplex::isFixed:
     148      break;
     149    case ClpSimplex::isFree:
     150    case ClpSimplex::superBasic:
     151      if (fabs(value)>bestFreeDj) {
     152        bestFreeDj = fabs(value);
     153        bestFreeSequence = iSequence;
     154      }
     155      break;
     156    case ClpSimplex::atUpperBound:
     157      if (value>bestDj) {
     158        bestDj = value;
     159        bestSequence = iSequence;
     160      }
     161      break;
     162    case ClpSimplex::atLowerBound:
     163      if (value<-bestDj) {
     164        bestDj = -value;
     165        bestSequence = iSequence;
    167166      }
    168167    }
  • trunk/ClpPrimalColumnSteepest.cpp

    r114 r137  
    204204  if (anyUpdates>0) {
    205205    model_->factorization()->updateColumnTranspose(spareRow2,updates);
     206   
    206207    // put row of tableau in rowArray and columnArray
    207208    model_->clpMatrix()->transposeTimes(model_,-1.0,
     
    229230        value -= updateBy[iSequence];
    230231        reducedCost[iSequence] = value;
    231         if (model_->fixed(iSequence+addSequence))
    232           continue;
    233232        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    234233
     
    237236        case ClpSimplex::basic:
    238237          infeasible_->zero(iSequence+addSequence);
     238        case ClpSimplex::isFixed:
    239239          break;
    240240        case ClpSimplex::isFree:
     
    290290  int sequenceOut = model_->sequenceOut();
    291291  if (sequenceOut>=0) {
    292     if (!model_->fixed(sequenceOut)) {
    293       ClpSimplex::Status status = model_->getStatus(sequenceOut);
    294       double value = model_->reducedCost(sequenceOut);
    295      
    296       switch(status) {
    297 
    298       case ClpSimplex::basic:
    299         break;
    300       case ClpSimplex::isFree:
    301       case ClpSimplex::superBasic:
    302         if (fabs(value)>1.0e2*tolerance) {
    303           // we are going to bias towards free (but only if reasonable)
    304           value *= FREE_BIAS;
    305           // store square in list
    306           if (infeas[sequenceOut])
    307             infeas[sequenceOut] = value*value; // already there
    308           else
    309             infeasible_->quickAdd(sequenceOut,value*value);
    310         } else {
    311           infeasible_->zero(sequenceOut);
    312         }
    313         break;
    314       case ClpSimplex::atUpperBound:
    315         if (value>tolerance) {
    316           // store square in list
    317           if (infeas[sequenceOut])
    318             infeas[sequenceOut] = value*value; // already there
    319           else
    320             infeasible_->quickAdd(sequenceOut,value*value);
    321         } else {
    322           infeasible_->zero(sequenceOut);
    323         }
    324         break;
    325       case ClpSimplex::atLowerBound:
    326         if (value<-tolerance) {
    327           // store square in list
    328           if (infeas[sequenceOut])
    329             infeas[sequenceOut] = value*value; // already there
    330           else
    331             infeasible_->quickAdd(sequenceOut,value*value);
    332         } else {
    333           infeasible_->zero(sequenceOut);
    334         }
     292    ClpSimplex::Status status = model_->getStatus(sequenceOut);
     293    double value = model_->reducedCost(sequenceOut);
     294   
     295    switch(status) {
     296     
     297    case ClpSimplex::basic:
     298    case ClpSimplex::isFixed:
     299      break;
     300    case ClpSimplex::isFree:
     301    case ClpSimplex::superBasic:
     302      if (fabs(value)>1.0e2*tolerance) {
     303        // we are going to bias towards free (but only if reasonable)
     304        value *= FREE_BIAS;
     305        // store square in list
     306        if (infeas[sequenceOut])
     307          infeas[sequenceOut] = value*value; // already there
     308        else
     309          infeasible_->quickAdd(sequenceOut,value*value);
     310      } else {
     311        infeasible_->zero(sequenceOut);
     312      }
     313      break;
     314    case ClpSimplex::atUpperBound:
     315      if (value>tolerance) {
     316        // store square in list
     317        if (infeas[sequenceOut])
     318          infeas[sequenceOut] = value*value; // already there
     319        else
     320          infeasible_->quickAdd(sequenceOut,value*value);
     321      } else {
     322        infeasible_->zero(sequenceOut);
     323      }
     324      break;
     325    case ClpSimplex::atLowerBound:
     326      if (value<-tolerance) {
     327        // store square in list
     328        if (infeas[sequenceOut])
     329          infeas[sequenceOut] = value*value; // already there
     330        else
     331          infeasible_->quickAdd(sequenceOut,value*value);
     332      } else {
     333        infeasible_->zero(sequenceOut);
    335334      }
    336335    }
     
    354353    for (iSequence=0;iSequence<number;iSequence++) {
    355354      double value = reducedCost[iSequence];
    356       if (model_->fixed(iSequence+addSequence))
    357         continue;
    358355      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    359356     
     
    362359      case ClpSimplex::basic:
    363360        infeasible_->zero(iSequence+addSequence);
     361      case ClpSimplex::isFixed:
    364362        break;
    365363      case ClpSimplex::isFree:
     
    508506    spareColumn2->setNumElements(0);
    509507#ifdef SOME_DEBUG_1
    510 #if 1
    511508    // check for accuracy
    512509    int iCheck=229;
    513510    printf("weight for iCheck is %g\n",weights_[iCheck]);
    514511    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
    515         !model_->fixed(iCheck))
     512        !model_->getStatus(iCheck)!=isFixed)
    516513      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
    517 #else
    518     // check for accuracy
    519     number = updates->getNumElements();
    520     index = updates->getIndices();
    521     for (j=0;j<number;j++) {
    522       if (model_->getStatus(index[j])!=ClpSimplex::basic&&
    523           !model_->fixed(index[j]))
    524         checkAccuracy(index[j],1.0e-1,updates,spareRow2);
    525     }
    526 #endif
    527514#endif
    528515    updates->setNumElements(0);
     
    697684     
    698685    for (iSequence=0;iSequence<number;iSequence++) {
    699       if (model_->fixed(iSequence))
    700         continue;
    701686      double value = reducedCost[iSequence];
    702687      ClpSimplex::Status status = model_->getStatus(iSequence);
     
    705690
    706691      case ClpSimplex::basic:
     692      case ClpSimplex::isFixed:
    707693        break;
    708694      case ClpSimplex::isFree:
     
    938924    for (iSequence=0;iSequence<number;iSequence++) {
    939925      weights_[iSequence]=1.0+ADD_ONE;
    940       if (model_->fixed(iSequence))
    941         continue;
    942       if (model_->getStatus(iSequence)!=ClpSimplex::basic) {
     926      if (model_->getStatus(iSequence)!=ClpSimplex::basic&&
     927          model_->getStatus(iSequence)!=ClpSimplex::isFixed) {
    943928        model_->unpack(alternateWeights_,iSequence);
    944929        double value=ADD_ONE;
  • trunk/ClpSimplex.cpp

    r133 r137  
    141141ClpSimplex::gutsOfSolution ( const double * rowActivities,
    142142                             const double * columnActivities,
     143                             double * givenDuals,
     144                             const double * givenPrimals,
    143145                             bool valuesPass)
    144146{
     
    164166  // do work
    165167  computePrimals(rowActivities, columnActivities);
     168  // If necessary - override results
     169  if (givenPrimals) {
     170    memcpy(columnActivityWork_,givenPrimals,numberColumns_*sizeof(double));
     171    memset(rowActivityWork_,0,numberRows_*sizeof(double));
     172    times(-1.0,columnActivityWork_,rowActivityWork_);
     173  }
    166174  double objectiveModification = 0.0;
    167175  if (algorithm_>0&&nonLinearCost_!=NULL) {
     
    222230  }
    223231
    224   computeDuals();
     232  computeDuals(givenDuals);
    225233
    226234  // now check solutions
     
    338346// now dual side
    339347void
    340 ClpSimplex::computeDuals()
    341 {
    342   double slackValue = factorization_->slackValue();
     348ClpSimplex::computeDuals(double * givenDjs)
     349{
    343350  //work space
    344351  CoinIndexedVector  * workSpace = rowArray_[0];
     
    352359  workSpace->checkClear();
    353360#endif
    354   for (iRow=0;iRow<numberRows_;iRow++) {
    355     int iPivot=pivotVariable_[iRow];
    356     array[iRow]=cost_[iPivot];
     361  if (!givenDjs) {
     362    for (iRow=0;iRow<numberRows_;iRow++) {
     363      int iPivot=pivotVariable_[iRow];
     364      array[iRow]=cost_[iPivot];
     365    }
     366  } else {
     367    // dual values pass - djs may not be zero
     368    for (iRow=0;iRow<numberRows_;iRow++) {
     369      int iPivot=pivotVariable_[iRow];
     370      // make sure zero if done
     371      if (!pivoted(iPivot))
     372        givenDjs[iPivot]=0.0;
     373      array[iRow]=cost_[iPivot]-givenDjs[iPivot];
     374    }
    357375  }
    358376  ClpDisjointCopyN ( array, numberRows_ , save);
     
    369387    ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    370388    transposeTimes(-1.0,array,reducedCostWork_);
    371     for (iRow=0;iRow<numberRows_;iRow++) {
    372       int iPivot=pivotVariable_[iRow];
    373       if (iPivot>=numberColumns_) {
    374         // slack
    375         //work[iRow] += slackValue*array[iPivot-numberColumns_];
    376         //work[iRow] += rowObjectiveWork_[iPivot-numberColumns_];
    377         work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
    378         - slackValue*array[iPivot-numberColumns_];
    379       } else {
    380         // column
     389    if (!givenDjs) {
     390      for (iRow=0;iRow<numberRows_;iRow++) {
     391        int iPivot=pivotVariable_[iRow];
     392        if (iPivot>=numberColumns_) {
     393          // slack
     394          work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
     395            + array[iPivot-numberColumns_];
     396        } else {
     397          // column
    381398        work[iRow] = reducedCostWork_[iPivot];
    382       }
    383       if (fabs(work[iRow])>largestDualError_) {
    384         largestDualError_=fabs(work[iRow]);
     399        }
     400        if (fabs(work[iRow])>largestDualError_) {
     401          largestDualError_=fabs(work[iRow]);
     402        }
     403      }
     404    } else {
     405      for (iRow=0;iRow<numberRows_;iRow++) {
     406        int iPivot=pivotVariable_[iRow];
     407        if (iPivot>=numberColumns_) {
     408          // slack
     409          work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
     410            + array[iPivot-numberColumns_]-givenDjs[iPivot];
     411        } else {
     412          // column
     413        work[iRow] = reducedCostWork_[iPivot]- givenDjs[iPivot];
     414        }
     415        if (fabs(work[iRow])>largestDualError_) {
     416          largestDualError_=fabs(work[iRow]);
     417          assert (largestDualError_<1.0e-7);
     418        }
    385419      }
    386420    }
     
    410444        work[iRow]=0.0;
    411445      }
     446    } else {
     447      break;
    412448    }
    413449  }
     
    418454    double value = array[iRow];
    419455    dual_[iRow]=value;
    420     value = rowObjectiveWork_[iRow] - value*slackValue;
     456    value += rowObjectiveWork_[iRow];
    421457    rowReducedCost_[iRow]=value;
    422458  }
    423459  ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    424460  transposeTimes(-1.0,dual_,reducedCostWork_);
     461  // If necessary - override results
     462  if (givenDjs) {
     463    // restore accurate duals
     464    memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double));
     465  }
     466
    425467  delete [] previous;
    426468  delete [] array;
     
    437479    createRim(7+8+16);
    438480    // do work
    439     gutsOfSolution ( rowActivities, columnActivities);
     481    gutsOfSolution ( rowActivities, columnActivities,NULL,NULL);
    440482    // release extra memory
    441483    deleteRim(false);
     
    553595            }
    554596            break;
     597          case ClpSimplex::isFixed:
    555598          case atLowerBound:
    556599            rowActivityWork_[iRow]=rowLowerWork_[iRow];
     
    635678            break;
    636679          case atLowerBound:
     680          case ClpSimplex::isFixed:
    637681            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    638682            if (columnActivityWork_[iColumn]<-largeValue_) {
     
    805849    }
    806850    numberRefinements_=1;
     851    // set fixed if they are
     852    for (iRow=0;iRow<numberRows_;iRow++) {
     853      if (getRowStatus(iRow)!=basic ) {
     854        if (rowLowerWork_[iRow]==rowUpperWork_[iRow]) {
     855          rowActivityWork_[iRow]=rowLowerWork_[iRow];
     856          setRowStatus(iRow,isFixed);
     857        }
     858      }
     859    }
     860    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     861      if (getColumnStatus(iColumn)!=basic ) {
     862        if (columnLowerWork_[iColumn]==columnUpperWork_[iColumn]) {
     863          columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
     864          setColumnStatus(iColumn,isFixed);
     865        }
     866      }
     867    }
    807868  }
    808869  int status=-99;
     
    9421003
    9431004  // sparse methods
    944   if (factorization_->sparseThreshold()) {
     1005  //if (factorization_->sparseThreshold()) {
    9451006    // get default value
    9461007    factorization_->sparseThreshold(0);
    9471008    factorization_->goSparse();
    948   }
     1009    //}
    9491010 
    9501011  return status;
     
    9831044    //assert( getStatus(sequenceOut_)== basic);
    9841045    setStatus(sequenceIn_,basic);
    985     if (directionOut_>0) {
    986       // going to lower
    987       setStatus(sequenceOut_,atLowerBound);
     1046    if (upper_[sequenceOut_]-lower_[sequenceOut_]>0) {
     1047      if (directionOut_>0) {
     1048        // going to lower
     1049        setStatus(sequenceOut_,atLowerBound);
     1050      } else {
     1051        // going to upper
     1052        setStatus(sequenceOut_,atUpperBound);
     1053      }
    9881054    } else {
    989       // going to upper
    990       setStatus(sequenceOut_,atUpperBound);
     1055      // fixed
     1056      setStatus(sequenceOut_,isFixed);
    9911057    }
    9921058    solution_[sequenceOut_]=valueOut_;
     
    14451511    if (infeasibility>largestSolutionError_)
    14461512      largestSolutionError_=infeasibility;
    1447     if (rowLowerWork_[iRow]!=rowUpperWork_[iRow])
    1448       clearFixed(iRow+numberColumns_);
    1449     else
    1450       setFixed(iRow+numberColumns_);
    14511513  }
    14521514  solution = columnActivityWork_;
     
    14731535    if (infeasibility>largestSolutionError_)
    14741536      largestSolutionError_=infeasibility;
    1475     if (columnUpperWork_[iColumn]-columnLowerWork_[iColumn]>primalTolerance_)
    1476       clearFixed(iColumn);
    1477     else
    1478       setFixed(iColumn);
    14791537  }
    14801538}
     
    22882346// dual
    22892347#include "ClpSimplexDual.hpp"
    2290 int ClpSimplex::dual ( )
    2291 {
    2292   return ((ClpSimplexDual *) this)->dual();
     2348int ClpSimplex::dual (int ifValuesPass )
     2349{
     2350  return ((ClpSimplexDual *) this)->dual(ifValuesPass);
    22932351}
    22942352// primal
     
    33363394            case basic:
    33373395              thisWay=0;
     3396            case ClpSimplex::isFixed:
    33383397              break;
    33393398            case isFree:
     
    35713630                    case basic:
    35723631                      thisWay=0;
     3632                    case isFixed:
    35733633                      break;
    35743634                    case isFree:
     
    36193679             
    36203680            case basic:
     3681            case ClpSimplex::isFixed:
    36213682              break;
    36223683            case isFree:
     
    38593920  } else {
    38603921    // just for now - recompute anyway
    3861     gutsOfSolution(rowActivityWork_,columnActivityWork_);
     3922    gutsOfSolution(rowActivityWork_,columnActivityWork_,NULL,NULL);
    38623923  }
    38633924  return returnCode;
     
    39303991    factorization_->slackValue(-1.0);
    39313992    factorization_->zeroTolerance(1.0e-13);
     3993    // Switch off dense
     3994    int saveThreshold = factorization_->denseThreshold();
     3995    factorization_->setDenseThreshold(0);
    39323996    numberIterations_=0;
    39333997
     
    39604024      // for this we need clean basis so it is after factorize
    39614025      numberThrownOut=gutsOfSolution(rowActivityWork_,columnActivityWork_,
     4026                                     NULL,NULL,
    39624027                                     ifValuesPass!=0);
    39634028      totalNumberThrownOut+= numberThrownOut;
     
    39694034        <<totalNumberThrownOut
    39704035        <<CoinMessageEol;
     4036    // Switch back dense
     4037    factorization_->setDenseThreshold(saveThreshold);
    39714038   
    39724039    problemStatus_ = -1;
     
    40244091  // put back original costs and then check
    40254092  createRim(4);
    4026   gutsOfSolution(rowActivityWork_, columnActivityWork_);
     4093  gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    40274094  return (primalFeasible()&&dualFeasible());
    40284095}
  • trunk/ClpSimplexDual.cpp

    r130 r137  
    107107#include <stdio.h>
    108108#include <iostream>
     109//#define CLP_DEBUG 1
    109110#define CHECK_DJ
    110111// dual
    111 int ClpSimplexDual::dual ( )
     112int ClpSimplexDual::dual (int ifValuesPass )
    112113{
    113114
     
    202203  ClpDataSave data = saveData();
    203204
     205  // If values pass then save given duals round check solution
     206  double * saveDuals = NULL;
     207  if (ifValuesPass) {
     208    saveDuals = new double [numberRows_+numberColumns_];
     209    memcpy(saveDuals,dual_,numberRows_*sizeof(double));
     210  }
     211
    204212  // sanity check
    205213  // initialize - no values pass and algorithm_ is -1
     
    212220    // looks okay
    213221   
     222    // If values pass then scale pi
     223    if (ifValuesPass) {
     224      int i;
     225      if (scalingFlag_>0) {
     226        for (i=0;i<numberRows_;i++) {
     227          dual_[i] = saveDuals[i]/rowScale_[i];
     228        }
     229      } else {
     230        memcpy(dual_,saveDuals,numberRows_*sizeof(double));
     231      }
     232      // now create my duals
     233      for (i=0;i<numberRows_;i++) {
     234        // slack
     235        double value = dual_[i];
     236        value += rowObjectiveWork_[i];
     237        saveDuals[i+numberColumns_]=value;
     238      }
     239      ClpDisjointCopyN(objectiveWork_,numberColumns_,saveDuals);
     240      transposeTimes(-1.0,dual_,saveDuals);
     241      memcpy(dj_,saveDuals,(numberColumns_+numberRows_)*sizeof(double));
     242      // set up possible ones
     243      for (i=0;i<numberRows_+numberColumns_;i++)
     244        clearPivoted(i);
     245      int iRow;
     246      for (iRow=0;iRow<numberRows_;iRow++) {
     247        int iPivot=pivotVariable_[iRow];
     248        if (fabs(saveDuals[iPivot])>dualTolerance_) {
     249          if (getStatus(iPivot)!=isFree)
     250            setPivoted(iPivot);
     251        }
     252      }
     253    }
     254
    214255    double objectiveChange;
    215256    numberFake_ =0; // Number of variables at fake bounds
     
    217258   
    218259    int lastCleaned=0; // last time objective or bounds cleaned up
    219    
    220260
    221261    // Progress indicator
     
    255295#endif
    256296      // may factorize, checks if problem finished
    257       statusOfProblemInDual(lastCleaned,factorType,progress);
     297      statusOfProblemInDual(lastCleaned,factorType,progress,saveDuals);
     298      // If values pass then do easy ones on first time
     299      if (ifValuesPass&&!numberIterations_) {
     300        doEasyOnesInValuesPass(saveDuals);
     301      }
    258302     
    259303      // Say good factorization
     
    266310     
    267311      // Do iterations
    268       whileIterating();
     312      whileIterating(saveDuals);
    269313    }
    270314  }
     
    274318  // clean up
    275319  finish();
     320  delete [] saveDuals;
    276321
    277322  // Restore any saved stuff
     
    288333 */
    289334int
    290 ClpSimplexDual::whileIterating()
     335ClpSimplexDual::whileIterating(double * & givenDuals)
    291336{
    292337#ifdef CLP_DEBUG
     
    301346  // status -3 to go to top without an invert
    302347  int returnCode = -1;
     348
     349  // compute average infeasibility for backward test
     350  double averagePrimalInfeasibility = sumPrimalInfeasibilities_/
     351    ((double ) (numberPrimalInfeasibilities_+1));
     352
     353  // If values pass then get list of candidates
     354  int * candidateList = NULL;
     355  int numberCandidates = 0;
     356#ifdef CLP_DEBUG
     357  bool wasInValuesPass= (givenDuals!=NULL);
     358#endif
     359  int candidate=-1;
     360  if (givenDuals) {
     361    candidateList = new int[numberRows_];
     362    // move reduced costs across
     363    memcpy(dj_,givenDuals,(numberRows_+numberColumns_)*sizeof(double));
     364    int iRow;
     365    for (iRow=0;iRow<numberRows_;iRow++) {
     366      int iPivot=pivotVariable_[iRow];
     367      if (flagged(iPivot))
     368        continue;
     369      if (fabs(dj_[iPivot])>dualTolerance_) {
     370        if (pivoted(iPivot))
     371          candidateList[numberCandidates++]=iRow;
     372      } else {
     373        clearPivoted(iPivot);
     374      }
     375    }
     376    // and set first candidate
     377    if (!numberCandidates) {
     378      delete [] candidateList;
     379      delete [] givenDuals;
     380      givenDuals=NULL;
     381      candidateList=NULL;
     382      int iRow;
     383      for (iRow=0;iRow<numberRows_;iRow++) {
     384        int iPivot=pivotVariable_[iRow];
     385        clearPivoted(iPivot);
     386      }
     387    }
     388  }
    303389  while (problemStatus_==-1) {
     390#ifdef CLP_DEBUG
     391    if (givenDuals) {
     392      double value5=0.0;
     393      int i;
     394      for (i=0;i<numberRows_+numberColumns_;i++) {
     395        if (dj_[i]<-1.0e-6)
     396          value5 += dj_[i]*upper_[i];
     397        else if (dj_[i] >1.0e-6)
     398          value5 += dj_[i]*lower_[i];
     399      }
     400      printf("Values objective Value %g\n",value5);
     401    }
     402    if ((handler_->logLevel()&32)&&wasInValuesPass) {
     403      double value5=0.0;
     404      int i;
     405      for (i=0;i<numberRows_+numberColumns_;i++) {
     406        if (dj_[i]<-1.0e-6)
     407          value5 += dj_[i]*upper_[i];
     408        else if (dj_[i] >1.0e-6)
     409          value5 += dj_[i]*lower_[i];
     410      }
     411      printf("Values objective Value %g\n",value5);
     412      {
     413        double tolerance = 1.0e-7;
     414        int i;
     415        for (i=0;i<numberRows_+numberColumns_;i++) {
     416          int iSequence = i;
     417          double oldValue;
     418         
     419          switch(getStatus(iSequence)) {
     420           
     421          case basic:
     422          case ClpSimplex::isFixed:
     423            break;
     424          case isFree:
     425          case superBasic:
     426            abort();
     427            break;
     428          case atUpperBound:
     429            oldValue = dj_[iSequence];
     430            assert (oldValue<=tolerance);
     431            assert (fabs(solution_[iSequence]-upper_[iSequence])<1.0e-7);
     432            break;
     433          case atLowerBound:
     434            oldValue = dj_[iSequence];
     435            assert (oldValue>=-tolerance);
     436            assert (fabs(solution_[iSequence]-lower_[iSequence])<1.0e-7);
     437            break;
     438          }
     439        }
     440      }
     441    }
     442#endif
    304443#ifdef CLP_DEBUG
    305444    {
     
    326465      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
    327466      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
    328       gutsOfSolution(rowActivityWork_,columnActivityWork_);
     467      gutsOfSolution(rowActivityWork_,columnActivityWork_,NULL,NULL);
    329468      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
    330469             numberIterations_,
     
    344483#endif
    345484#if 0
    346     {
     485    if (!factorization_->pivots()){
    347486      int iPivot;
    348487      double * array = rowArray_[3]->denseVector();
     
    352491        unpack(rowArray_[3],iSequence);
    353492        factorization_->updateColumn(rowArray_[2],rowArray_[3],false);
    354         assert (array[iPivot]==1.0);
     493        assert (fabs(array[iPivot]-1.0)<1.0e-4);
    355494        array[iPivot]=0.0;
    356495        for (i=0;i<numberRows_;i++)
    357           assert (array[i]==0.0);
     496          assert (fabs(array[i])<1.0e-4);
    358497        rowArray_[3]->clear();
    359498      }
     
    382521          break;
    383522        case atLowerBound:
     523        case ClpSimplex::isFixed:
    384524          assert (fabs(value-lowerValue)<=primalTolerance_) ;
    385525          break;
     
    393533    // choose row to go out
    394534    // dualRow will go to virtual row pivot choice algorithm
    395     dualRow();
     535    // make sure values pass off if it should be
     536    if (numberCandidates)
     537      candidate = candidateList[--numberCandidates];
     538    else
     539      candidate=-1;
     540    dualRow(candidate);
    396541    if (pivotRow_>=0) {
    397542      // we found a pivot row
     
    401546      // check accuracy of weights
    402547      dualRowPivot_->checkAccuracy();
     548      // Get good size for pivot
     549      double acceptablePivot=1.0e-7;
     550      if (factorization_->pivots()>10)
     551        acceptablePivot=1.0e-5; // if we have iterated be more strict
     552      else if (factorization_->pivots()>5)
     553        acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
    403554      // get sign for finding row of tableau
    404       rowArray_[0]->insert(pivotRow_,directionOut_);
    405       factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
    406       // put row of tableau in rowArray[0] and columnArray[0]
    407       matrix_->transposeTimes(this,-1.0,
     555      if (candidate<0) {
     556        // normal iteration
     557        rowArray_[0]->insert(pivotRow_,directionOut_);
     558        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
     559        // put row of tableau in rowArray[0] and columnArray[0]
     560        matrix_->transposeTimes(this,-1.0,
    408561                              rowArray_[0],columnArray_[1],columnArray_[0]);
    409       // rowArray has pi equivalent
    410       // do ratio test
    411       dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
    412                  rowArray_[3]);
     562        // do ratio test for normal iteration
     563        dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
     564                 rowArray_[3],acceptablePivot);
     565#if 0
     566        if (sequenceIn_<0) {
     567          // no obvious candidate - see if other way possible
     568          // This will only happen if dual bound too low (or actually infeasible)
     569          if (upperOut_-lowerOut_<1.0e15) {
     570            if (directionOut_==1) {
     571              directionOut_ = -1;
     572              dualOut_ = -(valueOut_ - upperOut_);
     573            } else {
     574              directionOut_ = 1;
     575              dualOut_ = -(lowerOut_ - valueOut_);
     576            }
     577            rowArray_[3]->clear();
     578            columnArray_[1]->clear();
     579            int i;
     580            double * work;
     581            int number;
     582            int * which;
     583            work = rowArray_[0]->denseVector();
     584            number = rowArray_[0]->getNumElements();
     585            which = rowArray_[0]->getIndices();
     586            for (i=0;i<number;i++) {
     587              int iSequence = which[i];
     588              work[iSequence] = - work[iSequence];
     589            }
     590            work = columnArray_[0]->denseVector();
     591            number = columnArray_[0]->getNumElements();
     592            which = columnArray_[0]->getIndices();
     593            for (i=0;i<number;i++) {
     594              int iSequence = which[i];
     595              work[iSequence] = - work[iSequence];
     596            }
     597            // and do ratio test again
     598            dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
     599                       rowArray_[3],acceptablePivot);
     600            if (sequenceIn_<0&&dualBound_<1.0e12)
     601              dualBound_ *= 1.1;
     602          }
     603        }
     604#endif
     605      } else {
     606        rowArray_[0]->insert(pivotRow_,directionOut_);
     607        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
     608        // put row of tableau in rowArray[0] and columnArray[0]
     609        matrix_->transposeTimes(this,-1.0,
     610                              rowArray_[0],columnArray_[1],columnArray_[0]);
     611        acceptablePivot *= 10.0;
     612        // do ratio test
     613        checkPossibleValuesMove(rowArray_[0],columnArray_[0],
     614                                            acceptablePivot);
     615       
     616        // recompute true dualOut_
     617        if (directionOut_<0) {
     618          dualOut_ = valueOut_ - upperOut_;
     619        } else {
     620          dualOut_ = lowerOut_ - valueOut_;
     621        }
     622        // check what happened if was values pass
     623        // may want to move part way i.e. movement
     624        bool normalIteration = (sequenceIn_!=sequenceOut_);
     625
     626        clearPivoted(sequenceOut_);  // make sure won't be done again
     627        // see if end of values pass
     628        if (!numberCandidates) {
     629          int iRow;
     630          delete [] candidateList;
     631          delete [] givenDuals;
     632          candidate=-2; // -2 signals end
     633          givenDuals=NULL;
     634          handler_->message(CLP_END_VALUES_PASS,messages_)
     635            <<numberIterations_;
     636          candidateList=NULL;
     637          for (iRow=0;iRow<numberRows_;iRow++) {
     638            int iPivot=pivotVariable_[iRow];
     639            //assert (fabs(dj_[iPivot]),1.0e-5);
     640            clearPivoted(iPivot);
     641          }
     642        }
     643        if (!normalIteration) {
     644          double objectiveChange=0.0;
     645          updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
     646                            objectiveChange);
     647          if (candidate==-2)
     648            problemStatus_=-2;
     649          continue; // skip rest of iteration
     650        } else {
     651          // recompute dualOut_
     652          if (directionOut_<0) {
     653            dualOut_ = valueOut_ - upperOut_;
     654          } else {
     655            dualOut_ = lowerOut_ - valueOut_;
     656          }
     657        }
     658      }
    413659      if (sequenceIn_>=0) {
    414660        // normal iteration
     
    473719        // which will change basic solution
    474720        if (nswapped) {
    475 #ifdef CLP_DEBUG
    476           if ((handler_->logLevel()&16))
    477             printf("old dualOut_ %g, v %g, l %g, u %g - new ",
    478                    dualOut_,valueOut_,lowerOut_,upperOut_);
    479           double oldOut=dualOut_;
    480 #endif
     721          assert (candidate<0);
    481722          factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
    482723          dualRowPivot_->updatePrimalSolution(rowArray_[2],
     
    490731            dualOut_ = lowerOut_ - valueOut_;
    491732          }
    492 #ifdef CLP_DEBUG
    493           if ((handler_->logLevel()&16))
    494             printf("%g\n",dualOut_);
    495           assert(dualOut_<=oldOut);
    496 #endif
    497           if(dualOut_<-10.0e-8&&factorization_->pivots()&&
     733          if(dualOut_<-max(1.0e-12*averagePrimalInfeasibility,1.0e-8)
     734             &&factorization_->pivots()&&
    498735             getStatus(sequenceIn_)!=isFree) {
    499736            // going backwards - factorize
     
    581818        }
    582819        // update primal solution
    583         if (theta_<0.0) {
     820        if (theta_<0.0&&candidate==-1) {
    584821#ifdef CLP_DEBUG
    585822          if (handler_->logLevel()&32)
     
    611848        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
    612849        // outgoing
     850        // set dj to zero unless values pass
    613851        if (directionOut_>0) {
    614852          valueOut_ = lowerOut_;
    615           //setStatus(sequenceOut_,atLowerBound);
    616           dj_[sequenceOut_]=theta_;
     853          if (candidate==-1)
     854            dj_[sequenceOut_] = theta_;
    617855        } else {
    618856          valueOut_ = upperOut_;
    619           //setStatus(sequenceOut_,atUpperBound);
    620           dj_[sequenceOut_]=-theta_;
     857          if (candidate==-1)
     858            dj_[sequenceOut_] = -theta_;
    621859        }
    622860        solution_[sequenceOut_]=valueOut_;
     
    629867          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
    630868#endif
    631         if (whatNext==1) {
     869        if (whatNext==1||candidate==-2) {
    632870          problemStatus_ =-2; // refactorize
    633871        } else if (whatNext==2) {
     
    696934      break;
    697935    }
     936  }
     937  if (givenDuals) {
     938    memcpy(givenDuals,dj_,(numberRows_+numberColumns_)*sizeof(double));
     939    // get rid of any values pass array
     940    delete [] candidateList;
    698941  }
    699942  return returnCode;
     
    713956  outputArray->clear();
    714957 
    715   double * work;
    716   int number;
    717   int * which;
    718958 
    719959  int numberInfeasibilities=0;
     
    729969 
    730970  double changeObj=0.0;
    731  
    732   int iSection;
    733  
    734   for (iSection=0;iSection<2;iSection++) {
     971
     972  // Coding is very similar but we can save a bit by splitting
     973  // Do rows
     974  if (!fullRecompute) {
    735975    int i;
    736     double * solution = solutionRegion(iSection);
    737     double * reducedCost = djRegion(iSection);
    738     double * lower = lowerRegion(iSection);
    739     double * upper = upperRegion(iSection);
    740     double * cost = costRegion(iSection);
    741     int addSequence;
    742     if (!iSection) {
    743       addSequence = numberColumns_;
    744       work = rowArray->denseVector();
    745       number = rowArray->getNumElements();
    746       which = rowArray->getIndices();
    747     } else {
    748       // set number of infeasibilities in row array
    749       addSequence=0;
    750       numberRowInfeasibilities=numberInfeasibilities;
    751       rowArray->setNumElements(numberInfeasibilities);
    752       numberInfeasibilities=0;
    753       work = columnArray->denseVector();
    754       number = columnArray->getNumElements();
    755       which = columnArray->getIndices();
    756     }
     976    double * reducedCost = djRegion(0);
     977    double * lower = lowerRegion(0);
     978    double * upper = upperRegion(0);
     979    double * cost = costRegion(0);
     980    double * work;
     981    int number;
     982    int * which;
     983    work = rowArray->denseVector();
     984    number = rowArray->getNumElements();
     985    which = rowArray->getIndices();
     986    for (i=0;i<number;i++) {
     987      int iSequence = which[i];
     988      double alphaI = work[iSequence];
     989      double value = reducedCost[iSequence]-theta*alphaI;
     990      work[iSequence]=0.0;
     991      reducedCost[iSequence]=value;
     992     
     993      Status status = getStatus(iSequence+numberColumns_);
     994      // more likely to be at upper bound ?
     995      if (status==atUpperBound) {
     996       
     997        if (value>tolerance) {
     998          // to lower bound (if swap)
     999          // put back alpha
     1000          which[numberInfeasibilities++]=iSequence;
     1001          work[iSequence]=alphaI;
     1002          double movement = lower[iSequence]-upper[iSequence];
     1003#ifdef CLP_DEBUG
     1004          if ((handler_->logLevel()&32))
     1005            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     1006                   0,iSequence,value,alphaI,movement);
     1007#endif
     1008          FakeBound bound = getFakeBound(iSequence+numberColumns_);
     1009          changeObj += movement*cost[iSequence];
     1010          if (bound==ClpSimplexDual::bothFake||
     1011              bound==ClpSimplexDual::lowerFake)
     1012            numberAtFake++;
     1013          outputArray->quickAdd(iSequence,-movement);
     1014        }
     1015      } else if (status==atLowerBound) {
     1016       
     1017        if (value<-tolerance) {
     1018          // to upper bound
     1019          // put back alpha
     1020          which[numberInfeasibilities++]=iSequence;
     1021          work[iSequence]=alphaI;
     1022          double movement = upper[iSequence] - lower[iSequence];
     1023#ifdef CLP_DEBUG
     1024          if ((handler_->logLevel()&32))
     1025            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     1026                   0,iSequence,value,alphaI,movement);
     1027#endif
     1028          FakeBound bound = getFakeBound(iSequence+numberColumns_);
     1029          changeObj += movement*cost[iSequence];
     1030          if (bound==ClpSimplexDual::bothFake||
     1031              bound==ClpSimplexDual::upperFake)
     1032            numberAtFake++;
     1033          outputArray->quickAdd(iSequence,-movement);
     1034        }
     1035      }
     1036    }
     1037  } else  {
     1038    int i;
     1039    double * solution = solutionRegion(0);
     1040    double * reducedCost = djRegion(0);
     1041    double * lower = lowerRegion(0);
     1042    double * upper = upperRegion(0);
     1043    double * cost = costRegion(0);
     1044    double * work;
     1045    int number;
     1046    int * which;
     1047    work = rowArray->denseVector();
     1048    number = rowArray->getNumElements();
     1049    which = rowArray->getIndices();
     1050    for (i=0;i<number;i++) {
     1051      int iSequence = which[i];
     1052      double alphaI = work[iSequence];
     1053      double value = reducedCost[iSequence]-theta*alphaI;
     1054      work[iSequence]=0.0;
     1055      reducedCost[iSequence]=value;
     1056     
     1057      Status status = getStatus(iSequence+numberColumns_);
     1058      // more likely to be at upper bound ?
     1059      if (status==atUpperBound) {
     1060        double movement=0.0;
     1061        FakeBound bound = getFakeBound(iSequence+numberColumns_);
     1062       
     1063        if (value>tolerance) {
     1064          // to lower bound (if swap)
     1065          // put back alpha
     1066          which[numberInfeasibilities++]=iSequence;
     1067          work[iSequence]=alphaI;
     1068          movement = lower[iSequence]-upper[iSequence];
     1069#ifdef CLP_DEBUG
     1070          if ((handler_->logLevel()&32))
     1071            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     1072                   0,iSequence,value,alphaI,movement);
     1073#endif
     1074          changeObj += movement*cost[iSequence];
     1075          if (bound==ClpSimplexDual::bothFake||
     1076              bound==ClpSimplexDual::lowerFake)
     1077            numberAtFake++;
     1078          outputArray->quickAdd(iSequence,-movement);
     1079        } else {
     1080          // at correct bound
     1081          if (bound==ClpSimplexDual::bothFake||
     1082              bound==ClpSimplexDual::upperFake) {
     1083            // but flip if dj would allow
     1084            if (bound==ClpSimplexDual::upperFake&&
     1085                value>=-tolerance) {
     1086              movement = lower[iSequence]-upper[iSequence];
     1087              setStatus(iSequence+numberColumns_,atLowerBound);
     1088              solution[iSequence] = lower[iSequence];
     1089              changeObj += movement*cost[iSequence];
     1090            } else {
     1091              numberAtFake++;
     1092            }
     1093          }
     1094        }
     1095      } else if (status==atLowerBound) {
     1096        double movement=0.0;
     1097        FakeBound bound = getFakeBound(iSequence+numberColumns_);
     1098       
     1099        if (value<-tolerance) {
     1100          // to upper bound
     1101          // put back alpha
     1102          which[numberInfeasibilities++]=iSequence;
     1103          work[iSequence]=alphaI;
     1104          movement = upper[iSequence] - lower[iSequence];
     1105#ifdef CLP_DEBUG
     1106          if ((handler_->logLevel()&32))
     1107            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     1108                   0,iSequence,value,alphaI,movement);
     1109#endif
     1110          changeObj += movement*cost[iSequence];
     1111          if (bound==ClpSimplexDual::bothFake||
     1112              bound==ClpSimplexDual::upperFake)
     1113            numberAtFake++;
     1114          outputArray->quickAdd(iSequence,-movement);
     1115        } else {
     1116          // at correct bound
     1117          if (bound==ClpSimplexDual::bothFake||
     1118              bound==ClpSimplexDual::lowerFake) {
     1119            // but flip if dj would allow
     1120            if (bound==ClpSimplexDual::lowerFake&&
     1121                value<=tolerance) {
     1122              movement = upper[iSequence] - lower[iSequence];
     1123              setStatus(iSequence+numberColumns_,atUpperBound);
     1124              solution[iSequence] = upper[iSequence];
     1125              changeObj += movement*cost[iSequence];
     1126            } else {
     1127              numberAtFake++;
     1128            }
     1129          }
     1130        }
     1131      }
     1132    }
     1133  }
     1134
     1135  // Do columns
     1136  if (!fullRecompute) {
     1137    int i;
     1138    double * reducedCost = djRegion(1);
     1139    double * lower = lowerRegion(1);
     1140    double * upper = upperRegion(1);
     1141    double * cost = costRegion(1);
     1142    double * work;
     1143    int number;
     1144    int * which;
     1145    // set number of infeasibilities in row array
     1146    numberRowInfeasibilities=numberInfeasibilities;
     1147    rowArray->setNumElements(numberInfeasibilities);
     1148    numberInfeasibilities=0;
     1149    work = columnArray->denseVector();
     1150    number = columnArray->getNumElements();
     1151    which = columnArray->getIndices();
    7571152   
    7581153    for (i=0;i<number;i++) {
     
    7631158      reducedCost[iSequence]=value;
    7641159     
    765       if (!fixed(iSequence+addSequence)) {
     1160      Status status = getStatus(iSequence);
     1161      if (status==atLowerBound) {
    7661162        double movement=0.0;
    767         FakeBound bound = getFakeBound(iSequence+addSequence);
    768         Status status = getStatus(iSequence+addSequence);
    769 
    770         switch(status) {
    771          
    772         case basic:
    773         case superBasic:
    774           break;
    775         case isFree:
    776           if (fabs(value)>tolerance) {
    777 #ifdef CLP_DEBUG
    778             if (handler_->logLevel()&32)
    779               printf("%d %d, free has dj of %g, alpha %g\n",
    780                      iSection,iSequence,value,alphaI);
    781 #endif
    782           }
    783           break;
    784         case atUpperBound:
    785           if (value>tolerance) {
    786             // to lower bound (if swap)
    787             // put back alpha
    788             which[numberInfeasibilities++]=iSequence;
    789             work[iSequence]=alphaI;
    790             movement = lower[iSequence]-upper[iSequence];
    791 #ifdef CLP_DEBUG
    792             if ((handler_->logLevel()&32))
    793               printf("%d %d, new dj %g, alpha %g, movement %g\n",
    794                    iSection,iSequence,value,alphaI,movement);
    795 #endif
    796             changeObj += movement*cost[iSequence];
    797             if (bound==ClpSimplexDual::bothFake||
    798                 bound==ClpSimplexDual::lowerFake)
     1163        FakeBound bound = getFakeBound(iSequence);
     1164
     1165        if (value<-tolerance) {
     1166          // to upper bound
     1167          // put back alpha
     1168          which[numberInfeasibilities++]=iSequence;
     1169          work[iSequence]=alphaI;
     1170          movement = upper[iSequence] - lower[iSequence];
     1171#ifdef CLP_DEBUG
     1172          if ((handler_->logLevel()&32))
     1173            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     1174                   1,iSequence,value,alphaI,movement);
     1175#endif
     1176          changeObj += movement*cost[iSequence];
     1177          if (bound==ClpSimplexDual::bothFake||
     1178              bound==ClpSimplexDual::upperFake)
     1179            numberAtFake++;
     1180          matrix_->add(this,outputArray,iSequence,movement);
     1181        }
     1182      } else if (status==atUpperBound) {
     1183        double movement=0.0;
     1184        FakeBound bound = getFakeBound(iSequence);
     1185
     1186        if (value>tolerance) {
     1187          // to lower bound (if swap)
     1188          // put back alpha
     1189          which[numberInfeasibilities++]=iSequence;
     1190          work[iSequence]=alphaI;
     1191          movement = lower[iSequence]-upper[iSequence];
     1192#ifdef CLP_DEBUG
     1193          if ((handler_->logLevel()&32))
     1194            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     1195                   1,iSequence,value,alphaI,movement);
     1196#endif
     1197          changeObj += movement*cost[iSequence];
     1198          if (bound==ClpSimplexDual::bothFake||
     1199              bound==ClpSimplexDual::lowerFake)
     1200            numberAtFake++;
     1201          matrix_->add(this,outputArray,iSequence,movement);
     1202        }
     1203      }
     1204    }
     1205  } else {
     1206    int i;
     1207    double * solution = solutionRegion(1);
     1208    double * reducedCost = djRegion(1);
     1209    double * lower = lowerRegion(1);
     1210    double * upper = upperRegion(1);
     1211    double * cost = costRegion(1);
     1212    double * work;
     1213    int number;
     1214    int * which;
     1215    // set number of infeasibilities in row array
     1216    numberRowInfeasibilities=numberInfeasibilities;
     1217    rowArray->setNumElements(numberInfeasibilities);
     1218    numberInfeasibilities=0;
     1219    work = columnArray->denseVector();
     1220    number = columnArray->getNumElements();
     1221    which = columnArray->getIndices();
     1222   
     1223    for (i=0;i<number;i++) {
     1224      int iSequence = which[i];
     1225      double alphaI = work[iSequence];
     1226      double value = reducedCost[iSequence]-theta*alphaI;
     1227      work[iSequence]=0.0;
     1228      reducedCost[iSequence]=value;
     1229     
     1230      Status status = getStatus(iSequence);
     1231      if (status==atLowerBound) {
     1232        double movement=0.0;
     1233        FakeBound bound = getFakeBound(iSequence);
     1234
     1235        if (value<-tolerance) {
     1236          // to upper bound
     1237          // put back alpha
     1238          which[numberInfeasibilities++]=iSequence;
     1239          work[iSequence]=alphaI;
     1240          movement = upper[iSequence] - lower[iSequence];
     1241#ifdef CLP_DEBUG
     1242          if ((handler_->logLevel()&32))
     1243            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     1244                   1,iSequence,value,alphaI,movement);
     1245#endif
     1246          changeObj += movement*cost[iSequence];
     1247          if (bound==ClpSimplexDual::bothFake||
     1248              bound==ClpSimplexDual::upperFake)
     1249            numberAtFake++;
     1250          matrix_->add(this,outputArray,iSequence,movement);
     1251        } else {
     1252          // at correct bound
     1253          if (bound==ClpSimplexDual::bothFake||
     1254              bound==ClpSimplexDual::lowerFake) {
     1255            // but flip if dj would allow
     1256            if (bound==ClpSimplexDual::lowerFake&&
     1257                value<=tolerance) {
     1258              movement = upper[iSequence] - lower[iSequence];
     1259              setStatus(iSequence,atUpperBound);
     1260              solution[iSequence] = upper[iSequence];
     1261              changeObj += movement*cost[iSequence];
     1262            } else {
    7991263              numberAtFake++;
    800           } else if (fullRecompute) {
    801             // at correct bound
    802             if (bound==ClpSimplexDual::bothFake||
    803                 bound==ClpSimplexDual::upperFake) {
    804               // but flip if dj would allow
    805               if (bound==ClpSimplexDual::upperFake&&
    806                   value>=-tolerance) {
    807                 movement = lower[iSequence]-upper[iSequence];
    808                 setStatus(iSequence+addSequence,atLowerBound);
    809                 solution[iSequence] = lower[iSequence];
    810                 changeObj += movement*cost[iSequence];
    811               } else {
    812                 numberAtFake++;
    813               }
    8141264            }
    8151265          }
    816           break;
    817         case atLowerBound:
    818           if (value<-tolerance) {
    819             // to upper bound
    820             // put back alpha
    821             which[numberInfeasibilities++]=iSequence;
    822             work[iSequence]=alphaI;
    823             movement = upper[iSequence] - lower[iSequence];
    824 #ifdef CLP_DEBUG
    825             if ((handler_->logLevel()&32))
    826               printf("%d %d, new dj %g, alpha %g, movement %g\n",
    827                    iSection,iSequence,value,alphaI,movement);
    828 #endif
    829             changeObj += movement*cost[iSequence];
    830             if (bound==ClpSimplexDual::bothFake||
    831                 bound==ClpSimplexDual::upperFake)
     1266        }
     1267      } else if (status==atUpperBound) {
     1268        double movement=0.0;
     1269        FakeBound bound = getFakeBound(iSequence);
     1270
     1271        if (value>tolerance) {
     1272          // to lower bound (if swap)
     1273          // put back alpha
     1274          which[numberInfeasibilities++]=iSequence;
     1275          work[iSequence]=alphaI;
     1276          movement = lower[iSequence]-upper[iSequence];
     1277#ifdef CLP_DEBUG
     1278          if ((handler_->logLevel()&32))
     1279            printf("%d %d, new dj %g, alpha %g, movement %g\n",
     1280                   1,iSequence,value,alphaI,movement);
     1281#endif
     1282          changeObj += movement*cost[iSequence];
     1283          if (bound==ClpSimplexDual::bothFake||
     1284              bound==ClpSimplexDual::lowerFake)
     1285            numberAtFake++;
     1286          matrix_->add(this,outputArray,iSequence,movement);
     1287        } else {
     1288          // at correct bound
     1289          if (bound==ClpSimplexDual::bothFake||
     1290              bound==ClpSimplexDual::upperFake) {
     1291            // but flip if dj would allow
     1292            if (bound==ClpSimplexDual::upperFake&&
     1293                value>=-tolerance) {
     1294              movement = lower[iSequence]-upper[iSequence];
     1295              setStatus(iSequence,atLowerBound);
     1296              solution[iSequence] = lower[iSequence];
     1297              changeObj += movement*cost[iSequence];
     1298            } else {
    8321299              numberAtFake++;
    833           } else if (fullRecompute) {
    834             // at correct bound
    835             if (bound==ClpSimplexDual::bothFake||
    836                 bound==ClpSimplexDual::lowerFake) {
    837               // but flip if dj would allow
    838               if (bound==ClpSimplexDual::lowerFake&&
    839                   value<=tolerance) {
    840                 movement = upper[iSequence] - lower[iSequence];
    841                 setStatus(iSequence+addSequence,atUpperBound);
    842                 solution[iSequence] = upper[iSequence];
    843                 changeObj += movement*cost[iSequence];
    844               } else {
    845                 numberAtFake++;
    846               }
    847             }
    848           }
    849           break;
    850         }
    851         if (!fullRecompute) {
    852           if (movement) {
    853             if (!iSection) {
    854               // row (sign ?)
    855               outputArray->quickAdd(iSequence,-movement);
    856             } else {
    857               matrix_->add(this,outputArray,iSequence,movement);
    8581300            }
    8591301          }
     
    8841326*/
    8851327void
    886 ClpSimplexDual::dualRow()
     1328ClpSimplexDual::dualRow(int alreadyChosen)
    8871329{
    8881330  // get pivot row using whichever method it is
    889   // first see if any free variables and put them in basis
    890   int nextFree = nextSuperBasic();
    8911331  int chosenRow=-1;
    892   //nextFree=-1; //off
    893   if (nextFree>=0) {
    894     // unpack vector and find a good pivot
    895     unpack(rowArray_[1],nextFree);
    896     factorization_->updateColumn(rowArray_[2],rowArray_[1],false);
    897 
    898     double * work=rowArray_[1]->denseVector();
    899     int number=rowArray_[1]->getNumElements();
    900     int * which=rowArray_[1]->getIndices();
    901     double bestFeasibleAlpha=0.0;
    902     int bestFeasibleRow=-1;
    903     double bestInfeasibleAlpha=0.0;
    904     int bestInfeasibleRow=-1;
    905     int i;
    906 
    907     for (i=0;i<number;i++) {
    908       int iRow = which[i];
    909       double alpha = fabs(work[iRow]);
    910       if (alpha>1.0e-3) {
    911         int iSequence=pivotVariable_[iRow];
    912         double value = solution_[iSequence];
    913         double lower = lower_[iSequence];
    914         double upper = upper_[iSequence];
    915         double infeasibility=0.0;
    916         if (value>upper)
    917           infeasibility = value-upper;
    918         else if (value<lower)
    919           infeasibility = lower-value;
    920         if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
    921           if (!flagged(iSequence)) {
    922             bestInfeasibleAlpha = infeasibility*alpha;
    923             bestInfeasibleRow=iRow;
    924           }
    925         }
    926         if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
    927           bestFeasibleAlpha = alpha;
    928           bestFeasibleRow=iRow;
    929         }
    930       }
    931     }
    932     if (bestInfeasibleRow>=0)
    933       chosenRow = bestInfeasibleRow;
    934     else if (bestFeasibleAlpha>1.0e-2)
    935       chosenRow = bestFeasibleRow;
    936     if (chosenRow>=0)
    937       pivotRow_=chosenRow;
    938     rowArray_[1]->clear();
    939   }
     1332  if (alreadyChosen<0) {
     1333    // first see if any free variables and put them in basis
     1334    int nextFree = nextSuperBasic();
     1335    //nextFree=-1; //off
     1336    if (nextFree>=0) {
     1337      // unpack vector and find a good pivot
     1338      unpack(rowArray_[1],nextFree);
     1339      factorization_->updateColumn(rowArray_[2],rowArray_[1],false);
     1340     
     1341      double * work=rowArray_[1]->denseVector();
     1342      int number=rowArray_[1]->getNumElements();
     1343      int * which=rowArray_[1]->getIndices();
     1344      double bestFeasibleAlpha=0.0;
     1345      int bestFeasibleRow=-1;
     1346      double bestInfeasibleAlpha=0.0;
     1347      int bestInfeasibleRow=-1;
     1348      int i;
     1349     
     1350      for (i=0;i<number;i++) {
     1351        int iRow = which[i];
     1352        double alpha = fabs(work[iRow]);
     1353        if (alpha>1.0e-3) {
     1354          int iSequence=pivotVariable_[iRow];
     1355          double value = solution_[iSequence];
     1356          double lower = lower_[iSequence];
     1357          double upper = upper_[iSequence];
     1358          double infeasibility=0.0;
     1359          if (value>upper)
     1360            infeasibility = value-upper;
     1361          else if (value<lower)
     1362            infeasibility = lower-value;
     1363          if (infeasibility*alpha>bestInfeasibleAlpha&&alpha>1.0e-1) {
     1364            if (!flagged(iSequence)) {
     1365              bestInfeasibleAlpha = infeasibility*alpha;
     1366              bestInfeasibleRow=iRow;
     1367            }
     1368          }
     1369          if (alpha>bestFeasibleAlpha&&(lower>-1.0e20||upper<1.0e20)) {
     1370            bestFeasibleAlpha = alpha;
     1371            bestFeasibleRow=iRow;
     1372          }
     1373        }
     1374      }
     1375      if (bestInfeasibleRow>=0)
     1376        chosenRow = bestInfeasibleRow;
     1377      else if (bestFeasibleAlpha>1.0e-2)
     1378        chosenRow = bestFeasibleRow;
     1379      if (chosenRow>=0)
     1380        pivotRow_=chosenRow;
     1381      rowArray_[1]->clear();
     1382    }
     1383  } else {
     1384    // in values pass
     1385    chosenRow=alreadyChosen;
     1386    pivotRow_=chosenRow;
     1387  }
    9401388  if (chosenRow<0)
    9411389    pivotRow_=dualRowPivot_->pivotRow();
     
    9461394    lowerOut_ = lower_[sequenceOut_];
    9471395    upperOut_ = upper_[sequenceOut_];
    948     // if we have problems we could try other way and hope we get a
    949     // zero pivot?
    950     if (valueOut_>upperOut_) {
    951       directionOut_ = -1;
    952       dualOut_ = valueOut_ - upperOut_;
    953     } else if (valueOut_<lowerOut_) {
    954       directionOut_ = 1;
    955       dualOut_ = lowerOut_ - valueOut_;
    956     } else {
    957       // odd (could be free) - it's feasible - go to nearest
    958       if (valueOut_-lowerOut_<upperOut_-valueOut_) {
     1396    if (alreadyChosen<0) {
     1397      // if we have problems we could try other way and hope we get a
     1398      // zero pivot?
     1399      if (valueOut_>upperOut_) {
     1400        directionOut_ = -1;
     1401        dualOut_ = valueOut_ - upperOut_;
     1402      } else if (valueOut_<lowerOut_) {
    9591403        directionOut_ = 1;
    9601404        dualOut_ = lowerOut_ - valueOut_;
    9611405      } else {
     1406        // odd (could be free) - it's feasible - go to nearest
     1407        if (valueOut_-lowerOut_<upperOut_-valueOut_) {
     1408          directionOut_ = 1;
     1409          dualOut_ = lowerOut_ - valueOut_;
     1410        } else {
     1411          directionOut_ = -1;
     1412          dualOut_ = valueOut_ - upperOut_;
     1413        }
     1414      }
     1415#ifdef CLP_DEBUG
     1416      assert(dualOut_>=0.0);
     1417#endif
     1418    } else {
     1419      // in values pass so just use sign of dj
     1420      // We don't want to go through any barriers so set dualOut low
     1421      // free variables will never be here
     1422      dualOut_ = 1.0e-6;
     1423      if (dj_[sequenceOut_]>0.0) {
     1424        // this will give a -1 in pivot row (as slacks are -1.0)
     1425        directionOut_ = 1;
     1426      } else {
    9621427        directionOut_ = -1;
    963         dualOut_ = valueOut_ - upperOut_;
    964       }
    965     }
    966 #ifdef CLP_DEBUG
    967     assert(dualOut_>=0.0);
    968 #endif
     1428      }
     1429    }
    9691430  }
    9701431  return ;
     
    10001461
    10011462      case basic:
     1463      case ClpSimplex::isFixed:
    10021464        break;
    10031465      case isFree:
     
    11201582void
    11211583ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
    1122                        CoinIndexedVector * columnArray,
    1123                        CoinIndexedVector * spareArray,
    1124                        CoinIndexedVector * spareArray2)
     1584                           CoinIndexedVector * columnArray,
     1585                           CoinIndexedVector * spareArray,
     1586                           CoinIndexedVector * spareArray2,
     1587                           double acceptablePivot)
    11251588{
    11261589  double * work;
     
    11281591  int * which;
    11291592  double * reducedCost;
    1130 
    11311593  int iSection;
    11321594
     
    11361598 
    11371599  double totalThru=0.0; // for when variables flip
    1138   double acceptablePivot=1.0e-7;
    1139   if (factorization_->pivots()>10)
    1140     acceptablePivot=1.0e-5; // if we have iterated be more strict
    1141   else if (factorization_->pivots()>5)
    1142     acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
     1600
    11431601  double bestEverPivot=acceptablePivot;
    11441602  int lastSequence = -1;
     
    12091667  upperTheta = 1.0e31;
    12101668  double freePivot = acceptablePivot;
     1669
    12111670  for (iSection=0;iSection<2;iSection++) {
    12121671
     
    12331692    for (i=0;i<number;i++) {
    12341693      int iSequence = which[i];
    1235       double alpha = work[iSequence];
    1236       if (fixed(iSequence+addSequence)||!alpha)
    1237         continue; // skip fixed ones or (zeroed out)
    1238       double oldValue = reducedCost[iSequence];
    1239       double value = oldValue-tentativeTheta*alpha;
    1240       int keep = 0;
     1694      double alpha;
     1695      double oldValue;
     1696      double value;
     1697      bool keep;
    12411698
    12421699      switch(getStatus(iSequence+addSequence)) {
    12431700         
    12441701      case basic:
     1702      case ClpSimplex::isFixed:
    12451703        break;
    12461704      case isFree:
    12471705      case superBasic:
     1706        alpha = work[iSequence];
     1707        oldValue = reducedCost[iSequence];
    12481708        if (oldValue>dualTolerance_) {
    1249           keep = 2;
     1709          keep = true;
    12501710        } else if (oldValue<-dualTolerance_) {
    1251           keep = 2;
     1711          keep = true;
    12521712        } else {
    12531713          if (fabs(alpha)>10.0*acceptablePivot)
    1254             keep = 2;
    1255 
    1256 
    1257         }
    1258         break;
    1259       case atUpperBound:
    1260         assert (oldValue<=dualTolerance_*1.0001);
    1261         if (value>newTolerance) {
    1262           keep = 1;
    1263           value = oldValue-upperTheta*alpha;
    1264           if (value>newTolerance && -alpha>=acceptablePivot)
    1265             upperTheta = (oldValue-newTolerance)/alpha;
    1266         }
    1267         break;
    1268       case atLowerBound:
    1269         assert (oldValue>=-dualTolerance_*1.0001);
    1270         if (value<-newTolerance) {
    1271           keep = 1;
    1272           value = oldValue-upperTheta*alpha;
    1273           if (value<-newTolerance && alpha>=acceptablePivot)
    1274             upperTheta = (oldValue+newTolerance)/alpha;
    1275         }
    1276         break;
    1277       }
    1278       if (keep) {
    1279         if (keep==2) {
     1714            keep = true;
     1715          else
     1716            keep=false;
     1717        }
     1718        if (keep) {
    12801719          // free - choose largest
    12811720          if (fabs(alpha)>freePivot) {
     
    12841723            theta_=oldValue/alpha;
    12851724          }
    1286         } else {
     1725        }
     1726        break;
     1727      case atUpperBound:
     1728        alpha = work[iSequence];
     1729        oldValue = reducedCost[iSequence];
     1730        value = oldValue-tentativeTheta*alpha;
     1731        assert (oldValue<=dualTolerance_*1.0001);
     1732        if (value>newTolerance) {
     1733          value = oldValue-upperTheta*alpha;
     1734          if (value>newTolerance && -alpha>=acceptablePivot)
     1735            upperTheta = (oldValue-newTolerance)/alpha;
    12871736          // add to list
    12881737          spare[numberRemaining]=alpha;
    12891738          index[numberRemaining++]=iSequence+addSequence;
    12901739        }
     1740        break;
     1741      case atLowerBound:
     1742        alpha = work[iSequence];
     1743        oldValue = reducedCost[iSequence];
     1744        value = oldValue-tentativeTheta*alpha;
     1745        assert (oldValue>=-dualTolerance_*1.0001);
     1746        if (value<-newTolerance) {
     1747          value = oldValue-upperTheta*alpha;
     1748          if (value<-newTolerance && alpha>=acceptablePivot)
     1749            upperTheta = (oldValue+newTolerance)/alpha;
     1750          // add to list
     1751          spare[numberRemaining]=alpha;
     1752          index[numberRemaining++]=iSequence+addSequence;
     1753        }
     1754        break;
    12911755      }
    12921756    }
     
    12961760
    12971761  if (!numberRemaining&&sequenceIn_<0)
    1298     return; // Looks infeasible
     1762    return ; // Looks infeasible
    12991763
    13001764  if (sequenceIn_>=0) {
     
    17902254void
    17912255ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
    1792                                ClpSimplexProgress &progress)
     2256                                      ClpSimplexProgress &progress,
     2257                                      double * givenDuals)
    17932258{
    17942259  if (type==2) {
     
    18052270  double changeCost;
    18062271  bool unflagVariables = true;
    1807 
    18082272  if (problemStatus_>-3||factorization_->pivots()) {
    18092273    // factorize
     
    18522316      }
    18532317    }
    1854     problemStatus_=-3;
     2318    if (problemStatus_!=-4||factorization_->pivots()>10)
     2319      problemStatus_=-3;
    18552320  }
    18562321  // at this stage status is -3 or -4 if looks infeasible
    18572322  // get primal and dual solutions
    1858   gutsOfSolution(rowActivityWork_,columnActivityWork_);
     2323  gutsOfSolution(rowActivityWork_,columnActivityWork_,givenDuals,NULL);
    18592324  // Check if looping
    1860   int loop = progress.looping();
     2325  int loop;
     2326  if (!givenDuals)
     2327    loop = progress.looping();
     2328  else
     2329    loop=-1;
    18612330  int situationChanged=0;
    18622331  if (loop>=0) {
     
    18702339  } else if (loop<-1) {
    18712340    // something may have changed
    1872     gutsOfSolution(rowActivityWork_,columnActivityWork_);
     2341    gutsOfSolution(rowActivityWork_,columnActivityWork_,NULL,NULL);
    18732342    situationChanged=1;
    18742343  }
     
    20762545        createRim(4);
    20772546        // make sure duals are current
    2078         computeDuals();
     2547        computeDuals(givenDuals);
    20792548        // put back bounds as they were if was optimal
    20802549        if (doOriginalTolerance==2) {
     
    20972566                          0.0,objectiveChange);
    20982567        // for now - recompute all
    2099         gutsOfSolution(rowActivityWork_,columnActivityWork_);
     2568        gutsOfSolution(rowActivityWork_,columnActivityWork_,NULL,NULL);
    21002569        //assert(numberDualInfeasibilitiesWithoutFree_==0);
    21012570        if (numberDualInfeasibilities_||situationChanged==2) {
     
    21972666      case isFree:
    21982667      case superBasic:
     2668      case ClpSimplex::isFixed:
    21992669        break;
    22002670      case atUpperBound:
     
    26463116  // for dual we will change bounds using dualBound_
    26473117  // for this we need clean basis so it is after factorize
    2648   gutsOfSolution(rowActivityWork_,columnActivityWork_);
     3118  gutsOfSolution(rowActivityWork_,columnActivityWork_,NULL,NULL);
    26493119  numberFake_ =0; // Number of variables at fake bounds
    26503120  changeBounds(true,NULL,objectiveChange);
     
    27023172    // may factorize, checks if problem finished
    27033173    // should be able to speed this up on first time
    2704     statusOfProblemInDual(lastCleaned,factorType,progress);
     3174    statusOfProblemInDual(lastCleaned,factorType,progress,NULL);
    27053175
    27063176    // Say good factorization
     
    27093179    // Do iterations
    27103180    if (problemStatus_<0) {
    2711       returnCode = whileIterating();
     3181      double * givenPi=NULL;
     3182      returnCode = whileIterating(givenPi);
    27123183      if (!alwaysFinish&&returnCode<1) {
    27133184        double limit = 0.0;
     
    27493220    case isFree:
    27503221    case superBasic:
     3222    case ClpSimplex::isFixed:
    27513223      break;
    27523224    case atUpperBound:
     
    27753247  return 0;
    27763248}
     3249/*
     3250   Row array has row part of pivot row
     3251   Column array has column part.
     3252   This is used in dual values pass
     3253*/
     3254int
     3255ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
     3256                                        CoinIndexedVector * columnArray,
     3257                                        double acceptablePivot)
     3258{
     3259  double * work;
     3260  int number;
     3261  int * which;
     3262  int iSection;
     3263
     3264  double tolerance = dualTolerance_*1.001;
     3265
     3266  double thetaDown = 1.0e31;
     3267  double changeDown ;
     3268  double thetaUp = 1.0e31;
     3269  double bestAlphaDown = acceptablePivot*0.99999;
     3270  double bestAlphaUp = acceptablePivot*0.99999;
     3271  int sequenceDown =-1;
     3272  int sequenceUp = sequenceOut_;
     3273
     3274  double djBasic = dj_[sequenceOut_];
     3275  if (djBasic>0.0) {
     3276    // basic at lower bound so directionOut_ 1 and -1 in pivot row
     3277    // dj will go to zero on other way
     3278    thetaUp = djBasic;
     3279    changeDown = -lower_[sequenceOut_];
     3280  } else {
     3281    // basic at upper bound so directionOut_ -1 and 1 in pivot row
     3282    // dj will go to zero on other way
     3283    thetaUp = -djBasic;
     3284    changeDown = upper_[sequenceOut_];
     3285  }
     3286  bestAlphaUp = 1.0;
     3287  int addSequence;
     3288
     3289  for (iSection=0;iSection<2;iSection++) {
     3290
     3291    int i;
     3292    if (!iSection) {
     3293      work = rowArray->denseVector();
     3294      number = rowArray->getNumElements();
     3295      which = rowArray->getIndices();
     3296      addSequence = numberColumns_;
     3297    } else {
     3298      work = columnArray->denseVector();
     3299      number = columnArray->getNumElements();
     3300      which = columnArray->getIndices();
     3301      addSequence = 0;
     3302    }
     3303   
     3304    for (i=0;i<number;i++) {
     3305      int iSequence = which[i];
     3306      int iSequence2 = iSequence + addSequence;
     3307      double alpha;
     3308      double oldValue;
     3309      double value;
     3310
     3311      switch(getStatus(iSequence2)) {
     3312         
     3313      case basic:
     3314        break;
     3315      case ClpSimplex::isFixed:
     3316        alpha = work[iSequence];
     3317        changeDown += alpha*upper_[iSequence2];
     3318        break;
     3319      case isFree:
     3320      case superBasic:
     3321        alpha = work[iSequence];
     3322        // dj must be effectively zero as dual feasible
     3323        if (fabs(alpha)>bestAlphaUp) {
     3324          thetaDown = 0.0;
     3325          thetaUp = 0.0;
     3326          bestAlphaDown = fabs(alpha);
     3327          bestAlphaUp = bestAlphaUp;
     3328          sequenceDown =iSequence2;
     3329          sequenceUp = sequenceDown;
     3330        }
     3331        break;
     3332      case atUpperBound:
     3333        alpha = work[iSequence];
     3334        oldValue = dj_[iSequence2];
     3335        changeDown += alpha*upper_[iSequence2];
     3336        if (alpha>=acceptablePivot) {
     3337          // might do other way
     3338          value = oldValue+thetaUp*alpha;
     3339          if (value>-tolerance) {
     3340            if (value>tolerance||fabs(alpha)>bestAlphaUp) {
     3341              thetaUp = -oldValue/alpha;
     3342              bestAlphaUp = fabs(alpha);
     3343              sequenceUp = iSequence2;
     3344            }
     3345          }
     3346        } else if (alpha<=-acceptablePivot) {
     3347          // might do this way
     3348          value = oldValue-thetaDown*alpha;
     3349          if (value>-tolerance) {
     3350            if (value>tolerance||fabs(alpha)>bestAlphaDown) {
     3351              thetaDown = oldValue/alpha;
     3352              bestAlphaDown = fabs(alpha);
     3353              sequenceDown = iSequence2;
     3354            }
     3355          }
     3356        }
     3357        break;
     3358      case atLowerBound:
     3359        alpha = work[iSequence];
     3360        oldValue = dj_[iSequence2];
     3361        changeDown += alpha*lower_[iSequence2];
     3362        if (alpha<=-acceptablePivot) {
     3363          // might do other way
     3364          value = oldValue+thetaUp*alpha;
     3365          if (value<tolerance) {
     3366            if (value<-tolerance||fabs(alpha)>bestAlphaUp) {
     3367              thetaUp = -oldValue/alpha;
     3368              bestAlphaUp = fabs(alpha);
     3369              sequenceUp = iSequence2;
     3370            }
     3371          }
     3372        } else if (alpha>=acceptablePivot) {
     3373          // might do this way
     3374          value = oldValue-thetaDown*alpha;
     3375          if (value<tolerance) {
     3376            if (value<-tolerance||fabs(alpha)>bestAlphaDown) {
     3377              thetaDown = oldValue/alpha;
     3378              bestAlphaDown = fabs(alpha);
     3379              sequenceDown = iSequence2;
     3380            }
     3381          }
     3382        }
     3383        break;
     3384      }
     3385    }
     3386  }
     3387  int returnCode;
     3388  thetaUp *= -1.0;
     3389  double changeUp = -thetaUp * changeDown;
     3390  changeDown = -thetaDown*changeDown;
     3391  // if small try other way
     3392  double alphaUp=0.0;
     3393  double alphaDown=0.0;
     3394  if (sequenceDown>=0) {
     3395    if (sequenceDown<numberColumns_)
     3396      alphaDown = columnArray->denseVector()[sequenceDown];
     3397    else
     3398      alphaDown = rowArray->denseVector()[sequenceDown-numberColumns_];
     3399  }
     3400  if (sequenceUp>=0) {
     3401    if (sequenceUp<numberColumns_)
     3402      alphaUp = columnArray->denseVector()[sequenceUp];
     3403    else
     3404      alphaUp = rowArray->denseVector()[sequenceUp-numberColumns_];
     3405  }
     3406  if (max(fabs(thetaDown),fabs(thetaUp))<1.0e-8) {
     3407    // largest
     3408    if (fabs(alphaDown)<fabs(alphaUp)) {
     3409      sequenceDown =-1;
     3410    }
     3411  }
     3412  // choose
     3413  if (changeDown>changeUp&&sequenceDown>=0) {
     3414    theta_ = thetaDown;
     3415    sequenceIn_ = sequenceDown;
     3416#ifdef CLP_DEBUG
     3417    if ((handler_->logLevel()&32))
     3418      printf("predicted way - dirout %d, change %g,%g theta %g\n",
     3419             directionOut_,changeDown,changeUp,theta_);
     3420#endif
     3421    returnCode = 0;
     3422  } else {
     3423    theta_ = thetaUp;
     3424    sequenceIn_ = sequenceUp;
     3425    if (sequenceIn_!=sequenceOut_) {
     3426#ifdef CLP_DEBUG
     3427      if ((handler_->logLevel()&32))
     3428        printf("opposite way - dirout %d, change %g,%g theta %g\n",
     3429               directionOut_,changeDown,changeUp,theta_);
     3430#endif
     3431      returnCode = 0;
     3432    } else {
     3433#ifdef CLP_DEBUG
     3434      if ((handler_->logLevel()&32))
     3435        printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
     3436               directionOut_,changeDown,changeUp,theta_);
     3437#endif
     3438      returnCode = 1;
     3439    }
     3440  }
     3441  if (sequenceIn_<numberColumns_)
     3442    alpha_ = columnArray->denseVector()[sequenceIn_];
     3443  else
     3444    alpha_ = rowArray->denseVector()[sequenceIn_-numberColumns_];
     3445  return returnCode;
     3446}
     3447/*
     3448   This sees if we can move duals in dual values pass.
     3449   This is done before any pivoting
     3450*/
     3451void ClpSimplexDual::doEasyOnesInValuesPass(double * dj)
     3452{
     3453  // Get column copy
     3454  CoinPackedMatrix * columnCopy = matrix();
     3455  // Get a row copy in standard format
     3456  CoinPackedMatrix copy;
     3457  copy.reverseOrderedCopyOf(*columnCopy);
     3458  // get matrix data pointers
     3459  const int * column = copy.getIndices();
     3460  const CoinBigIndex * rowStart = copy.getVectorStarts();
     3461  const int * rowLength = copy.getVectorLengths();
     3462  const double * elementByRow = copy.getElements();
     3463  double tolerance = dualTolerance_*1.001;
     3464
     3465  int iRow;
     3466#ifdef CLP_DEBUG
     3467  {
     3468    double value5=0.0;
     3469    int i;
     3470    for (i=0;i<numberRows_+numberColumns_;i++) {
     3471      if (dj[i]<-1.0e-6)
     3472        value5 += dj[i]*upper_[i];
     3473      else if (dj[i] >1.0e-6)
     3474        value5 += dj[i]*lower_[i];
     3475    }
     3476    printf("Values objective Value before %g\n",value5);
     3477  }
     3478#endif
     3479  // for scaled row
     3480  double * scaled=NULL;
     3481  if (rowScale_)
     3482    scaled = new double[numberColumns_];
     3483  for (iRow=0;iRow<numberRows_;iRow++) {
     3484
     3485    int iSequence = iRow + numberColumns_;
     3486    double djBasic = dj[iSequence];
     3487    if (getRowStatus(iRow)==basic&&fabs(djBasic)>tolerance) {
     3488
     3489      double changeUp ;
     3490      // always -1 in pivot row
     3491      if (djBasic>0.0) {
     3492        // basic at lower bound
     3493        changeUp = -lower_[iSequence];
     3494      } else {
     3495        // basic at upper bound
     3496        changeUp = upper_[iSequence];
     3497      }
     3498      bool canMove=true;
     3499      int i;
     3500      const double * thisElements = elementByRow + rowStart[iRow];
     3501      const int * thisIndices = column+rowStart[iRow];
     3502      if (rowScale_) {
     3503        // scale row
     3504        double scale = rowScale_[iRow];
     3505        for (i = 0;i<rowLength[iRow];i++) {
     3506          int iColumn = thisIndices[i];
     3507          double alpha = thisElements[i];
     3508          scaled[i] = scale*alpha*columnScale_[iColumn];
     3509        }
     3510        thisElements = scaled;
     3511      }
     3512      for (i = 0;i<rowLength[iRow];i++) {
     3513        int iColumn = thisIndices[i];
     3514        double alpha = thisElements[i];
     3515        double oldValue = dj[iColumn];;
     3516        double value;
     3517
     3518        switch(getStatus(iColumn)) {
     3519         
     3520        case basic:
     3521          if (dj[iColumn]<-tolerance&&
     3522              fabs(solution_[iColumn]-upper_[iColumn])<1.0e-8) {
     3523            // at ub
     3524            changeUp += alpha*upper_[iColumn];
     3525            // might do other way
     3526            value = oldValue+djBasic*alpha;
     3527            if (value>tolerance)
     3528              canMove=false;
     3529          } else if (dj[iColumn]>tolerance&&
     3530              fabs(solution_[iColumn]-lower_[iColumn])<1.0e-8) {
     3531            changeUp += alpha*lower_[iColumn];
     3532            // might do other way
     3533            value = oldValue+djBasic*alpha;
     3534            if (value<-tolerance)
     3535              canMove=false;
     3536          } else {
     3537            canMove=false;
     3538          }
     3539          break;
     3540        case ClpSimplex::isFixed:
     3541          changeUp += alpha*upper_[iColumn];
     3542          break;
     3543        case isFree:
     3544        case superBasic:
     3545          canMove=false;
     3546        break;
     3547      case atUpperBound:
     3548        changeUp += alpha*upper_[iColumn];
     3549        // might do other way
     3550        value = oldValue+djBasic*alpha;
     3551        if (value>tolerance)
     3552          canMove=false;
     3553        break;
     3554      case atLowerBound:
     3555        changeUp += alpha*lower_[iColumn];
     3556        // might do other way
     3557        value = oldValue+djBasic*alpha;
     3558        if (value<-tolerance)
     3559          canMove=false;
     3560        break;
     3561        }
     3562      }
     3563      if (canMove) {
     3564        if (changeUp*djBasic>1.0e-12||fabs(changeUp)<1.0e-8) {
     3565          // move
     3566          for (i = 0;i<rowLength[iRow];i++) {
     3567            int iColumn = thisIndices[i];
     3568            double alpha = thisElements[i];
     3569            dj[iColumn] += djBasic * alpha;
     3570          }
     3571          dj[iSequence]=0.0;
     3572#ifdef CLP_DEBUG
     3573          {
     3574            double value5=0.0;
     3575            int i;
     3576            for (i=0;i<numberRows_+numberColumns_;i++) {
     3577              if (dj[i]<-1.0e-6)
     3578                value5 += dj[i]*upper_[i];
     3579              else if (dj[i] >1.0e-6)
     3580                value5 += dj[i]*lower_[i];
     3581            }
     3582            printf("Values objective Value after row %d old dj %g %g\n",
     3583                   iRow,djBasic,value5);
     3584          }
     3585#endif
     3586        }
     3587      }
     3588    }     
     3589  }
     3590  delete [] scaled;
     3591}
  • trunk/ClpSimplexPrimal.cpp

    r132 r137  
    254254    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
    255255    // and get good feasible duals
    256     computeDuals();
     256    computeDuals(NULL);
    257257  }
    258258  // clean up
     
    308308      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
    309309      createRim(7);
    310       gutsOfSolution(rowActivityWork_,columnActivityWork_);
     310      gutsOfSolution(rowActivityWork_,columnActivityWork_,NULL,NULL);
    311311      printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
    312312             numberIterations_,
     
    336336      if (sequenceIn<0) {
    337337        // end of values pass - initialize weights etc
     338        handler_->message(CLP_END_VALUES_PASS,messages_)
     339          <<numberIterations_;
    338340        primalColumnPivot_->saveWeights(this,5);
    339341        problemStatus_=-2; // factorize now
     
    444446  // put back original costs and then check
    445447  createRim(4);
    446   gutsOfSolution(rowActivityWork_, columnActivityWork_);
     448  gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    447449  // Check if looping
    448450  int loop = progress.looping();
     
    452454  } else if (loop<-1) {
    453455    // something may have changed
    454     gutsOfSolution(rowActivityWork_, columnActivityWork_);
     456    gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    455457  }
    456458  // really for free variables in
     
    495497      createRim(4);
    496498      nonLinearCost_->checkInfeasibilities(true);
    497       gutsOfSolution(rowActivityWork_, columnActivityWork_);
     499      gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    498500
    499501      infeasibilityCost_=1.0e100;
     
    512514        for (i=0;i<numberRows_+numberColumns_;i++)
    513515          cost_[i] *= 1.0e-100;
    514         gutsOfSolution(rowActivityWork_, columnActivityWork_);
     516        gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    515517        nonLinearCost_=nonLinear;
    516518        infeasibilityCost_=saveWeight;
     
    530532          createRim(4);
    531533          nonLinearCost_->checkInfeasibilities(true);
    532           gutsOfSolution(rowActivityWork_, columnActivityWork_);
     534          gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    533535          // so will exit
    534536          infeasibilityCost_=1.0e30;
     
    544546          createRim(4);
    545547          nonLinearCost_->checkInfeasibilities(true);
    546           gutsOfSolution(rowActivityWork_, columnActivityWork_);
     548          gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    547549          problemStatus_=-1; //continue
    548550        } else {
     
    576578          createRim(4);
    577579          nonLinearCost_->checkInfeasibilities(true);
    578           gutsOfSolution(rowActivityWork_, columnActivityWork_);
     580          gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    579581          problemStatus_ = -1;
    580582        } else {
     
    607609          // put back original costs and then check
    608610          createRim(4);
    609           gutsOfSolution(rowActivityWork_, columnActivityWork_);
     611          gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    610612          problemStatus_=-1; //continue
    611613        } else {
     
    12171219  // move non basic variables to new bounds
    12181220  nonLinearCost_->checkInfeasibilities(true);
    1219   gutsOfSolution(rowActivityWork_, columnActivityWork_);
     1221  gutsOfSolution(rowActivityWork_, columnActivityWork_,NULL,NULL);
    12201222}
    12211223// Unflag all variables and return number unflagged
  • trunk/Makefile.Clp

    r131 r137  
    55# highest level of optimization the compiler supports. If want something in
    66# between then specify the exact level you want, e.g., -O1 or -O2
    7 #OptLevel := -g
     7OptLevel := -g
    88OptLevel := -O3
    99
  • trunk/include/ClpMessage.hpp

    r102 r137  
    7272  CLP_EMPTY_PROBLEM,
    7373  CLP_CRASH,
     74  CLP_END_VALUES_PASS,
    7475  CLP_DUMMY_END
    7576};
  • trunk/include/ClpSimplex.hpp

    r136 r137  
    5858    atLowerBound = 0x03,
    5959    superBasic = 0x04,
     60    isFixed = 0x05
    6061  };
    6162
     
    134135  //@{
    135136  /** Dual algorithm - see ClpSimplexDual.hpp for method */
    136   int dual();
     137  int dual(int ifValuesPass=0);
    137138  /** Primal algorithm - see ClpSimplexPrimal.hpp for method */
    138139  int primal(int ifValuesPass=0);
     
    323324  /// Factorizes using current basis. For external use
    324325  int factorize();
    325   /// Computes duals from scratch
    326   void computeDuals();
     326  /** Computes duals from scratch. If givenDjs then
     327      allows for nonzero basic djs */
     328  void computeDuals(double * givenDjs);
    327329  /// Computes primals from scratch
    328330  void computePrimals (  const double * rowActivities,
     
    447449  /**@name protected methods */
    448450  //@{
    449   /// May change basis and then returns number changed
     451  /** May change basis and then returns number changed.
     452      Computation of solutions may be overriden by given pi and solution
     453  */
    450454  int gutsOfSolution ( const double * rowActivities,
    451455                       const double * columnActivities,
     456                       double * givenDuals,
     457                       const double * givenPrimals,
    452458                       bool valuesPass=false);
    453459  /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
     
    601607  inline Status getColumnStatus(int sequence) const
    602608  {return static_cast<Status> (status_[sequence]&7);};
    603   inline void setFixed( int sequence)
     609  inline void setPivoted( int sequence)
    604610  { status_[sequence] |= 32;};
    605   inline void clearFixed( int sequence)
     611  inline void clearPivoted( int sequence)
    606612  { status_[sequence] &= ~32; };
    607   inline bool fixed(int sequence) const
     613  inline bool pivoted(int sequence) const
    608614  {return (((status_[sequence]>>5)&1)!=0);};
    609615  inline void setFlagged( int sequence)
  • trunk/include/ClpSimplexDual.hpp

    r109 r137  
    114114  */
    115115
    116   int dual();
     116  int dual(int ifValuesPass);
    117117  /** For strong branching.  On input lower and upper are new bounds
    118118      while on output they are change in objective function values
     
    139139      +1 looks infeasible
    140140      +3 max iterations
     141
     142      If givenPi not NULL then in values pass
    141143   */
    142   int whileIterating();
     144  int whileIterating(double * & givenPi);
    143145  /** The duals are updated by the given arrays.
    144146      Returns number of infeasibilities.
     
    169171      For speed, we may need to go to a bucket approach when many
    170172      variables are being flipped
     173
    171174  */
    172175  void dualColumn(CoinIndexedVector * rowArray,
    173176                  CoinIndexedVector * columnArray,
    174177                  CoinIndexedVector * spareArray,
    175                   CoinIndexedVector * spareArray2);
     178                  CoinIndexedVector * spareArray2,
     179                  double accpetablePivot);
     180  /**
     181      Row array has row part of pivot row
     182      Column array has column part.
     183      This sees what is best thing to do in dual values pass
     184      Returns 0 if theta_ move will put basic variable out to bound,
     185      1 if can change dual on chosen row and leave variable in basis
     186  */
     187  int checkPossibleValuesMove(CoinIndexedVector * rowArray,
     188                               CoinIndexedVector * columnArray,
     189                               double acceptablePivot);
     190  /**
     191      This sees if we can move duals in dual values pass.
     192      This is done before any pivoting
     193  */
     194  void doEasyOnesInValuesPass(double * givenReducedCosts);
    176195  /**
    177196      Chooses dual pivot row
     
    179198      and will have this (with square of infeasibility) when steepest
    180199      For easy problems we can just choose one of the first rows we look at
    181   */
    182   void dualRow();
     200     
     201      If alreadyChosen >=0 then in values pass and that row has been
     202      selected
     203  */
     204  void dualRow(int alreadyChosen);
    183205  /** Checks if any fake bounds active - if so returns number and modifies
    184206      updatedDualBound_ and everything.
     
    210232  */
    211233  void statusOfProblemInDual(int & lastCleaned, int type,
    212                              ClpSimplexProgress & progress);
     234                             ClpSimplexProgress & progress,
     235                             double * givenDjs);
    213236  /// Perturbs problem (method depends on perturbation())
    214237  void perturb();
Note: See TracChangeset for help on using the changeset viewer.