Changeset 180 for branches


Ignore:
Timestamp:
Jul 18, 2003 3:17:33 PM (16 years ago)
Author:
forrest
Message:

new stuff

Location:
branches/pre
Files:
49 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpDualRowDantzig.cpp

    r63 r180  
    7676  return chosenRow;
    7777}
     78// Returns pivot alpha
     79double
     80ClpDualRowDantzig::updateWeights(CoinIndexedVector * input,
     81                                  CoinIndexedVector * spare,
     82                                  CoinIndexedVector * updatedColumn)
     83{
     84  // pivot element
     85  double alpha=0.0;
     86  // look at updated column
     87  double * work = updatedColumn->denseVector();
     88  int number = updatedColumn->getNumElements();
     89  int * which = updatedColumn->getIndices();
     90  int i;
     91  int pivotRow = model_->pivotRow();
     92
     93  if (updatedColumn->packedMode()) {
     94    for (i =0; i < number; i++) {
     95      int iRow = which[i];
     96      if (iRow==pivotRow) {
     97        alpha = work[i];
     98        break;
     99      }
     100    }
     101  } else {
     102    alpha = work[pivotRow];
     103  }
     104  return alpha;
     105}
    78106 
    79107/* Updates primal solution (and maybe list of candidates)
     
    92120  double changeObj=0.0;
    93121  const int * pivotVariable = model_->pivotVariable();
    94   for (i=0;i<number;i++) {
    95     int iRow=which[i];
    96     int iPivot=pivotVariable[iRow];
    97     double & value = model_->solutionAddress(iPivot);
    98     double cost = model_->cost(iPivot);
    99     double change = primalRatio*work[iRow];
    100     value -= change;
    101     changeObj -= change*cost;
    102     work[iRow]=0.0;
     122  if (primalUpdate->packedMode()) {
     123    for (i=0;i<number;i++) {
     124      int iRow=which[i];
     125      int iPivot=pivotVariable[iRow];
     126      double & value = model_->solutionAddress(iPivot);
     127      double cost = model_->cost(iPivot);
     128      double change = primalRatio*work[i];
     129      value -= change;
     130      changeObj -= change*cost;
     131      work[i]=0.0;
     132    }
     133  } else {
     134    for (i=0;i<number;i++) {
     135      int iRow=which[i];
     136      int iPivot=pivotVariable[iRow];
     137      double & value = model_->solutionAddress(iPivot);
     138      double cost = model_->cost(iPivot);
     139      double change = primalRatio*work[iRow];
     140      value -= change;
     141      changeObj -= change*cost;
     142      work[iRow]=0.0;
     143    }
    103144  }
    104145  primalUpdate->setNumElements(0);
  • branches/pre/ClpDualRowPivot.cpp

    r63 r180  
    6060{
    6161}
    62 
    63 void
    64 ClpDualRowPivot::updateWeights(CoinIndexedVector * input,
    65                                CoinIndexedVector * spare,
    66                                CoinIndexedVector * updatedColumn)
    67 {
    68 }
    6962void
    7063ClpDualRowPivot::unrollWeights()
  • branches/pre/ClpDualRowSteepest.cpp

    r137 r180  
    135135  // allow tolerance at least slightly bigger than standard
    136136  tolerance = tolerance +  error;
     137  // But cap
     138  tolerance = min(1000.0,tolerance);
    137139  tolerance *= tolerance; // as we are using squares
    138140  double * solution = model_->solutionRegion();
     
    156158      if (iPivot<numberColumns)
    157159        value *= COLUMN_BIAS; // bias towards columns
     160k
    158161#endif
    159162      // store square in list
     
    181184    number = infeasible_->getNumElements();
    182185  }
    183   for (i=0;i<number;i++) {
    184     iRow = index[i];
    185     double value = infeas[iRow];
    186     if (value>largest*weights_[iRow]&&value>tolerance) {
    187       // make last pivot row last resort choice
    188       if (iRow==lastPivotRow) {
    189         if (value*1.0e-10<largest*weights_[iRow])
    190           continue;
    191         else
    192           value *= 1.0e-10;
    193       }
    194       int iSequence = pivotVariable[iRow];
    195       if (!model_->flagged(iSequence)) {
    196         //#define CLP_DEBUG 1
     186  if(model_->numberIterations()<model_->lastBadIteration()+200) {
     187    // we can't really trust infeasibilities if there is dual error
     188    double checkTolerance = 1.0e-8;
     189    if (!model_->factorization()->pivots())
     190      checkTolerance = 1.0e-6;
     191    if (model_->largestPrimalError()>checkTolerance)
     192      tolerance *= model_->largestPrimalError()/checkTolerance;
     193  }
     194  int numberWanted;
     195  if (mode_<2)
     196    numberWanted = number+1;
     197  else
     198    numberWanted = max(2000,number/8);
     199  int iPass;
     200  // Setup two passes
     201  int start[4];
     202  start[1]=number;
     203  start[2]=0;
     204  double dstart = ((double) number) * CoinDrand48();
     205  start[0]=(int) dstart;
     206  start[3]=start[0];
     207  for (iPass=0;iPass<2;iPass++) {
     208    int end = start[2*iPass+1];
     209    for (i=start[2*iPass];i<end;i++) {
     210      iRow = index[i];
     211      double value = infeas[iRow];
     212      if (value>tolerance) {
     213        if (value>largest*weights_[iRow]) {
     214          // make last pivot row last resort choice
     215          if (iRow==lastPivotRow) {
     216            if (value*1.0e-10<largest*weights_[iRow])
     217              continue;
     218            else
     219              value *= 1.0e-10;
     220          }
     221          int iSequence = pivotVariable[iRow];
     222          if (!model_->flagged(iSequence)) {
     223            //#define CLP_DEBUG 1
    197224#ifdef CLP_DEBUG
    198         double value2=0.0;
    199         if (solution[iSequence]>upper[iSequence]+tolerance)
    200           value2=solution[iSequence]-upper[iSequence];
    201         else if (solution[iSequence]<lower[iSequence]-tolerance)
    202           value2=solution[iSequence]-lower[iSequence];
    203         assert(fabs(value2*value2-infeas[iRow])<1.0e-8*min(value2*value2,infeas[iRow]));
    204 #endif
    205         if (solution[iSequence]>upper[iSequence]+tolerance||
    206             solution[iSequence]<lower[iSequence]-tolerance) {
    207           chosenRow=iRow;
    208           largest=value/weights_[iRow];
     225            double value2=0.0;
     226            if (solution[iSequence]>upper[iSequence]+tolerance)
     227              value2=solution[iSequence]-upper[iSequence];
     228            else if (solution[iSequence]<lower[iSequence]-tolerance)
     229              value2=solution[iSequence]-lower[iSequence];
     230            assert(fabs(value2*value2-infeas[iRow])<1.0e-8*min(value2*value2,infeas[iRow]));
     231#endif
     232            if (solution[iSequence]>upper[iSequence]+tolerance||
     233                solution[iSequence]<lower[iSequence]-tolerance) {
     234              chosenRow=iRow;
     235              largest=value/weights_[iRow];
     236            }
     237          }
    209238        }
    210       }
    211     }
     239        numberWanted--;
     240        if (!numberWanted)
     241          break;
     242      }
     243    }
     244    if (!numberWanted)
     245      break;
    212246  }
    213247  return chosenRow;
    214248}
    215249#define TRY_NORM 1.0e-4
    216 // Updates weights
    217 void
     250// Updates weights and returns pivot alpha
     251double
    218252ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
    219253                                  CoinIndexedVector * spare,
    220254                                  CoinIndexedVector * updatedColumn)
    221255{
    222   // clear other region
    223   alternateWeights_->clear();
    224   double norm = 0.0;
    225   double * work = input->denseVector();
    226   int number = input->getNumElements();
    227   int * which = input->getIndices();
    228   double * work2 = spare->denseVector();
    229   int * which2 = spare->getIndices();
    230   double * work3 = alternateWeights_->denseVector();
    231   int * which3 = alternateWeights_->getIndices();
    232   int i;
    233256#if CLP_DEBUG>2
    234257  // Very expensive debug
     
    261284  }
    262285#endif
    263   for (i=0;i<number;i++) {
    264     int iRow = which[i];
    265     double value = work[iRow];
    266     norm += value*value;
    267     work2[iRow]=value;
    268     which2[i]=iRow;
    269   }
    270   spare->setNumElements(number);
    271   // ftran
    272   model_->factorization()->updateColumn(alternateWeights_,spare);
    273   // alternateWeights_ should still be empty
    274   int pivotRow = model_->pivotRow();
     286  assert (input->packedMode());
     287  assert (updatedColumn->packedMode());
     288  double alpha=0.0;
     289  if (!model_->factorization()->networkBasis()) {
     290    // clear other region
     291    alternateWeights_->clear();
     292    double norm = 0.0;
     293    int i;
     294    double * work = input->denseVector();
     295    int numberNonZero = input->getNumElements();
     296    int * which = input->getIndices();
     297    double * work2 = spare->denseVector();
     298    int * which2 = spare->getIndices();
     299    // ftran
     300    //permute and move indices into index array
     301    //also compute norm
     302    //int *regionIndex = alternateWeights_->getIndices (  );
     303    const int *permute = model_->factorization()->permute();
     304    //double * region = alternateWeights_->denseVector();
     305    for ( i = 0; i < numberNonZero; i ++ ) {
     306      int iRow = which[i];
     307      double value = work[i];
     308      norm += value*value;
     309      iRow = permute[iRow];
     310      work2[iRow] = value;
     311      which2[i] = iRow;
     312    }
     313    spare->setNumElements ( numberNonZero );
     314    model_->factorization()->updateColumn(spare,spare,true);
     315    // permute back
     316    numberNonZero = spare->getNumElements();
     317    // alternateWeights_ should still be empty
     318    int pivotRow = model_->pivotRow();
     319    // could re-initialize here (could be expensive)
     320    norm /= model_->alpha() * model_->alpha();
     321   
     322    assert(norm);
     323    // pivot element
     324    alpha=0.0;
     325    double multiplier = 2.0 / model_->alpha();
     326    // look at updated column
     327    work = updatedColumn->denseVector();
     328    numberNonZero = updatedColumn->getNumElements();
     329    which = updatedColumn->getIndices();
     330   
     331    int nSave=0;
     332    double * work3 = alternateWeights_->denseVector();
     333    int * which3 = alternateWeights_->getIndices();
     334    const int * pivotColumn = model_->factorization()->pivotColumn();
     335    for (i =0; i < numberNonZero; i++) {
     336      int iRow = which[i];
     337      double theta = work[i];
     338      if (iRow==pivotRow)
     339        alpha = theta;
     340      double devex = weights_[iRow];
     341      work3[nSave]=devex; // save old
     342      which3[nSave++]=iRow;
     343      // transform to match spare
     344      int jRow = pivotColumn[iRow];
     345      double value = work2[jRow];
     346      devex +=  theta * (theta*norm+value * multiplier);
     347      if (devex < TRY_NORM)
     348        devex = TRY_NORM;
     349      weights_[iRow]=devex;
     350    }
     351    assert (alpha);
     352    alternateWeights_->setPackedMode(true);
     353    alternateWeights_->setNumElements(nSave);
     354    if (norm < TRY_NORM)
     355      norm = TRY_NORM;
     356    weights_[pivotRow] = norm;
     357    spare->clear();
    275358#ifdef CLP_DEBUG
    276   if ( model_->logLevel (  ) >4  &&
    277        fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm))
    278     printf("on row %d, true weight %g, old %g\n",
    279            pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
    280 #endif
    281   // could re-initialize here (could be expensive)
    282   norm /= model_->alpha() * model_->alpha();
    283 
    284   assert(norm);
    285   if (norm < TRY_NORM)
    286     norm = TRY_NORM;
    287   if (norm != 0.) {
     359    spare->checkClear();
     360#endif
     361  } else {
     362    // clear other region
     363    alternateWeights_->clear();
     364    double norm = 0.0;
     365    int i;
     366    double * work = input->denseVector();
     367    int number = input->getNumElements();
     368    int * which = input->getIndices();
     369    double * work2 = spare->denseVector();
     370    int * which2 = spare->getIndices();
     371    for (i=0;i<number;i++) {
     372      int iRow = which[i];
     373      double value = work[i];
     374      norm += value*value;
     375      work2[iRow]=value;
     376      which2[i]=iRow;
     377    }
     378    spare->setNumElements(number);
     379    // ftran
     380    model_->factorization()->updateColumn(alternateWeights_,spare);
     381    // alternateWeights_ should still be empty
     382    int pivotRow = model_->pivotRow();
     383#ifdef CLP_DEBUG
     384    if ( model_->logLevel (  ) >4  &&
     385         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm))
     386      printf("on row %d, true weight %g, old %g\n",
     387             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
     388#endif
     389    // could re-initialize here (could be expensive)
     390    norm /= model_->alpha() * model_->alpha();
     391   
     392    assert(norm);
     393    //if (norm < TRY_NORM)
     394    //norm = TRY_NORM;
     395    // pivot element
     396    alpha=0.0;
    288397    double multiplier = 2.0 / model_->alpha();
    289398    // look at updated column
     
    291400    number = updatedColumn->getNumElements();
    292401    which = updatedColumn->getIndices();
    293 
     402   
    294403    int nSave=0;
    295 
     404    double * work3 = alternateWeights_->denseVector();
     405    int * which3 = alternateWeights_->getIndices();
    296406    for (i =0; i < number; i++) {
    297407      int iRow = which[i];
    298       double theta = work[iRow];
    299       if (theta) {
    300         double devex = weights_[iRow];
    301         work3[iRow]=devex; // save old
    302         which3[nSave++]=iRow;
    303         double value = work2[iRow];
    304         devex +=  theta * (theta*norm+value * multiplier);
    305         if (devex < TRY_NORM)
    306           devex = TRY_NORM;
    307         weights_[iRow]=devex;
    308       }
    309     }
    310 #ifdef CLP_DEBUG
    311     assert(work3[pivotRow]&&work[pivotRow]);
    312 #endif
     408      double theta = work[i];
     409      if (iRow==pivotRow)
     410        alpha = theta;
     411      double devex = weights_[iRow];
     412      work3[nSave]=devex; // save old
     413      which3[nSave++]=iRow;
     414      double value = work2[iRow];
     415      devex +=  theta * (theta*norm+value * multiplier);
     416      if (devex < TRY_NORM)
     417        devex = TRY_NORM;
     418      weights_[iRow]=devex;
     419    }
     420    assert (alpha);
     421    alternateWeights_->setPackedMode(true);
    313422    alternateWeights_->setNumElements(nSave);
    314423    if (norm < TRY_NORM)
    315424      norm = TRY_NORM;
    316425    weights_[pivotRow] = norm;
    317   }
    318   spare->clear();
     426    spare->clear();
     427  }
    319428#ifdef CLP_DEBUG
    320429  spare->checkClear();
    321430#endif
     431  return alpha;
    322432}
    323433 
     
    345455  int numberColumns = model_->numberColumns();
    346456#endif
    347   for (i=0;i<number;i++) {
    348     int iRow=which[i];
    349     int iPivot=pivotVariable[iRow];
    350     double value = solution[iPivot];
    351     double cost = model_->cost(iPivot);
    352     double change = primalRatio*work[iRow];
    353     value -= change;
    354     changeObj -= change*cost;
    355     solution[iPivot] = value;
    356     double lower = model_->lower(iPivot);
    357     double upper = model_->upper(iPivot);
    358     // But if pivot row then use value of incoming
    359     // Although it is safer to recompute before next selection
    360     // in case something odd happens
    361     if (iRow==pivotRow) {
    362       iPivot = model_->sequenceIn();
    363       lower = model_->lower(iPivot);
    364       upper = model_->upper(iPivot);
    365       value = model_->valueIncomingDual();
    366     }
    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
     457  if (primalUpdate->packedMode()) {
     458    for (i=0;i<number;i++) {
     459      int iRow=which[i];
     460      int iPivot=pivotVariable[iRow];
     461      double value = solution[iPivot];
     462      double cost = model_->cost(iPivot);
     463      double change = primalRatio*work[i];
     464      work[i]=0.0;
     465      value -= change;
     466      changeObj -= change*cost;
     467      solution[iPivot] = value;
     468      double lower = model_->lower(iPivot);
     469      double upper = model_->upper(iPivot);
     470      // But if pivot row then use value of incoming
     471      // Although it is safer to recompute before next selection
     472      // in case something odd happens
     473      if (iRow==pivotRow) {
     474        iPivot = model_->sequenceIn();
     475        lower = model_->lower(iPivot);
     476        upper = model_->upper(iPivot);
     477        value = model_->valueIncomingDual();
     478      }
     479      if (value<lower-tolerance) {
     480        value -= lower;
     481        value *= value;
     482#ifdef COLUMN_BIAS
     483        if (iPivot<numberColumns)
     484          value *= COLUMN_BIAS; // bias towards columns
    373485#endif
    374486#ifdef FIXED_BIAS
    375       if (lower==upper)
    376         value *= FIXED_BIAS; // bias towards taking out fixed variables
    377 #endif
    378       // store square in list
    379       if (infeas[iRow])
    380         infeas[iRow] = value; // already there
    381       else
    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
     487        if (lower==upper)
     488          value *= FIXED_BIAS; // bias towards taking out fixed variables
     489#endif
     490        // store square in list
     491        if (infeas[iRow])
     492          infeas[iRow] = value; // already there
     493        else
     494          infeasible_->quickAdd(iRow,value);
     495      } else if (value>upper+tolerance) {
     496        value -= upper;
     497        value *= value;
     498#ifdef COLUMN_BIAS
     499        if (iPivot<numberColumns)
     500          value *= COLUMN_BIAS; // bias towards columns
    389501#endif
    390502#ifdef FIXED_BIAS
    391       if (lower==upper)
    392         value *= FIXED_BIAS; // bias towards taking out fixed variables
    393 #endif
    394       // store square in list
    395       if (infeas[iRow])
    396         infeas[iRow] = value; // already there
    397       else
    398         infeasible_->quickAdd(iRow,value);
    399     } else {
    400       // feasible - was it infeasible - if so set tiny
    401       if (infeas[iRow])
    402         infeas[iRow] = 1.0e-100;
    403     }
    404     work[iRow]=0.0;
     503        if (lower==upper)
     504          value *= FIXED_BIAS; // bias towards taking out fixed variables
     505#endif
     506        // store square in list
     507        if (infeas[iRow])
     508          infeas[iRow] = value; // already there
     509        else
     510          infeasible_->quickAdd(iRow,value);
     511      } else {
     512        // feasible - was it infeasible - if so set tiny
     513        if (infeas[iRow])
     514          infeas[iRow] = 1.0e-100;
     515      }
     516    }
     517  } else {
     518    for (i=0;i<number;i++) {
     519      int iRow=which[i];
     520      int iPivot=pivotVariable[iRow];
     521      double value = solution[iPivot];
     522      double cost = model_->cost(iPivot);
     523      double change = primalRatio*work[iRow];
     524      value -= change;
     525      changeObj -= change*cost;
     526      solution[iPivot] = value;
     527      double lower = model_->lower(iPivot);
     528      double upper = model_->upper(iPivot);
     529      // But if pivot row then use value of incoming
     530      // Although it is safer to recompute before next selection
     531      // in case something odd happens
     532      if (iRow==pivotRow) {
     533        iPivot = model_->sequenceIn();
     534        lower = model_->lower(iPivot);
     535        upper = model_->upper(iPivot);
     536        value = model_->valueIncomingDual();
     537      }
     538      if (value<lower-tolerance) {
     539        value -= lower;
     540        value *= value;
     541#ifdef COLUMN_BIAS
     542        if (iPivot<numberColumns)
     543          value *= COLUMN_BIAS; // bias towards columns
     544#endif
     545#ifdef FIXED_BIAS
     546        if (lower==upper)
     547          value *= FIXED_BIAS; // bias towards taking out fixed variables
     548#endif
     549        // store square in list
     550        if (infeas[iRow])
     551          infeas[iRow] = value; // already there
     552        else
     553          infeasible_->quickAdd(iRow,value);
     554      } else if (value>upper+tolerance) {
     555        value -= upper;
     556        value *= value;
     557#ifdef COLUMN_BIAS
     558        if (iPivot<numberColumns)
     559          value *= COLUMN_BIAS; // bias towards columns
     560#endif
     561#ifdef FIXED_BIAS
     562        if (lower==upper)
     563          value *= FIXED_BIAS; // bias towards taking out fixed variables
     564#endif
     565        // store square in list
     566        if (infeas[iRow])
     567          infeas[iRow] = value; // already there
     568        else
     569          infeasible_->quickAdd(iRow,value);
     570      } else {
     571        // feasible - was it infeasible - if so set tiny
     572        if (infeas[iRow])
     573          infeas[iRow] = 1.0e-100;
     574      }
     575      work[iRow]=0.0;
     576    }
    405577  }
    406578  primalUpdate->setNumElements(0);
     
    442614      alternateWeights_->reserve(numberRows+
    443615                                 model_->factorization()->maximumPivots());
    444       if (!mode_||mode==5) {
     616      if (!mode_||mode_==2||mode==5) {
    445617        // initialize to 1.0 (can we do better?)
    446618        for (i=0;i<numberRows;i++) {
     
    455627        for (i=0;i<numberRows;i++) {
    456628          double value=0.0;
    457           array[i] = 1.0;
     629          array[0] = 1.0;
    458630          which[0] = i;
    459631          alternateWeights_->setNumElements(1);
     632          alternateWeights_->setPackedMode(true);
    460633          model_->factorization()->updateColumnTranspose(temp,
    461634                                                         alternateWeights_);
     
    463636          int j;
    464637          for (j=0;j<number;j++) {
    465             int iRow=which[j];
    466             value+=array[iRow]*array[iRow];
    467             array[iRow]=0.0;
     638            value+=array[j]*array[j];
     639            array[j]=0.0;
    468640          }
    469641          alternateWeights_->setNumElements(0);
     
    566738  int * which = alternateWeights_->getIndices();
    567739  int i;
    568   for (i=0;i<number;i++) {
    569     int iRow = which[i];
    570     weights_[iRow]=saved[iRow];
    571     saved[iRow]=0.0;
     740  if (alternateWeights_->packedMode()) {
     741    for (i=0;i<number;i++) {
     742      int iRow = which[i];
     743      weights_[iRow]=saved[i];
     744      saved[i]=0.0;
     745    }
     746  } else {
     747    for (i=0;i<number;i++) {
     748      int iRow = which[i];
     749      weights_[iRow]=saved[iRow];
     750      saved[iRow]=0.0;
     751    }
    572752  }
    573753  alternateWeights_->setNumElements(0);
     
    598778  state_ =-1;
    599779}
     780// Returns true if would not find any row
     781bool
     782ClpDualRowSteepest::looksOptimal() const
     783{
     784  int iRow;
     785  const int * pivotVariable = model_->pivotVariable();
     786  double tolerance=model_->currentPrimalTolerance();
     787  // we can't really trust infeasibilities if there is primal error
     788  // this coding has to mimic coding in checkPrimalSolution
     789  double error = min(1.0e-3,model_->largestPrimalError());
     790  // allow tolerance at least slightly bigger than standard
     791  tolerance = tolerance +  error;
     792  // But cap
     793  tolerance = min(1000.0,tolerance);
     794  int numberRows = model_->numberRows();
     795  int numberInfeasible=0;
     796  for (iRow=0;iRow<numberRows;iRow++) {
     797    int iPivot=pivotVariable[iRow];
     798    double value = model_->solution(iPivot);
     799    double lower = model_->lower(iPivot);
     800    double upper = model_->upper(iPivot);
     801    if (value<lower-tolerance) {
     802      numberInfeasible++;
     803    } else if (value>upper+tolerance) {
     804      numberInfeasible++;
     805    }
     806  }
     807  return (numberInfeasible==0);
     808}
    600809
  • branches/pre/ClpFactorization.cpp

    r138 r180  
    7373}
    7474int
    75 ClpFactorization::factorize ( const ClpSimplex * model,
    76                               const ClpMatrixBase * matrix,
    77                               int numberRows, int numberColumns,
    78                               int rowIsBasic[], int columnIsBasic[] ,
    79                               double areaFactor )
    80 {
     75ClpFactorization::factorize ( ClpSimplex * model,
     76                              int solveType, bool valuesPass)
     77{
     78  const ClpMatrixBase * matrix = model->clpMatrix();
     79  int numberRows = model->numberRows();
     80  int numberColumns = model->numberColumns();
     81  // If too many compressions increase area
     82  if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
     83    areaFactor_ *= 1.1;
     84  }
    8185  if (!networkBasis_||doCheck) {
    82     // see if matrix a network
    83     ClpNetworkMatrix* networkMatrix =
    84       dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
    85     // If network - still allow ordinary factorization first time for laziness
    86     int saveMaximumPivots = maximumPivots();
    87     delete networkBasis_;
    88     networkBasis_ = NULL;
    89     if (networkMatrix&&!doCheck)
    90       maximumPivots(1);
    91     // maybe for speed will be better to leave as many regions as possible
    92     gutsOfDestructor();
    93     gutsOfInitialize(2);
    94     if (areaFactor)
    95       areaFactor_ = areaFactor;
    96     int numberBasic = 0;
    97     CoinBigIndex numberElements=0;
    98     int numberRowBasic=0;
    99    
    100     // compute how much in basis
    101    
    102     int i;
    103    
    104     for (i=0;i<numberRows;i++) {
    105       if (rowIsBasic[i]>=0)
    106         numberRowBasic++;
    107     }
    108    
    109     numberBasic = numberRowBasic;
    110     for (i=0;i<numberColumns;i++) {
    111       if (columnIsBasic[i]>=0) {
    112         numberBasic++;
     86    status_=-99;
     87    int * pivotVariable = model->pivotVariable();
     88    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
     89    while (status_<-98) {
     90     
     91      int i;
     92      int numberBasic=0;
     93      int numberRowBasic;
     94      for (i=0;i<numberRows;i++) {
     95        if (model->getRowStatus(i) == ClpSimplex::basic)
     96          pivotVariable[numberBasic++]=i;
    11397      }
    114     }
    115     numberElements += matrix->numberInBasis(columnIsBasic);
    116     if ( numberBasic > numberRows ) {
    117     return -2; // say too many in basis
    118     }
    119     numberElements = 3 * numberBasic + 3 * numberElements + 10000;
    120     getAreas ( numberRows, numberBasic, numberElements,
    121                2 * numberElements );
    122     //fill
    123     //copy
    124     numberBasic=0;
    125     numberElements=0;
    126     for (i=0;i<numberRows;i++) {
    127       if (rowIsBasic[i]>=0) {
    128         indexRowU_[numberElements]=i;
    129         indexColumnU_[numberElements]=numberBasic;
    130         elementU_[numberElements++]=slackValue_;
    131         numberBasic++;
     98      numberRowBasic=numberBasic;
     99      for (i=0;i<numberColumns;i++) {
     100        if (model->getColumnStatus(i) == ClpSimplex::basic)
     101          pivotVariable[numberBasic++]=i;
    132102      }
    133     }
    134     numberElements +=matrix->fillBasis(model, columnIsBasic, numberBasic,
    135                                        indexRowU_+numberElements,
    136                                        indexColumnU_+numberElements,
    137                                        elementU_+numberElements);
    138     lengthU_ = numberElements;
    139    
    140     preProcess ( 0 );
    141     factor (  );
    142     numberBasic=0;
    143     if (status_ == 0) {
    144       int * permuteBack = permuteBack_;
    145       int * back = pivotColumnBack_;
    146       for (i=0;i<numberRows;i++) {
    147         if (rowIsBasic[i]>=0) {
    148           rowIsBasic[i]=permuteBack[back[numberBasic++]];
     103      assert (numberBasic<=numberRows);
     104      // see if matrix a network
     105      ClpNetworkMatrix* networkMatrix =
     106        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
     107      // If network - still allow ordinary factorization first time for laziness
     108      int saveMaximumPivots = maximumPivots();
     109      delete networkBasis_;
     110      networkBasis_ = NULL;
     111      if (networkMatrix&&!doCheck)
     112        maximumPivots(1);
     113      while (status_==-99) {
     114        // maybe for speed will be better to leave as many regions as possible
     115        gutsOfDestructor();
     116        gutsOfInitialize(2);
     117        CoinBigIndex numberElements=numberRowBasic;
     118
     119        // compute how much in basis
     120
     121        int i;
     122
     123        numberElements +=matrix->fillBasis(model,
     124                                           pivotVariable+numberRowBasic,
     125                                           numberRowBasic,
     126                                           numberBasic-numberRowBasic,
     127                                           NULL,NULL,NULL);
     128
     129        numberElements = 3 * numberBasic + 3 * numberElements + 10000;
     130        getAreas ( numberRows, numberBasic, numberElements,
     131                   2 * numberElements );
     132        //fill
     133        //copy
     134        numberElements=numberRowBasic;
     135        for (i=0;i<numberRowBasic;i++) {
     136          int iRow = pivotVariable[i];
     137          indexRowU_[i]=iRow;
     138          indexColumnU_[i]=i;
     139          elementU_[i]=slackValue_;
     140        }
     141        numberElements +=matrix->fillBasis(model,
     142                                           pivotVariable+numberRowBasic,
     143                                           numberRowBasic,
     144                                           numberBasic-numberRowBasic,
     145                                           indexRowU_+numberElements,
     146                                           indexColumnU_+numberElements,
     147                                           elementU_+numberElements);
     148        lengthU_ = numberElements;
     149
     150        preProcess ( 0 );
     151        factor (  );
     152        if (status_==-99) {
     153          // get more memory
     154          areaFactor(2.0*areaFactor());
    149155        }
    150156      }
    151       for (i=0;i<numberColumns;i++) {
    152         if (columnIsBasic[i]>=0) {
    153           columnIsBasic[i]=permuteBack[back[numberBasic++]];
    154         }
    155       }
    156       if (increasingRows_>1) {
     157      // If we get here status is 0 or -1
     158      if (status_ == 0) {
     159        int * permuteBack = permuteBack_;
     160        int * back = pivotColumnBack_;
     161        int * pivotTemp = pivotColumn_;
     162        ClpDisjointCopyN ( pivotVariable, numberRows_ , pivotTemp  );
     163        // Redo pivot order
     164        for (i=0;i<numberRowBasic;i++) {
     165          int k = pivotTemp[i];
     166          // so rowIsBasic[k] would be permuteBack[back[i]]
     167          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
     168        }
     169        for (;i<numberRows;i++) {
     170          int k = pivotTemp[i];
     171          // so rowIsBasic[k] would be permuteBack[back[i]]
     172          pivotVariable[permuteBack[back[i]]]=k;
     173        }
    157174        // Set up permutation vector
    158         if (increasingRows_<3) {
    159           // these arrays start off as copies of permute
    160           // (and we could use permute_ instead of pivotColumn (not back though))
    161           ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
    162           ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
    163         }
    164       } else {
    165         // Set up permutation vector
    166         // (we could use permute_ instead of pivotColumn (not back though))
    167         for (i=0;i<numberRows_;i++) {
    168           int k=pivotColumn_[i];
    169           pivotColumn_[i]=pivotColumnBack_[i];
    170           pivotColumnBack_[i]=k;
    171         }
    172       }
    173     } else if (status_ == -1) {
    174       // mark as basic or non basic
    175       for (i=0;i<numberRows_;i++) {
    176         if (rowIsBasic[i]>=0) {
    177           if (pivotColumn_[numberBasic]>=0)
    178             rowIsBasic[i]=pivotColumn_[numberBasic];
    179           else
    180             rowIsBasic[i]=-1;
    181           numberBasic++;
    182         }
    183       }
    184       for (i=0;i<numberColumns;i++) {
    185         if (columnIsBasic[i]>=0) {
    186           if (pivotColumn_[numberBasic]>=0)
    187             columnIsBasic[i]=pivotColumn_[numberBasic];
    188           else
    189             columnIsBasic[i]=-1;
    190           numberBasic++;
    191         }
    192       }
    193     }
    194     if (networkMatrix) {
    195       maximumPivots(saveMaximumPivots);
    196       if (!status_) {
    197         // create network factorization
    198         if (doCheck)
    199           delete networkBasis_; // temp
    200         networkBasis_ = new ClpNetworkBasis(model,numberRows_,
    201                                             pivotRegion_,
    202                                             permuteBack_,
    203                                             startColumnU_,
    204                                             numberInColumn_,
    205                                             indexRowU_,
    206                                             elementU_);
    207         // kill off arrays in ordinary factorization
    208         if (!doCheck)
    209           gutsOfDestructor();
    210       }
    211     }
    212     if (!status_&&!networkBasis_) {
    213       // See if worth going sparse and when
    214       if (numberFtranCounts_>100) {
    215         ftranAverageAfterL_ = max(ftranCountAfterL_/ftranCountInput_,1.0);
    216         ftranAverageAfterR_ = max(ftranCountAfterR_/ftranCountAfterL_,1.0);
    217         ftranAverageAfterU_ = max(ftranCountAfterU_/ftranCountAfterR_,1.0);
    218         if (btranCountInput_) {
    219           btranAverageAfterU_ = max(btranCountAfterU_/btranCountInput_,1.0);
    220           btranAverageAfterR_ = max(btranCountAfterR_/btranCountAfterU_,1.0);
    221           btranAverageAfterL_ = max(btranCountAfterL_/btranCountAfterR_,1.0);
     175        // these arrays start off as copies of permute
     176        // (and we could use permute_ instead of pivotColumn (not back though))
     177        ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
     178        ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
     179        if (networkMatrix) {
     180          maximumPivots(saveMaximumPivots);
     181          // create network factorization
     182          if (doCheck)
     183            delete networkBasis_; // temp
     184          networkBasis_ = new ClpNetworkBasis(model,numberRows_,
     185                                              pivotRegion_,
     186                                              permuteBack_,
     187                                              startColumnU_,
     188                                              numberInColumn_,
     189                                              indexRowU_,
     190                                              elementU_);
     191          // kill off arrays in ordinary factorization
     192          if (!doCheck) {
     193            gutsOfDestructor();
     194#if 0
     195            // but put back permute arrays so odd things will work
     196            int numberRows = model->numberRows();
     197            pivotColumnBack_ = new int [numberRows];
     198            permute_ = new int [numberRows];
     199            int i;
     200            for (i=0;i<numberRows;i++) {
     201              pivotColumnBack_[i]=i;
     202              permute_[i]=i;
     203            }
     204#endif
     205          }
    222206        } else {
    223           // odd - we have not done any btrans
    224           btranAverageAfterU_ = 1.0;
    225           btranAverageAfterR_ = 1.0;
    226           btranAverageAfterL_ = 1.0;
    227         }
    228       }
    229       // scale back
    230 
    231       ftranCountInput_ *= 0.8;
    232       ftranCountAfterL_ *= 0.8;
    233       ftranCountAfterR_ *= 0.8;
    234       ftranCountAfterU_ *= 0.8;
    235       btranCountInput_ *= 0.8;
    236       btranCountAfterU_ *= 0.8;
    237       btranCountAfterR_ *= 0.8;
    238       btranCountAfterL_ *= 0.8;
     207          // See if worth going sparse and when
     208          if (numberFtranCounts_>100) {
     209            ftranAverageAfterL_ = max(ftranCountAfterL_/ftranCountInput_,1.0);
     210            ftranAverageAfterR_ = max(ftranCountAfterR_/ftranCountAfterL_,1.0);
     211            ftranAverageAfterU_ = max(ftranCountAfterU_/ftranCountAfterR_,1.0);
     212            if (btranCountInput_) {
     213              btranAverageAfterU_ = max(btranCountAfterU_/btranCountInput_,1.0);
     214              btranAverageAfterR_ = max(btranCountAfterR_/btranCountAfterU_,1.0);
     215              btranAverageAfterL_ = max(btranCountAfterL_/btranCountAfterR_,1.0);
     216            } else {
     217              // odd - we have not done any btrans
     218              btranAverageAfterU_ = 1.0;
     219              btranAverageAfterR_ = 1.0;
     220              btranAverageAfterL_ = 1.0;
     221            }
     222          }
     223          // scale back
     224         
     225          ftranCountInput_ *= 0.8;
     226          ftranCountAfterL_ *= 0.8;
     227          ftranCountAfterR_ *= 0.8;
     228          ftranCountAfterU_ *= 0.8;
     229          btranCountInput_ *= 0.8;
     230          btranCountAfterU_ *= 0.8;
     231          btranCountAfterR_ *= 0.8;
     232          btranCountAfterL_ *= 0.8;
     233        }
     234      } else if (status_ == -1&&(solveType==0||solveType==2)) {
     235        // This needs redoing as it was merged coding - does not need array
     236        int numberTotal = numberRows+numberColumns;
     237        int * isBasic = new int [numberTotal];
     238        int * rowIsBasic = isBasic+numberColumns;
     239        int * columnIsBasic = isBasic;
     240        for (i=0;i<numberTotal;i++)
     241          isBasic[i]=-1;
     242        for (i=0;i<numberRowBasic;i++) {
     243          int iRow = pivotVariable[i];
     244          rowIsBasic[iRow]=1;
     245        }
     246        for (;i<numberBasic;i++) {
     247          int iColumn = pivotVariable[i];
     248          columnIsBasic[iColumn]=1;
     249        }
     250        numberBasic=0;
     251        for (i=0;i<numberRows;i++)
     252          pivotVariable[i]=-1;
     253        // mark as basic or non basic
     254        for (i=0;i<numberRows;i++) {
     255          if (rowIsBasic[i]>=0) {
     256            if (pivotColumn_[numberBasic]>=0)
     257              rowIsBasic[i]=pivotColumn_[numberBasic];
     258            else
     259              rowIsBasic[i]=-1;
     260            numberBasic++;
     261          }
     262        }
     263        for (i=0;i<numberColumns;i++) {
     264          if (columnIsBasic[i]>=0) {
     265            if (pivotColumn_[numberBasic]>=0)
     266              columnIsBasic[i]=pivotColumn_[numberBasic];
     267            else
     268              columnIsBasic[i]=-1;
     269            numberBasic++;
     270          }
     271        }
     272        // leave pivotVariable in useful form for cleaning basis
     273        int * pivotVariable = model->pivotVariable();
     274        for (i=0;i<numberRows;i++) {
     275          pivotVariable[i]=-1;
     276        }
     277        for (i=0;i<numberRows;i++) {
     278          if (model->getRowStatus(i) == ClpSimplex::basic) {
     279            int iPivot = rowIsBasic[i];
     280            if (iPivot>=0)
     281              pivotVariable[iPivot]=i+numberColumns;
     282          }
     283        }
     284        for (i=0;i<numberColumns;i++) {
     285          if (model->getColumnStatus(i) == ClpSimplex::basic) {
     286            int iPivot = columnIsBasic[i];
     287            if (iPivot>=0)
     288              pivotVariable[iPivot]=i;
     289          }
     290        }
     291        delete [] isBasic;
     292        double * columnLower = model->lowerRegion();
     293        double * columnUpper = model->upperRegion();
     294        double * columnActivity = model->solutionRegion();
     295        double * rowLower = model->lowerRegion(0);
     296        double * rowUpper = model->upperRegion(0);
     297        double * rowActivity = model->solutionRegion(0);
     298        //redo basis - first take ALL columns out
     299        int iColumn;
     300        double largeValue = model->largeValue();
     301        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     302          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
     303            // take out
     304            if (!valuesPass) {
     305              double lower=columnLower[iColumn];
     306              double upper=columnUpper[iColumn];
     307              double value=columnActivity[iColumn];
     308              if (lower>-largeValue||upper<largeValue) {
     309                if (fabs(value-lower)<fabs(value-upper)) {
     310                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
     311                  columnActivity[iColumn]=lower;
     312                } else {
     313                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
     314                  columnActivity[iColumn]=upper;
     315                }
     316              } else {
     317                model->setColumnStatus(iColumn,ClpSimplex::isFree);
     318              }
     319            } else {
     320              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
     321            }
     322          }
     323        }
     324        int iRow;
     325        for (iRow=0;iRow<numberRows;iRow++) {
     326          int iSequence=pivotVariable[iRow];
     327          if (iSequence>=0) {
     328            // basic
     329            if (iSequence>=numberColumns) {
     330              // slack in - leave
     331              //assert (iSequence-numberColumns==iRow);
     332            } else {
     333              // put back structural
     334              model->setColumnStatus(iSequence,ClpSimplex::basic);
     335            }
     336          } else {
     337            // put in slack
     338            model->setRowStatus(iRow,ClpSimplex::basic);
     339          }
     340        }
     341        // signal repeat
     342        status_=-99;
     343        // set fixed if they are
     344        for (iRow=0;iRow<numberRows;iRow++) {
     345          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
     346            if (rowLower[iRow]==rowUpper[iRow]) {
     347              rowActivity[iRow]=rowLower[iRow];
     348              model->setRowStatus(iRow,ClpSimplex::isFixed);
     349            }
     350          }
     351        }
     352        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     353          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
     354            if (columnLower[iColumn]==columnUpper[iColumn]) {
     355              columnActivity[iColumn]=columnLower[iColumn];
     356              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
     357            }
     358          }
     359        }
     360      }
    239361    }
    240362  } else {
     
    285407   region1 starts as zero and is zero at end */
    286408int
    287 ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
    288                                  CoinIndexedVector * regionSparse2,
    289                                  bool FTUpdate)
     409ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
     410                                   CoinIndexedVector * regionSparse2)
    290411{
    291412#ifdef CLP_DEBUG
     
    294415  if (!networkBasis_) {
    295416    collectStatistics_ = true;
    296     return CoinFactorization::updateColumn(regionSparse,
    297                                            regionSparse2,
    298                                            FTUpdate);
     417    int returnValue= CoinFactorization::updateColumnFT(regionSparse,
     418                                                       regionSparse2);
    299419    collectStatistics_ = false;
     420    return returnValue;
    300421  } else {
    301422#ifdef CHECK_NETWORK
    302423    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    303     int returnCode = CoinFactorization::updateColumn(regionSparse,
    304                                                      regionSparse2, FTUpdate);
     424    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
     425                                                       regionSparse2);
    305426    networkBasis_->updateColumn(regionSparse,save);
    306427    int i;
     
    315436    return returnCode;
    316437#else
    317     return networkBasis_->updateColumn(regionSparse,regionSparse2);
    318 #endif
    319   }
    320 }
    321 /* Updates one column (FTRAN) to/from array
     438    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
     439    return 1;
     440#endif
     441  }
     442}
     443/* Updates one column (FTRAN) from region2
    322444   number returned is negative if no room
    323    ** For large problems you should ALWAYS know where the nonzeros
    324    are, so please try and migrate to previous method after you
    325    have got code working using this simple method - thank you!
    326    (the only exception is if you know input is dense e.g. rhs)
    327    region starts as zero and is zero at end */
     445   region1 starts as zero and is zero at end */
    328446int
    329447ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
    330                         double array[] ) const
    331 {
     448                                 CoinIndexedVector * regionSparse2,
     449                                 bool noPermute) const
     450{
     451#ifdef CLP_DEBUG
     452  regionSparse->checkClear();
     453#endif
    332454  if (!networkBasis_) {
    333     return CoinFactorization::updateColumn(regionSparse,
    334                                            array);
     455    collectStatistics_ = true;
     456    int returnValue= CoinFactorization::updateColumn(regionSparse,
     457                                                     regionSparse2,
     458                                                     noPermute);
     459    collectStatistics_ = false;
     460    return returnValue;
    335461  } else {
    336462#ifdef CHECK_NETWORK
    337     double * save = new double [numberRows_+1];
    338     memcpy(save,array,(numberRows_+1)*sizeof(double));
     463    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    339464    int returnCode = CoinFactorization::updateColumn(regionSparse,
    340                                                      array);
    341     networkBasis_->updateColumn(regionSparse, save);
     465                                                     regionSparse2,
     466                                                     noPermute);
     467    networkBasis_->updateColumn(regionSparse,save);
    342468    int i;
    343     for (i=0;i<numberRows_;i++)
    344       assert (fabs(save[i]-array[i])<1.0e-8*(1.0+fabs(array[i])));
    345     delete [] save;
     469    double * array = regionSparse2->denseVector();
     470    double * array2 = save->denseVector();
     471    for (i=0;i<numberRows_;i++) {
     472      double value1 = array[i];
     473      double value2 = array2[i];
     474      assert (value1==value2);
     475    }
     476    delete save;
    346477    return returnCode;
    347478#else
    348     return networkBasis_->updateColumn(regionSparse, array);
    349 #endif
    350   }
    351 }
    352 /* Updates one column transpose (BTRAN)
    353    For large problems you should ALWAYS know where the nonzeros
    354    are, so please try and migrate to previous method after you
    355    have got code working using this simple method - thank you!
    356    (the only exception is if you know input is dense e.g. dense objective)
    357    returns number of nonzeros */
    358 int
    359 ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
    360                                           double array[] ) const
    361 {
    362   if (!networkBasis_) {
    363     return CoinFactorization::updateColumnTranspose(regionSparse,
    364                                                     array);
    365   } else {
    366 #ifdef CHECK_NETWORK
    367     double * save = new double [numberRows_+1];
    368     memcpy(save,array,(numberRows_+1)*sizeof(double));
    369     int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
    370                                                               array);
    371     networkBasis_->updateColumnTranspose(regionSparse, save);
    372     int i;
    373     for (i=0;i<numberRows_;i++)
    374       assert (fabs(save[i]-array[i])<1.0e-8*(1.0+fabs(array[i])));
    375     delete [] save;
    376     return returnCode;
    377 #else
    378     return networkBasis_->updateColumnTranspose(regionSparse, array);
     479    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
     480    return 1;
    379481#endif
    380482  }
  • branches/pre/ClpLinearObjective.cpp

    r124 r180  
    5353  }
    5454}
     55/* Subset constructor.  Duplicates are allowed
     56   and order is as given.
     57*/
     58ClpLinearObjective::ClpLinearObjective (const ClpLinearObjective &rhs,
     59                                        int numberColumns,
     60                                        const int * whichColumn)
     61: ClpObjective(rhs)
     62{
     63  objective_=NULL;
     64  numberColumns_=0;
     65  if (numberColumns>0) {
     66    // check valid lists
     67    int numberBad=0;
     68    int i;
     69    for (i=0;i<numberColumns;i++)
     70      if (whichColumn[i]<0||whichColumn[i]>=rhs.numberColumns_)
     71        numberBad++;
     72    if (numberBad)
     73      throw CoinError("bad column list", "subset constructor",
     74                      "ClpLinearObjective");
     75    numberColumns_ = numberColumns;
     76    objective_ = new double[numberColumns_];
     77    for (i=0;i<numberColumns_;i++)
     78      objective_[i]=rhs.objective_[whichColumn[i]];
     79  }
     80}
     81 
    5582
    5683//-------------------------------------------------------------------
     
    95122{
    96123  return new ClpLinearObjective(*this);
     124}
     125/* Subset clone.  Duplicates are allowed
     126   and order is as given.
     127*/
     128ClpObjective *
     129ClpLinearObjective::subsetClone (int numberColumns,
     130                           const int * whichColumns) const
     131{
     132  return new ClpLinearObjective(*this, numberColumns, whichColumns);
    97133}
    98134// Resize objective
  • branches/pre/ClpMatrixBase.cpp

    r63 r180  
    6969  abort();
    7070}
     71/* Subset clone (without gaps).  Duplicates are allowed
     72   and order is as given.
     73   Derived classes need not provide this as it may not always make
     74   sense */
     75ClpMatrixBase *
     76ClpMatrixBase::subsetClone (
     77                            int numberRows, const int * whichRows,
     78                            int numberColumns, const int * whichColumns) const
     79 
    7180
     81{
     82  std::cerr<<"subsetClone not supported - ClpMatrixBase"<<std::endl;
     83  abort();
     84}
  • branches/pre/ClpMessage.cpp

    r161 r180  
    2626  {CLP_SIMPLEX_FLAG,10,3,"Flagging variable %c%d"},
    2727  {CLP_SIMPLEX_GIVINGUP,11,2,"Stopping as close enough"},
    28   {CLP_DUAL_CHECKB,12,3,"Looking infeasible checking bounds with %g"},
     28  {CLP_DUAL_CHECKB,12,2,"New dual bound of %g"},
    2929  {CLP_DUAL_ORIGINAL,13,3,"Going back to original objective"},
    3030  {CLP_SIMPLEX_PERTURB,14,1,"Perturbing problem by %g"},
     
    4646  {CLP_DUAL_CHECK,106,4,"Btran alpha %g, ftran alpha %g"},
    4747  {CLP_PRIMAL_DJ,107,4,"Btran dj %g, ftran dj %g"},
    48   {CLP_PRESOLVE_COLINFEAS,501,2,"Problem is infeasible due to column %d, %g %g"},
    49   {CLP_PRESOLVE_ROWINFEAS,502,2,"Problem is infeasible due to row %d, %g %g"},
    50   {CLP_PRESOLVE_COLUMNBOUNDA,503,2,"Problem looks unbounded above due to column %d, %g %g"},
    51   {CLP_PRESOLVE_COLUMNBOUNDB,504,2,"Problem looks unbounded below due to column %d, %g %g"},
    52   {CLP_PRESOLVE_NONOPTIMAL,505,1,"Problem not optimal, resolve after postsolve"},
    53   {CLP_PRESOLVE_STATS,506,1,"Presolve %d (%d) rows, %d (%d) columns and %d (%d) elements"},
    54   {CLP_PRESOLVE_INFEAS,507,0,"Presolve determined that the problem was infeasible with tolerance of %g"},
    55   {CLP_PRESOLVE_UNBOUND,508,0,"Presolve thinks problem is unbounded"},
    56   {CLP_PRESOLVE_INFEASUNBOUND,509,0,"Presolve thinks problem is infeasible AND unbounded???"},
    57   {CLP_PRESOLVE_INTEGERMODS,510,1,"Presolve is modifying %d integer bounds and re-presolving"},
    58   {CLP_PRESOLVE_POSTSOLVE,511,0,"After Postsolve, objective %g, infeasibilities - dual %g (%d), primal %g (%d)"},
    59   {CLP_PRESOLVE_NEEDS_CLEANING,512,0,"Presolved model was optimal, full model needs cleaning up"},
    6048  {CLP_PACKEDSCALE_INITIAL,1001,2,"Initial range of elements is %g to %g"},
    6149  {CLP_PACKEDSCALE_WHILE,1002,3,"Range of elements is %g to %g"},
  • branches/pre/ClpModel.cpp

    r169 r180  
    99#include <iostream>
    1010
    11 #include <time.h>
    1211
    1312#include "CoinPragma.hpp"
    14 #ifndef _MSC_VER
    15 #include <sys/times.h>
    16 #include <sys/resource.h>
    17 #include <unistd.h>
    18 #endif
    1913
    2014#include "CoinHelperFunctions.hpp"
     
    2721#include "ClpLinearObjective.hpp"
    2822
    29 static double cpuTime()
    30 {
    31   double cpu_temp;
    32 #if defined(_MSC_VER)
    33   unsigned int ticksnow;        /* clock_t is same as int */
    34  
    35   ticksnow = (unsigned int)clock();
    36  
    37   cpu_temp = (double)((double)ticksnow/CLOCKS_PER_SEC);
    38 #else
    39   struct rusage usage;
    40   getrusage(RUSAGE_SELF,&usage);
    41   cpu_temp = usage.ru_utime.tv_sec;
    42   cpu_temp += 1.0e-6*((double) usage.ru_utime.tv_usec);
    43 #endif
    44   return cpu_temp;
    45 }
    4623//#############################################################################
    4724
     
    339316  strParam_[ClpProbName] = rhs.strParam_[ClpProbName];
    340317
     318  optimizationDirection_ = rhs.optimizationDirection_;
    341319  objectiveValue_=rhs.objectiveValue_;
    342320  smallElement_ = rhs.smallElement_;
     
    344322  solveType_ = rhs.solveType_;
    345323  problemStatus_ = rhs.problemStatus_;
     324  numberRows_ = rhs.numberRows_;
     325  numberColumns_ = rhs.numberColumns_;
    346326  if (trueCopy) {
    347327    lengthNames_ = rhs.lengthNames_;
     
    527507  case ClpMaxSeconds:
    528508    if(value>=0)
    529       value += cpuTime();
     509      value += CoinCpuTime();
    530510    else
    531511      value = -1.0;
     
    733713                              number, which, newSize);
    734714  matrix_->deleteRows(number,which);
     715  //matrix_->removeGaps();
    735716  // status
    736717  if (status_) {
     
    770751                              number, which, newSize);
    771752  matrix_->deleteCols(number,which);
     753  //matrix_->removeGaps();
    772754  if (quadraticObjective_) {
    773755    quadraticObjective_->deleteCols(number,which);
     
    787769    status_ = temp;
    788770  }
     771  integerType_ = deleteChar(integerType_,numberColumns_,
     772                            number, which, newSize,true);
    789773  numberColumns_=newSize;
    790774  // set state back to unknown
     
    796780  rowNames_ = std::vector<std::string> ();
    797781  columnNames_ = std::vector<std::string> ();
    798   integerType_ = deleteChar(integerType_,numberColumns_,
    799                             number, which, newSize,true);
    800782}
    801783// Add rows
     
    815797      rows[iRow] =
    816798        new CoinPackedVector(rowStarts[iRow+1]-iStart,
     799                             columns+iStart,elements+iStart);
     800    }
     801    addRows(number, rowLower, rowUpper,
     802            rows);
     803    for (iRow=0;iRow<number;iRow++)
     804      delete rows[iRow];
     805    delete [] rows;
     806  }
     807}
     808// Add rows
     809void
     810ClpModel::addRows(int number, const double * rowLower,
     811                  const double * rowUpper,
     812                  const int * rowStarts,
     813                  const int * rowLengths, const int * columns,
     814                  const double * elements)
     815{
     816  // Create a list of CoinPackedVectors
     817  if (number) {
     818    CoinPackedVectorBase ** rows=
     819      new CoinPackedVectorBase * [number];
     820    int iRow;
     821    for (iRow=0;iRow<number;iRow++) {
     822      int iStart = rowStarts[iRow];
     823      rows[iRow] =
     824        new CoinPackedVector(rowLengths[iRow],
    817825                             columns+iStart,elements+iStart);
    818826    }
     
    886894      columns[iColumn] =
    887895        new CoinPackedVector(columnStarts[iColumn+1]-iStart,
     896                             rows+iStart,elements+iStart);
     897    }
     898    addColumns(number, columnLower, columnUpper,
     899               objIn, columns);
     900    for (iColumn=0;iColumn<number;iColumn++)
     901      delete columns[iColumn];
     902    delete [] columns;
     903
     904  }
     905}
     906// Add columns
     907void
     908ClpModel::addColumns(int number, const double * columnLower,
     909                     const double * columnUpper,
     910                     const double * objIn,
     911                     const int * columnStarts,
     912                     const int * columnLengths, const int * rows,
     913                     const double * elements)
     914{
     915  // Create a list of CoinPackedVectors
     916  if (number) {
     917    CoinPackedVectorBase ** columns=
     918      new CoinPackedVectorBase * [number];
     919    int iColumn;
     920    for (iColumn=0;iColumn<number;iColumn++) {
     921      int iStart = columnStarts[iColumn];
     922      columns[iColumn] =
     923        new CoinPackedVector(columnLengths[iColumn],
    888924                             rows+iStart,elements+iStart);
    889925    }
     
    9791015{
    9801016  if(value>=0)
    981     dblParam_[ClpMaxSeconds]=value+cpuTime();
     1017    dblParam_[ClpMaxSeconds]=value+CoinCpuTime();
    9821018  else
    9831019    dblParam_[ClpMaxSeconds]=-1.0;
     
    9891025  bool hitMax= (numberIterations_>=maximumIterations());
    9901026  if (dblParam_[ClpMaxSeconds]>=0.0&&!hitMax)
    991     hitMax = (cpuTime()>=dblParam_[ClpMaxSeconds]);
     1027    hitMax = (CoinCpuTime()>=dblParam_[ClpMaxSeconds]);
    9921028  return hitMax;
    9931029}
     
    10301066  }
    10311067  CoinMpsIO m;
    1032   double time1 = cpuTime(),time2;
     1068  double time1 = CoinCpuTime(),time2;
    10331069  int status=m.readMps(fileName,"");
    10341070  if (!status||ignoreErrors) {
     
    10681104    }
    10691105    setDblParam(ClpObjOffset,m.objectiveOffset());
    1070     time2 = cpuTime();
     1106    time2 = CoinCpuTime();
    10711107    handler_->message(CLP_IMPORT_RESULT,messages_)
    10721108      <<fileName
     
    11881224  delete quadraticObjective_;
    11891225  quadraticObjective_ = NULL;
     1226}
     1227// Returns resized array and updates size
     1228double * whichDouble(double * array , int number, const int * which)
     1229{
     1230  double * newArray=NULL;
     1231  if (array&&number) {
     1232    int i ;
     1233    newArray = new double[number];
     1234    for (i=0;i<number;i++)
     1235      newArray[i]=array[which[i]];
     1236  }
     1237  return newArray;
     1238}
     1239char * whichChar(char * array , int number, const int * which)
     1240{
     1241  char * newArray=NULL;
     1242  if (array&&number) {
     1243    int i ;
     1244    newArray = new char[number];
     1245    for (i=0;i<number;i++)
     1246      newArray[i]=array[which[i]];
     1247  }
     1248  return newArray;
     1249}
     1250unsigned char * whichUnsignedChar(unsigned char * array ,
     1251                                  int number, const int * which)
     1252{
     1253  unsigned char * newArray=NULL;
     1254  if (array&&number) {
     1255    int i ;
     1256    newArray = new unsigned char[number];
     1257    for (i=0;i<number;i++)
     1258      newArray[i]=array[which[i]];
     1259  }
     1260  return newArray;
     1261}
     1262// Subproblem constructor
     1263ClpModel::ClpModel ( const ClpModel * rhs,
     1264                     int numberRows, const int * whichRow,
     1265                     int numberColumns, const int * whichColumn,
     1266                     bool dropNames, bool dropIntegers)
     1267{
     1268#if 0
     1269  // Could be recoded to be faster and take less memory
     1270  // and to allow re-ordering and duplicates etc
     1271  gutsOfCopy(*rhs,true);
     1272  int numberRowsWhole = rhs->numberRows();
     1273  int * delRow = new int[numberRowsWhole];
     1274  memset(delRow,0,numberRowsWhole*sizeof(int));
     1275  int i;
     1276  for (i=0;i<numberRows;i++) {
     1277    int iRow=whichRow[i];
     1278    assert (iRow>=0&&iRow<numberRowsWhole);
     1279    delRow[iRow]=1;
     1280  }
     1281  numberRows=0;
     1282  for (i=0;i<numberRowsWhole;i++) {
     1283    if (delRow[i]==0)
     1284      delRow[numberRows++]=i;
     1285  }
     1286  deleteRows(numberRows,delRow);
     1287  delete [] delRow;
     1288  int numberColumnsWhole = rhs->numberColumns();
     1289  int * delColumn = new int[numberColumnsWhole];
     1290  memset(delColumn,0,numberColumnsWhole*sizeof(int));
     1291  for (i=0;i<numberColumns;i++) {
     1292    int iColumn=whichColumn[i];
     1293    assert (iColumn>=0&&iColumn<numberColumnsWhole);
     1294    delColumn[iColumn]=1;
     1295  }
     1296  numberColumns=0;
     1297  for (i=0;i<numberColumnsWhole;i++) {
     1298    if (delColumn[i]==0)
     1299      delColumn[numberColumns++]=i;
     1300  }
     1301  deleteColumns(numberColumns,delColumn);
     1302  delete [] delColumn;
     1303#else
     1304  defaultHandler_ = rhs->defaultHandler_;
     1305  if (defaultHandler_)
     1306    handler_ = new CoinMessageHandler(*rhs->handler_);
     1307   else
     1308    handler_ = rhs->handler_;
     1309  messages_ = rhs->messages_;
     1310  intParam_[ClpMaxNumIteration] = rhs->intParam_[ClpMaxNumIteration];
     1311  intParam_[ClpMaxNumIterationHotStart] =
     1312    rhs->intParam_[ClpMaxNumIterationHotStart];
     1313
     1314  dblParam_[ClpDualObjectiveLimit] = rhs->dblParam_[ClpDualObjectiveLimit];
     1315  dblParam_[ClpPrimalObjectiveLimit] = rhs->dblParam_[ClpPrimalObjectiveLimit];
     1316  dblParam_[ClpDualTolerance] = rhs->dblParam_[ClpDualTolerance];
     1317  dblParam_[ClpPrimalTolerance] = rhs->dblParam_[ClpPrimalTolerance];
     1318  dblParam_[ClpObjOffset] = rhs->dblParam_[ClpObjOffset];
     1319  dblParam_[ClpMaxSeconds] = rhs->dblParam_[ClpMaxSeconds];
     1320
     1321  strParam_[ClpProbName] = rhs->strParam_[ClpProbName];
     1322
     1323  optimizationDirection_ = rhs->optimizationDirection_;
     1324  objectiveValue_=rhs->objectiveValue_;
     1325  smallElement_ = rhs->smallElement_;
     1326  numberIterations_ = rhs->numberIterations_;
     1327  solveType_ = rhs->solveType_;
     1328  problemStatus_ = rhs->problemStatus_;
     1329  // check valid lists
     1330  int numberBad=0;
     1331  int i;
     1332  for (i=0;i<numberRows;i++)
     1333    if (whichRow[i]<0||whichRow[i]>=rhs->numberRows_)
     1334      numberBad++;
     1335  if (numberBad)
     1336    throw CoinError("bad row list", "subproblem constructor", "ClpModel");
     1337  numberBad=0;
     1338  for (i=0;i<numberColumns;i++)
     1339    if (whichColumn[i]<0||whichColumn[i]>=rhs->numberColumns_)
     1340      numberBad++;
     1341  if (numberBad)
     1342    throw CoinError("bad column list", "subproblem constructor", "ClpModel");
     1343  numberRows_ = numberRows;
     1344  numberColumns_ = numberColumns;
     1345  if (!dropNames) {
     1346    unsigned int maxLength=0;
     1347    int iRow;
     1348    rowNames_.reserve(numberRows_);
     1349    for (iRow=0;iRow<numberRows_;iRow++) {
     1350      rowNames_[iRow] = rhs->rowNames_[whichRow[iRow]];
     1351      maxLength = max(maxLength,(unsigned int) strlen(rowNames_[iRow].c_str()));
     1352    }
     1353    int iColumn;
     1354    columnNames_.reserve(numberColumns_);
     1355    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1356      columnNames_[iColumn] = rhs->columnNames_[whichColumn[iColumn]];
     1357      maxLength = max(maxLength,(unsigned int) strlen(columnNames_[iColumn].c_str()));
     1358    }
     1359    lengthNames_=(int) maxLength;
     1360  } else {
     1361    lengthNames_ = 0;
     1362    rowNames_ = std::vector<std::string> ();
     1363    columnNames_ = std::vector<std::string> ();
     1364  }
     1365  if (rhs->integerType_&&!dropIntegers) {
     1366    integerType_ = whichChar(rhs->integerType_,numberColumns,whichColumn);
     1367  } else {
     1368    integerType_ = NULL;
     1369  }
     1370  if (rhs->rowActivity_) {
     1371    rowActivity_=whichDouble(rhs->rowActivity_,numberRows,whichRow);
     1372    dual_=whichDouble(rhs->dual_,numberRows,whichRow);
     1373    columnActivity_=whichDouble(rhs->columnActivity_,numberColumns,
     1374                                whichColumn);
     1375    reducedCost_=whichDouble(rhs->reducedCost_,numberColumns,
     1376                                whichColumn);
     1377  } else {
     1378    rowActivity_=NULL;
     1379    columnActivity_=NULL;
     1380    dual_=NULL;
     1381    reducedCost_=NULL;
     1382  }
     1383  rowLower_=whichDouble(rhs->rowLower_,numberRows,whichRow);
     1384  rowUpper_=whichDouble(rhs->rowUpper_,numberRows,whichRow);
     1385  columnLower_=whichDouble(rhs->columnLower_,numberColumns,whichColumn);
     1386  columnUpper_=whichDouble(rhs->columnUpper_,numberColumns,whichColumn);
     1387  if (rhs->objective_)
     1388    objective_  = rhs->objective_->subsetClone(numberColumns,whichColumn);
     1389  else
     1390    objective_ = NULL;
     1391  rowObjective_=whichDouble(rhs->rowObjective_,numberRows,whichRow);
     1392  // status has to be done in two stages
     1393  status_ = new unsigned char[numberColumns_+numberRows_];
     1394  unsigned char * rowStatus = whichUnsignedChar(rhs->status_+rhs->numberColumns_,
     1395                                                numberRows_,whichRow);
     1396  unsigned char * columnStatus = whichUnsignedChar(rhs->status_,
     1397                                                numberColumns_,whichColumn);
     1398  memcpy(status_+numberColumns_,rowStatus,numberRows_);
     1399  delete [] rowStatus;
     1400  memcpy(status_,columnStatus,numberColumns_);
     1401  delete [] columnStatus;
     1402  ray_ = NULL;
     1403  if (problemStatus_==1)
     1404    ray_ = whichDouble (rhs->ray_,numberRows,whichRow);
     1405  else if (problemStatus_==2)
     1406    ray_ = whichDouble (rhs->ray_,numberColumns,whichColumn);
     1407  if (rhs->rowCopy_) {
     1408    rowCopy_ = rhs->rowCopy_->subsetClone(numberRows,whichRow,
     1409                                          numberColumns,whichColumn);
     1410  } else {
     1411    rowCopy_=NULL;
     1412  }
     1413  if (rhs->quadraticObjective_) {
     1414    quadraticObjective_ = new CoinPackedMatrix(*rhs->quadraticObjective_,
     1415                                               numberColumns,whichColumn,
     1416                                               numberColumns,whichColumn);
     1417  } else {
     1418    quadraticObjective_=NULL;
     1419  }
     1420  matrix_=NULL;
     1421  if (rhs->matrix_) {
     1422    matrix_ = rhs->matrix_->subsetClone(numberRows,whichRow,
     1423                                        numberColumns,whichColumn);
     1424  }
     1425#endif
    11901426}
    11911427//#############################################################################
  • branches/pre/ClpNetworkBasis.cpp

    r131 r180  
    561561
    562562/* Updates one column (FTRAN) from region2 */
    563 int
     563double
    564564ClpNetworkBasis::updateColumn (  CoinIndexedVector * regionSparse,
    565                                  CoinIndexedVector * regionSparse2)
     565                                 CoinIndexedVector * regionSparse2,
     566                                 int pivotRow)
    566567{
    567568  regionSparse->clear (  );
     
    573574  int i;
    574575  bool doTwo = (numberNonZero==2);
    575   int i0;
    576   int i1;
     576  int i0=-1;
     577  int i1=-1;
    577578  if (doTwo) {
    578579    i0 = regionIndex2[0];
    579580    i1 = regionIndex2[1];
    580     doTwo =  (region2[i0]*region2[i1]<0.0);
    581   }
    582   if (doTwo) {
    583     // If just +- 1 then could go backwards on depth until join
    584     region[i0] = region2[i0];
    585     region2[i0]=0.0;
    586     region[i1] = region2[i1];
    587     region2[i1]=0.0;
    588     int iDepth0 = depth_[i0];
    589     int iDepth1 = depth_[i1];
    590     if (iDepth1>iDepth0) {
    591       int temp = i0;
    592       i0 = i1;
    593       i1 = temp;
    594       temp = iDepth0;
    595       iDepth0 = iDepth1;
    596       iDepth1 = temp;
    597     }
    598     numberNonZero=0;
    599     while (iDepth0>iDepth1) {
    600       double pivotValue = region[i0];
    601       // put back now ?
    602       int iBack = permuteBack_[i0];
    603       regionIndex2[numberNonZero++] = iBack;
    604       int otherRow = parent_[i0];
    605       region2[iBack] = pivotValue*sign_[i0];
    606       region[i0] =0.0;
    607       region[otherRow] += pivotValue;
    608       iDepth0--;
    609       i0 = otherRow;
    610     }
    611     while (i0!=i1) {
    612       double pivotValue = region[i0];
    613       // put back now ?
    614       int iBack = permuteBack_[i0];
    615       regionIndex2[numberNonZero++] = iBack;
    616       int otherRow = parent_[i0];
    617       region2[iBack] = pivotValue*sign_[i0];
    618       region[i0] =0.0;
    619       region[otherRow] += pivotValue;
    620       i0 = otherRow;
    621       double pivotValue1 = region[i1];
    622       // put back now ?
    623       int iBack1 = permuteBack_[i1];
    624       regionIndex2[numberNonZero++] = iBack1;
    625       int otherRow1 = parent_[i1];
    626       region2[iBack1] = pivotValue1*sign_[i1];
    627       region[i1] =0.0;
    628       region[otherRow1] += pivotValue1;
    629       i1 = otherRow1;
    630     }
    631   } else {
    632     // set up linked lists at each depth
    633     // stack2 is start, stack is next
    634     int greatestDepth=-1;
    635     //mark_[numberRows_]=1;
    636     for (i=0;i<numberNonZero;i++) {
    637       int j = regionIndex2[i];
    638       double value = region2[j];
    639       region2[j]=0.0;
    640       region[j]=value;
    641       regionIndex[i]=j;
    642       int iDepth = depth_[j];
    643       if (iDepth>greatestDepth)
    644         greatestDepth = iDepth;
    645       // and back until marked
    646       while (!mark_[j]) {
    647         int iNext = stack2_[iDepth];
    648         stack2_[iDepth]=j;
    649         stack_[j]=iNext;
    650         mark_[j]=1;
    651         iDepth--;
    652         j=parent_[j];
    653       }
    654     }
    655     numberNonZero=0;
    656     for (;greatestDepth>=0; greatestDepth--) {
    657       int iPivot = stack2_[greatestDepth];
    658       stack2_[greatestDepth]=-1;
    659       while (iPivot>=0) {
    660         mark_[iPivot]=0;
    661         double pivotValue = region[iPivot];
    662         if (pivotValue) {
     581  }
     582  double returnValue=0.0;
     583  bool packed = regionSparse2->packedMode();
     584  if (packed) {
     585    if (doTwo&&region2[0]*region2[1]<0.0) {
     586      region[i0] = region2[0];
     587      region2[0]=0.0;
     588      region[i1] = region2[1];
     589      region2[1]=0.0;
     590      int iDepth0 = depth_[i0];
     591      int iDepth1 = depth_[i1];
     592      if (iDepth1>iDepth0) {
     593        int temp = i0;
     594        i0 = i1;
     595        i1 = temp;
     596        temp = iDepth0;
     597        iDepth0 = iDepth1;
     598        iDepth1 = temp;
     599      }
     600      numberNonZero=0;
     601      if (pivotRow<0) {
     602        while (iDepth0>iDepth1) {
     603          double pivotValue = region[i0];
    663604          // put back now ?
    664           int iBack = permuteBack_[iPivot];
     605          int iBack = permuteBack_[i0];
     606          region2[numberNonZero] = pivotValue*sign_[i0];
    665607          regionIndex2[numberNonZero++] = iBack;
    666           int otherRow = parent_[iPivot];
    667           region2[iBack] = pivotValue*sign_[iPivot];
    668           region[iPivot] =0.0;
     608          int otherRow = parent_[i0];
     609          region[i0] =0.0;
    669610          region[otherRow] += pivotValue;
    670         }
    671         iPivot = stack_[iPivot];
    672       }
    673     }
     611          iDepth0--;
     612          i0 = otherRow;
     613        }
     614        while (i0!=i1) {
     615          double pivotValue = region[i0];
     616          // put back now ?
     617          int iBack = permuteBack_[i0];
     618          region2[numberNonZero] = pivotValue*sign_[i0];
     619          regionIndex2[numberNonZero++] = iBack;
     620          int otherRow = parent_[i0];
     621          region[i0] =0.0;
     622          region[otherRow] += pivotValue;
     623          i0 = otherRow;
     624          double pivotValue1 = region[i1];
     625          // put back now ?
     626          int iBack1 = permuteBack_[i1];
     627          region2[numberNonZero] = pivotValue1*sign_[i1];
     628          regionIndex2[numberNonZero++] = iBack1;
     629          int otherRow1 = parent_[i1];
     630          region[i1] =0.0;
     631          region[otherRow1] += pivotValue1;
     632          i1 = otherRow1;
     633        }
     634      } else {
     635        while (iDepth0>iDepth1) {
     636          double pivotValue = region[i0];
     637          // put back now ?
     638          int iBack = permuteBack_[i0];
     639          double value = pivotValue*sign_[i0];
     640          region2[numberNonZero] = value;
     641          regionIndex2[numberNonZero++] = iBack;
     642          if (iBack==pivotRow)
     643            returnValue = value;
     644          int otherRow = parent_[i0];
     645          region[i0] =0.0;
     646          region[otherRow] += pivotValue;
     647          iDepth0--;
     648          i0 = otherRow;
     649        }
     650        while (i0!=i1) {
     651          double pivotValue = region[i0];
     652          // put back now ?
     653          int iBack = permuteBack_[i0];
     654          double value = pivotValue*sign_[i0];
     655          region2[numberNonZero] = value;
     656          regionIndex2[numberNonZero++] = iBack;
     657          if (iBack==pivotRow)
     658            returnValue = value;
     659          int otherRow = parent_[i0];
     660          region[i0] =0.0;
     661          region[otherRow] += pivotValue;
     662          i0 = otherRow;
     663          double pivotValue1 = region[i1];
     664          // put back now ?
     665          int iBack1 = permuteBack_[i1];
     666          value = pivotValue1*sign_[i1];
     667          region2[numberNonZero] = value;
     668          regionIndex2[numberNonZero++] = iBack1;
     669          if (iBack1==pivotRow)
     670            returnValue = value;
     671          int otherRow1 = parent_[i1];
     672          region[i1] =0.0;
     673          region[otherRow1] += pivotValue1;
     674          i1 = otherRow1;
     675        }
     676      }
     677    } else {
     678      // set up linked lists at each depth
     679      // stack2 is start, stack is next
     680      int greatestDepth=-1;
     681      //mark_[numberRows_]=1;
     682      for (i=0;i<numberNonZero;i++) {
     683        int j = regionIndex2[i];
     684        double value = region2[j];
     685        region2[j]=0.0;
     686        region[j]=value;
     687        regionIndex[i]=j;
     688        int iDepth = depth_[j];
     689        if (iDepth>greatestDepth)
     690          greatestDepth = iDepth;
     691        // and back until marked
     692        while (!mark_[j]) {
     693          int iNext = stack2_[iDepth];
     694          stack2_[iDepth]=j;
     695          stack_[j]=iNext;
     696          mark_[j]=1;
     697          iDepth--;
     698          j=parent_[j];
     699        }
     700      }
     701      numberNonZero=0;
     702      if (pivotRow<0) {
     703        for (;greatestDepth>=0; greatestDepth--) {
     704          int iPivot = stack2_[greatestDepth];
     705          stack2_[greatestDepth]=-1;
     706          while (iPivot>=0) {
     707            mark_[iPivot]=0;
     708            double pivotValue = region[iPivot];
     709            if (pivotValue) {
     710              // put back now ?
     711              int iBack = permuteBack_[iPivot];
     712              region2[numberNonZero] = pivotValue*sign_[iPivot];
     713              regionIndex2[numberNonZero++] = iBack;
     714              int otherRow = parent_[iPivot];
     715              region[iPivot] =0.0;
     716            region[otherRow] += pivotValue;
     717            }
     718            iPivot = stack_[iPivot];
     719          }
     720        }
     721      } else {
     722        for (;greatestDepth>=0; greatestDepth--) {
     723          int iPivot = stack2_[greatestDepth];
     724          stack2_[greatestDepth]=-1;
     725          while (iPivot>=0) {
     726            mark_[iPivot]=0;
     727            double pivotValue = region[iPivot];
     728            if (pivotValue) {
     729              // put back now ?
     730              int iBack = permuteBack_[iPivot];
     731              double value = pivotValue*sign_[iPivot];
     732              region2[numberNonZero] = value;
     733              regionIndex2[numberNonZero++] = iBack;
     734              if (iBack==pivotRow)
     735                returnValue = value;
     736              int otherRow = parent_[iPivot];
     737              region[iPivot] =0.0;
     738              region[otherRow] += pivotValue;
     739            }
     740            iPivot = stack_[iPivot];
     741          }
     742        }
     743      }
     744    }
     745  } else {
     746    if (doTwo&&region2[i0]*region2[i1]<0.0) {
     747      // If just +- 1 then could go backwards on depth until join
     748      region[i0] = region2[i0];
     749      region2[i0]=0.0;
     750      region[i1] = region2[i1];
     751      region2[i1]=0.0;
     752      int iDepth0 = depth_[i0];
     753      int iDepth1 = depth_[i1];
     754      if (iDepth1>iDepth0) {
     755        int temp = i0;
     756        i0 = i1;
     757        i1 = temp;
     758        temp = iDepth0;
     759        iDepth0 = iDepth1;
     760        iDepth1 = temp;
     761      }
     762      numberNonZero=0;
     763      while (iDepth0>iDepth1) {
     764        double pivotValue = region[i0];
     765        // put back now ?
     766        int iBack = permuteBack_[i0];
     767        regionIndex2[numberNonZero++] = iBack;
     768        int otherRow = parent_[i0];
     769        region2[iBack] = pivotValue*sign_[i0];
     770        region[i0] =0.0;
     771        region[otherRow] += pivotValue;
     772        iDepth0--;
     773        i0 = otherRow;
     774      }
     775      while (i0!=i1) {
     776        double pivotValue = region[i0];
     777        // put back now ?
     778        int iBack = permuteBack_[i0];
     779        regionIndex2[numberNonZero++] = iBack;
     780        int otherRow = parent_[i0];
     781        region2[iBack] = pivotValue*sign_[i0];
     782        region[i0] =0.0;
     783        region[otherRow] += pivotValue;
     784        i0 = otherRow;
     785        double pivotValue1 = region[i1];
     786        // put back now ?
     787        int iBack1 = permuteBack_[i1];
     788        regionIndex2[numberNonZero++] = iBack1;
     789        int otherRow1 = parent_[i1];
     790        region2[iBack1] = pivotValue1*sign_[i1];
     791        region[i1] =0.0;
     792        region[otherRow1] += pivotValue1;
     793        i1 = otherRow1;
     794      }
     795    } else {
     796      // set up linked lists at each depth
     797      // stack2 is start, stack is next
     798      int greatestDepth=-1;
     799      //mark_[numberRows_]=1;
     800      for (i=0;i<numberNonZero;i++) {
     801        int j = regionIndex2[i];
     802        double value = region2[j];
     803        region2[j]=0.0;
     804        region[j]=value;
     805        regionIndex[i]=j;
     806        int iDepth = depth_[j];
     807        if (iDepth>greatestDepth)
     808          greatestDepth = iDepth;
     809        // and back until marked
     810        while (!mark_[j]) {
     811          int iNext = stack2_[iDepth];
     812          stack2_[iDepth]=j;
     813          stack_[j]=iNext;
     814          mark_[j]=1;
     815          iDepth--;
     816          j=parent_[j];
     817        }
     818      }
     819      numberNonZero=0;
     820      for (;greatestDepth>=0; greatestDepth--) {
     821        int iPivot = stack2_[greatestDepth];
     822        stack2_[greatestDepth]=-1;
     823        while (iPivot>=0) {
     824          mark_[iPivot]=0;
     825          double pivotValue = region[iPivot];
     826          if (pivotValue) {
     827            // put back now ?
     828            int iBack = permuteBack_[iPivot];
     829            regionIndex2[numberNonZero++] = iBack;
     830            int otherRow = parent_[iPivot];
     831            region2[iBack] = pivotValue*sign_[iPivot];
     832            region[iPivot] =0.0;
     833            region[otherRow] += pivotValue;
     834          }
     835          iPivot = stack_[iPivot];
     836        }
     837      }
     838    }
     839    if (pivotRow>=0)
     840      returnValue = region2[pivotRow];
    674841  }
    675842  region[numberRows_]=0.0;
     
    684851 }
    685852#endif
    686   return numberNonZero;
     853  return returnValue;
    687854}
    688855/* Updates one column (FTRAN) to/from array
     
    832999  int i;
    8331000  int numberNonZero=0;
    834   for (i=0;i<numberNonZero2;i++) {
    835     int k = regionIndex2[i];
    836     int j = permute_[k];
    837     double value = region2[k];
    838     region2[k]=0.0;
    839     region[j]=value;
    840     mark_[j]=1;
    841     regionIndex[numberNonZero++]=j;
    842   }
    843   // copy back
    844   // set up linked lists at each depth
    845   // stack2 is start, stack is next
    846   int greatestDepth=-1;
    847   int smallestDepth=numberRows_;
    848   //mark_[numberRows_]=1;
    849   for (i=0;i<numberNonZero2;i++) {
    850     int j = regionIndex[i];
    851     double value = region[j];
    852     region[j]=0.0;
    853     region2[j]=value;
    854     regionIndex2[i]=j;
    855     // add in
    856     int iDepth = depth_[j];
    857     smallestDepth = min(iDepth,smallestDepth) ;
    858     greatestDepth = max(iDepth,greatestDepth) ;
    859     int jNext = stack2_[iDepth];
    860     stack2_[iDepth]=j;
    861     stack_[j]=jNext;
    862     // and put all descendants on list
    863     int iChild = descendant_[j];
    864     while (iChild>=0) {
    865       if (!mark_[iChild]) {
    866         regionIndex2[numberNonZero++] = iChild;
    867         mark_[iChild]=1;
    868       }
    869       iChild = rightSibling_[iChild];
    870     }
    871   }
    872   for (;i<numberNonZero;i++) {
    873     int j = regionIndex2[i];
    874     // add in
    875     int iDepth = depth_[j];
    876     smallestDepth = min(iDepth,smallestDepth) ;
    877     greatestDepth = max(iDepth,greatestDepth) ;
    878     int jNext = stack2_[iDepth];
    879     stack2_[iDepth]=j;
    880     stack_[j]=jNext;
    881     // and put all descendants on list
    882     int iChild = descendant_[j];
    883     while (iChild>=0) {
    884       if (!mark_[iChild]) {
    885         regionIndex2[numberNonZero++] = iChild;
    886         mark_[iChild]=1;
    887       }
    888       iChild = rightSibling_[iChild];
    889     }
    890   }
    891   numberNonZero2=0;
    892   region2[numberRows_]=0.0;
    893   int iDepth;
    894   for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
    895     int iPivot = stack2_[iDepth];
    896     stack2_[iDepth]=-1;
    897     while (iPivot>=0) {
    898       mark_[iPivot]=0;
    899       double pivotValue = region2[iPivot];
    900       int otherRow = parent_[iPivot];
    901       double otherValue = region2[otherRow];
    902       pivotValue = sign_[iPivot]*pivotValue+otherValue;
    903       region2[iPivot]=pivotValue;
    904       if (pivotValue)
    905         regionIndex2[numberNonZero2++]=iPivot;
    906       iPivot = stack_[iPivot];
     1001  bool packed = regionSparse2->packedMode();
     1002  if (packed) {
     1003    for (i=0;i<numberNonZero2;i++) {
     1004      int k = regionIndex2[i];
     1005      int j = permute_[k];
     1006      double value = region2[i];
     1007      region2[i]=0.0;
     1008      region[j]=value;
     1009      mark_[j]=1;
     1010      regionIndex[numberNonZero++]=j;
     1011    }
     1012    // set up linked lists at each depth
     1013    // stack2 is start, stack is next
     1014    int greatestDepth=-1;
     1015    int smallestDepth=numberRows_;
     1016    //mark_[numberRows_]=1;
     1017    for (i=0;i<numberNonZero2;i++) {
     1018      int j = regionIndex[i];
     1019      regionIndex2[i]=j;
     1020      // add in
     1021      int iDepth = depth_[j];
     1022      smallestDepth = min(iDepth,smallestDepth) ;
     1023      greatestDepth = max(iDepth,greatestDepth) ;
     1024      int jNext = stack2_[iDepth];
     1025      stack2_[iDepth]=j;
     1026      stack_[j]=jNext;
     1027      // and put all descendants on list
     1028      int iChild = descendant_[j];
     1029      while (iChild>=0) {
     1030        if (!mark_[iChild]) {
     1031          regionIndex2[numberNonZero++] = iChild;
     1032          mark_[iChild]=1;
     1033        }
     1034        iChild = rightSibling_[iChild];
     1035      }
     1036    }
     1037    for (;i<numberNonZero;i++) {
     1038      int j = regionIndex2[i];
     1039      // add in
     1040      int iDepth = depth_[j];
     1041      smallestDepth = min(iDepth,smallestDepth) ;
     1042      greatestDepth = max(iDepth,greatestDepth) ;
     1043      int jNext = stack2_[iDepth];
     1044      stack2_[iDepth]=j;
     1045      stack_[j]=jNext;
     1046      // and put all descendants on list
     1047      int iChild = descendant_[j];
     1048      while (iChild>=0) {
     1049        if (!mark_[iChild]) {
     1050          regionIndex2[numberNonZero++] = iChild;
     1051          mark_[iChild]=1;
     1052        }
     1053        iChild = rightSibling_[iChild];
     1054      }
     1055    }
     1056    numberNonZero2=0;
     1057    region[numberRows_]=0.0;
     1058    int iDepth;
     1059    for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
     1060      int iPivot = stack2_[iDepth];
     1061      stack2_[iDepth]=-1;
     1062      while (iPivot>=0) {
     1063        mark_[iPivot]=0;
     1064        double pivotValue = region[iPivot];
     1065        int otherRow = parent_[iPivot];
     1066        double otherValue = region[otherRow];
     1067        pivotValue = sign_[iPivot]*pivotValue+otherValue;
     1068        region[iPivot]=pivotValue;
     1069        if (pivotValue) {
     1070          region2[numberNonZero2]=pivotValue;
     1071          regionIndex2[numberNonZero2++]=iPivot;
     1072        }
     1073        iPivot = stack_[iPivot];
     1074      }
     1075    }
     1076    // zero out
     1077    for (i=0;i<numberNonZero2;i++) {
     1078      int k = regionIndex2[i];
     1079      region[k]=0.0;
     1080    }
     1081  } else {
     1082    for (i=0;i<numberNonZero2;i++) {
     1083      int k = regionIndex2[i];
     1084      int j = permute_[k];
     1085      double value = region2[k];
     1086      region2[k]=0.0;
     1087      region[j]=value;
     1088      mark_[j]=1;
     1089      regionIndex[numberNonZero++]=j;
     1090    }
     1091    // copy back
     1092    // set up linked lists at each depth
     1093    // stack2 is start, stack is next
     1094    int greatestDepth=-1;
     1095    int smallestDepth=numberRows_;
     1096    //mark_[numberRows_]=1;
     1097    for (i=0;i<numberNonZero2;i++) {
     1098      int j = regionIndex[i];
     1099      double value = region[j];
     1100      region[j]=0.0;
     1101      region2[j]=value;
     1102      regionIndex2[i]=j;
     1103      // add in
     1104      int iDepth = depth_[j];
     1105      smallestDepth = min(iDepth,smallestDepth) ;
     1106      greatestDepth = max(iDepth,greatestDepth) ;
     1107      int jNext = stack2_[iDepth];
     1108      stack2_[iDepth]=j;
     1109      stack_[j]=jNext;
     1110      // and put all descendants on list
     1111      int iChild = descendant_[j];
     1112      while (iChild>=0) {
     1113        if (!mark_[iChild]) {
     1114          regionIndex2[numberNonZero++] = iChild;
     1115          mark_[iChild]=1;
     1116        }
     1117        iChild = rightSibling_[iChild];
     1118      }
     1119    }
     1120    for (;i<numberNonZero;i++) {
     1121      int j = regionIndex2[i];
     1122      // add in
     1123      int iDepth = depth_[j];
     1124      smallestDepth = min(iDepth,smallestDepth) ;
     1125      greatestDepth = max(iDepth,greatestDepth) ;
     1126      int jNext = stack2_[iDepth];
     1127      stack2_[iDepth]=j;
     1128      stack_[j]=jNext;
     1129      // and put all descendants on list
     1130      int iChild = descendant_[j];
     1131      while (iChild>=0) {
     1132        if (!mark_[iChild]) {
     1133          regionIndex2[numberNonZero++] = iChild;
     1134          mark_[iChild]=1;
     1135        }
     1136        iChild = rightSibling_[iChild];
     1137      }
     1138    }
     1139    numberNonZero2=0;
     1140    region2[numberRows_]=0.0;
     1141    int iDepth;
     1142    for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
     1143      int iPivot = stack2_[iDepth];
     1144      stack2_[iDepth]=-1;
     1145      while (iPivot>=0) {
     1146        mark_[iPivot]=0;
     1147        double pivotValue = region2[iPivot];
     1148        int otherRow = parent_[iPivot];
     1149        double otherValue = region2[otherRow];
     1150        pivotValue = sign_[iPivot]*pivotValue+otherValue;
     1151        region2[iPivot]=pivotValue;
     1152        if (pivotValue)
     1153          regionIndex2[numberNonZero2++]=iPivot;
     1154        iPivot = stack_[iPivot];
     1155      }
    9071156    }
    9081157  }
  • branches/pre/ClpNetworkMatrix.cpp

    r124 r180  
    369369  ClpPlusMinusOneMatrix* rowCopy =
    370370    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
     371  bool packed = rowArray->packedMode();
    371372  if (numberInRowArray>0.3*numberRows||!rowCopy) {
    372373    // do by column
    373374    int iColumn;
    374     double * markVector = y->denseVector(); // probably empty
     375    assert (!y->getNumElements());
    375376    CoinBigIndex j=0;
    376     if (trueNetwork_) {
    377       for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    378         double value = markVector[iColumn];
    379         markVector[iColumn]=0.0;
    380         int iRowM = indices_[j];
    381         int iRowP = indices_[j+1];
    382         value -= scalar*pi[iRowM];
    383         value += scalar*pi[iRowP];
    384         if (fabs(value)>zeroTolerance) {
    385           index[numberNonZero++]=iColumn;
    386           array[iColumn]=value;
    387         }
     377    if (packed) {
     378      // need to expand pi into y
     379      assert(y->capacity()>=numberRows);
     380      double * piOld = pi;
     381      pi = y->denseVector();
     382      const int * whichRow = rowArray->getIndices();
     383      int i;
     384      // modify pi so can collapse to one loop
     385      for (i=0;i<numberInRowArray;i++) {
     386        int iRow = whichRow[i];
     387        pi[iRow]=scalar*piOld[i];
     388      }
     389      if (trueNetwork_) {
     390        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     391          double value = 0.0;
     392          int iRowM = indices_[j];
     393          int iRowP = indices_[j+1];
     394          value -= pi[iRowM];
     395          value += pi[iRowP];
     396          if (fabs(value)>zeroTolerance) {
     397            array[numberNonZero]=value;
     398            index[numberNonZero++]=iColumn;
     399          }
     400        }
     401      } else {
     402        // skip negative rows
     403        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     404          double value = 0.0;
     405          int iRowM = indices_[j];
     406          int iRowP = indices_[j+1];
     407          if (iRowM>=0)
     408            value -= pi[iRowM];
     409          if (iRowP>=0)
     410            value += pi[iRowP];
     411          if (fabs(value)>zeroTolerance) {
     412            array[numberNonZero]=value;
     413            index[numberNonZero++]=iColumn;
     414          }
     415        }
     416      }
     417      for (i=0;i<numberInRowArray;i++) {
     418        int iRow = whichRow[i];
     419        pi[iRow]=0.0;
    388420      }
    389421    } else {
    390       // skip negative rows
    391       for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
    392         double value = markVector[iColumn];
    393         markVector[iColumn]=0.0;
    394         int iRowM = indices_[j];
    395         int iRowP = indices_[j+1];
    396         if (iRowM>=0)
     422      if (trueNetwork_) {
     423        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     424          double value = 0.0;
     425          int iRowM = indices_[j];
     426          int iRowP = indices_[j+1];
    397427          value -= scalar*pi[iRowM];
    398         if (iRowP>=0)
    399428          value += scalar*pi[iRowP];
    400         if (fabs(value)>zeroTolerance) {
    401           index[numberNonZero++]=iColumn;
    402           array[iColumn]=value;
     429          if (fabs(value)>zeroTolerance) {
     430            index[numberNonZero++]=iColumn;
     431            array[iColumn]=value;
     432          }
     433        }
     434      } else {
     435        // skip negative rows
     436        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     437          double value = 0.0;
     438          int iRowM = indices_[j];
     439          int iRowP = indices_[j+1];
     440          if (iRowM>=0)
     441            value -= scalar*pi[iRowM];
     442          if (iRowP>=0)
     443            value += scalar*pi[iRowP];
     444          if (fabs(value)>zeroTolerance) {
     445            index[numberNonZero++]=iColumn;
     446            array[iColumn]=value;
     447          }
    403448        }
    404449      }
    405450    }
    406451    columnArray->setNumElements(numberNonZero);
    407     y->setNumElements(0);
    408452  } else {
    409453    // do by row
     
    430474  int numberToDo = y->getNumElements();
    431475  const int * which = y->getIndices();
    432   if (trueNetwork_) {
    433     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    434       int iColumn = which[jColumn];
    435       double value = 0.0;
    436       CoinBigIndex j=iColumn<<1;
    437       int iRowM = indices_[j];
    438       int iRowP = indices_[j+1];
    439       value -= pi[iRowM];
    440       value += pi[iRowP];
    441       if (fabs(value)>zeroTolerance) {
    442         index[numberNonZero++]=iColumn;
    443         array[iColumn]=value;
    444       }
     476  bool packed = rowArray->packedMode();
     477  if (packed) {
     478    // Must line up with y
     479    // need to expand pi into y
     480    int numberInRowArray = rowArray->getNumElements();
     481    assert(y->capacity()>=model->numberRows());
     482    double * piOld = pi;
     483    pi = y->denseVector();
     484    const int * whichRow = rowArray->getIndices();
     485    int i;
     486    for (i=0;i<numberInRowArray;i++) {
     487      int iRow = whichRow[i];
     488      pi[iRow]=piOld[i];
     489    }
     490    if (trueNetwork_) {
     491      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     492        int iColumn = which[jColumn];
     493        double value = 0.0;
     494        CoinBigIndex j=iColumn<<1;
     495        int iRowM = indices_[j];
     496        int iRowP = indices_[j+1];
     497        value -= pi[iRowM];
     498        value += pi[iRowP];
     499        array[jColumn]=value;
     500      }
     501    } else {
     502      // skip negative rows
     503      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     504        int iColumn = which[jColumn];
     505        double value = 0.0;
     506        CoinBigIndex j=iColumn<<1;
     507        int iRowM = indices_[j];
     508        int iRowP = indices_[j+1];
     509        if (iRowM>=0)
     510          value -= pi[iRowM];
     511        if (iRowP>=0)
     512          value += pi[iRowP];
     513        array[jColumn]=value;
     514      }
     515    }
     516    for (i=0;i<numberInRowArray;i++) {
     517      int iRow = whichRow[i];
     518      pi[iRow]=0.0;
    445519    }
    446520  } else {
    447     // skip negative rows
    448     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    449       int iColumn = which[jColumn];
    450       double value = 0.0;
    451       CoinBigIndex j=iColumn<<1;
    452       int iRowM = indices_[j];
    453       int iRowP = indices_[j+1];
    454       if (iRowM>=0)
     521    if (trueNetwork_) {
     522      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     523        int iColumn = which[jColumn];
     524        double value = 0.0;
     525        CoinBigIndex j=iColumn<<1;
     526        int iRowM = indices_[j];
     527        int iRowP = indices_[j+1];
    455528        value -= pi[iRowM];
    456       if (iRowP>=0)
    457529        value += pi[iRowP];
    458       if (fabs(value)>zeroTolerance) {
    459         index[numberNonZero++]=iColumn;
    460         array[iColumn]=value;
     530        if (fabs(value)>zeroTolerance) {
     531          index[numberNonZero++]=iColumn;
     532          array[iColumn]=value;
     533        }
     534      }
     535    } else {
     536      // skip negative rows
     537      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     538        int iColumn = which[jColumn];
     539        double value = 0.0;
     540        CoinBigIndex j=iColumn<<1;
     541        int iRowM = indices_[j];
     542        int iRowP = indices_[j+1];
     543        if (iRowM>=0)
     544          value -= pi[iRowM];
     545        if (iRowP>=0)
     546          value += pi[iRowP];
     547        if (fabs(value)>zeroTolerance) {
     548          index[numberNonZero++]=iColumn;
     549          array[iColumn]=value;
     550        }
    461551      }
    462552    }
     
    542632  return numberElements;
    543633}
     634/* If element NULL returns number of elements in column part of basis,
     635   If not NULL fills in as well */
     636CoinBigIndex
     637ClpNetworkMatrix::fillBasis(const ClpSimplex * model,
     638                                 const int * whichColumn,
     639                                 int numberBasic,
     640                                 int numberColumnBasic,
     641                                 int * indexRowU, int * indexColumnU,
     642                                 double * elementU) const
     643{
     644  int i;
     645  CoinBigIndex numberElements=0;
     646  if (elementU!=NULL) {
     647    if (trueNetwork_) {
     648      for (i=0;i<numberColumnBasic;i++) {
     649        int iColumn = whichColumn[i];
     650        CoinBigIndex j=iColumn<<1;
     651        int iRowM = indices_[j];
     652        int iRowP = indices_[j+1];
     653        indexRowU[numberElements]=iRowM;
     654        indexColumnU[numberElements]=numberBasic;
     655        elementU[numberElements]=-1.0;
     656        indexRowU[numberElements+1]=iRowP;
     657        indexColumnU[numberElements+1]=numberBasic;
     658        elementU[numberElements+1]=1.0;
     659        numberElements+=2;
     660        numberBasic++;
     661      }
     662    } else {
     663      for (i=0;i<numberColumnBasic;i++) {
     664        int iColumn = whichColumn[i];
     665        CoinBigIndex j=iColumn<<1;
     666        int iRowM = indices_[j];
     667        int iRowP = indices_[j+1];
     668        if (iRowM>=0) {
     669          indexRowU[numberElements]=iRowM;
     670          indexColumnU[numberElements]=numberBasic;
     671          elementU[numberElements++]=-1.0;
     672        }
     673        if (iRowP>=0) {
     674          indexRowU[numberElements]=iRowP;
     675          indexColumnU[numberElements]=numberBasic;
     676          elementU[numberElements++]=1.0;
     677        }
     678        numberBasic++;
     679      }
     680    }
     681  } else {
     682    if (trueNetwork_) {
     683      numberElements = 2*numberColumnBasic;
     684    } else {
     685      for (i=0;i<numberColumnBasic;i++) {
     686        int iColumn = whichColumn[i];
     687        CoinBigIndex j=iColumn<<1;
     688        int iRowM = indices_[j];
     689        int iRowP = indices_[j+1];
     690        if (iRowM>=0)
     691          numberElements ++;
     692        if (iRowP>=0)
     693          numberElements ++;
     694      }
     695    }
     696  }
     697  return numberElements;
     698}
    544699/* Unpacks a column into an CoinIndexedvector
    545       Note that model is NOT const.  Bounds and objective could
    546       be modified if doing column generation */
     700 */
    547701void
    548702ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
     
    557711    rowArray->add(iRowP,1.0);
    558712}
     713/* Unpacks a column into an CoinIndexedvector
     714** in packed foramt
     715Note that model is NOT const.  Bounds and objective could
     716be modified if doing column generation (just for this variable) */
     717void
     718ClpNetworkMatrix::unpackPacked(ClpSimplex * model,
     719                            CoinIndexedVector * rowArray,
     720                            int iColumn) const
     721{
     722  int * index = rowArray->getIndices();
     723  double * array = rowArray->denseVector();
     724  int number = 0;
     725  CoinBigIndex j=iColumn<<1;
     726  int iRowM = indices_[j];
     727  int iRowP = indices_[j+1];
     728  if (iRowM>=0) {
     729    array[number]=-1.0;
     730    index[number++]=iRowM;
     731  }
     732  if (iRowP>=0) {
     733    array[number]=1.0;
     734    index[number++]=iRowP;
     735  }
     736  rowArray->setNumElements(number);
     737  rowArray->setPackedMode(true);
     738}
    559739/* Adds multiple of a column into an CoinIndexedvector
    560740      You can use quickAdd to add to vector */
  • branches/pre/ClpNonLinearCost.cpp

    r170 r180  
    364364// purpose which will come to me later
    365365void
    366 ClpNonLinearCost::checkInfeasibilities(bool toNearest)
     366ClpNonLinearCost::checkInfeasibilities(double oldTolerance)
    367367{
    368368  numberInfeasibilities_=0;
     
    372372  sumInfeasibilities_ = 0.0;
    373373  double primalTolerance = model_->currentPrimalTolerance();
    374  
    375374  int iSequence;
    376375  double * solution = model_->solutionRegion();
     
    378377  double * lower = model_->lowerRegion();
    379378  double * cost = model_->costRegion();
     379  bool toNearest = oldTolerance<=0.0;
    380380   
    381381  // nonbasic should be at a valid bound
     
    454454      if (!toNearest) {
    455455        // With increasing tolerances - we may be at wrong place
    456         if (fabs(value-upperValue)>primalTolerance*1.0001) {
    457           assert(fabs(value-lowerValue)<=primalTolerance*1.0001);
     456        if (fabs(value-upperValue)>oldTolerance*1.0001) {
     457          assert(fabs(value-lowerValue)<=oldTolerance*1.0001);
     458          if  (fabs(value-lowerValue)>primalTolerance)
     459            solution[iSequence]=lowerValue;
    458460          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
     461        } else if  (fabs(value-upperValue)>primalTolerance) {
     462          solution[iSequence]=upperValue;
    459463        }
    460464      } else {
     
    470474      if (!toNearest) {
    471475        // With increasing tolerances - we may be at wrong place
    472         if (fabs(value-lowerValue)>primalTolerance*1.0001) {
    473           assert(fabs(value-upperValue)<=primalTolerance*1.0001);
    474           model_->setStatus(iSequence,ClpSimplex::atUpperBound);
     476        if (fabs(value-lowerValue)>oldTolerance*1.0001) {
     477          assert(fabs(value-upperValue)<=oldTolerance*1.0001);
     478          if  (fabs(value-upperValue)>primalTolerance)
     479            solution[iSequence]=upperValue;
     480          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
     481        } else if  (fabs(value-lowerValue)>primalTolerance) {
     482          solution[iSequence]=lowerValue;
    475483        }
    476484      } else {
  • branches/pre/ClpObjective.cpp

    r124 r180  
    4747  return *this;
    4848}
     49/* Subset clone.  Duplicates are allowed
     50   and order is as given.
     51*/
     52ClpObjective *
     53ClpObjective::subsetClone (int numberColumns,
     54                           const int * whichColumns) const
     55{
     56  std::cerr<<"subsetClone not supported - ClpObjective"<<std::endl;
     57  abort();
     58}
     59
  • branches/pre/ClpPackedMatrix.cpp

    r145 r180  
    9090{
    9191  return new ClpPackedMatrix(*this);
     92}
     93/* Subset clone (without gaps).  Duplicates are allowed
     94   and order is as given */
     95ClpMatrixBase *
     96ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
     97                              int numberColumns,
     98                              const int * whichColumns) const
     99{
     100  return new ClpPackedMatrix(*this, numberRows, whichRows,
     101                                   numberColumns, whichColumns);
     102}
     103/* Subset constructor (without gaps).  Duplicates are allowed
     104   and order is as given */
     105ClpPackedMatrix::ClpPackedMatrix (
     106                       const ClpPackedMatrix & rhs,
     107                       int numberRows, const int * whichRows,
     108                       int numberColumns, const int * whichColumns)
     109: ClpMatrixBase(rhs)
     110{
     111  matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
     112                                 numberColumns,whichColumns);
     113  zeroElements_ = rhs.zeroElements_;
     114}
     115ClpPackedMatrix::ClpPackedMatrix (
     116                       const CoinPackedMatrix & rhs,
     117                       int numberRows, const int * whichRows,
     118                       int numberColumns, const int * whichColumns)
     119: ClpMatrixBase()
     120{
     121  matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
     122                                 numberColumns,whichColumns);
     123  zeroElements_ = false;
     124  setType(1);
    92125}
    93126
     
    227260  ClpPackedMatrix* rowCopy =
    228261    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
     262  bool packed = rowArray->packedMode();
    229263  if (numberInRowArray>0.333*numberRows||!rowCopy) {
    230264    // do by column
     
    238272    int numberColumns = model->numberColumns();
    239273    if (!y->getNumElements()) {
    240       if (!rowScale) {
    241         if (scalar==-1.0) {
     274      if (packed) {
     275        // need to expand pi into y
     276        assert(y->capacity()>=numberRows);
     277        double * piOld = pi;
     278        pi = y->denseVector();
     279        const int * whichRow = rowArray->getIndices();
     280        int i;
     281        if (!rowScale) {
     282          // modify pi so can collapse to one loop
     283          for (i=0;i<numberInRowArray;i++) {
     284            int iRow = whichRow[i];
     285            pi[iRow]=scalar*piOld[i];
     286          }
    242287          for (iColumn=0;iColumn<numberColumns;iColumn++) {
    243288            double value = 0.0;
     
    249294            }
    250295            if (fabs(value)>zeroTolerance) {
     296              array[numberNonZero]=value;
    251297              index[numberNonZero++]=iColumn;
    252               array[iColumn]=-value;
    253             }
    254           }
    255         } else if (scalar==1.0) {
    256           for (iColumn=0;iColumn<numberColumns;iColumn++) {
    257             double value = 0.0;
    258             CoinBigIndex j;
    259             for (j=columnStart[iColumn];
    260                  j<columnStart[iColumn]+columnLength[iColumn];j++) {
    261               int iRow = row[j];
    262               value += pi[iRow]*elementByColumn[j];
    263             }
    264             if (fabs(value)>zeroTolerance) {
    265               index[numberNonZero++]=iColumn;
    266               array[iColumn]=value;
    267298            }
    268299          }
    269300        } else {
    270           for (iColumn=0;iColumn<numberColumns;iColumn++) {
    271             double value = 0.0;
    272             CoinBigIndex j;
    273             for (j=columnStart[iColumn];
    274                  j<columnStart[iColumn]+columnLength[iColumn];j++) {
    275               int iRow = row[j];
    276               value += pi[iRow]*elementByColumn[j];
    277             }
    278             value *= scalar;
    279             if (fabs(value)>zeroTolerance) {
    280               index[numberNonZero++]=iColumn;
    281               array[iColumn]=value;
    282             }
    283           }
    284         }
    285       } else {
    286         // scaled
    287         if (scalar==-1.0) {
     301          // scaled
     302          // modify pi so can collapse to one loop
     303          for (i=0;i<numberInRowArray;i++) {
     304            int iRow = whichRow[i];
     305            pi[iRow]=scalar*piOld[i]*rowScale[iRow];
     306          }
    288307          for (iColumn=0;iColumn<numberColumns;iColumn++) {
    289308            double value = 0.0;
     
    293312                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
    294313              int iRow = row[j];
    295               value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     314              value += pi[iRow]*elementByColumn[j];
    296315            }
    297316            value *= columnScale[iColumn];
    298317            if (fabs(value)>zeroTolerance) {
     318              array[numberNonZero]=value;
    299319              index[numberNonZero++]=iColumn;
    300               array[iColumn]=-value;
    301             }
    302           }
    303         } else if (scalar==1.0) {
    304           for (iColumn=0;iColumn<numberColumns;iColumn++) {
    305             double value = 0.0;
    306             CoinBigIndex j;
    307             const double * columnScale = model->columnScale();
    308             for (j=columnStart[iColumn];
    309                  j<columnStart[iColumn]+columnLength[iColumn];j++) {
    310               int iRow = row[j];
    311               value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    312             }
    313             value *= columnScale[iColumn];
    314             if (fabs(value)>zeroTolerance) {
    315               index[numberNonZero++]=iColumn;
    316               array[iColumn]=value;
     320            }
     321          }
     322        }
     323        // zero out
     324        for (i=0;i<numberInRowArray;i++) {
     325          int iRow = whichRow[i];
     326          pi[iRow]=0.0;
     327        }
     328      } else {
     329        if (!rowScale) {
     330          if (scalar==-1.0) {
     331            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     332              double value = 0.0;
     333              CoinBigIndex j;
     334              for (j=columnStart[iColumn];
     335                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     336                int iRow = row[j];
     337                value += pi[iRow]*elementByColumn[j];
     338              }
     339              if (fabs(value)>zeroTolerance) {
     340                index[numberNonZero++]=iColumn;
     341                array[iColumn]=-value;
     342              }
     343            }
     344          } else if (scalar==1.0) {
     345            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     346              double value = 0.0;
     347              CoinBigIndex j;
     348              for (j=columnStart[iColumn];
     349                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     350                int iRow = row[j];
     351                value += pi[iRow]*elementByColumn[j];
     352              }
     353              if (fabs(value)>zeroTolerance) {
     354                index[numberNonZero++]=iColumn;
     355                array[iColumn]=value;
     356              }
     357            }
     358          } else {
     359            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     360              double value = 0.0;
     361              CoinBigIndex j;
     362              for (j=columnStart[iColumn];
     363                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     364                int iRow = row[j];
     365                value += pi[iRow]*elementByColumn[j];
     366              }
     367              value *= scalar;
     368              if (fabs(value)>zeroTolerance) {
     369                index[numberNonZero++]=iColumn;
     370                array[iColumn]=value;
     371              }
    317372            }
    318373          }
    319374        } else {
    320           for (iColumn=0;iColumn<numberColumns;iColumn++) {
    321             double value = 0.0;
    322             CoinBigIndex j;
    323             const double * columnScale = model->columnScale();
    324             for (j=columnStart[iColumn];
    325                  j<columnStart[iColumn]+columnLength[iColumn];j++) {
    326               int iRow = row[j];
    327               value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    328             }
    329             value *= scalar*columnScale[iColumn];
    330             if (fabs(value)>zeroTolerance) {
    331               index[numberNonZero++]=iColumn;
    332               array[iColumn]=value;
     375          // scaled
     376          if (scalar==-1.0) {
     377            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     378              double value = 0.0;
     379              CoinBigIndex j;
     380              const double * columnScale = model->columnScale();
     381              for (j=columnStart[iColumn];
     382                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     383                int iRow = row[j];
     384                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     385              }
     386              value *= columnScale[iColumn];
     387              if (fabs(value)>zeroTolerance) {
     388                index[numberNonZero++]=iColumn;
     389                array[iColumn]=-value;
     390              }
     391            }
     392          } else if (scalar==1.0) {
     393            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     394              double value = 0.0;
     395              CoinBigIndex j;
     396              const double * columnScale = model->columnScale();
     397              for (j=columnStart[iColumn];
     398                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     399                int iRow = row[j];
     400                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     401              }
     402              value *= columnScale[iColumn];
     403              if (fabs(value)>zeroTolerance) {
     404                index[numberNonZero++]=iColumn;
     405                array[iColumn]=value;
     406              }
     407            }
     408          } else {
     409            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     410              double value = 0.0;
     411              CoinBigIndex j;
     412              const double * columnScale = model->columnScale();
     413              for (j=columnStart[iColumn];
     414                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     415                int iRow = row[j];
     416                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     417              }
     418              value *= scalar*columnScale[iColumn];
     419              if (fabs(value)>zeroTolerance) {
     420                index[numberNonZero++]=iColumn;
     421                array[iColumn]=value;
     422              }
    333423            }
    334424          }
     
    336426      }
    337427    } else {
     428      assert(!packed);
    338429      double * markVector = y->denseVector(); // not empty
    339430      if (!rowScale) {
     
    381472    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
    382473  }
     474  if (packed)
     475    columnArray->setPackedMode(true);
     476  if (0) {
     477    columnArray->checkClean();
     478    int numberNonZero=columnArray->getNumElements();;
     479    int * index = columnArray->getIndices();
     480    double * array = columnArray->denseVector();
     481    int i;
     482    for (i=0;i<numberNonZero;i++) {
     483      int j=index[i];
     484      double value;
     485      if (packed)
     486        value=array[i];
     487      else
     488        value=array[j];
     489      printf("Ti %d %d %g\n",i,j,value);
     490    }
     491  }
    383492}
    384493/* Return <code>x * A + y</code> in <code>z</code>.
     
    402511  const double * element = getElements();
    403512  const int * whichRow = rowArray->getIndices();
     513  bool packed = rowArray->packedMode();
    404514  if (numberInRowArray>2||y->getNumElements()) {
    405515    // do by rows
    406516    // ** Row copy is already scaled
    407517    int iRow;
    408     double * markVector = y->denseVector(); // probably empty .. but
    409518    int * mark = y->getIndices();
    410519    int numberOriginal=y->getNumElements();
    411520    int i;
    412     for (i=0;i<numberOriginal;i++) {
    413       int iColumn = mark[i];
    414       index[i]=iColumn;
    415       array[iColumn]=markVector[iColumn];
    416       markVector[iColumn]=0.0;
    417     }
    418     numberNonZero=numberOriginal;
    419     // and set up mark as char array
    420     char * marked = (char *) markVector;
    421     for (i=0;i<numberOriginal;i++) {
    422       int iColumn = index[i];
    423       marked[iColumn]=0;
    424     }
    425 
    426     for (i=0;i<numberInRowArray;i++) {
    427       iRow = whichRow[i];
    428       double value = pi[iRow]*scalar;
    429       CoinBigIndex j;
     521    if (packed) {
     522      assert(!numberOriginal);
     523      numberNonZero=0;
     524      // and set up mark as char array
     525      char * marked = (char *) (index+columnArray->capacity());
     526      double * array2 = y->denseVector();
     527#ifdef CLP_DEBUG
     528      int numberColumns = model->numberColumns();
     529      for (i=0;i<numberColumns;i++) {
     530        assert(!marked[i]);
     531        assert(!array2[i]);
     532      }
     533#endif     
     534      for (i=0;i<numberInRowArray;i++) {
     535        iRow = whichRow[i];
     536        double value = pi[i]*scalar;
     537        CoinBigIndex j;
     538        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     539          int iColumn = column[j];
     540          if (!marked[iColumn]) {
     541            marked[iColumn]=1;
     542            index[numberNonZero++]=iColumn;
     543          }
     544          array2[iColumn] += value*element[j];
     545        }
     546      }
     547      // get rid of tiny values and zero out marked
     548      numberOriginal=numberNonZero;
     549      numberNonZero=0;
     550      for (i=0;i<numberOriginal;i++) {
     551        int iColumn = index[i];
     552        if (marked[iColumn]) {
     553          double value = array2[iColumn];
     554          array2[iColumn]=0.0;
     555          marked[iColumn]=0;
     556          if (fabs(value)>zeroTolerance) {
     557            array[numberNonZero]=value;
     558            index[numberNonZero++]=iColumn;
     559          }
     560        }
     561      }
     562    } else {
     563      double * markVector = y->denseVector(); // probably empty .. but
     564      for (i=0;i<numberOriginal;i++) {
     565        int iColumn = mark[i];
     566        index[i]=iColumn;
     567        array[iColumn]=markVector[iColumn];
     568        markVector[iColumn]=0.0;
     569      }
     570      numberNonZero=numberOriginal;
     571      // and set up mark as char array
     572      char * marked = (char *) markVector;
     573      for (i=0;i<numberOriginal;i++) {
     574        int iColumn = index[i];
     575        marked[iColumn]=0;
     576      }
     577     
     578      for (i=0;i<numberInRowArray;i++) {
     579        iRow = whichRow[i];
     580        double value = pi[iRow]*scalar;
     581        CoinBigIndex j;
     582        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     583          int iColumn = column[j];
     584          if (!marked[iColumn]) {
     585            marked[iColumn]=1;
     586            index[numberNonZero++]=iColumn;
     587          }
     588          array[iColumn] += value*element[j];
     589        }
     590      }
     591      // get rid of tiny values and zero out marked
     592      numberOriginal=numberNonZero;
     593      numberNonZero=0;
     594      for (i=0;i<numberOriginal;i++) {
     595        int iColumn = index[i];
     596        marked[iColumn]=0;
     597        if (fabs(array[iColumn])>zeroTolerance) {
     598          index[numberNonZero++]=iColumn;
     599        } else {
     600          array[iColumn]=0.0;
     601        }
     602      }
     603    }
     604  } else if (numberInRowArray==2) {
     605    // do by rows when two rows
     606    int numberOriginal;
     607    int i;
     608    CoinBigIndex j;
     609    numberNonZero=0;
     610
     611    double value;
     612    if (packed) {
     613      int iRow0 = whichRow[0];
     614      int iRow1 = whichRow[1];
     615      double pi0 = pi[0];
     616      double pi1 = pi[1];
     617      if (rowStart[iRow0+1]-rowStart[iRow0]>
     618          rowStart[iRow1+1]-rowStart[iRow1]) {
     619        // do one with fewer first
     620        iRow0=iRow1;
     621        iRow1=whichRow[0];
     622        pi0=pi1;
     623        pi1=pi[0];
     624      }
     625      // and set up mark as char array
     626      char * marked = (char *) (index+columnArray->capacity());
     627      int * lookup = y->getIndices();
     628      value = pi0*scalar;
     629      for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
     630        int iColumn = column[j];
     631        double value2 = value*element[j];
     632        array[numberNonZero] = value2;
     633        marked[iColumn]=1;
     634        lookup[iColumn]=numberNonZero;
     635        index[numberNonZero++]=iColumn;
     636      }
     637      numberOriginal = numberNonZero;
     638      value = pi1*scalar;
     639      for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
     640        int iColumn = column[j];
     641        double value2 = value*element[j];
     642        // I am assuming no zeros in matrix
     643        if (marked[iColumn]) {
     644          int iLookup = lookup[iColumn];
     645          array[iLookup] += value2;
     646        } else {
     647          if (fabs(value2)>zeroTolerance) {
     648            array[numberNonZero] = value2;
     649            index[numberNonZero++]=iColumn;
     650          }
     651        }
     652      }
     653      // get rid of tiny values and zero out marked
     654      int nDelete=0;
     655      for (i=0;i<numberOriginal;i++) {
     656        int iColumn = index[i];
     657        marked[iColumn]=0;
     658        if (fabs(array[i])<=zeroTolerance)
     659          nDelete++;
     660      }
     661      if (nDelete) {
     662        numberOriginal=numberNonZero;
     663        numberNonZero=0;
     664        for (i=0;i<numberOriginal;i++) {
     665          int iColumn = index[i];
     666          double value = array[i];
     667          array[i]=0.0;
     668          if (fabs(value)>zeroTolerance) {
     669            array[numberNonZero]=value;
     670            index[numberNonZero++]=iColumn;
     671          }
     672        }
     673      }
     674    } else {
     675      int iRow = whichRow[0];
     676      value = pi[iRow]*scalar;
    430677      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    431678        int iColumn = column[j];
    432         if (!marked[iColumn]) {
    433           marked[iColumn]=1;
     679        double value2 = value*element[j];
     680        index[numberNonZero++]=iColumn;
     681        array[iColumn] = value2;
     682      }
     683      iRow = whichRow[1];
     684      value = pi[iRow]*scalar;
     685      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     686        int iColumn = column[j];
     687        double value2 = value*element[j];
     688        // I am assuming no zeros in matrix
     689        if (array[iColumn])
     690          value2 += array[iColumn];
     691        else
    434692          index[numberNonZero++]=iColumn;
    435         }
    436         array[iColumn] += value*element[j];
    437       }
    438     }
    439     // get rid of tiny values and zero out marked
    440     numberOriginal=numberNonZero;
    441     numberNonZero=0;
    442     for (i=0;i<numberOriginal;i++) {
    443       int iColumn = index[i];
    444       marked[iColumn]=0;
    445       if (fabs(array[iColumn])>zeroTolerance) {
    446         index[numberNonZero++]=iColumn;
    447       } else {
    448         array[iColumn]=0.0;
    449       }
    450     }
    451   } else if (numberInRowArray==2) {
    452     // do by rows when two rows
    453     int iRow;
    454     int numberOriginal;
    455     int i;
    456     numberNonZero=0;
    457 
    458     double value;
    459     iRow = whichRow[0];
    460     value = pi[iRow]*scalar;
    461     CoinBigIndex j;
    462     for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    463       int iColumn = column[j];
    464       double value2 = value*element[j];
    465       index[numberNonZero++]=iColumn;
    466       array[iColumn] = value2;
    467     }
    468     iRow = whichRow[1];
    469     value = pi[iRow]*scalar;
    470     for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    471       int iColumn = column[j];
    472       double value2 = value*element[j];
    473       // I am assuming no zeros in matrix
    474       if (array[iColumn])
    475         value2 += array[iColumn];
    476       else
    477         index[numberNonZero++]=iColumn;
    478       array[iColumn] = value2;
    479     }
    480     // get rid of tiny values and zero out marked
    481     numberOriginal=numberNonZero;
    482     numberNonZero=0;
    483     for (i=0;i<numberOriginal;i++) {
    484       int iColumn = index[i];
    485       if (fabs(array[iColumn])>zeroTolerance) {
    486         index[numberNonZero++]=iColumn;
    487       } else {
    488         array[iColumn]=0.0;
     693        array[iColumn] = value2;
     694      }
     695      // get rid of tiny values and zero out marked
     696      numberOriginal=numberNonZero;
     697      numberNonZero=0;
     698      for (i=0;i<numberOriginal;i++) {
     699        int iColumn = index[i];
     700        if (fabs(array[iColumn])>zeroTolerance) {
     701          index[numberNonZero++]=iColumn;
     702        } else {
     703          array[iColumn]=0.0;
     704        }
    489705      }
    490706    }
     
    493709    int iRow=rowArray->getIndices()[0];
    494710    numberNonZero=0;
    495     double value = pi[iRow]*scalar;
    496711    CoinBigIndex j;
    497     for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    498       int iColumn = column[j];
    499       double value2 = value*element[j];
    500       if (fabs(value2)>zeroTolerance) {
    501         index[numberNonZero++]=iColumn;
    502         array[iColumn] = value2;
     712    if (packed) {
     713      double value = pi[0]*scalar;
     714      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     715        int iColumn = column[j];
     716        double value2 = value*element[j];
     717        if (fabs(value2)>zeroTolerance) {
     718          array[numberNonZero] = value2;
     719          index[numberNonZero++]=iColumn;
     720        }
     721      }
     722    } else {
     723      double value = pi[iRow]*scalar;
     724      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     725        int iColumn = column[j];
     726        double value2 = value*element[j];
     727        if (fabs(value2)>zeroTolerance) {
     728          index[numberNonZero++]=iColumn;
     729          array[iColumn] = value2;
     730        }
    503731      }
    504732    }
     
    532760  int numberToDo = y->getNumElements();
    533761  const int * which = y->getIndices();
    534   if (!rowScale) {
    535     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    536       int iColumn = which[jColumn];
    537       double value = 0.0;
    538       CoinBigIndex j;
    539       for (j=columnStart[iColumn];
    540            j<columnStart[iColumn]+columnLength[iColumn];j++) {
    541         int iRow = row[j];
    542         value += pi[iRow]*elementByColumn[j];
    543       }
    544       if (fabs(value)>zeroTolerance) {
    545         index[numberNonZero++]=iColumn;
    546         array[iColumn]=value;
    547       }
     762  bool packed = rowArray->packedMode();
     763  if (packed) {
     764    // need to expand pi into y
     765    int numberInRowArray = rowArray->getNumElements();
     766    assert(y->capacity()>=model->numberRows());
     767    double * piOld = pi;
     768    pi = y->denseVector();
     769    const int * whichRow = rowArray->getIndices();
     770    int i;
     771    // Do NOT squash small elements - must line up with  y
     772    if (!rowScale) {
     773      for (i=0;i<numberInRowArray;i++) {
     774        int iRow = whichRow[i];
     775        pi[iRow]=piOld[i];
     776      }
     777      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     778        int iColumn = which[jColumn];
     779        double value = 0.0;
     780        CoinBigIndex j;
     781        for (j=columnStart[iColumn];
     782             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     783          int iRow = row[j];
     784          value += pi[iRow]*elementByColumn[j];
     785        }
     786        array[jColumn]=value;
     787      }
     788    } else {
     789      // scaled
     790      for (i=0;i<numberInRowArray;i++) {
     791        int iRow = whichRow[i];
     792        pi[iRow]=rowScale[iRow]*piOld[i];
     793      }
     794      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     795        int iColumn = which[jColumn];
     796        double value = 0.0;
     797        CoinBigIndex j;
     798        const double * columnScale = model->columnScale();
     799        for (j=columnStart[iColumn];
     800             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     801          int iRow = row[j];
     802          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     803        }
     804        value *= columnScale[iColumn];
     805        array[jColumn]=value;
     806      }
     807    }
     808    for (i=0;i<numberInRowArray;i++) {
     809      int iRow = whichRow[i];
     810      pi[iRow]=0.0;
    548811    }
    549812  } else {
    550     // scaled
    551     for (jColumn=0;jColumn<numberToDo;jColumn++) {
    552       int iColumn = which[jColumn];
    553       double value = 0.0;
    554       CoinBigIndex j;
    555       const double * columnScale = model->columnScale();
    556       for (j=columnStart[iColumn];
    557            j<columnStart[iColumn]+columnLength[iColumn];j++) {
    558         int iRow = row[j];
    559         value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
    560       }
    561       value *= columnScale[iColumn];
    562       if (fabs(value)>zeroTolerance) {
    563         index[numberNonZero++]=iColumn;
    564         array[iColumn]=value;
     813    if (!rowScale) {
     814      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     815        int iColumn = which[jColumn];
     816        double value = 0.0;
     817        CoinBigIndex j;
     818        for (j=columnStart[iColumn];
     819             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     820          int iRow = row[j];
     821          value += pi[iRow]*elementByColumn[j];
     822        }
     823        if (fabs(value)>zeroTolerance) {
     824          index[numberNonZero++]=iColumn;
     825          array[iColumn]=value;
     826        }
     827      }
     828    } else {
     829      // scaled
     830      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     831        int iColumn = which[jColumn];
     832        double value = 0.0;
     833        CoinBigIndex j;
     834        const double * columnScale = model->columnScale();
     835        for (j=columnStart[iColumn];
     836             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     837          int iRow = row[j];
     838          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     839        }
     840        value *= columnScale[iColumn];
     841        if (fabs(value)>zeroTolerance) {
     842          index[numberNonZero++]=iColumn;
     843          array[iColumn]=value;
     844        }
    565845      }
    566846    }
     
    681961        }
    682962      }
     963    }
     964  }
     965  return numberElements;
     966}
     967/* If element NULL returns number of elements in column part of basis,
     968   If not NULL fills in as well */
     969CoinBigIndex
     970ClpPackedMatrix::fillBasis(const ClpSimplex * model,
     971                           const int * whichColumn,
     972                           int numberBasic,
     973                           int numberColumnBasic,
     974                           int * indexRowU, int * indexColumnU,
     975                           double * elementU) const
     976{
     977  const int * columnLength = matrix_->getVectorLengths();
     978  int i;
     979  CoinBigIndex numberElements=0;
     980  if (elementU!=NULL) {
     981    // fill
     982    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     983    const double * rowScale = model->rowScale();
     984    const int * row = matrix_->getIndices();
     985    const double * elementByColumn = matrix_->getElements();
     986    if (!zeroElements_) {
     987      if (!rowScale) {
     988        // no scaling
     989        for (i=0;i<numberColumnBasic;i++) {
     990          int iColumn = whichColumn[i];
     991          CoinBigIndex j;
     992          for (j=columnStart[iColumn];
     993               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     994            indexRowU[numberElements]=row[j];
     995            indexColumnU[numberElements]=numberBasic;
     996            elementU[numberElements++]=elementByColumn[j];
     997          }
     998          numberBasic++;
     999        }
     1000      } else {
     1001        // scaling
     1002        const double * columnScale = model->columnScale();
     1003        for (i=0;i<numberColumnBasic;i++) {
     1004          int iColumn = whichColumn[i];
     1005          CoinBigIndex j;
     1006          double scale = columnScale[iColumn];
     1007          for (j=columnStart[iColumn];
     1008               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1009            int iRow = row[j];
     1010            indexRowU[numberElements]=iRow;
     1011            indexColumnU[numberElements]=numberBasic;
     1012            elementU[numberElements++]=
     1013              elementByColumn[j]*scale*rowScale[iRow];
     1014          }
     1015          numberBasic++;
     1016        }
     1017      }
     1018    } else {
     1019      // there are zero elements so need to look more closely
     1020      if (!rowScale) {
     1021        // no scaling
     1022        for (i=0;i<numberColumnBasic;i++) {
     1023          int iColumn = whichColumn[i];
     1024          CoinBigIndex j;
     1025          for (j=columnStart[iColumn];
     1026               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1027            double value = elementByColumn[j];
     1028            if (value) {
     1029              indexRowU[numberElements]=row[j];
     1030              indexColumnU[numberElements]=numberBasic;
     1031              elementU[numberElements++]=value;
     1032            }
     1033          }
     1034          numberBasic++;
     1035        }
     1036      } else {
     1037        // scaling
     1038        const double * columnScale = model->columnScale();
     1039        for (i=0;i<numberColumnBasic;i++) {
     1040          int iColumn = whichColumn[i];
     1041          CoinBigIndex j;
     1042          double scale = columnScale[iColumn];
     1043          for (j=columnStart[iColumn];
     1044               j<columnStart[iColumn]+columnLength[i];j++) {
     1045            double value = elementByColumn[j];
     1046            if (value) {
     1047              int iRow = row[j];
     1048              indexRowU[numberElements]=iRow;
     1049              indexColumnU[numberElements]=numberBasic;
     1050              elementU[numberElements++]=value*scale*rowScale[iRow];
     1051            }
     1052          }
     1053          numberBasic++;
     1054        }
     1055      }
     1056    }
     1057  } else {
     1058    // just count - can be over so ignore zero problem
     1059    for (i=0;i<numberColumnBasic;i++) {
     1060      int iColumn = whichColumn[i];
     1061      numberElements += columnLength[iColumn];
    6831062    }
    6841063  }
     
    8931272}
    8941273/* Unpacks a column into an CoinIndexedvector
    895       Note that model is NOT const.  Bounds and objective could
    896       be modified if doing column generation */
     1274 */
    8971275void
    8981276ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
     
    9181296      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
    9191297    }
     1298  }
     1299}
     1300/* Unpacks a column into an CoinIndexedvector
     1301** in packed foramt
     1302Note that model is NOT const.  Bounds and objective could
     1303be modified if doing column generation (just for this variable) */
     1304void
     1305ClpPackedMatrix::unpackPacked(ClpSimplex * model,
     1306                            CoinIndexedVector * rowArray,
     1307                            int iColumn) const
     1308{
     1309  const double * rowScale = model->rowScale();
     1310  const int * row = matrix_->getIndices();
     1311  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1312  const int * columnLength = matrix_->getVectorLengths();
     1313  const double * elementByColumn = matrix_->getElements();
     1314  CoinBigIndex i;
     1315  if (!rowScale) {
     1316    int j=columnStart[iColumn];
     1317    rowArray->createPacked(columnLength[iColumn],
     1318                           row+j,elementByColumn+j);
     1319  } else {
     1320    // apply scaling
     1321    double scale = model->columnScale()[iColumn];
     1322    int * index = rowArray->getIndices();
     1323    double * array = rowArray->denseVector();
     1324    int number = 0;
     1325    for (i=columnStart[iColumn];
     1326         i<columnStart[iColumn]+columnLength[iColumn];i++) {
     1327      int iRow = row[i];
     1328      array[number]=elementByColumn[i]*scale*rowScale[iRow];
     1329      index[number++]=iRow;
     1330    }
     1331    rowArray->setNumElements(number);
     1332    rowArray->setPackedMode(true);
    9201333  }
    9211334}
  • branches/pre/ClpPlusMinusOneMatrix.cpp

    r124 r180  
    5151  columnOrdered_=rhs.columnOrdered_;
    5252  if (numberColumns_) {
    53     indices_ = new int [ 2*numberColumns_];
    54     memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
     53    int numberElements = rhs.startPositive_[numberColumns_];
     54    indices_ = new int [ numberElements];
     55    memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
    5556    startPositive_ = new int [ numberColumns_+1];
    5657    memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
     
    165166    columnOrdered_=rhs.columnOrdered_;
    166167    if (numberColumns_) {
    167       indices_ = new int [ 2*numberColumns_];
    168       memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
     168      int numberElements = rhs.startPositive_[numberColumns_];
     169      indices_ = new int [ numberElements];
     170      memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
    169171      startPositive_ = new int [ numberColumns_+1];
    170172      memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
     
    182184  return new ClpPlusMinusOneMatrix(*this);
    183185}
     186/* Subset clone (without gaps).  Duplicates are allowed
     187   and order is as given */
     188ClpMatrixBase *
     189ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
     190                              int numberColumns,
     191                              const int * whichColumns) const
     192{
     193  return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
     194                                   numberColumns, whichColumns);
     195}
     196/* Subset constructor (without gaps).  Duplicates are allowed
     197   and order is as given */
     198ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
     199                       const ClpPlusMinusOneMatrix & rhs,
     200                       int numberRows, const int * whichRow,
     201                       int numberColumns, const int * whichColumn)
     202: ClpMatrixBase(rhs)
     203
     204  elements_ = NULL;
     205  startPositive_ = NULL;
     206  startNegative_ = NULL;
     207  lengths_=NULL;
     208  indices_=NULL;
     209  numberRows_=0;
     210  numberColumns_=0;
     211  columnOrdered_=rhs.columnOrdered_;
     212  if (numberRows<=0||numberColumns<=0) {
     213    startPositive_ = new int[1];
     214    startPositive_[0] = 0;
     215  } else {
     216    numberColumns_ = numberColumns;
     217    numberRows_ = numberRows;
     218    const int * index1 = rhs.indices_;
     219    int * startPositive1 = rhs.startPositive_;
     220
     221    int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
     222    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
     223    int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
     224    int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
     225    // Throw exception if rhs empty
     226    if (numberMajor1 <= 0 || numberMinor1 <= 0)
     227      throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
     228    // Array to say if an old row is in new copy
     229    int * newRow = new int [numberMinor1];
     230    int iRow;
     231    for (iRow=0;iRow<numberMinor1;iRow++)
     232      newRow[iRow] = -1;
     233    // and array for duplicating rows
     234    int * duplicateRow = new int [numberMinor];
     235    int numberBad=0;
     236    for (iRow=0;iRow<numberMinor;iRow++) {
     237      duplicateRow[iRow] = -1;
     238      int kRow = whichRow[iRow];
     239      if (kRow>=0  && kRow < numberMinor1) {
     240        if (newRow[kRow]<0) {
     241          // first time
     242          newRow[kRow]=iRow;
     243        } else {
     244          // duplicate
     245          int lastRow = newRow[kRow];
     246          newRow[kRow]=iRow;
     247          duplicateRow[iRow] = lastRow;
     248        }
     249      } else {
     250        // bad row
     251        numberBad++;
     252      }
     253    }
     254
     255    if (numberBad)
     256      throw CoinError("bad minor entries",
     257                      "subset constructor", "ClpPlusMinusOneMatrix");
     258    // now get size and check columns
     259    int size = 0;
     260    int iColumn;
     261    numberBad=0;
     262    for (iColumn=0;iColumn<numberMajor;iColumn++) {
     263      int kColumn = whichColumn[iColumn];
     264      if (kColumn>=0  && kColumn <numberMajor1) {
     265        int i;
     266        for (i=startPositive1[kColumn];i<startPositive1[kColumn+1];i++) {
     267          int kRow = index1[i];
     268          kRow = newRow[kRow];
     269          while (kRow>=0) {
     270            size++;
     271            kRow = duplicateRow[kRow];
     272          }
     273        }
     274      } else {
     275        // bad column
     276        numberBad++;
     277      }
     278    }
     279    if (numberBad)
     280      throw CoinError("bad major entries",
     281                      "subset constructor", "ClpPlusMinusOneMatrix");
     282    // now create arrays
     283    startPositive_ = new int [numberMajor+1];
     284    startNegative_ = new int [numberMajor];
     285    indices_ = new int[size];
     286    // and fill them
     287    size = 0;
     288    startPositive_[0]=0;
     289    int * startNegative1 = rhs.startNegative_;
     290    for (iColumn=0;iColumn<numberMajor;iColumn++) {
     291      int kColumn = whichColumn[iColumn];
     292      int i;
     293      for (i=startPositive1[kColumn];i<startNegative1[kColumn];i++) {
     294        int kRow = index1[i];
     295        kRow = newRow[kRow];
     296        while (kRow>=0) {
     297          indices_[size++] = kRow;
     298          kRow = duplicateRow[kRow];
     299        }
     300      }
     301      startNegative_[iColumn] = size;
     302      for (;i<startPositive1[kColumn+1];i++) {
     303        int kRow = index1[i];
     304        kRow = newRow[kRow];
     305        while (kRow>=0) {
     306          indices_[size++] = kRow;
     307          kRow = duplicateRow[kRow];
     308        }
     309      }
     310      startPositive_[iColumn+1] = size;
     311    }
     312    delete [] newRow;
     313    delete [] duplicateRow;
     314  }
     315}
     316
    184317
    185318/* Returns a new matrix in reverse order without gaps */
     
    323456  double zeroTolerance = model->factorization()->zeroTolerance();
    324457  int numberRows = model->numberRows();
     458  bool packed = rowArray->packedMode();
    325459  ClpPlusMinusOneMatrix* rowCopy =
    326460    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    327461  if (numberInRowArray>0.3*numberRows||!rowCopy) {
     462    assert (!y->getNumElements());
    328463    // do by column
     464    // Need to expand if packed mode
    329465    int iColumn;
    330     double * markVector = y->denseVector(); // probably empty
    331466    CoinBigIndex j=0;
    332467    assert (columnOrdered_);
    333     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    334       double value2 = 0.0;
    335       for (;j<startNegative_[iColumn];j++) {
    336         int iRow = indices_[j];
    337         value2 += pi[iRow];
    338       }
    339       for (;j<startPositive_[iColumn+1];j++) {
    340         int iRow = indices_[j];
    341         value2 -= pi[iRow];
    342       }
    343       double value = markVector[iColumn];
    344       markVector[iColumn]=0.0;
    345       value += scalar*value2;
    346       if (fabs(value)>zeroTolerance) {
    347         index[numberNonZero++]=iColumn;
    348         array[iColumn]=value;
     468    if (packed) {
     469      // need to expand pi into y
     470      assert(y->capacity()>=numberRows);
     471      double * piOld = pi;
     472      pi = y->denseVector();
     473      const int * whichRow = rowArray->getIndices();
     474      int i;
     475      // modify pi so can collapse to one loop
     476      for (i=0;i<numberInRowArray;i++) {
     477        int iRow = whichRow[i];
     478        pi[iRow]=scalar*piOld[i];
     479      }
     480      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     481        double value = 0.0;
     482        for (;j<startNegative_[iColumn];j++) {
     483          int iRow = indices_[j];
     484          value += pi[iRow];
     485        }
     486        for (;j<startPositive_[iColumn+1];j++) {
     487          int iRow = indices_[j];
     488          value -= pi[iRow];
     489        }
     490        if (fabs(value)>zeroTolerance) {
     491          array[numberNonZero]=value;
     492          index[numberNonZero++]=iColumn;
     493        }
     494      }
     495      for (i=0;i<numberInRowArray;i++) {
     496        int iRow = whichRow[i];
     497        pi[iRow]=0.0;
     498      }
     499    } else {
     500      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     501        double value = 0.0;
     502        for (;j<startNegative_[iColumn];j++) {
     503          int iRow = indices_[j];
     504          value += pi[iRow];
     505        }
     506        for (;j<startPositive_[iColumn+1];j++) {
     507          int iRow = indices_[j];
     508          value -= pi[iRow];
     509        }
     510        value += scalar;
     511        if (fabs(value)>zeroTolerance) {
     512          index[numberNonZero++]=iColumn;
     513          array[iColumn]=value;
     514        }
    349515      }
    350516    }
    351517    columnArray->setNumElements(numberNonZero);
    352     y->setNumElements(0);
    353518  } else {
    354519    // do by row
     
    376541  const CoinBigIndex * startNegative = startNegative_;
    377542  const int * whichRow = rowArray->getIndices();
     543  bool packed = rowArray->packedMode();
    378544  if (numberInRowArray>2||y->getNumElements()) {
    379545    // do by rows
     
    383549    int numberOriginal=y->getNumElements();
    384550    int i;
    385     for (i=0;i<numberOriginal;i++) {
    386       int iColumn = mark[i];
    387       index[i]=iColumn;
    388       array[iColumn]=markVector[iColumn];
    389       markVector[iColumn]=0.0;
    390     }
    391     numberNonZero=numberOriginal;
    392     // and set up mark as char array
    393     char * marked = (char *) markVector;
    394     for (i=0;i<numberOriginal;i++) {
    395       int iColumn = index[i];
    396       marked[iColumn]=0;
    397     }
    398     for (i=0;i<numberInRowArray;i++) {
    399       iRow = whichRow[i];
    400       double value = pi[iRow]*scalar;
    401       CoinBigIndex j;
    402       for (j=startPositive[iRow];j<startNegative[iRow];j++) {
    403         int iColumn = column[j];
    404         if (!marked[iColumn]) {
    405           marked[iColumn]=1;
     551    if (packed) {
     552      assert(!numberOriginal);
     553      numberNonZero=0;
     554      // and set up mark as char array
     555      char * marked = (char *) (index+columnArray->capacity());
     556      double * array2 = y->denseVector();
     557#ifdef CLP_DEBUG
     558      int numberColumns = model->numberColumns();
     559      for (i=0;i<numberColumns;i++) {
     560        assert(!marked[i]);
     561        assert(!array2[i]);
     562      }
     563#endif
     564      for (i=0;i<numberInRowArray;i++) {
     565        iRow = whichRow[i];
     566        double value = pi[i]*scalar;
     567        CoinBigIndex j;
     568        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     569          int iColumn = column[j];
     570          if (!marked[iColumn]) {
     571            marked[iColumn]=1;
     572            index[numberNonZero++]=iColumn;
     573          }
     574          array2[iColumn] += value;
     575        }
     576        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     577          int iColumn = column[j];
     578          if (!marked[iColumn]) {
     579            marked[iColumn]=1;
     580            index[numberNonZero++]=iColumn;
     581          }
     582          array2[iColumn] -= value;
     583        }
     584      }
     585      // get rid of tiny values and zero out marked
     586      numberOriginal=numberNonZero;
     587      numberNonZero=0;
     588      for (i=0;i<numberOriginal;i++) {
     589        int iColumn = index[i];
     590        if (marked[iColumn]) {
     591          double value = array2[iColumn];
     592          array2[iColumn]=0.0;
     593          marked[iColumn]=0;
     594          if (fabs(value)>zeroTolerance) {
     595            array[numberNonZero]=value;
     596            index[numberNonZero++]=iColumn;
     597          }
     598        }
     599      }
     600    } else {     
     601      for (i=0;i<numberOriginal;i++) {
     602        int iColumn = mark[i];
     603        index[i]=iColumn;
     604        array[iColumn]=markVector[iColumn];
     605        markVector[iColumn]=0.0;
     606      }
     607      numberNonZero=numberOriginal;
     608      // and set up mark as char array
     609      char * marked = (char *) markVector;
     610      for (i=0;i<numberOriginal;i++) {
     611        int iColumn = index[i];
     612        marked[iColumn]=0;
     613      }
     614      for (i=0;i<numberInRowArray;i++) {
     615        iRow = whichRow[i];
     616        double value = pi[iRow]*scalar;
     617        CoinBigIndex j;
     618        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     619          int iColumn = column[j];
     620          if (!marked[iColumn]) {
     621            marked[iColumn]=1;
     622            index[numberNonZero++]=iColumn;
     623          }
     624          array[iColumn] += value;
     625        }
     626        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     627          int iColumn = column[j];
     628          if (!marked[iColumn]) {
     629            marked[iColumn]=1;
     630            index[numberNonZero++]=iColumn;
     631          }
     632          array[iColumn] -= value;
     633        }
     634      }
     635      // get rid of tiny values and zero out marked
     636      numberOriginal=numberNonZero;
     637      numberNonZero=0;
     638      for (i=0;i<numberOriginal;i++) {
     639        int iColumn = index[i];
     640        marked[iColumn]=0;
     641        if (fabs(array[iColumn])>zeroTolerance) {
    406642          index[numberNonZero++]=iColumn;
    407         }
    408         array[iColumn] += value;
    409       }
    410       for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
    411         int iColumn = column[j];
    412         if (!marked[iColumn]) {
    413           marked[iColumn]=1;
    414           index[numberNonZero++]=iColumn;
    415         }
    416         array[iColumn] -= value;
    417       }
    418     }
    419     // get rid of tiny values and zero out marked
    420     numberOriginal=numberNonZero;
    421     numberNonZero=0;
    422     for (i=0;i<numberOriginal;i++) {
    423       int iColumn = index[i];
    424       marked[iColumn]=0;
    425       if (fabs(array[iColumn])>zeroTolerance) {
    426         index[numberNonZero++]=iColumn;
    427       } else {
    428         array[iColumn]=0.0;
     643        } else {
     644          array[iColumn]=0.0;
     645        }
    429646      }
    430647    }
    431648  } else if (numberInRowArray==2) {
    432     // do by rows when two rows (do longer first)
     649    /* do by rows when two rows (do longer first when not packed
     650       and shorter first if packed */
    433651    int iRow0 = whichRow[0];
    434652    int iRow1 = whichRow[1];
    435     if (startPositive[iRow0+1]-startPositive[iRow0]<
    436         startPositive[iRow1+1]-startPositive[iRow1]) {
    437       int temp = iRow0;
    438       iRow0 = iRow1;
    439       iRow1 = temp;
    440     }
    441     int numberOriginal;
    442     int i;
    443     numberNonZero=0;
    444     double value;
    445     value = pi[iRow0]*scalar;
    446     CoinBigIndex j;
    447     for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
    448       int iColumn = column[j];
    449       index[numberNonZero++]=iColumn;
    450       array[iColumn] = value;
    451     }
    452     for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
    453       int iColumn = column[j];
    454       index[numberNonZero++]=iColumn;
    455       array[iColumn] = -value;
    456     }
    457     value = pi[iRow1]*scalar;
    458     for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
    459       int iColumn = column[j];
    460       double value2= array[iColumn];
    461       if (value2) {
    462         value2 += value;
    463       } else {
    464         value2 = value;
     653    int j;
     654    if (packed) {
     655      double pi0 = pi[0];
     656      double pi1 = pi[1];
     657      if (startPositive[iRow0+1]-startPositive[iRow0]>
     658          startPositive[iRow1+1]-startPositive[iRow1]) {
     659        int temp = iRow0;
     660        iRow0 = iRow1;
     661        iRow1 = temp;
     662        pi0=pi1;
     663        pi1=pi[0];
     664      }
     665      // and set up mark as char array
     666      char * marked = (char *) (index+columnArray->capacity());
     667      int * lookup = y->getIndices();
     668      double value = pi0*scalar;
     669      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
     670        int iColumn = column[j];
     671        array[numberNonZero] = value;
     672        marked[iColumn]=1;
     673        lookup[iColumn]=numberNonZero;
    465674        index[numberNonZero++]=iColumn;
    466675      }
    467       array[iColumn] = value2;
    468     }
    469     for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
    470       int iColumn = column[j];
    471       double value2= array[iColumn];
    472       if (value2) {
    473         value2 -= value;
    474       } else {
    475         value2 = -value;
     676      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
     677        int iColumn = column[j];
     678        array[numberNonZero] = -value;
     679        marked[iColumn]=1;
     680        lookup[iColumn]=numberNonZero;
    476681        index[numberNonZero++]=iColumn;
    477682      }
    478       array[iColumn] = value2;
    479     }
    480     // get rid of tiny values and zero out marked
    481     numberOriginal=numberNonZero;
    482     numberNonZero=0;
    483     for (i=0;i<numberOriginal;i++) {
    484       int iColumn = index[i];
    485       if (fabs(array[iColumn])>zeroTolerance) {
     683      int numberOriginal = numberNonZero;
     684      value = pi1*scalar;
     685      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
     686        int iColumn = column[j];
     687        if (marked[iColumn]) {
     688          int iLookup = lookup[iColumn];
     689          array[iLookup] += value;
     690        } else {
     691          if (fabs(value)>zeroTolerance) {
     692            array[numberNonZero] = value;
     693            index[numberNonZero++]=iColumn;
     694          }
     695        }
     696      }
     697      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
     698        int iColumn = column[j];
     699        if (marked[iColumn]) {
     700          int iLookup = lookup[iColumn];
     701          array[iLookup] -= value;
     702        } else {
     703          if (fabs(value)>zeroTolerance) {
     704            array[numberNonZero] = -value;
     705            index[numberNonZero++]=iColumn;
     706          }
     707        }
     708      }
     709      // get rid of tiny values and zero out marked
     710      int nDelete=0;
     711      for (j=0;j<numberOriginal;j++) {
     712        int iColumn = index[j];
     713        marked[iColumn]=0;
     714        if (fabs(array[j])<=zeroTolerance)
     715          nDelete++;
     716      }
     717      if (nDelete) {
     718        numberOriginal=numberNonZero;
     719        numberNonZero=0;
     720        for (j=0;j<numberOriginal;j++) {
     721          int iColumn = index[j];
     722          double value = array[j];
     723          array[j]=0.0;
     724          if (fabs(value)>zeroTolerance) {
     725            array[numberNonZero]=value;
     726            index[numberNonZero++]=iColumn;
     727          }
     728        }
     729      }
     730    } else {
     731      if (startPositive[iRow0+1]-startPositive[iRow0]<
     732          startPositive[iRow1+1]-startPositive[iRow1]) {
     733        int temp = iRow0;
     734        iRow0 = iRow1;
     735        iRow1 = temp;
     736      }
     737      int numberOriginal;
     738      int i;
     739      numberNonZero=0;
     740      double value;
     741      value = pi[iRow0]*scalar;
     742      CoinBigIndex j;
     743      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
     744        int iColumn = column[j];
    486745        index[numberNonZero++]=iColumn;
    487       } else {
    488         array[iColumn]=0.0;
     746        array[iColumn] = value;
     747      }
     748      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
     749        int iColumn = column[j];
     750        index[numberNonZero++]=iColumn;
     751        array[iColumn] = -value;
     752      }
     753      value = pi[iRow1]*scalar;
     754      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
     755        int iColumn = column[j];
     756        double value2= array[iColumn];
     757        if (value2) {
     758          value2 += value;
     759        } else {
     760          value2 = value;
     761          index[numberNonZero++]=iColumn;
     762        }
     763        array[iColumn] = value2;
     764      }
     765      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
     766        int iColumn = column[j];
     767        double value2= array[iColumn];
     768        if (value2) {
     769          value2 -= value;
     770        } else {
     771          value2 = -value;
     772          index[numberNonZero++]=iColumn;
     773        }
     774        array[iColumn] = value2;
     775      }
     776      // get rid of tiny values and zero out marked
     777      numberOriginal=numberNonZero;
     778      numberNonZero=0;
     779      for (i=0;i<numberOriginal;i++) {
     780        int iColumn = index[i];
     781        if (fabs(array[iColumn])>zeroTolerance) {
     782          index[numberNonZero++]=iColumn;
     783        } else {
     784          array[iColumn]=0.0;
     785        }
    489786      }
    490787    }
     
    495792    double value;
    496793    iRow = whichRow[0];
    497     value = pi[iRow]*scalar;
    498     if (fabs(value)>zeroTolerance) {
    499       CoinBigIndex j;
    500       for (j=startPositive[iRow];j<startNegative[iRow];j++) {
    501         int iColumn = column[j];
    502         index[numberNonZero++]=iColumn;
    503         array[iColumn] = value;
    504       }
    505       for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
    506         int iColumn = column[j];
    507         index[numberNonZero++]=iColumn;
    508         array[iColumn] = -value;
     794    CoinBigIndex j;
     795    if (packed) {
     796      value = pi[0]*scalar;
     797      if (fabs(value)>zeroTolerance) {
     798        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     799          int iColumn = column[j];
     800          array[numberNonZero] = value;
     801          index[numberNonZero++]=iColumn;
     802        }
     803        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     804          int iColumn = column[j];
     805          array[numberNonZero] = -value;
     806          index[numberNonZero++]=iColumn;
     807        }
     808      }
     809    } else {
     810      value = pi[iRow]*scalar;
     811      if (fabs(value)>zeroTolerance) {
     812        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     813          int iColumn = column[j];
     814          array[iColumn] = value;
     815          index[numberNonZero++]=iColumn;
     816        }
     817        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     818          int iColumn = column[j];
     819          array[iColumn] = -value;
     820          index[numberNonZero++]=iColumn;
     821        }
    509822      }
    510823    }
     
    532845  int numberToDo = y->getNumElements();
    533846  const int * which = y->getIndices();
    534   for (jColumn=0;jColumn<numberToDo;jColumn++) {
    535     int iColumn = which[jColumn];
    536     double value = 0.0;
    537     CoinBigIndex j=startPositive_[iColumn];
    538     for (;j<startNegative_[iColumn];j++) {
    539       int iRow = indices_[j];
    540       value += pi[iRow];
    541     }
    542     for (;j<startPositive_[iColumn+1];j++) {
    543       int iRow = indices_[j];
    544       value -= pi[iRow];
    545     }
    546     if (fabs(value)>zeroTolerance) {
    547       index[numberNonZero++]=iColumn;
    548       array[iColumn]=value;
     847  bool packed = rowArray->packedMode();
     848  if (packed) {
     849    // need to expand pi into y
     850    int numberInRowArray = rowArray->getNumElements();
     851    assert(y->capacity()>=model->numberRows());
     852    double * piOld = pi;
     853    pi = y->denseVector();
     854    const int * whichRow = rowArray->getIndices();
     855    int i;
     856    for (i=0;i<numberInRowArray;i++) {
     857      int iRow = whichRow[i];
     858      pi[iRow]=piOld[i];
     859    }
     860    // Must line up with y
     861    for (jColumn=0;jColumn<numberToDo;jColumn++) {
     862      int iColumn = which[jColumn];
     863      double value = 0.0;
     864      CoinBigIndex j=startPositive_[iColumn];
     865      for (;j<startNegative_[iColumn];j++) {
     866        int iRow = indices_[j];
     867        value += pi[iRow];
     868      }
     869      for (;j<startPositive_[iColumn+1];j++) {
     870        int iRow = indices_[j];
     871        value -= pi[iRow];
     872      }
     873      array[jColumn]=value;
     874    }
     875    for (i=0;i<numberInRowArray;i++) {
     876      int iRow = whichRow[i];
     877      pi[iRow]=0.0;
     878    }
     879  } else {
     880    for (jColumn=0;jColumn<numberToDo;jColumn++) {
     881      int iColumn = which[jColumn];
     882      double value = 0.0;
     883      CoinBigIndex j=startPositive_[iColumn];
     884      for (;j<startNegative_[iColumn];j++) {
     885        int iRow = indices_[j];
     886        value += pi[iRow];
     887      }
     888      for (;j<startPositive_[iColumn+1];j++) {
     889        int iRow = indices_[j];
     890        value -= pi[iRow];
     891      }
     892      if (fabs(value)>zeroTolerance) {
     893        index[numberNonZero++]=iColumn;
     894        array[iColumn]=value;
     895      }
    549896    }
    550897  }
     
    598945  return numberElements;
    599946}
     947/* If element NULL returns number of elements in column part of basis,
     948   If not NULL fills in as well */
     949CoinBigIndex
     950ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
     951                                 const int * whichColumn,
     952                                 int numberBasic,
     953                                 int numberColumnBasic,
     954                                 int * indexRowU, int * indexColumnU,
     955                                 double * elementU) const
     956{
     957  int i;
     958  CoinBigIndex numberElements=0;
     959  if (elementU!=NULL) {
     960    assert (columnOrdered_);
     961    for (i=0;i<numberColumnBasic;i++) {
     962      int iColumn = whichColumn[i];
     963      CoinBigIndex j=startPositive_[iColumn];
     964      for (;j<startNegative_[iColumn];j++) {
     965        int iRow = indices_[j];
     966        indexRowU[numberElements]=iRow;
     967        indexColumnU[numberElements]=numberBasic;
     968        elementU[numberElements++]=1.0;
     969      }
     970      for (;j<startPositive_[iColumn+1];j++) {
     971        int iRow = indices_[j];
     972        indexRowU[numberElements]=iRow;
     973        indexColumnU[numberElements]=numberBasic;
     974        elementU[numberElements++]=-1.0;
     975      }
     976      numberBasic++;
     977    }
     978  } else {
     979    for (i=0;i<numberColumnBasic;i++) {
     980      int iColumn = whichColumn[i];
     981      numberElements += startPositive_[iColumn+1]-startPositive_[iColumn];
     982    }
     983  }
     984  return numberElements;
     985}
    600986/* Unpacks a column into an CoinIndexedvector
    601       Note that model is NOT const.  Bounds and objective could
    602       be modified if doing column generation */
     987 */
    603988void
    604989ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
     
    6151000    rowArray->add(iRow,-1.0);
    6161001  }
     1002}
     1003/* Unpacks a column into an CoinIndexedvector
     1004** in packed foramt
     1005Note that model is NOT const.  Bounds and objective could
     1006be modified if doing column generation (just for this variable) */
     1007void
     1008ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * model,
     1009                            CoinIndexedVector * rowArray,
     1010                            int iColumn) const
     1011{
     1012  int * index = rowArray->getIndices();
     1013  double * array = rowArray->denseVector();
     1014  int number = 0;
     1015  CoinBigIndex j=startPositive_[iColumn];
     1016  for (;j<startNegative_[iColumn];j++) {
     1017    int iRow = indices_[j];
     1018    array[number]=1.0;
     1019    index[number++]=iRow;
     1020  }
     1021  for (;j<startPositive_[iColumn+1];j++) {
     1022    int iRow = indices_[j];
     1023    array[number]=-1.0;
     1024    index[number++]=iRow;
     1025  }
     1026  rowArray->setNumElements(number);
     1027  rowArray->setPackedMode(true);
    6171028}
    6181029/* Adds multiple of a column into an CoinIndexedvector
  • branches/pre/ClpPrimalColumnDantzig.cpp

    r137 r180  
    139139
    140140  for (iSequence=0;iSequence<number;iSequence++) {
    141     double value = reducedCost[iSequence];
    142     ClpSimplex::Status status = model_->getStatus(iSequence);
    143    
    144     switch(status) {
     141    // check flagged variable
     142    if (!model_->flagged(iSequence)) {
     143      double value = reducedCost[iSequence];
     144      ClpSimplex::Status status = model_->getStatus(iSequence);
    145145     
    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;
     146      switch(status) {
     147
     148      case ClpSimplex::basic:
     149      case ClpSimplex::isFixed:
     150        break;
     151      case ClpSimplex::isFree:
     152      case ClpSimplex::superBasic:
     153        if (fabs(value)>bestFreeDj) {
     154          bestFreeDj = fabs(value);
     155          bestFreeSequence = iSequence;
     156        }
     157        break;
     158      case ClpSimplex::atUpperBound:
     159        if (value>bestDj) {
     160          bestDj = value;
     161          bestSequence = iSequence;
     162        }
     163        break;
     164      case ClpSimplex::atLowerBound:
     165        if (value<-bestDj) {
     166          bestDj = -value;
     167          bestSequence = iSequence;
     168        }
    166169      }
    167170    }
  • branches/pre/ClpPrimalColumnSteepest.cpp

    r162 r180  
    1515
    1616// bias for free variables
    17 
    1817#define FREE_BIAS 1.0e1
     18// Acceptance criteria for free variables
     19#define FREE_ACCEPT 1.0e2
    1920//#############################################################################
    2021// Constructors / Destructor / Assignment
     
    227228      }
    228229      if (!model_->nonLinearCost()->lookBothWays()) {
    229        
     230
    230231        for (j=0;j<number;j++) {
    231232          int iSequence = index[j];
     
    243244          case ClpSimplex::isFree:
    244245          case ClpSimplex::superBasic:
    245             if (fabs(value)>1.0e2*tolerance) {
     246            if (fabs(value)>FREE_ACCEPT*tolerance) {
    246247              // we are going to bias towards free (but only if reasonable)
    247248              value *= FREE_BIAS;
     
    296297          case ClpSimplex::isFree:
    297298          case ClpSimplex::superBasic:
    298             if (fabs(value)>1.0e2*tolerance) {
     299            if (fabs(value)>FREE_ACCEPT*tolerance) {
    299300              // we are going to bias towards free (but only if reasonable)
    300301              value *= FREE_BIAS;
     
    377378    case ClpSimplex::isFree:
    378379    case ClpSimplex::superBasic:
    379       if (fabs(value)>1.0e2*tolerance) {
     380      if (fabs(value)>FREE_ACCEPT*tolerance) {
    380381        // we are going to bias towards free (but only if reasonable)
    381382        value *= FREE_BIAS;
     
    432433        double value = reducedCost[iSequence];
    433434        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    434        
     435
    435436        switch(status) {
    436437         
     
    481482        double value = reducedCost[iSequence];
    482483        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    483        
     484
    484485        switch(status) {
    485486         
     
    681682    if (model_->largestDualError()>checkTolerance)
    682683      tolerance *= model_->largestDualError()/checkTolerance;
    683   }
     684    // But cap
     685    tolerance = min(1000.0,tolerance);
     686  }
     687#ifdef CLP_DEBUG
     688  if (model_->numberDualInfeasibilities()==1)
     689    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
     690           model_->largestDualError(),model_,model_->messageHandler(),
     691           number);
     692#endif
    684693  // stop last one coming immediately
    685694  double saveOutInfeasibility=0.0;
     
    693702    double value = infeas[iSequence];
    694703    double weight = weights_[iSequence];
    695     /*if (model_->numberIterations()%100==0)
    696       printf("%d inf %g wt %g\n",iSequence,value,weight);*/
    697704    //weight=1.0;
    698705    if (value>bestDj*weight&&value>tolerance) {
     
    761768      initializeWeights();
    762769      // create saved weights
     770      delete [] savedWeights_;
    763771      savedWeights_ = new double[numberRows+numberColumns];
    764772      memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
     
    828836    int number = model_->numberRows() + model_->numberColumns();
    829837    int iSequence;
    830 
     838   
    831839    double * reducedCost = model_->djRegion();
    832840     
     
    835843        double value = reducedCost[iSequence];
    836844        ClpSimplex::Status status = model_->getStatus(iSequence);
    837        
     845
    838846        switch(status) {
    839847         
     
    843851        case ClpSimplex::isFree:
    844852        case ClpSimplex::superBasic:
    845           if (fabs(value)>1.0e2*tolerance) {
     853          if (fabs(value)>FREE_ACCEPT*tolerance) {
    846854            // we are going to bias towards free (but only if reasonable)
    847855            value *= FREE_BIAS;
     
    867875        double value = reducedCost[iSequence];
    868876        ClpSimplex::Status status = model_->getStatus(iSequence);
    869        
     877
    870878        switch(status) {
    871879         
     
    875883        case ClpSimplex::isFree:
    876884        case ClpSimplex::superBasic:
    877           if (fabs(value)>1.0e2*tolerance) {
     885          if (fabs(value)>FREE_ACCEPT*tolerance) {
    878886            // we are going to bias towards free (but only if reasonable)
    879887            value *= FREE_BIAS;
     
    11581166  devex_ = 0.0;
    11591167}
     1168// Returns true if would not find any column
     1169bool
     1170ClpPrimalColumnSteepest::looksOptimal() const
     1171{
     1172  //**** THIS MUST MATCH the action coding above
     1173  double tolerance=model_->currentDualTolerance();
     1174  // we can't really trust infeasibilities if there is dual error
     1175  // this coding has to mimic coding in checkDualSolution
     1176  double error = min(1.0e-3,model_->largestDualError());
     1177  // allow tolerance at least slightly bigger than standard
     1178  tolerance = tolerance +  error;
     1179  if(model_->numberIterations()<model_->lastBadIteration()+200) {
     1180    // we can't really trust infeasibilities if there is dual error
     1181    double checkTolerance = 1.0e-8;
     1182    if (!model_->factorization()->pivots())
     1183      checkTolerance = 1.0e-6;
     1184    if (model_->largestDualError()>checkTolerance)
     1185      tolerance *= model_->largestDualError()/checkTolerance;
     1186    // But cap
     1187    tolerance = min(1000.0,tolerance);
     1188  }
     1189  int number = model_->numberRows() + model_->numberColumns();
     1190  int iSequence;
     1191 
     1192  double * reducedCost = model_->djRegion();
     1193  int numberInfeasible=0;
     1194  if (!model_->nonLinearCost()->lookBothWays()) {
     1195    for (iSequence=0;iSequence<number;iSequence++) {
     1196      double value = reducedCost[iSequence];
     1197      ClpSimplex::Status status = model_->getStatus(iSequence);
     1198     
     1199      switch(status) {
     1200
     1201      case ClpSimplex::basic:
     1202      case ClpSimplex::isFixed:
     1203        break;
     1204      case ClpSimplex::isFree:
     1205      case ClpSimplex::superBasic:
     1206        if (fabs(value)>FREE_ACCEPT*tolerance)
     1207          numberInfeasible++;
     1208        break;
     1209      case ClpSimplex::atUpperBound:
     1210        if (value>tolerance)
     1211          numberInfeasible++;
     1212        break;
     1213      case ClpSimplex::atLowerBound:
     1214        if (value<-tolerance)
     1215          numberInfeasible++;
     1216      }
     1217    }
     1218  } else {
     1219    ClpNonLinearCost * nonLinear = model_->nonLinearCost();
     1220    // can go both ways
     1221    for (iSequence=0;iSequence<number;iSequence++) {
     1222      double value = reducedCost[iSequence];
     1223      ClpSimplex::Status status = model_->getStatus(iSequence);
     1224     
     1225      switch(status) {
     1226
     1227      case ClpSimplex::basic:
     1228      case ClpSimplex::isFixed:
     1229        break;
     1230      case ClpSimplex::isFree:
     1231      case ClpSimplex::superBasic:
     1232        if (fabs(value)>FREE_ACCEPT*tolerance)
     1233          numberInfeasible++;
     1234        break;
     1235      case ClpSimplex::atUpperBound:
     1236        if (value>tolerance) {
     1237          numberInfeasible++;
     1238        } else {
     1239          // look other way - change up should be negative
     1240          value -= nonLinear->changeUpInCost(iSequence);
     1241          if (value<-tolerance)
     1242            numberInfeasible++;
     1243        }
     1244        break;
     1245      case ClpSimplex::atLowerBound:
     1246        if (value<-tolerance) {
     1247          numberInfeasible++;
     1248        } else {
     1249          // look other way - change down should be positive
     1250          value -= nonLinear->changeDownInCost(iSequence);
     1251          if (value>tolerance)
     1252            numberInfeasible++;
     1253        }
     1254      }
     1255    }
     1256  }
     1257  return numberInfeasible==0;
     1258}
  • branches/pre/ClpSimplex.cpp

    r171 r180  
    104104  numberFake_(0),
    105105  progressFlag_(0),
    106   firstFree_(-1)
    107 
     106  firstFree_(-1),
     107  incomingInfeasibility_(1.0),
     108  allowedInfeasibility_(10.0),
     109  progress_(NULL)
    108110{
    109111  int i;
     
    125127}
    126128
     129// Subproblem constructor
     130ClpSimplex::ClpSimplex ( const ClpModel * rhs,
     131                     int numberRows, const int * whichRow,
     132                     int numberColumns, const int * whichColumn,
     133                     bool dropNames, bool dropIntegers)
     134  : ClpModel(rhs, numberRows, whichRow,
     135             numberColumns,whichColumn,dropNames,dropIntegers),
     136  columnPrimalInfeasibility_(0.0),
     137  rowPrimalInfeasibility_(0.0),
     138  columnPrimalSequence_(-2),
     139  rowPrimalSequence_(-2),
     140  columnDualInfeasibility_(0.0),
     141  rowDualInfeasibility_(0.0),
     142  columnDualSequence_(-2),
     143  rowDualSequence_(-2),
     144  primalToleranceToGetOptimal_(-1.0),
     145  remainingDualInfeasibility_(0.0),
     146  largeValue_(1.0e15),
     147  largestPrimalError_(0.0),
     148  largestDualError_(0.0),
     149  largestSolutionError_(0.0),
     150  dualBound_(1.0e10),
     151  alpha_(0.0),
     152  theta_(0.0),
     153  lowerIn_(0.0),
     154  valueIn_(0.0),
     155  upperIn_(0.0),
     156  dualIn_(0.0),
     157  lowerOut_(-1),
     158  valueOut_(-1),
     159  upperOut_(-1),
     160  dualOut_(-1),
     161  dualTolerance_(0.0),
     162  primalTolerance_(0.0),
     163  sumDualInfeasibilities_(0.0),
     164  sumPrimalInfeasibilities_(0.0),
     165  infeasibilityCost_(1.0e10),
     166  sumOfRelaxedDualInfeasibilities_(0.0),
     167  sumOfRelaxedPrimalInfeasibilities_(0.0),
     168  lower_(NULL),
     169  rowLowerWork_(NULL),
     170  columnLowerWork_(NULL),
     171  upper_(NULL),
     172  rowUpperWork_(NULL),
     173  columnUpperWork_(NULL),
     174  cost_(NULL),
     175  rowObjectiveWork_(NULL),
     176  objectiveWork_(NULL),
     177  sequenceIn_(-1),
     178  directionIn_(-1),
     179  sequenceOut_(-1),
     180  directionOut_(-1),
     181  pivotRow_(-1),
     182  lastGoodIteration_(-100),
     183  dj_(NULL),
     184  rowReducedCost_(NULL),
     185  reducedCostWork_(NULL),
     186  solution_(NULL),
     187  rowActivityWork_(NULL),
     188  columnActivityWork_(NULL),
     189  numberDualInfeasibilities_(0),
     190  numberDualInfeasibilitiesWithoutFree_(0),
     191  numberPrimalInfeasibilities_(100),
     192  numberRefinements_(0),
     193  pivotVariable_(NULL),
     194  factorization_(NULL),
     195  rowScale_(NULL),
     196  savedSolution_(NULL),
     197  columnScale_(NULL),
     198  scalingFlag_(1),
     199  numberTimesOptimal_(0),
     200  changeMade_(1),
     201  algorithm_(0),
     202  forceFactorization_(-1),
     203  perturbation_(100),
     204  nonLinearCost_(NULL),
     205  specialOptions_(0),
     206  lastBadIteration_(-999999),
     207  numberFake_(0),
     208  progressFlag_(0),
     209  firstFree_(-1),
     210  incomingInfeasibility_(1.0),
     211  allowedInfeasibility_(10.0),
     212  progress_(NULL)
     213{
     214  int i;
     215  for (i=0;i<6;i++) {
     216    rowArray_[i]=NULL;
     217    columnArray_[i]=NULL;
     218  }
     219  saveStatus_=NULL;
     220  // get an empty factorization so we can set tolerances etc
     221  factorization_ = new ClpFactorization();
     222  // say Steepest pricing
     223  dualRowPivot_ = new ClpDualRowSteepest();
     224  // say Steepest pricing
     225  primalColumnPivot_ = new ClpPrimalColumnSteepest();
     226  solveType_=1; // say simplex based life form
     227 
     228}
    127229
    128230//-----------------------------------------------------------------------------
     
    140242}
    141243int
    142 ClpSimplex::gutsOfSolution ( const double * rowActivities,
    143                              const double * columnActivities,
    144                              double * givenDuals,
     244ClpSimplex::gutsOfSolution ( double * givenDuals,
    145245                             const double * givenPrimals,
    146246                             bool valuesPass)
     
    155255    assert(nonLinearCost_);
    156256    int iRow;
    157     checkPrimalSolution( rowActivities, columnActivities);
     257    checkPrimalSolution( rowActivityWork_, columnActivityWork_);
    158258    // get correct bounds on all variables
    159     nonLinearCost_->checkInfeasibilities();
     259    nonLinearCost_->checkInfeasibilities(primalTolerance_);
    160260    oldValue = nonLinearCost_->largestInfeasibility();
    161261    save = new double[numberRows_];
     
    166266  }
    167267  // do work
    168   computePrimals(rowActivities, columnActivities);
     268  computePrimals(rowActivityWork_, columnActivityWork_);
    169269  // If necessary - override results
    170270  if (givenPrimals) {
     
    177277    // primal algorithm
    178278    // get correct bounds on all variables
    179     nonLinearCost_->checkInfeasibilities();
     279    // If  4 bit set - Force outgoing variables to exact bound (primal)
     280    if ((specialOptions_&4)==0)
     281      nonLinearCost_->checkInfeasibilities(primalTolerance_);
     282    else
     283      nonLinearCost_->checkInfeasibilities(0.0);
    180284    objectiveModification += nonLinearCost_->changeInCost();
    181285    if (nonLinearCost_->numberInfeasibilities())
     
    191295#endif
    192296    int numberOut=0;
    193     if (oldValue<1.0
    194         &&nonLinearCost_->largestInfeasibility()>1.1*oldValue+1.0e-4||
     297    if (oldValue<incomingInfeasibility_
     298        &&nonLinearCost_->largestInfeasibility()>
     299        max(incomingInfeasibility_,allowedInfeasibility_)||
    195300        largestPrimalError_>1.0e-3) {
     301      printf("Original largest infeas %g, now %g, primalError %g\n",
     302             oldValue,nonLinearCost_->largestInfeasibility(),
     303             largestPrimalError_);
    196304      // throw out up to 1000 structurals
    197305      int iRow;
     
    200308      for (iRow=0;iRow<numberRows_;iRow++) {
    201309        int iPivot=pivotVariable_[iRow];
    202         double difference = fabs(solution_[iPivot]=save[iRow]);
     310        double difference = fabs(solution_[iPivot]-save[iRow]);
    203311        solution_[iPivot]=save[iRow];
    204312        save[iRow]=difference;
     
    234342
    235343  // now check solutions
    236   checkPrimalSolution( rowActivities, columnActivities);
     344  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
    237345  objectiveValue_ += objectiveModification;
    238346  checkDualSolution();
     
    243351      <<largestDualError_
    244352      <<CoinMessageEol;
     353  // Switch off false values pass indicator
     354  if (!valuesPass&&algorithm_>0)
     355    firstFree_ = -1;
    245356  return 0;
    246357}
     
    253364  CoinIndexedVector  * workSpace = rowArray_[0];
    254365
    255   double * array = new double [numberRows_+1]; // +1 for network
    256   double * save = new double [numberRows_+1];
    257   double * previous = new double [numberRows_+1];
     366  CoinIndexedVector arrayVector;
     367  arrayVector.reserve(numberRows_+1);
     368  CoinIndexedVector previousVector;
     369  previousVector.reserve(numberRows_+1);
    258370
    259371  // accumulate non basic stuff
    260   ClpFillN(array,numberRows_,0.0);
    261372
    262373  int iRow;
     
    272383    solution_[iPivot] = 0.0;
    273384  }
     385  double * array = arrayVector.denseVector();
    274386  times(-1.0,columnActivityWork_,array);
    275387
     388  int * index = arrayVector.getIndices();
     389  int number=0;
    276390  for (iRow=0;iRow<numberRows_;iRow++) {
    277     array[iRow] += rowActivityWork_[iRow];
    278   }
     391    double value = array[iRow] + rowActivityWork_[iRow];
     392    if (value) {
     393      array[iRow]=value;
     394      index[number++]=iRow;
     395    } else {
     396      array[iRow]=0.0;
     397    }
     398  }
     399
     400  arrayVector.setNumElements(number);
    279401
    280402  // Ftran adjusted RHS and iterate to improve accuracy
     
    282404  int iRefine;
    283405  double * work = workSpace->denseVector();
    284   factorization_->updateColumn(workSpace,array);
    285 
     406  CoinIndexedVector * thisVector = &arrayVector;
     407  CoinIndexedVector * lastVector = &previousVector;
     408  factorization_->updateColumn(workSpace,thisVector);
     409  bool goodSolution=true;
    286410  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
    287     // save in case we want to restore
    288     ClpDisjointCopyN ( array, numberRows_ , save);
    289 
     411
     412    int numberIn = thisVector->getNumElements();
     413    int * indexIn = thisVector->getIndices();
     414    double * arrayIn = thisVector->denseVector();
     415    // put solution in correct place
     416    int j;
     417    for (j=0;j<numberIn;j++) {
     418      iRow = indexIn[j];
     419      int iPivot=pivotVariable_[iRow];
     420      solution_[iPivot] = arrayIn[iRow];
     421    }
     422    // check Ax == b  (for all)
     423    times(-1.0,columnActivityWork_,work);
     424    largestPrimalError_=0.0;
     425    double multiplier = 131072.0;
     426    for (iRow=0;iRow<numberRows_;iRow++) {
     427      double value = work[iRow] + rowActivityWork_[iRow];
     428      work[iRow] = value*multiplier;
     429      if (fabs(value)>largestPrimalError_) {
     430        largestPrimalError_=fabs(value);
     431      }
     432    }
     433    if (largestPrimalError_>=lastError) {
     434      // restore
     435      CoinIndexedVector * temp = thisVector;
     436      thisVector = lastVector;
     437      lastVector=temp;
     438      goodSolution=false;
     439      break;
     440    }
     441    if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-10) {
     442      // try and make better
     443      // save this
     444      CoinIndexedVector * temp = thisVector;
     445      thisVector = lastVector;
     446      lastVector=temp;
     447      int * indexOut = thisVector->getIndices();
     448      int number=0;
     449      array = thisVector->denseVector();
     450      thisVector->clear();
     451      for (iRow=0;iRow<numberRows_;iRow++) {
     452        double value = work[iRow];
     453        if (value) {
     454          array[iRow]=value;
     455          indexOut[number++]=iRow;
     456          work[iRow]=0.0;
     457        }
     458      }
     459      thisVector->setNumElements(number);
     460      lastError=largestPrimalError_;
     461      factorization_->updateColumn(workSpace,thisVector);
     462      multiplier = 1.0/multiplier;
     463      double * previous = lastVector->denseVector();
     464      number=0;
     465      for (iRow=0;iRow<numberRows_;iRow++) {
     466        double value = previous[iRow] + multiplier*array[iRow];
     467        if (value) {
     468          array[iRow]=value;
     469          indexOut[number++]=iRow;
     470        } else {
     471          array[iRow]=0.0;
     472        }
     473      }
     474      thisVector->setNumElements(number);
     475    } else {
     476      break;
     477    }
     478  }
     479
     480  // solution as accurate as we are going to get
     481  ClpFillN(work,numberRows_,0.0);
     482  if (!goodSolution) {
     483    array = thisVector->denseVector();
    290484    // put solution in correct place
    291485    for (iRow=0;iRow<numberRows_;iRow++) {
     
    293487      solution_[iPivot] = array[iRow];
    294488    }
    295     // check Ax == b  (for all)
    296     times(-1.0,columnActivityWork_,work);
    297     for (iRow=0;iRow<numberRows_;iRow++) {
    298       work[iRow] += rowActivityWork_[iRow];
    299     }
    300 
    301     largestPrimalError_=0.0;
    302     for (iRow=0;iRow<numberRows_;iRow++) {
    303       if (fabs(work[iRow])>largestPrimalError_) {
    304         largestPrimalError_=fabs(work[iRow]);
    305       }
    306       //work[iRow] -= save[iRow];
    307     }
    308     if (largestPrimalError_>=lastError) {
    309       // restore
    310       double * temp = array;
    311       array = previous;
    312       previous=temp;
    313       break;
    314     }
    315     if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-12) {
    316       // try and make better
    317       // save this
    318       double * temp = array;
    319       array = previous;
    320       previous=temp;
    321       double multiplier = 131072.0;
    322       for (iRow=0;iRow<numberRows_;iRow++) {
    323         array[iRow] = multiplier*work[iRow];
    324         work[iRow]=0.0;
    325       }
    326       lastError=largestPrimalError_;
    327       factorization_->updateColumn(workSpace,array);
    328       multiplier = 1.0/multiplier;
    329       for (iRow=0;iRow<numberRows_;iRow++) {
    330         array[iRow] = previous[iRow] + multiplier*array[iRow];
    331         work[iRow]=0.0;
    332       }
    333     }
    334   }
    335 
    336   // solution as accurate as we are going to get
    337   ClpFillN(work,numberRows_,0.0);
    338   // put solution in correct place
    339   for (iRow=0;iRow<numberRows_;iRow++) {
    340     int iPivot=pivotVariable_[iRow];
    341     solution_[iPivot] = array[iRow];
    342   }
    343   delete [] previous;
    344   delete [] array;
    345   delete [] save;
     489  }
    346490}
    347491// now dual side
     
    352496  CoinIndexedVector  * workSpace = rowArray_[0];
    353497
    354   double * array = new double [numberRows_+1]; // +1 for network
    355   double * save = new double [numberRows_+1];
    356   double * previous = new double [numberRows_+1];
     498  CoinIndexedVector arrayVector;
     499  arrayVector.reserve(numberRows_+1);
     500  CoinIndexedVector previousVector;
     501  previousVector.reserve(numberRows_+1);
     502
    357503
    358504  int iRow;
     
    360506  workSpace->checkClear();
    361507#endif
     508  double * array = arrayVector.denseVector();
     509  int * index = arrayVector.getIndices();
     510  int number=0;
    362511  if (!givenDjs) {
    363512    for (iRow=0;iRow<numberRows_;iRow++) {
    364513      int iPivot=pivotVariable_[iRow];
    365       array[iRow]=cost_[iPivot];
     514      double value = cost_[iPivot];
     515      if (value) {
     516        array[iRow]=value;
     517        index[number++]=iRow;
     518      }
    366519    }
    367520  } else {
     
    372525      if (!pivoted(iPivot))
    373526        givenDjs[iPivot]=0.0;
    374       array[iRow]=cost_[iPivot]-givenDjs[iPivot];
    375     }
    376   }
    377   ClpDisjointCopyN ( array, numberRows_ , save);
     527      double value =cost_[iPivot]-givenDjs[iPivot];
     528      if (value) {
     529        array[iRow]=value;
     530        index[number++]=iRow;
     531      }
     532    }
     533  }
     534  arrayVector.setNumElements(number);
    378535
    379536  // Btran basic costs and get as accurate as possible
     
    381538  int iRefine;
    382539  double * work = workSpace->denseVector();
    383   factorization_->updateColumnTranspose(workSpace,array);
     540  CoinIndexedVector * thisVector = &arrayVector;
     541  CoinIndexedVector * lastVector = &previousVector;
     542  factorization_->updateColumnTranspose(workSpace,thisVector);
     543
    384544  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
    385545    // check basic reduced costs zero
     
    391551      for (iRow=0;iRow<numberRows_;iRow++) {
    392552        int iPivot=pivotVariable_[iRow];
     553        double value;
    393554        if (iPivot>=numberColumns_) {
    394555          // slack
    395           work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
     556          value = rowObjectiveWork_[iPivot-numberColumns_]
    396557            + array[iPivot-numberColumns_];
    397558        } else {
    398559          // column
    399         work[iRow] = reducedCostWork_[iPivot];
    400         }
    401         if (fabs(work[iRow])>largestDualError_) {
    402           largestDualError_=fabs(work[iRow]);
     560          value = reducedCostWork_[iPivot];
     561        }
     562        work[iRow]=value;
     563        if (fabs(value)>largestDualError_) {
     564          largestDualError_=fabs(value);
    403565        }
    404566      }
     
    424586    if (largestDualError_>=lastError) {
    425587      // restore
    426       double * temp = array;
    427       array = previous;
    428       previous=temp;
     588      CoinIndexedVector * temp = thisVector;
     589      thisVector = lastVector;
     590      lastVector=temp;
    429591      break;
    430592    }
    431     if (iRefine<numberRefinements_&&largestDualError_>1.0e-20
     593    if (iRefine<numberRefinements_&&largestDualError_>1.0e-10
    432594        &&!givenDjs) {
    433595      // try and make better
    434596      // save this
    435       double * temp = array;
    436       array = previous;
    437       previous=temp;
     597      CoinIndexedVector * temp = thisVector;
     598      thisVector = lastVector;
     599      lastVector=temp;
     600      int * indexOut = thisVector->getIndices();
     601      int number=0;
     602      array = thisVector->denseVector();
     603      thisVector->clear();
    438604      double multiplier = 131072.0;
    439605      for (iRow=0;iRow<numberRows_;iRow++) {
    440         array[iRow] = multiplier*work[iRow];
     606        double value = multiplier*work[iRow];
     607        if (value) {
     608          array[iRow]=value;
     609          indexOut[number++]=iRow;
     610          work[iRow]=0.0;
     611        }
    441612        work[iRow]=0.0;
    442613      }
     614      thisVector->setNumElements(number);
    443615      lastError=largestDualError_;
    444       factorization_->updateColumnTranspose(workSpace,array);
     616      factorization_->updateColumnTranspose(workSpace,thisVector);
    445617      multiplier = 1.0/multiplier;
     618      double * previous = lastVector->denseVector();
     619      number=0;
    446620      for (iRow=0;iRow<numberRows_;iRow++) {
    447         array[iRow] = previous[iRow] + multiplier*array[iRow];
    448         work[iRow]=0.0;
    449       }
     621        double value = previous[iRow] + multiplier*array[iRow];
     622        if (value) {
     623          array[iRow]=value;
     624          indexOut[number++]=iRow;
     625        } else {
     626          array[iRow]=0.0;
     627        }
     628      }
     629      thisVector->setNumElements(number);
    450630    } else {
    451631      break;
     
    454634  ClpFillN(work,numberRows_,0.0);
    455635  // now look at dual solution
     636  array = thisVector->denseVector();
    456637  for (iRow=0;iRow<numberRows_;iRow++) {
    457638    // slack
     
    469650  }
    470651
    471   delete [] previous;
    472   delete [] array;
    473   delete [] save;
    474652}
    475653/* Given an existing factorization computes and checks
     
    483661    createRim(7+8+16+32);
    484662    // do work
    485     gutsOfSolution ( rowActivities, columnActivities,NULL,NULL);
     663    gutsOfSolution ( NULL,NULL);
    486664    // release extra memory
    487665    deleteRim(0);
     
    519697/* Factorizes using current basis. 
    520698      solveType - 1 iterating, 0 initial, -1 external
     699      - 2 then iterating but can throw out of basis
    521700      If 10 added then in primal values pass
    522701*/
     
    528707  int iRow,iColumn;
    529708  int totalSlacks=numberRows_;
     709  if (!status_)
     710    createStatus();
    530711
    531712  bool valuesPass=false;
     
    535716  }
    536717#ifdef CLP_DEBUG
    537   if (solveType==1) {
     718  if (solveType>0) {
    538719    int numberFreeIn=0,numberFreeOut=0;
    539720    double biggestDj=0.0;
    540721    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    541722      switch(getColumnStatus(iColumn)) {
    542        
     723
    543724      case basic:
    544725        if (columnLower_[iColumn]<-largeValue_
     
    560741  }
    561742#endif
    562   if (!solveType) {
     743  if (solveType<=0) {
     744    // Make sure everything is clean
     745    for (iRow=0;iRow<numberRows_;iRow++) {
     746      if(getRowStatus(iRow)==isFixed) {
     747        // double check fixed
     748        if (rowUpperWork_[iRow]>rowLowerWork_[iRow])
     749          setRowStatus(iRow,atLowerBound);
     750      } else if (getRowStatus(iRow)==isFree) {
     751        // may not be free after all
     752        if (rowLowerWork_[iRow]>-largeValue_||rowUpperWork_[iRow]<largeValue_)
     753          setRowStatus(iRow,superBasic);
     754      }
     755    }
     756    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     757      if(getColumnStatus(iColumn)==isFixed) {
     758        // double check fixed
     759        if (columnUpperWork_[iColumn]>columnLowerWork_[iColumn])
     760          setColumnStatus(iColumn,atLowerBound);
     761      } else if (getColumnStatus(iColumn)==isFree) {
     762        // may not be free after all
     763        if (columnLowerWork_[iColumn]>-largeValue_||columnUpperWork_[iColumn]<largeValue_)
     764          setColumnStatus(iColumn,superBasic);
     765      }
     766    }
    563767    if (!valuesPass) {
    564768      // not values pass so set to bounds
     
    609813            break;
    610814          case isFree:
    611             // may not be free after all
    612             if (rowLowerWork_[iRow]<-largeValue_&&rowUpperWork_[iRow]>largeValue_)
    613815              break;
    614816            // not really free - fall through to superbasic
     
    682884            }
    683885            break;
     886          case isFixed:
    684887          case atLowerBound:
    685           case ClpSimplex::isFixed:
    686888            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    687889            if (columnActivityWork_[iColumn]<-largeValue_) {
     
    696898            break;
    697899          case isFree:
    698             // may not be free after all
    699             if (columnLowerWork_[iColumn]<-largeValue_&&columnUpperWork_[iColumn]>largeValue_)
    700900              break;
    701901            // not really free - fall through to superbasic
     
    8811081    }
    8821082  }
     1083#if 0
    8831084  int status=-99;
    8841085  int * rowIsBasic = new int[numberRows_];
     
    9071108    assert (numberBasic<=numberRows_);
    9081109    while (status==-99) {
    909       status =  factorization_->factorize(this,matrix_,
     1110      status =  factorization_->factorizeOld(this,matrix_,
    9101111                                          numberRows_,numberColumns_,
    9111112                                          rowIsBasic, columnIsBasic,
     
    9171118    }
    9181119    if (!status) {
    919   /// See whether to redo pivot order
    920       if (factorization_->needToReorder()||!numberIterations_) {
     1120      // See whether to redo pivot order
     1121      if (factorization_->needToReorder()||
     1122          progress_->lastIterationNumber(0)<0) {
    9211123        // do pivot information
    9221124        for (i=0;i<numberRows_;i++) {
     
    10181220  delete [] rowIsBasic;
    10191221  delete [] columnIsBasic;
     1222#else
     1223  int status = factorization_->factorize(this, solveType,valuesPass);
     1224#endif
    10201225  if (status) {
    10211226    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
     
    11081313  if (hitMaximumIterations())
    11091314    return 2;
     1315#if 0
     1316  // check for small cycles
     1317  int cycle=progress_->cycle(sequenceIn_,sequenceOut_,
     1318                            directionIn_,directionOut_);
     1319  if (cycle>0) {
     1320    printf("Cycle of %d\n",cycle);
     1321    if (factorization_->pivots()>cycle)
     1322      forceFactorization_=cycle-1;
     1323    else
     1324      forceFactorization_=1;
     1325    return 1;
     1326  }
     1327#endif
    11101328  // only time to re-factorize if one before real time
    11111329  // this is so user won't be surprised that maximumPivots has exact meaning
     
    11161334        factorization_->pivots()==forceFactorization_) {
    11171335      // relax
    1118       forceFactorization_ = (1+3*forceFactorization_)/2;
     1336      forceFactorization_ = (3+5*forceFactorization_)/4;
    11191337      if (forceFactorization_>factorization_->maximumPivots())
    11201338        forceFactorization_ = -1; //off
     
    12021420  numberFake_(0),
    12031421  progressFlag_(0),
    1204   firstFree_(-1)
     1422  firstFree_(-1),
     1423  incomingInfeasibility_(1.0),
     1424  allowedInfeasibility_(10.0),
     1425  progress_(NULL)
    12051426{
    12061427  int i;
     
    12961517  numberFake_(0),
    12971518  progressFlag_(0),
    1298   firstFree_(-1)
     1519  firstFree_(-1),
     1520  incomingInfeasibility_(1.0),
     1521  allowedInfeasibility_(10.0),
     1522  progress_(NULL)
    12991523{
    13001524  int i;
     
    13081532  // say Steepest pricing
    13091533  dualRowPivot_ = new ClpDualRowSteepest();
    1310   // say Dantzig pricing
    1311   primalColumnPivot_ = new ClpPrimalColumnDantzig();
     1534  // say Steepest pricing
     1535  primalColumnPivot_ = new ClpPrimalColumnSteepest();
    13121536  solveType_=1; // say simplex based life form
    13131537 
     
    13361560  rowUpperWork_ = upper_+numberColumns_;
    13371561  columnUpperWork_ = upper_;
     1562  //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
    13381563  cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows_));
    13391564  objectiveWork_ = cost_;
     
    14291654  progressFlag_ = rhs.progressFlag_;
    14301655  firstFree_ = rhs.firstFree_;
     1656  incomingInfeasibility_ = rhs.incomingInfeasibility_;
     1657  allowedInfeasibility_ = rhs.allowedInfeasibility_;
     1658  if (rhs.progress_)
     1659    progress_ = new ClpSimplexProgress(*rhs.progress_);
     1660  else
     1661    progress_=NULL;
    14311662  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    14321663  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     
    14921723    delete primalColumnPivot_;
    14931724    primalColumnPivot_ = NULL;
     1725    delete progress_;
     1726    progress_=NULL;
    14941727  } else {
    14951728    // delete any size information in methods
     
    15121745
    15131746  objectiveValue_ = 0.0;
    1514   firstFree_ = -1;
    15151747  // now look at primal solution
    15161748  columnPrimalInfeasibility_=0.0;
     
    15901822  rowDualInfeasibility_=0.0;
    15911823  rowDualSequence_=-1;
    1592   int firstFree = -1;
     1824  int firstFreePrimal = -1;
     1825  int firstFreeDual = -1;
     1826  int numberSuperBasicWithDj=0;
    15931827  primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
    15941828                                   columnPrimalInfeasibility_);
     
    16021836
    16031837  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1604     if (getColumnStatus(iColumn) != basic) {
     1838    if (getColumnStatus(iColumn) != basic&&!flagged(iColumn)) {
    16051839      // not basic
    16061840      double distanceUp = columnUpperWork_[iColumn]-
     
    16101844      if (distanceUp>primalTolerance_) {
    16111845        double value = reducedCostWork_[iColumn];
    1612         if (!flagged(iColumn)) {
    1613           // Check if "free"
    1614           if (firstFree<0&&distanceDown>primalTolerance_) {
    1615             if (algorithm_>0) {
    1616               // primal
    1617               firstFree = iColumn;
    1618             } else if (fabs(value)>1.0e2*relaxedTolerance) {
    1619               // dual with dj
    1620               firstFree = iColumn;
    1621             }
     1846        // Check if "free"
     1847        if (distanceDown>primalTolerance_) {
     1848          if (fabs(value)>1.0e2*relaxedTolerance) {
     1849            numberSuperBasicWithDj++;
     1850            if (firstFreeDual<0)
     1851              firstFreeDual = iColumn;
    16221852          }
     1853          if (firstFreePrimal<0)
     1854            firstFreePrimal = iColumn;
    16231855        }
    16241856        // should not be negative
     
    16831915  }
    16841916  for (iRow=0;iRow<numberRows_;iRow++) {
    1685     if (getRowStatus(iRow) != basic) {
     1917    if (getRowStatus(iRow) != basic&&!flagged(iRow+numberColumns_)) {
    16861918      // not basic
    16871919      double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow];
     
    16891921      if (distanceUp>primalTolerance_) {
    16901922        double value = rowReducedCost_[iRow];
    1691         if (!flagged(iRow+numberColumns_)) {
    1692           // Check if "free"
    1693           if (firstFree<0&&distanceDown>primalTolerance_) {
    1694             if (algorithm_>0) {
    1695               // primal
    1696               firstFree = iRow+numberColumns_;
    1697             } else if (fabs(value)>1.0e2*relaxedTolerance) {
    1698               // dual with dj
    1699               firstFree = iRow+numberColumns_;
    1700             }
     1923        // Check if "free"
     1924        if (distanceDown>primalTolerance_) {
     1925          if (fabs(value)>1.0e2*relaxedTolerance) {
     1926            numberSuperBasicWithDj++;
     1927            if (firstFreeDual<0)
     1928              firstFreeDual = iRow+numberColumns_;
    17011929          }
     1930          if (firstFreePrimal<0)
     1931            firstFreePrimal = iRow+numberColumns_;
    17021932        }
    17031933        // should not be negative
     
    17511981    }
    17521982  }
    1753   if (firstFree>=0) {
    1754     if (algorithm_<0||!numberIterations_)
    1755       firstFree_=firstFree;
     1983  if (algorithm_<0&&firstFreeDual>=0) {
     1984    // dual
     1985    firstFree_ = firstFreeDual;
     1986  } else if (numberSuperBasicWithDj||
     1987             (progress_&&progress_->lastIterationNumber(0)<=0)) {
     1988    firstFree_=firstFreePrimal;
    17561989  }
    17571990}
     
    17832016  }
    17842017}
     2018/*
     2019  Unpacks one column of the matrix into indexed array
     2020*/
     2021void
     2022ClpSimplex::unpackPacked(CoinIndexedVector * rowArray)
     2023{
     2024  rowArray->clear();
     2025  if (sequenceIn_>=numberColumns_) {
     2026    //slack
     2027    int * index = rowArray->getIndices();
     2028    double * array = rowArray->denseVector();
     2029    array[0]=-1.0;
     2030    index[0]=sequenceIn_-numberColumns_;
     2031    rowArray->setNumElements(1);
     2032    rowArray->setPackedMode(true);
     2033  } else {
     2034    // column
     2035    matrix_->unpackPacked(this,rowArray,sequenceIn_);
     2036  }
     2037}
     2038void
     2039ClpSimplex::unpackPacked(CoinIndexedVector * rowArray,int sequence)
     2040{
     2041  rowArray->clear();
     2042  if (sequence>=numberColumns_) {
     2043    //slack
     2044    int * index = rowArray->getIndices();
     2045    double * array = rowArray->denseVector();
     2046    array[0]=-1.0;
     2047    index[0]=sequence-numberColumns_;
     2048    rowArray->setNumElements(1);
     2049    rowArray->setPackedMode(true);
     2050  } else {
     2051    // column
     2052    matrix_->unpackPacked(this,rowArray,sequence);
     2053  }
     2054}
    17852055bool
    17862056ClpSimplex::createRim(int what,bool makeRowCopy)
    17872057{
    17882058  bool goodMatrix=true;
     2059  int saveLevel=handler_->logLevel();
     2060  if (problemStatus_==10)
     2061    handler_->setLogLevel(0); // switch off messages
    17892062  int i;
    17902063  if ((what&(16+32))!=0) {
     
    18452118  if ((what&4)!=0) {
    18462119    delete [] cost_;
    1847     cost_ = new double[numberColumns_+numberRows_];
     2120    // extra copy with original costs
     2121    int nTotal = numberRows_+numberColumns_;
     2122    //cost_ = new double[2*nTotal];
     2123    cost_ = new double[nTotal];
    18482124    objectiveWork_ = cost_;
    18492125    rowObjectiveWork_ = cost_+numberColumns_;
     
    18532129    else
    18542130      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
     2131    // and initialize changes to zero
     2132    //memset(cost_+nTotal,0,nTotal*sizeof(double));
    18552133  }
    18562134  if ((what&1)!=0) {
     
    19152193      if (columnUpperWork_[i]<1.0e50)
    19162194        columnUpperWork_[i] *= multiplier;
     2195     
    19172196    }
    19182197    for (i=0;i<numberRows_;i++) {
     
    19712250    }   
    19722251  }
     2252  double primalTolerance=dblParam_[ClpPrimalTolerance];
     2253  if ((what&1)!=0) {
     2254    // fix any variables with tiny gaps
     2255    for (i=0;i<numberColumns_;i++) {
     2256      if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
     2257        if (columnLowerWork_[i]>=0.0) {
     2258          columnUpperWork_[i] = columnLowerWork_[i];
     2259        } else if (columnUpperWork_[i]<=0.0) {
     2260          columnLowerWork_[i] = columnUpperWork_[i];
     2261        } else {
     2262          columnUpperWork_[i] = 0.0;
     2263          columnLowerWork_[i] = 0.0;
     2264        }
     2265      }
     2266    }
     2267    for (i=0;i<numberRows_;i++) {
     2268      if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
     2269        if (rowLowerWork_[i]>=0.0) {
     2270          rowUpperWork_[i] = rowLowerWork_[i];
     2271        } else if (rowUpperWork_[i]<=0.0) {
     2272          rowLowerWork_[i] = rowUpperWork_[i];
     2273        } else {
     2274          rowUpperWork_[i] = 0.0;
     2275          rowLowerWork_[i] = 0.0;
     2276        }
     2277      }
     2278    }
     2279  }
     2280  if (problemStatus_==10) {
     2281    problemStatus_=-1;
     2282    handler_->setLogLevel(saveLevel); // switch back messages
     2283  }
    19732284  return goodMatrix;
    19742285}
     
    20372348    infeasibilityCost_=value;
    20382349}
    2039 void ClpSimplex::setnumberRefinements( int value)
     2350void ClpSimplex::setNumberRefinements( int value)
    20402351{
    20412352  if (value>=0&&value<10)
     
    21792490    }
    21802491  }
    2181 #define MAXPASS 10 
     2492#define MAXPASS 10
    21822493
    21832494  // Loop round seeing if we can tighten bounds
    21842495  // Would be faster to have a stack of possible rows
    21852496  // and we put altered rows back on stack
    2186 
    2187   while(numberChanged) {
     2497  int numberCheck=-1;
     2498  while(numberChanged>numberCheck) {
    21882499
    21892500    numberChanged = 0; // Bounds tightened this pass
     
    22052516        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
    22062517        CoinBigIndex j;
    2207 
    22082518        // Compute possible lower and upper ranges
    2209 
     2519     
    22102520        for (j = rStart; j < rEnd; ++j) {
    22112521          double value=element[j];
     
    22132523          if (value > 0.0) {
    22142524            if (columnUpper_[iColumn] >= large) {
    2215               maximumUp = 1e31;
    22162525              ++infiniteUpper;
    22172526            } else {
     
    22192528            }
    22202529            if (columnLower_[iColumn] <= -large) {
    2221               maximumDown = -1e31;
    22222530              ++infiniteLower;
    22232531            } else {
     
    22262534          } else if (value<0.0) {
    22272535            if (columnUpper_[iColumn] >= large) {
    2228               maximumDown = -1e31;
    22292536              ++infiniteLower;
    22302537            } else {
     
    22322539            }
    22332540            if (columnLower_[iColumn] <= -large) {
    2234               maximumUp = 1e31;
    22352541              ++infiniteUpper;
    22362542            } else {
     
    22392545          }
    22402546        }
    2241         if (maximumUp <= rowUpper_[iRow] + tolerance &&
    2242             maximumDown >= rowLower_[iRow] - tolerance) {
    2243 
     2547        double maxUp = maximumUp+infiniteUpper*1.0e31;
     2548        double maxDown = maximumDown-infiniteLower*1.0e31;
     2549        if (maxUp <= rowUpper_[iRow] + tolerance &&
     2550            maxDown >= rowLower_[iRow] - tolerance) {
     2551         
    22442552          // Row is redundant - make totally free
    2245           rowLower_[iRow]=-COIN_DBL_MAX;
    2246           rowUpper_[iRow]=COIN_DBL_MAX;
     2553          // NO - can't do this for postsolve
     2554          // rowLower_[iRow]=-COIN_DBL_MAX;
     2555          // rowUpper_[iRow]=COIN_DBL_MAX;
     2556          //printf("Redundant row in presolveX %d\n",iRow);
    22472557
    22482558        } else {
    2249           if (maximumUp < rowLower_[iRow] -tolerance ||
    2250               maximumDown > rowUpper_[iRow]+tolerance) {
     2559          if (maxUp < rowLower_[iRow] -tolerance ||
     2560              maxDown > rowUpper_[iRow]+tolerance) {
    22512561            // problem is infeasible - exit at once
    22522562            numberInfeasible++;
    2253             if (handler_->logLevel()>3)
    2254               printf("a %g %g %g %g\n",maximumUp,rowLower_[iRow],
    2255                      maximumDown,rowUpper_[iRow]);
    22562563            break;
    22572564          }
    2258 
    2259           if (infiniteUpper == 0 && rowLower_[iRow] > -large) {
    2260             for (j = rStart; j < rEnd; ++j) {
    2261               double value=element[j];
    2262               iColumn = column[j];
    2263               if (value > 0.0) {
    2264                 if (columnUpper_[iColumn] < large) {
    2265                   newBound = columnUpper_[iColumn] +
    2266                     (rowLower_[iRow] - maximumUp) / value;
    2267                   if (newBound > columnLower_[iColumn] + 1.0e-12) {
    2268                     // Tighten the lower bound
    2269 
    2270                     columnLower_[iColumn] = newBound;
    2271                     ++numberChanged;
    2272 
    2273                     // check infeasible (relaxed)
    2274                     if (columnUpper_[iColumn] - columnLower_[iColumn] <
    2275                         -100.0*tolerance) {
    2276                       numberInfeasible++;
    2277                       if (handler_->logLevel()>3)
    2278                         printf("b %g %g \n",columnLower_[iColumn],columnUpper_[iColumn]);
    2279                     }
    2280                     infiniteLower=1; // skip looking at other bound
     2565          double lower = rowLower_[iRow];
     2566          double upper = rowUpper_[iRow];
     2567          for (j = rStart; j < rEnd; ++j) {
     2568            double value=element[j];
     2569            iColumn = column[j];
     2570            double nowLower = columnLower_[iColumn];
     2571            double nowUpper = columnUpper_[iColumn];
     2572            if (value > 0.0) {
     2573              // positive value
     2574              if (lower>-large) {
     2575                if (!infiniteUpper) {
     2576                  assert(nowUpper < large);
     2577                  newBound = nowUpper +
     2578                    (lower - maximumUp) / value;
     2579                } else if (infiniteUpper==1&&nowUpper>large) {
     2580                  newBound = (lower -maximumUp) / value;
     2581                } else {
     2582                  newBound = -COIN_DBL_MAX;
     2583                }
     2584                if (newBound > nowLower + 1.0e-12) {
     2585                  // Tighten the lower bound
     2586                  columnLower_[iColumn] = newBound;
     2587                  numberChanged++;
     2588                  // check infeasible (relaxed)
     2589                  if (nowUpper - newBound <
     2590                      -100.0*tolerance) {
     2591                    numberInfeasible++;
    22812592                  }
     2593                  // adjust
     2594                  double now;
     2595                  if (nowLower<-large) {
     2596                    now=0.0;
     2597                    infiniteLower--;
     2598                  } else {
     2599                    now = nowLower;
     2600                  }
     2601                  maximumDown += (newBound-now) * value;
     2602                  nowLower = newBound;
    22822603                }
    2283               } else {
    2284                 if (columnLower_[iColumn] > -large) {
    2285                   newBound = columnLower_[iColumn] +
    2286                     (rowLower_[iRow] - maximumUp) / value;
    2287                   if (newBound < columnUpper_[iColumn] - 1.0e-12) {
    2288                     // Tighten the upper bound
    2289 
    2290                     columnUpper_[iColumn] = newBound;
    2291                     ++numberChanged;
    2292 
    2293                     // check infeasible (relaxed)
    2294                     if (columnUpper_[iColumn] - columnLower_[iColumn] <
    2295                         -100.0*tolerance) {
    2296                       numberInfeasible++;
    2297                       if (handler_->logLevel()>3)
    2298                         printf("c %g %g \n",columnLower_[iColumn],columnUpper_[iColumn]);
    2299                     }
    2300                     infiniteLower=1; // skip looking at other bound
     2604              }
     2605              if (upper <large) {
     2606                if (!infiniteLower) {
     2607                  assert(nowLower >- large);
     2608                  newBound = nowLower +
     2609                    (upper - maximumDown) / value;
     2610                } else if (infiniteLower==1&&nowLower<-large) {
     2611                  newBound =   (upper - maximumDown) / value;
     2612                } else {
     2613                  newBound = COIN_DBL_MAX;
     2614                }
     2615                if (newBound < nowUpper - 1.0e-12) {
     2616                  // Tighten the upper bound
     2617                  columnUpper_[iColumn] = newBound;
     2618                  numberChanged++;
     2619                  // check infeasible (relaxed)
     2620                  if (newBound - nowLower <
     2621                      -100.0*tolerance) {
     2622                    numberInfeasible++;
    23012623                  }
     2624                  // adjust
     2625                  double now;
     2626                  if (nowUpper>large) {
     2627                    now=0.0;
     2628                    infiniteUpper--;
     2629                  } else {
     2630                    now = nowUpper;
     2631                  }
     2632                  maximumUp += (newBound-now) * value;
     2633                  nowUpper = newBound;
     2634                }
     2635              }
     2636            } else {
     2637              // negative value
     2638              if (lower>-large) {
     2639                if (!infiniteUpper) {
     2640                  assert(nowLower >- large);
     2641                  newBound = nowLower +
     2642                    (lower - maximumUp) / value;
     2643                } else if (infiniteUpper==1&&nowLower<-large) {
     2644                  newBound = (lower -maximumUp) / value;
     2645                } else {
     2646                  newBound = COIN_DBL_MAX;
     2647                }
     2648                if (newBound < nowUpper - 1.0e-12) {
     2649                  // Tighten the upper bound
     2650                  columnUpper_[iColumn] = newBound;
     2651                  numberChanged++;
     2652                  // check infeasible (relaxed)
     2653                  if (newBound - nowLower <
     2654                      -100.0*tolerance) {
     2655                    numberInfeasible++;
     2656                  }
     2657                  // adjust
     2658                  double now;
     2659                  if (nowUpper>large) {
     2660                    now=0.0;
     2661                    infiniteLower--;
     2662                  } else {
     2663                    now = nowUpper;
     2664                  }
     2665                  maximumDown += (newBound-now) * value;
     2666                  nowUpper = newBound;
     2667                }
     2668              }
     2669              if (upper <large) {
     2670                if (!infiniteLower) {
     2671                  assert(nowUpper < large);
     2672                  newBound = nowUpper +
     2673                    (upper - maximumDown) / value;
     2674                } else if (infiniteLower==1&&nowUpper>large) {
     2675                  newBound =   (upper - maximumDown) / value;
     2676                } else {
     2677                  newBound = -COIN_DBL_MAX;
     2678                }
     2679                if (newBound > nowLower + 1.0e-12) {
     2680                  // Tighten the lower bound
     2681                  columnLower_[iColumn] = newBound;
     2682                  numberChanged++;
     2683                  // check infeasible (relaxed)
     2684                  if (nowUpper - newBound <
     2685                      -100.0*tolerance) {
     2686                    numberInfeasible++;
     2687                  }
     2688                  // adjust
     2689                  double now;
     2690                  if (nowLower<-large) {
     2691                    now=0.0;
     2692                    infiniteUpper--;
     2693                  } else {
     2694                    now = nowLower;
     2695                  }
     2696                  maximumUp += (newBound-now) * value;
     2697                  nowLower = newBound;
    23022698                }