Changeset 225


Ignore:
Timestamp:
Oct 15, 2003 5:34:57 PM (16 years ago)
Author:
forrest
Message:

This should break everything

Location:
trunk
Files:
24 added
33 deleted
57 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpDualRowDantzig.cpp

    r63 r225  
    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);
  • trunk/ClpDualRowPivot.cpp

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

    r137 r225  
    2424    infeasible_(NULL),
    2525    alternateWeights_(NULL),
    26     savedWeights_(NULL)
     26    savedWeights_(NULL),
     27    dubiousWeights_(NULL)
    2728{
    2829  type_=2+64*mode;
     
    6162    savedWeights_=NULL;
    6263  }
     64  if (rhs.dubiousWeights_) {
     65    assert(model_);
     66    int number = model_->numberRows();
     67    dubiousWeights_= new int[number];
     68    ClpDisjointCopyN(rhs.dubiousWeights_,number,dubiousWeights_);
     69  } else {
     70    dubiousWeights_=NULL;
     71  }
    6372}
    6473
     
    6978{
    7079  delete [] weights_;
     80  delete [] dubiousWeights_;
    7181  delete infeasible_;
    7282  delete alternateWeights_;
     
    8696    model_ = rhs.model_;
    8797    delete [] weights_;
     98    delete [] dubiousWeights_;
    8899    delete infeasible_;
    89100    delete alternateWeights_;
     
    111122    } else {
    112123      savedWeights_=NULL;
     124    }
     125    if (rhs.dubiousWeights_) {
     126      assert(model_);
     127      int number = model_->numberRows();
     128      dubiousWeights_= new int[number];
     129      ClpDisjointCopyN(rhs.dubiousWeights_,number,dubiousWeights_);
     130    } else {
     131      dubiousWeights_=NULL;
    113132    }
    114133  }
     
    135154  // allow tolerance at least slightly bigger than standard
    136155  tolerance = tolerance +  error;
     156  // But cap
     157  tolerance = min(1000.0,tolerance);
    137158  tolerance *= tolerance; // as we are using squares
    138159  double * solution = model_->solutionRegion();
     
    156177      if (iPivot<numberColumns)
    157178        value *= COLUMN_BIAS; // bias towards columns
     179k
    158180#endif
    159181      // store square in list
     
    181203    number = infeasible_->getNumElements();
    182204  }
    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
     205  if(model_->numberIterations()<model_->lastBadIteration()+200) {
     206    // we can't really trust infeasibilities if there is dual error
     207    double checkTolerance = 1.0e-8;
     208    if (!model_->factorization()->pivots())
     209      checkTolerance = 1.0e-6;
     210    if (model_->largestPrimalError()>checkTolerance)
     211      tolerance *= model_->largestPrimalError()/checkTolerance;
     212  }
     213  int numberWanted;
     214  if (mode_<2 ) {
     215    numberWanted = number+1;
     216  } else if (mode_==2) {
     217    numberWanted = max(2000,number/8);
     218  } else {
     219    int numberElements = model_->factorization()->numberElements();
     220    double ratio = (double) numberElements/(double) model_->numberRows();
     221    numberWanted = max(2000,number/8);
     222    if (ratio<1.0) {
     223      numberWanted = max(2000,number/20);
     224    } else if (ratio>10.0) {
     225      ratio = number * (ratio/80.0);
     226      if (ratio>number)
     227        numberWanted=number+1;
     228      else
     229        numberWanted = max(2000,(int) ratio);
     230    }
     231  }
     232  int iPass;
     233  // Setup two passes
     234  int start[4];
     235  start[1]=number;
     236  start[2]=0;
     237  double dstart = ((double) number) * CoinDrand48();
     238  start[0]=(int) dstart;
     239  start[3]=start[0];
     240  //double largestWeight=0.0;
     241  //double smallestWeight=1.0e100;
     242  for (iPass=0;iPass<2;iPass++) {
     243    int end = start[2*iPass+1];
     244    for (i=start[2*iPass];i<end;i++) {
     245      iRow = index[i];
     246      double value = infeas[iRow];
     247      if (value>tolerance) {
     248        //#define OUT_EQ
     249#ifdef OUT_EQ
     250        {
     251          int iSequence = pivotVariable[iRow];
     252          if (upper[iSequence]==lower[iSequence])
     253            value *= 2.0;
     254        }
     255#endif
     256        double weight = weights_[iRow];
     257        //largestWeight = max(largestWeight,weight);
     258        //smallestWeight = min(smallestWeight,weight);
     259        //double dubious = dubiousWeights_[iRow];
     260        //weight *= dubious;
     261        //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
     262        if (value>largest*weight) {
     263          // make last pivot row last resort choice
     264          if (iRow==lastPivotRow) {
     265            if (value*1.0e-10<largest*weight)
     266              continue;
     267            else
     268              value *= 1.0e-10;
     269          }
     270          int iSequence = pivotVariable[iRow];
     271          if (!model_->flagged(iSequence)) {
     272            //#define CLP_DEBUG 1
    197273#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];
     274            double value2=0.0;
     275            if (solution[iSequence]>upper[iSequence]+tolerance)
     276              value2=solution[iSequence]-upper[iSequence];
     277            else if (solution[iSequence]<lower[iSequence]-tolerance)
     278              value2=solution[iSequence]-lower[iSequence];
     279            assert(fabs(value2*value2-infeas[iRow])<1.0e-8*min(value2*value2,infeas[iRow]));
     280#endif
     281            if (solution[iSequence]>upper[iSequence]+tolerance||
     282                solution[iSequence]<lower[iSequence]-tolerance) {
     283              chosenRow=iRow;
     284              largest=value/weight;
     285              //largestWeight = dubious;
     286            }
     287          } else {
     288            // just to make sure we don't exit before got something
     289            numberWanted++;
     290          }
    209291        }
    210       }
    211     }
    212   }
     292        numberWanted--;
     293        if (!numberWanted)
     294          break;
     295      }
     296    }
     297    if (!numberWanted)
     298      break;
     299  }
     300  //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
    213301  return chosenRow;
    214302}
    215303#define TRY_NORM 1.0e-4
    216 // Updates weights
    217 void
     304// Updates weights and returns pivot alpha
     305double
    218306ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
    219307                                  CoinIndexedVector * spare,
    220308                                  CoinIndexedVector * updatedColumn)
    221309{
    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;
    233310#if CLP_DEBUG>2
    234311  // Very expensive debug
     
    240317    double * array = alternateWeights_->denseVector();
    241318    int * which = alternateWeights_->getIndices();
     319    int i;
    242320    for (i=0;i<numberRows;i++) {
    243321      double value=0.0;
     
    257335      if (fabs(weights_[i]-value)>1.0e-4)
    258336        printf("%d old %g, true %g\n",i,weights_[i],value);
     337      //else
     338      //printf("%d matches %g\n",i,value);
    259339    }
    260340    delete temp;
    261341  }
    262342#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();
     343  assert (input->packedMode());
     344  assert (updatedColumn->packedMode());
     345  double alpha=0.0;
     346  if (!model_->factorization()->networkBasis()) {
     347    // clear other region
     348    alternateWeights_->clear();
     349    double norm = 0.0;
     350    int i;
     351    double * work = input->denseVector();
     352    int numberNonZero = input->getNumElements();
     353    int * which = input->getIndices();
     354    double * work2 = spare->denseVector();
     355    int * which2 = spare->getIndices();
     356    // ftran
     357    //permute and move indices into index array
     358    //also compute norm
     359    //int *regionIndex = alternateWeights_->getIndices (  );
     360    const int *permute = model_->factorization()->permute();
     361    //double * region = alternateWeights_->denseVector();
     362    for ( i = 0; i < numberNonZero; i ++ ) {
     363      int iRow = which[i];
     364      double value = work[i];
     365      norm += value*value;
     366      iRow = permute[iRow];
     367      work2[iRow] = value;
     368      which2[i] = iRow;
     369    }
     370    spare->setNumElements ( numberNonZero );
     371    // Only one array active as already permuted
     372    model_->factorization()->updateColumn(spare,spare,true);
     373    // permute back
     374    numberNonZero = spare->getNumElements();
     375    // alternateWeights_ should still be empty
     376    int pivotRow = model_->pivotRow();
    275377#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.) {
     378    if ( model_->logLevel (  ) >4  &&
     379         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm))
     380      printf("on row %d, true weight %g, old %g\n",
     381             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
     382#endif
     383    // could re-initialize here (could be expensive)
     384    norm /= model_->alpha() * model_->alpha();
     385   
     386    assert(norm);
     387    // pivot element
     388    alpha=0.0;
     389    double multiplier = 2.0 / model_->alpha();
     390    // look at updated column
     391    work = updatedColumn->denseVector();
     392    numberNonZero = updatedColumn->getNumElements();
     393    which = updatedColumn->getIndices();
     394   
     395    int nSave=0;
     396    double * work3 = alternateWeights_->denseVector();
     397    int * which3 = alternateWeights_->getIndices();
     398    const int * pivotColumn = model_->factorization()->pivotColumn();
     399    for (i =0; i < numberNonZero; i++) {
     400      int iRow = which[i];
     401      double theta = work[i];
     402      if (iRow==pivotRow)
     403        alpha = theta;
     404      double devex = weights_[iRow];
     405      work3[nSave]=devex; // save old
     406      which3[nSave++]=iRow;
     407      // transform to match spare
     408      int jRow = pivotColumn[iRow];
     409      double value = work2[jRow];
     410      devex +=  theta * (theta*norm+value * multiplier);
     411      if (devex < TRY_NORM)
     412        devex = TRY_NORM;
     413      weights_[iRow]=devex;
     414    }
     415    assert (alpha);
     416    alternateWeights_->setPackedMode(true);
     417    alternateWeights_->setNumElements(nSave);
     418    if (norm < TRY_NORM)
     419      norm = TRY_NORM;
     420    weights_[pivotRow] = norm;
     421    spare->clear();
     422#ifdef CLP_DEBUG
     423    spare->checkClear();
     424#endif
     425  } else {
     426    // clear other region
     427    alternateWeights_->clear();
     428    double norm = 0.0;
     429    int i;
     430    double * work = input->denseVector();
     431    int number = input->getNumElements();
     432    int * which = input->getIndices();
     433    double * work2 = spare->denseVector();
     434    int * which2 = spare->getIndices();
     435    for (i=0;i<number;i++) {
     436      int iRow = which[i];
     437      double value = work[i];
     438      norm += value*value;
     439      work2[iRow]=value;
     440      which2[i]=iRow;
     441    }
     442    spare->setNumElements(number);
     443    // ftran
     444    model_->factorization()->updateColumn(alternateWeights_,spare);
     445    // alternateWeights_ should still be empty
     446    int pivotRow = model_->pivotRow();
     447#ifdef CLP_DEBUG
     448    if ( model_->logLevel (  ) >4  &&
     449         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm))
     450      printf("on row %d, true weight %g, old %g\n",
     451             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
     452#endif
     453    // could re-initialize here (could be expensive)
     454    norm /= model_->alpha() * model_->alpha();
     455   
     456    assert(norm);
     457    //if (norm < TRY_NORM)
     458    //norm = TRY_NORM;
     459    // pivot element
     460    alpha=0.0;
    288461    double multiplier = 2.0 / model_->alpha();
    289462    // look at updated column
     
    291464    number = updatedColumn->getNumElements();
    292465    which = updatedColumn->getIndices();
    293 
     466   
    294467    int nSave=0;
    295 
     468    double * work3 = alternateWeights_->denseVector();
     469    int * which3 = alternateWeights_->getIndices();
    296470    for (i =0; i < number; i++) {
    297471      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
     472      double theta = work[i];
     473      if (iRow==pivotRow)
     474        alpha = theta;
     475      double devex = weights_[iRow];
     476      work3[nSave]=devex; // save old
     477      which3[nSave++]=iRow;
     478      double value = work2[iRow];
     479      devex +=  theta * (theta*norm+value * multiplier);
     480      if (devex < TRY_NORM)
     481        devex = TRY_NORM;
     482      weights_[iRow]=devex;
     483    }
     484    assert (alpha);
     485    alternateWeights_->setPackedMode(true);
    313486    alternateWeights_->setNumElements(nSave);
    314487    if (norm < TRY_NORM)
    315488      norm = TRY_NORM;
    316489    weights_[pivotRow] = norm;
    317   }
    318   spare->clear();
     490    spare->clear();
     491  }
    319492#ifdef CLP_DEBUG
    320493  spare->checkClear();
    321494#endif
     495  return alpha;
    322496}
    323497 
     
    345519  int numberColumns = model_->numberColumns();
    346520#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
     521  if (primalUpdate->packedMode()) {
     522    for (i=0;i<number;i++) {
     523      int iRow=which[i];
     524      int iPivot=pivotVariable[iRow];
     525      double value = solution[iPivot];
     526      double cost = model_->cost(iPivot);
     527      double change = primalRatio*work[i];
     528      work[i]=0.0;
     529      value -= change;
     530      changeObj -= change*cost;
     531      solution[iPivot] = value;
     532      double lower = model_->lower(iPivot);
     533      double upper = model_->upper(iPivot);
     534      // But if pivot row then use value of incoming
     535      // Although it is safer to recompute before next selection
     536      // in case something odd happens
     537      if (iRow==pivotRow) {
     538        iPivot = model_->sequenceIn();
     539        lower = model_->lower(iPivot);
     540        upper = model_->upper(iPivot);
     541        value = model_->valueIncomingDual();
     542      }
     543      if (value<lower-tolerance) {
     544        value -= lower;
     545        value *= value;
     546#ifdef COLUMN_BIAS
     547        if (iPivot<numberColumns)
     548          value *= COLUMN_BIAS; // bias towards columns
    373549#endif
    374550#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
     551        if (lower==upper)
     552          value *= FIXED_BIAS; // bias towards taking out fixed variables
     553#endif
     554        // store square in list
     555        if (infeas[iRow])
     556          infeas[iRow] = value; // already there
     557        else
     558          infeasible_->quickAdd(iRow,value);
     559      } else if (value>upper+tolerance) {
     560        value -= upper;
     561        value *= value;
     562#ifdef COLUMN_BIAS
     563        if (iPivot<numberColumns)
     564          value *= COLUMN_BIAS; // bias towards columns
    389565#endif
    390566#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;
     567        if (lower==upper)
     568          value *= FIXED_BIAS; // bias towards taking out fixed variables
     569#endif
     570        // store square in list
     571        if (infeas[iRow])
     572          infeas[iRow] = value; // already there
     573        else
     574          infeasible_->quickAdd(iRow,value);
     575      } else {
     576        // feasible - was it infeasible - if so set tiny
     577        if (infeas[iRow])
     578          infeas[iRow] = 1.0e-100;
     579      }
     580    }
     581  } else {
     582    for (i=0;i<number;i++) {
     583      int iRow=which[i];
     584      int iPivot=pivotVariable[iRow];
     585      double value = solution[iPivot];
     586      double cost = model_->cost(iPivot);
     587      double change = primalRatio*work[iRow];
     588      value -= change;
     589      changeObj -= change*cost;
     590      solution[iPivot] = value;
     591      double lower = model_->lower(iPivot);
     592      double upper = model_->upper(iPivot);
     593      // But if pivot row then use value of incoming
     594      // Although it is safer to recompute before next selection
     595      // in case something odd happens
     596      if (iRow==pivotRow) {
     597        iPivot = model_->sequenceIn();
     598        lower = model_->lower(iPivot);
     599        upper = model_->upper(iPivot);
     600        value = model_->valueIncomingDual();
     601      }
     602      if (value<lower-tolerance) {
     603        value -= lower;
     604        value *= value;
     605#ifdef COLUMN_BIAS
     606        if (iPivot<numberColumns)
     607          value *= COLUMN_BIAS; // bias towards columns
     608#endif
     609#ifdef FIXED_BIAS
     610        if (lower==upper)
     611          value *= FIXED_BIAS; // bias towards taking out fixed variables
     612#endif
     613        // store square in list
     614        if (infeas[iRow])
     615          infeas[iRow] = value; // already there
     616        else
     617          infeasible_->quickAdd(iRow,value);
     618      } else if (value>upper+tolerance) {
     619        value -= upper;
     620        value *= value;
     621#ifdef COLUMN_BIAS
     622        if (iPivot<numberColumns)
     623          value *= COLUMN_BIAS; // bias towards columns
     624#endif
     625#ifdef FIXED_BIAS
     626        if (lower==upper)
     627          value *= FIXED_BIAS; // bias towards taking out fixed variables
     628#endif
     629        // store square in list
     630        if (infeas[iRow])
     631          infeas[iRow] = value; // already there
     632        else
     633          infeasible_->quickAdd(iRow,value);
     634      } else {
     635        // feasible - was it infeasible - if so set tiny
     636        if (infeas[iRow])
     637          infeas[iRow] = 1.0e-100;
     638      }
     639      work[iRow]=0.0;
     640    }
    405641  }
    406642  primalUpdate->setNumElements(0);
     
    442678      alternateWeights_->reserve(numberRows+
    443679                                 model_->factorization()->maximumPivots());
    444       if (!mode_||mode==5) {
     680      if (!mode_||mode_==2||mode==5) {
    445681        // initialize to 1.0 (can we do better?)
    446682        for (i=0;i<numberRows;i++) {
     
    455691        for (i=0;i<numberRows;i++) {
    456692          double value=0.0;
    457           array[i] = 1.0;
     693          array[0] = 1.0;
    458694          which[0] = i;
    459695          alternateWeights_->setNumElements(1);
     696          alternateWeights_->setPackedMode(true);
    460697          model_->factorization()->updateColumnTranspose(temp,
    461698                                                         alternateWeights_);
     
    463700          int j;
    464701          for (j=0;j<number;j++) {
    465             int iRow=which[j];
    466             value+=array[iRow]*array[iRow];
    467             array[iRow]=0.0;
     702            value+=array[j]*array[j];
     703            array[j]=0.0;
    468704          }
    469705          alternateWeights_->setNumElements(0);
     
    519755  }
    520756  if (mode>=2) {
     757    // Get dubious weights
     758    //if (!dubiousWeights_)
     759    //dubiousWeights_=new int[numberRows];
     760    //model_->factorization()->getWeights(dubiousWeights_);
    521761    infeasible_->clear();
    522762    int iRow;
     
    566806  int * which = alternateWeights_->getIndices();
    567807  int i;
    568   for (i=0;i<number;i++) {
    569     int iRow = which[i];
    570     weights_[iRow]=saved[iRow];
    571     saved[iRow]=0.0;
     808  if (alternateWeights_->packedMode()) {
     809    for (i=0;i<number;i++) {
     810      int iRow = which[i];
     811      weights_[iRow]=saved[i];
     812      saved[i]=0.0;
     813    }
     814  } else {
     815    for (i=0;i<number;i++) {
     816      int iRow = which[i];
     817      weights_[iRow]=saved[iRow];
     818      saved[iRow]=0.0;
     819    }
    572820  }
    573821  alternateWeights_->setNumElements(0);
     
    590838  delete [] weights_;
    591839  weights_=NULL;
     840  delete [] dubiousWeights_;
     841  dubiousWeights_=NULL;
    592842  delete infeasible_;
    593843  infeasible_ = NULL;
     
    598848  state_ =-1;
    599849}
     850// Returns true if would not find any row
     851bool
     852ClpDualRowSteepest::looksOptimal() const
     853{
     854  int iRow;
     855  const int * pivotVariable = model_->pivotVariable();
     856  double tolerance=model_->currentPrimalTolerance();
     857  // we can't really trust infeasibilities if there is primal error
     858  // this coding has to mimic coding in checkPrimalSolution
     859  double error = min(1.0e-3,model_->largestPrimalError());
     860  // allow tolerance at least slightly bigger than standard
     861  tolerance = tolerance +  error;
     862  // But cap
     863  tolerance = min(1000.0,tolerance);
     864  int numberRows = model_->numberRows();
     865  int numberInfeasible=0;
     866  for (iRow=0;iRow<numberRows;iRow++) {
     867    int iPivot=pivotVariable[iRow];
     868    double value = model_->solution(iPivot);
     869    double lower = model_->lower(iPivot);
     870    double upper = model_->upper(iPivot);
     871    if (value<lower-tolerance) {
     872      numberInfeasible++;
     873    } else if (value>upper+tolerance) {
     874      numberInfeasible++;
     875    }
     876  }
     877  return (numberInfeasible==0);
     878}
    600879
  • trunk/ClpFactorization.cpp

    r138 r225  
    44#include "CoinPragma.hpp"
    55#include "ClpFactorization.hpp"
     6#include "ClpQuadraticObjective.hpp"
    67#include "CoinHelperFunctions.hpp"
    78#include "CoinIndexedVector.hpp"
     
    7374}
    7475int
    75 ClpFactorization::factorize ( const ClpSimplex * model,
    76                               const ClpMatrixBase * matrix,
    77                               int numberRows, int numberColumns,
    78                               int rowIsBasic[], int columnIsBasic[] ,
    79                               double areaFactor )
    80 {
     76ClpFactorization::factorize ( ClpSimplex * model,
     77                              int solveType, bool valuesPass)
     78{
     79  const ClpMatrixBase * matrix = model->clpMatrix();
     80  int numberRows = model->numberRows();
     81  int numberColumns = model->numberColumns();
     82  // If too many compressions increase area
     83  if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
     84    areaFactor_ *= 1.1;
     85  }
    8186  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++;
     87    status_=-99;
     88    int * pivotVariable = model->pivotVariable();
     89    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
     90    while (status_<-98) {
     91     
     92      int i;
     93      int numberBasic=0;
     94      int numberRowBasic;
     95      for (i=0;i<numberRows;i++) {
     96        if (model->getRowStatus(i) == ClpSimplex::basic)
     97          pivotVariable[numberBasic++]=i;
    11398      }
    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++;
     99      numberRowBasic=numberBasic;
     100      for (i=0;i<numberColumns;i++) {
     101        if (model->getColumnStatus(i) == ClpSimplex::basic)
     102          pivotVariable[numberBasic++]=i;
    132103      }
    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++]];
     104      assert (numberBasic<=numberRows);
     105      // see if matrix a network
     106      ClpNetworkMatrix* networkMatrix =
     107        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
     108      // If network - still allow ordinary factorization first time for laziness
     109      int saveMaximumPivots = maximumPivots();
     110      delete networkBasis_;
     111      networkBasis_ = NULL;
     112      if (networkMatrix&&!doCheck)
     113        maximumPivots(1);
     114      while (status_==-99) {
     115        // maybe for speed will be better to leave as many regions as possible
     116        gutsOfDestructor();
     117        gutsOfInitialize(2);
     118        CoinBigIndex numberElements=numberRowBasic;
     119
     120        // compute how much in basis
     121
     122        int i;
     123
     124        numberElements +=matrix->fillBasis(model,
     125                                           pivotVariable+numberRowBasic,
     126                                           numberRowBasic,
     127                                           numberBasic-numberRowBasic,
     128                                           NULL,NULL,NULL);
     129
     130        numberElements = 3 * numberBasic + 3 * numberElements + 10000;
     131        getAreas ( numberRows, numberBasic, numberElements,
     132                   2 * numberElements );
     133        //fill
     134        //copy
     135        numberElements=numberRowBasic;
     136        for (i=0;i<numberRowBasic;i++) {
     137          int iRow = pivotVariable[i];
     138          indexRowU_[i]=iRow;
     139          indexColumnU_[i]=i;
     140          elementU_[i]=slackValue_;
     141        }
     142        numberElements +=matrix->fillBasis(model,
     143                                           pivotVariable+numberRowBasic,
     144                                           numberRowBasic,
     145                                           numberBasic-numberRowBasic,
     146                                           indexRowU_+numberElements,
     147                                           indexColumnU_+numberElements,
     148                                           elementU_+numberElements);
     149        lengthU_ = numberElements;
     150
     151        preProcess ( 0 );
     152        factor (  );
     153        if (status_==-99) {
     154          // get more memory
     155          areaFactor(2.0*areaFactor());
    149156        }
    150157      }
    151       for (i=0;i<numberColumns;i++) {
    152         if (columnIsBasic[i]>=0) {
    153           columnIsBasic[i]=permuteBack[back[numberBasic++]];
    154         }
    155       }
    156       if (increasingRows_>1) {
     158      // If we get here status is 0 or -1
     159      if (status_ == 0) {
     160        int * permuteBack = permuteBack_;
     161        int * back = pivotColumnBack_;
     162        int * pivotTemp = pivotColumn_;
     163        ClpDisjointCopyN ( pivotVariable, numberRows_ , pivotTemp  );
     164        // Redo pivot order
     165        for (i=0;i<numberRowBasic;i++) {
     166          int k = pivotTemp[i];
     167          // so rowIsBasic[k] would be permuteBack[back[i]]
     168          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
     169        }
     170        for (;i<numberRows;i++) {
     171          int k = pivotTemp[i];
     172          // so rowIsBasic[k] would be permuteBack[back[i]]
     173          pivotVariable[permuteBack[back[i]]]=k;
     174        }
    157175        // 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);
     176        // these arrays start off as copies of permute
     177        // (and we could use permute_ instead of pivotColumn (not back though))
     178        ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
     179        ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
     180        if (networkMatrix) {
     181          maximumPivots(saveMaximumPivots);
     182          // create network factorization
     183          if (doCheck)
     184            delete networkBasis_; // temp
     185          networkBasis_ = new ClpNetworkBasis(model,numberRows_,
     186                                              pivotRegion_,
     187                                              permuteBack_,
     188                                              startColumnU_,
     189                                              numberInColumn_,
     190                                              indexRowU_,
     191                                              elementU_);
     192          // kill off arrays in ordinary factorization
     193          if (!doCheck) {
     194            gutsOfDestructor();
     195#if 0
     196            // but put back permute arrays so odd things will work
     197            int numberRows = model->numberRows();
     198            pivotColumnBack_ = new int [numberRows];
     199            permute_ = new int [numberRows];
     200            int i;
     201            for (i=0;i<numberRows;i++) {
     202              pivotColumnBack_[i]=i;
     203              permute_[i]=i;
     204            }
     205#endif
     206          }
    222207        } 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;
     208          // See if worth going sparse and when
     209          if (numberFtranCounts_>100) {
     210            ftranAverageAfterL_ = max(ftranCountAfterL_/ftranCountInput_,1.0);
     211            ftranAverageAfterR_ = max(ftranCountAfterR_/ftranCountAfterL_,1.0);
     212            ftranAverageAfterU_ = max(ftranCountAfterU_/ftranCountAfterR_,1.0);
     213            assert (ftranCountInput_&&ftranCountAfterL_&&ftranCountAfterR_);
     214            if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
     215              btranAverageAfterU_ = max(btranCountAfterU_/btranCountInput_,1.0);
     216              btranAverageAfterR_ = max(btranCountAfterR_/btranCountAfterU_,1.0);
     217              btranAverageAfterL_ = max(btranCountAfterL_/btranCountAfterR_,1.0);
     218            } else {
     219              // we have not done any useful btrans (values pass?)
     220              btranAverageAfterU_ = 1.0;
     221              btranAverageAfterR_ = 1.0;
     222              btranAverageAfterL_ = 1.0;
     223            }
     224          }
     225          // scale back
     226         
     227          ftranCountInput_ *= 0.8;
     228          ftranCountAfterL_ *= 0.8;
     229          ftranCountAfterR_ *= 0.8;
     230          ftranCountAfterU_ *= 0.8;
     231          btranCountInput_ *= 0.8;
     232          btranCountAfterU_ *= 0.8;
     233          btranCountAfterR_ *= 0.8;
     234          btranCountAfterL_ *= 0.8;
     235        }
     236      } else if (status_ == -1&&(solveType==0||solveType==2)) {
     237        // This needs redoing as it was merged coding - does not need array
     238        int numberTotal = numberRows+numberColumns;
     239        int * isBasic = new int [numberTotal];
     240        int * rowIsBasic = isBasic+numberColumns;
     241        int * columnIsBasic = isBasic;
     242        for (i=0;i<numberTotal;i++)
     243          isBasic[i]=-1;
     244        for (i=0;i<numberRowBasic;i++) {
     245          int iRow = pivotVariable[i];
     246          rowIsBasic[iRow]=1;
     247        }
     248        for (;i<numberBasic;i++) {
     249          int iColumn = pivotVariable[i];
     250          columnIsBasic[iColumn]=1;
     251        }
     252        numberBasic=0;
     253        for (i=0;i<numberRows;i++)
     254          pivotVariable[i]=-1;
     255        // mark as basic or non basic
     256        for (i=0;i<numberRows;i++) {
     257          if (rowIsBasic[i]>=0) {
     258            if (pivotColumn_[numberBasic]>=0)
     259              rowIsBasic[i]=pivotColumn_[numberBasic];
     260            else
     261              rowIsBasic[i]=-1;
     262            numberBasic++;
     263          }
     264        }
     265        for (i=0;i<numberColumns;i++) {
     266          if (columnIsBasic[i]>=0) {
     267            if (pivotColumn_[numberBasic]>=0)
     268              columnIsBasic[i]=pivotColumn_[numberBasic];
     269            else
     270              columnIsBasic[i]=-1;
     271            numberBasic++;
     272          }
     273        }
     274        // leave pivotVariable in useful form for cleaning basis
     275        int * pivotVariable = model->pivotVariable();
     276        for (i=0;i<numberRows;i++) {
     277          pivotVariable[i]=-1;
     278        }
     279        for (i=0;i<numberRows;i++) {
     280          if (model->getRowStatus(i) == ClpSimplex::basic) {
     281            int iPivot = rowIsBasic[i];
     282            if (iPivot>=0)
     283              pivotVariable[iPivot]=i+numberColumns;
     284          }
     285        }
     286        for (i=0;i<numberColumns;i++) {
     287          if (model->getColumnStatus(i) == ClpSimplex::basic) {
     288            int iPivot = columnIsBasic[i];
     289            if (iPivot>=0)
     290              pivotVariable[iPivot]=i;
     291          }
     292        }
     293        delete [] isBasic;
     294        double * columnLower = model->lowerRegion();
     295        double * columnUpper = model->upperRegion();
     296        double * columnActivity = model->solutionRegion();
     297        double * rowLower = model->lowerRegion(0);
     298        double * rowUpper = model->upperRegion(0);
     299        double * rowActivity = model->solutionRegion(0);
     300        //redo basis - first take ALL columns out
     301        int iColumn;
     302        double largeValue = model->largeValue();
     303        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     304          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
     305            // take out
     306            if (!valuesPass) {
     307              double lower=columnLower[iColumn];
     308              double upper=columnUpper[iColumn];
     309              double value=columnActivity[iColumn];
     310              if (lower>-largeValue||upper<largeValue) {
     311                if (fabs(value-lower)<fabs(value-upper)) {
     312                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
     313                  columnActivity[iColumn]=lower;
     314                } else {
     315                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
     316                  columnActivity[iColumn]=upper;
     317                }
     318              } else {
     319                model->setColumnStatus(iColumn,ClpSimplex::isFree);
     320              }
     321            } else {
     322              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
     323            }
     324          }
     325        }
     326        int iRow;
     327        for (iRow=0;iRow<numberRows;iRow++) {
     328          int iSequence=pivotVariable[iRow];
     329          if (iSequence>=0) {
     330            // basic
     331            if (iSequence>=numberColumns) {
     332              // slack in - leave
     333              //assert (iSequence-numberColumns==iRow);
     334            } else {
     335              // put back structural
     336              model->setColumnStatus(iSequence,ClpSimplex::basic);
     337            }
     338          } else {
     339            // put in slack
     340            model->setRowStatus(iRow,ClpSimplex::basic);
     341          }
     342        }
     343        // signal repeat
     344        status_=-99;
     345        // set fixed if they are
     346        for (iRow=0;iRow<numberRows;iRow++) {
     347          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
     348            if (rowLower[iRow]==rowUpper[iRow]) {
     349              rowActivity[iRow]=rowLower[iRow];
     350              model->setRowStatus(iRow,ClpSimplex::isFixed);
     351            }
     352          }
     353        }
     354        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     355          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
     356            if (columnLower[iColumn]==columnUpper[iColumn]) {
     357              columnActivity[iColumn]=columnLower[iColumn];
     358              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
     359            }
     360          }
     361        }
     362      }
    239363    }
    240364  } else {
     
    243367  }
    244368
     369  if (!status_) {
     370    // take out part if quadratic
     371    if (model->algorithm()==2) {
     372      ClpObjective * obj = model->objectiveAsObject();
     373      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
     374      assert (quadraticObj);
     375      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     376      int numberXColumns = quadratic->getNumCols();
     377      assert (numberXColumns<numberColumns);
     378      int base = numberColumns-numberXColumns;
     379      int * which = new int [numberXColumns];
     380      int * pivotVariable = model->pivotVariable();
     381      int * permute = pivotColumn();
     382      int i;
     383      int n=0;
     384      for (i=0;i<numberRows;i++) {
     385        int iSj = pivotVariable[i]-base;
     386        if (iSj>=0&&iSj<numberXColumns)
     387          which[n++]=permute[i];
     388      }
     389      if (n)
     390        emptyRows(n,which);
     391      delete [] which;
     392    }
     393  }
    245394  return status_;
    246395}
     
    285434   region1 starts as zero and is zero at end */
    286435int
    287 ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
    288                                  CoinIndexedVector * regionSparse2,
    289                                  bool FTUpdate)
     436ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
     437                                   CoinIndexedVector * regionSparse2)
    290438{
    291439#ifdef CLP_DEBUG
     
    294442  if (!networkBasis_) {
    295443    collectStatistics_ = true;
    296     return CoinFactorization::updateColumn(regionSparse,
    297                                            regionSparse2,
    298                                            FTUpdate);
     444    int returnValue= CoinFactorization::updateColumnFT(regionSparse,
     445                                                       regionSparse2);
    299446    collectStatistics_ = false;
     447    return returnValue;
    300448  } else {
    301449#ifdef CHECK_NETWORK
    302450    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    303     int returnCode = CoinFactorization::updateColumn(regionSparse,
    304                                                      regionSparse2, FTUpdate);
     451    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
     452                                                       regionSparse2);
    305453    networkBasis_->updateColumn(regionSparse,save);
    306454    int i;
     
    315463    return returnCode;
    316464#else
    317     return networkBasis_->updateColumn(regionSparse,regionSparse2);
    318 #endif
    319   }
    320 }
    321 /* Updates one column (FTRAN) to/from array
     465    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
     466    return 1;
     467#endif
     468  }
     469}
     470/* Updates one column (FTRAN) from region2
    322471   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 */
     472   region1 starts as zero and is zero at end */
    328473int
    329474ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
    330                         double array[] ) const
    331 {
     475                                 CoinIndexedVector * regionSparse2,
     476                                 bool noPermute) const
     477{
     478#ifdef CLP_DEBUG
     479  if (!noPermute)
     480    regionSparse->checkClear();
     481#endif
    332482  if (!networkBasis_) {
    333     return CoinFactorization::updateColumn(regionSparse,
    334                                            array);
     483    collectStatistics_ = true;
     484    int returnValue= CoinFactorization::updateColumn(regionSparse,
     485                                                     regionSparse2,
     486                                                     noPermute);
     487    collectStatistics_ = false;
     488    return returnValue;
    335489  } else {
    336490#ifdef CHECK_NETWORK
    337     double * save = new double [numberRows_+1];
    338     memcpy(save,array,(numberRows_+1)*sizeof(double));
     491    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
    339492    int returnCode = CoinFactorization::updateColumn(regionSparse,
    340                                                      array);
    341     networkBasis_->updateColumn(regionSparse, save);
     493                                                     regionSparse2,
     494                                                     noPermute);
     495    networkBasis_->updateColumn(regionSparse,save);
    342496    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;
     497    double * array = regionSparse2->denseVector();
     498    double * array2 = save->denseVector();
     499    for (i=0;i<numberRows_;i++) {
     500      double value1 = array[i];
     501      double value2 = array2[i];
     502      assert (value1==value2);
     503    }
     504    delete save;
    346505    return returnCode;
    347506#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);
     507    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
     508    return 1;
    379509#endif
    380510  }
  • trunk/ClpLinearObjective.cpp

    r124 r225  
    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//-------------------------------------------------------------------
     
    83110// Returns gradient
    84111double * 
    85 ClpLinearObjective::gradient(const double * solution, double & offset)
     112ClpLinearObjective::gradient(const double * solution, double & offset,bool refresh)
    86113{
    87114  offset=0.0;
     
    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
  • trunk/ClpMatrixBase.cpp

    r63 r225  
    77
    88#include "ClpMatrixBase.hpp"
     9#include "ClpSimplex.hpp"
    910
    1011//#############################################################################
     
    5657                     const double * columnScale) const
    5758{
    58   std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
    59   abort();
     59  if (rowScale) {
     60    std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
     61    abort();
     62  } else {
     63    times(scalar,x,y);
     64  }
    6065}
    6166// And for scaling - default aborts for when scaling not supported
     
    6671                                const double * columnScale) const
    6772{
    68   std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
     73  if (rowScale) {
     74    std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
     75    abort();
     76  } else {
     77    transposeTimes(scalar,x,y);
     78  }
     79}
     80/* Subset clone (without gaps).  Duplicates are allowed
     81   and order is as given.
     82   Derived classes need not provide this as it may not always make
     83   sense */
     84ClpMatrixBase *
     85ClpMatrixBase::subsetClone (
     86                            int numberRows, const int * whichRows,
     87                            int numberColumns, const int * whichColumns) const
     88 
     89
     90{
     91  std::cerr<<"subsetClone not supported - ClpMatrixBase"<<std::endl;
    6992  abort();
    7093}
    71 
     94/* Given positive integer weights for each row fills in sum of weights
     95   for each column (and slack).
     96   Returns weights vector
     97   Default returns vector of ones
     98*/
     99CoinBigIndex *
     100ClpMatrixBase::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
     101{
     102  int number = model->numberRows()+model->numberColumns();
     103  CoinBigIndex * weights = new CoinBigIndex[number];
     104  int i;
     105  for (i=0;i<number;i++)
     106    weights[i]=1;
     107  return weights;
     108}
     109// Append Columns
     110void
     111ClpMatrixBase::appendCols(int number, const CoinPackedVectorBase * const * columns)
     112{
     113  std::cerr<<"appendCols not supported - ClpMatrixBase"<<std::endl;
     114  abort();
     115}
     116// Append Rows
     117void
     118ClpMatrixBase::appendRows(int number, const CoinPackedVectorBase * const * rows)
     119{
     120  std::cerr<<"appendRows not supported - ClpMatrixBase"<<std::endl;
     121  abort();
     122}
     123/* Returns largest and smallest elements of both signs.
     124   Largest refers to largest absolute value.
     125*/
     126void
     127ClpMatrixBase::rangeOfElements(double & smallestNegative, double & largestNegative,
     128                       double & smallestPositive, double & largestPositive)
     129{
     130  smallestNegative=0.0;
     131  largestNegative=0.0;
     132  smallestPositive=0.0;
     133  largestPositive=0.0;
     134}
  • trunk/ClpMessage.cpp

    r161 r225  
    2222  {CLP_SIMPLEX_ACCURACY,6,3,"Primal error %g, dual error %g"},
    2323  {CLP_SIMPLEX_BADFACTOR,7,1,"Singular factorization of basis - status %d"},
    24   {CLP_SIMPLEX_BOUNDTIGHTEN,8,1,"Bounds were tightened %d times"},
     24  {CLP_SIMPLEX_BOUNDTIGHTEN,8,3,"Bounds were tightened %d times"},
    2525  {CLP_SIMPLEX_INFEASIBILITIES,9,1,"%d infeasibilities"},
    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"},
    30   {CLP_SIMPLEX_PERTURB,14,1,"Perturbing problem by %g"},
     30  {CLP_SIMPLEX_PERTURB,14,1,"Perturbing problem by %g %% of %g - largest nonzero change %g (%% %g) - largest zero change %g"},
    3131  {CLP_PRIMAL_ORIGINAL,15,2,"Going back to original tolerance"},
    3232  {CLP_PRIMAL_WEIGHT,16,2,"New infeasibility weight of %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"},
    6250  {CLP_PACKEDSCALE_FINAL,1003,2,"Final range of elements is %g to %g"},
    6351  {CLP_PACKEDSCALE_FORGET,1004,2,"Not bothering to scale as good enough"},
    64   {CLP_INITIALIZE_STEEP,1005,1,"Initializing steepest edge weights - old %g, new %g"},
     52  {CLP_INITIALIZE_STEEP,1005,3,"Initializing steepest edge weights - old %g, new %g"},
    6553  {CLP_UNABLE_OPEN,6001,0,"Unable to open file %s for reading"},
    6654  {CLP_BAD_BOUNDS,6002,0,"%d bad bound pairs or bad objectives were found - first at %c%d"},
     
    7361  {CLP_CRASH,28,1,"Crash put %d variables in basis, %d dual infeasibilities"},
    7462  {CLP_END_VALUES_PASS,29,1,"End of values pass after %d iterations"},
     63  {CLP_QUADRATIC_BOTH,108,32,"%s %d (%g) and %d (%g) both basic"},
     64  {CLP_QUADRATIC_PRIMAL_DETAILS,109,32,"coeff %g, %g, %g - dj %g - deriv zero at %g, sj at %g"},
     65  {CLP_IDIOT_ITERATION,30,1,"%d infeas %g, obj %g - mu %g, its %d, %d interior"},
     66  {CLP_INFEASIBLE,3003,1,"Analysis indicates model infeasible or unbounded"},
     67  {CLP_MATRIX_CHANGE,31,2,"Matrix can not be converted into %s"},
     68  {CLP_TIMING,32,1,"%s objective %.10g - %d iterations time %.2f2%?, Presolve %.2f%?, Idiot %.2f%?"},
     69  {CLP_INTERVAL_TIMING,33,2,"%s took %.2f seconds (total %.2f)"},
     70  {CLP_SPRINT,34,1,"Pass %d took %d iterations, objective %g, dual infeasibilities %g( %d)"},
    7571  {CLP_DUMMY_END,999999,0,""}
    7672};
  • trunk/ClpModel.cpp

    r169 r225  
    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"
     15#include "CoinTime.hpp"
    2116#include "ClpModel.hpp"
    2217#include "ClpPackedMatrix.hpp"
     
    2621#include "ClpMessage.hpp"
    2722#include "ClpLinearObjective.hpp"
    28 
    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 }
     23#include "ClpQuadraticObjective.hpp"
     24
    4625//#############################################################################
    4726
     
    6544  matrix_(NULL),
    6645  rowCopy_(NULL),
    67   quadraticObjective_(NULL),
    6846  ray_(NULL),
    6947  status_(NULL),
     
    7250  solveType_(0),
    7351  problemStatus_(-1),
     52  secondaryStatus_(0),
    7453  lengthNames_(0),
    7554  defaultHandler_(true),
     
    8968  strParam_[ClpProbName] = "ClpDefaultName";
    9069  handler_ = new CoinMessageHandler();
    91   handler_->setLogLevel(2);
     70  handler_->setLogLevel(1);
    9271  messages_ = ClpMessage();
    9372}
     
    129108  delete rowCopy_;
    130109  rowCopy_=NULL;
    131   delete quadraticObjective_;
    132   quadraticObjective_ = NULL;
    133110  delete [] ray_;
    134111  ray_ = NULL;
     
    149126    dblParam_[ClpDualTolerance]=value;
    150127}
    151 void ClpModel::setOptimizationDirection( int value)
    152 {
    153   if (value>=-1&&value<=1)
    154     optimizationDirection_=value;
     128void ClpModel::setOptimizationDirection( double value)
     129{
     130  optimizationDirection_=value;
    155131}
    156132void
     
    339315  strParam_[ClpProbName] = rhs.strParam_[ClpProbName];
    340316
     317  optimizationDirection_ = rhs.optimizationDirection_;
    341318  objectiveValue_=rhs.objectiveValue_;
    342319  smallElement_ = rhs.smallElement_;
     
    344321  solveType_ = rhs.solveType_;
    345322  problemStatus_ = rhs.problemStatus_;
     323  secondaryStatus_ = rhs.secondaryStatus_;
     324  numberRows_ = rhs.numberRows_;
     325  numberColumns_ = rhs.numberColumns_;
    346326  if (trueCopy) {
    347327    lengthNames_ = rhs.lengthNames_;
     
    384364    status_ = ClpCopyOfArray( rhs.status_,numberColumns_+numberRows_);
    385365    ray_ = NULL;
    386     if (problemStatus_==1)
     366    if (problemStatus_==1&&!secondaryStatus_)
    387367      ray_ = ClpCopyOfArray (rhs.ray_,numberRows_);
    388368    else if (problemStatus_==2)
     
    392372    } else {
    393373      rowCopy_=NULL;
    394     }
    395     if (rhs.quadraticObjective_) {
    396       quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
    397     } else {
    398       quadraticObjective_=NULL;
    399374    }
    400375    matrix_=NULL;
     
    415390    matrix_ = rhs.matrix_;
    416391    rowCopy_ = NULL;
    417     quadraticObjective_ = rhs.quadraticObjective_;
    418392    ray_ = rhs.ray_;
    419393    lengthNames_ = 0;
     
    449423  otherModel.numberIterations_ = numberIterations_;
    450424  otherModel.problemStatus_ = problemStatus_;
     425  otherModel.secondaryStatus_ = secondaryStatus_;
    451426  rowActivity_ = NULL;
    452427  columnActivity_ = NULL;
     
    461436  matrix_ = NULL;
    462437  rowCopy_ = NULL;
    463   quadraticObjective_=NULL,
    464438  delete [] otherModel.ray_;
    465439  otherModel.ray_ = ray_;
     
    527501  case ClpMaxSeconds:
    528502    if(value>=0)
    529       value += cpuTime();
     503      value += CoinCpuTime();
    530504    else
    531505      value = -1.0;
     
    678652    // set state back to unknown
    679653    problemStatus_ = -1;
     654    secondaryStatus_ = 0;
    680655    delete [] ray_;
    681656    ray_ = NULL;
     
    697672      which[i-newNumberColumns]=i;
    698673    matrix_->deleteCols(numberColumns_-newNumberColumns,which);
    699     if (quadraticObjective_) {
    700       quadraticObjective_->deleteCols(numberColumns_-newNumberColumns,which);
    701       quadraticObjective_->deleteRows(numberColumns_-newNumberColumns,which);
    702     }
    703674    delete [] which;
    704675  }
     
    733704                              number, which, newSize);
    734705  matrix_->deleteRows(number,which);
     706  //matrix_->removeGaps();
    735707  // status
    736708  if (status_) {
     
    748720  // set state back to unknown
    749721  problemStatus_ = -1;
     722  secondaryStatus_ = 0;
    750723  delete [] ray_;
    751724  ray_ = NULL;
     
    770743                              number, which, newSize);
    771744  matrix_->deleteCols(number,which);
    772   if (quadraticObjective_) {
    773     quadraticObjective_->deleteCols(number,which);
    774     quadraticObjective_->deleteRows(number,which);
    775   }
     745  //matrix_->removeGaps();
    776746  // status
    777747  if (status_) {
     
    787757    status_ = temp;
    788758  }
     759  integerType_ = deleteChar(integerType_,numberColumns_,
     760                            number, which, newSize,true);
    789761  numberColumns_=newSize;
    790762  // set state back to unknown
    791763  problemStatus_ = -1;
     764  secondaryStatus_ = 0;
    792765  delete [] ray_;
    793766  ray_ = NULL;
     
    796769  rowNames_ = std::vector<std::string> ();
    797770  columnNames_ = std::vector<std::string> ();
    798   integerType_ = deleteChar(integerType_,numberColumns_,
    799                             number, which, newSize,true);
    800771}
    801772// Add rows
     
    815786      rows[iRow] =
    816787        new CoinPackedVector(rowStarts[iRow+1]-iStart,
     788                             columns+iStart,elements+iStart);
     789    }
     790    addRows(number, rowLower, rowUpper,
     791            rows);
     792    for (iRow=0;iRow<number;iRow++)
     793      delete rows[iRow];
     794    delete [] rows;
     795  }
     796}
     797// Add rows
     798void
     799ClpModel::addRows(int number, const double * rowLower,
     800                  const double * rowUpper,
     801                  const int * rowStarts,
     802                  const int * rowLengths, const int * columns,
     803                  const double * elements)
     804{
     805  // Create a list of CoinPackedVectors
     806  if (number) {
     807    CoinPackedVectorBase ** rows=
     808      new CoinPackedVectorBase * [number];
     809    int iRow;
     810    for (iRow=0;iRow<number;iRow++) {
     811      int iStart = rowStarts[iRow];
     812      rows[iRow] =
     813        new CoinPackedVector(rowLengths[iRow],
    817814                             columns+iStart,elements+iStart);
    818815    }
     
    866863  if (!matrix_)
    867864    createEmptyMatrix();
    868   // Use matrix() to get round virtual problem
    869   matrix()->appendRows(number,rows);
     865  matrix_->appendRows(number,rows);
    870866}
    871867// Add columns
     
    886882      columns[iColumn] =
    887883        new CoinPackedVector(columnStarts[iColumn+1]-iStart,
     884                             rows+iStart,elements+iStart);
     885    }
     886    addColumns(number, columnLower, columnUpper,
     887               objIn, columns);
     888    for (iColumn=0;iColumn<number;iColumn++)
     889      delete columns[iColumn];
     890    delete [] columns;
     891
     892  }
     893}
     894// Add columns
     895void
     896ClpModel::addColumns(int number, const double * columnLower,
     897                     const double * columnUpper,
     898                     const double * objIn,
     899                     const int * columnStarts,
     900                     const int * columnLengths, const int * rows,
     901                     const double * elements)
     902{
     903  // Create a list of CoinPackedVectors
     904  if (number) {
     905    CoinPackedVectorBase ** columns=
     906      new CoinPackedVectorBase * [number];
     907    int iColumn;
     908    for (iColumn=0;iColumn<number;iColumn++) {
     909      int iStart = columnStarts[iColumn];
     910      columns[iColumn] =
     911        new CoinPackedVector(columnLengths[iColumn],
    888912                             rows+iStart,elements+iStart);
    889913    }
     
    949973  if (!matrix_)
    950974    createEmptyMatrix();
    951   // Use matrix() to get round virtual problem
    952   matrix()->appendCols(number,columns);
     975  matrix_->appendCols(number,columns);
    953976}
    954977// Infeasibility/unbounded ray (NULL returned if none/wrong)
     
    957980{
    958981  double * array = NULL;
    959   if (problemStatus_==1)
     982  if (problemStatus_==1&&!secondaryStatus_)
    960983    array = ClpCopyOfArray(ray_,numberRows_);
    961984  return array;
     
    9791002{
    9801003  if(value>=0)
    981     dblParam_[ClpMaxSeconds]=value+cpuTime();
     1004    dblParam_[ClpMaxSeconds]=value+CoinCpuTime();
    9821005  else
    9831006    dblParam_[ClpMaxSeconds]=-1.0;
     
    9891012  bool hitMax= (numberIterations_>=maximumIterations());
    9901013  if (dblParam_[ClpMaxSeconds]>=0.0&&!hitMax)
    991     hitMax = (cpuTime()>=dblParam_[ClpMaxSeconds]);
     1014    hitMax = (CoinCpuTime()>=dblParam_[ClpMaxSeconds]);
    9921015  return hitMax;
    9931016}
     
    10301053  }
    10311054  CoinMpsIO m;
    1032   double time1 = cpuTime(),time2;
     1055  bool savePrefix =m.messageHandler()->prefix();
     1056  m.messageHandler()->setPrefix(handler_->prefix());
     1057  double time1 = CoinCpuTime(),time2;
    10331058  int status=m.readMps(fileName,"");
     1059  m.messageHandler()->setPrefix(savePrefix);
    10341060  if (!status||ignoreErrors) {
    10351061    loadProblem(*m.getMatrixByCol(),
     
    10491075      unsigned int maxLength=0;
    10501076      int iRow;
     1077      rowNames_ = std::vector<std::string> ();
     1078      columnNames_ = std::vector<std::string> ();
    10511079      rowNames_.reserve(numberRows_);
    10521080      for (iRow=0;iRow<numberRows_;iRow++) {
     
    10681096    }
    10691097    setDblParam(ClpObjOffset,m.objectiveOffset());
    1070     time2 = cpuTime();
     1098    time2 = CoinCpuTime();
    10711099    handler_->message(CLP_IMPORT_RESULT,messages_)
    10721100      <<fileName
     
    10901118   
    10911119  const double obj = objectiveValue();
    1092   const int maxmin = optimizationDirection();
     1120  const double maxmin = optimizationDirection();
    10931121
    10941122  if (problemStatus_ == 0) // optimal
     
    11111139   
    11121140  const double obj = objectiveValue();
    1113   const int maxmin = optimizationDirection();
     1141  const double maxmin = optimizationDirection();
    11141142
    11151143  if (problemStatus_ == 0) // optimal
     
    11731201{
    11741202  assert (numberColumns==numberColumns_);
    1175   quadraticObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
    1176                                              start[numberColumns],element,column,start,NULL);
     1203  assert ((dynamic_cast< ClpLinearObjective*>(objective_)));
     1204  double offset;
     1205  ClpObjective * obj = new ClpQuadraticObjective(objective_->gradient(NULL,offset,false),numberColumns,
     1206                                             start,column,element);
     1207  delete objective_;
     1208  objective_ = obj;
     1209
    11771210}
    11781211void
     
    11801213{
    11811214  assert (matrix.getNumCols()==numberColumns_);
    1182   quadraticObjective_ = new CoinPackedMatrix(matrix);
     1215  assert ((dynamic_cast< ClpLinearObjective*>(objective_)));
     1216  double offset;
     1217  ClpQuadraticObjective * obj =
     1218    new ClpQuadraticObjective(objective_->gradient(NULL,offset,false),numberColumns_,
     1219                                                 NULL,NULL,NULL);
     1220  delete objective_;
     1221  objective_ = obj;
     1222  obj->loadQuadraticObjective(matrix);
    11831223}
    11841224// Get rid of quadratic objective
     
    11861226ClpModel::deleteQuadraticObjective()
    11871227{
    1188   delete quadraticObjective_;
    1189   quadraticObjective_ = NULL;
     1228  ClpQuadraticObjective * obj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     1229  if (obj)
     1230    obj->deleteQuadraticObjective();
     1231}
     1232void
     1233ClpModel::setObjective(ClpObjective * objective)
     1234{
     1235  delete objective_;
     1236  objective_=objective->clone();
     1237}
     1238// Returns resized array and updates size
     1239double * whichDouble(double * array , int number, const int * which)
     1240{
     1241  double * newArray=NULL;
     1242  if (array&&number) {
     1243    int i ;
     1244    newArray = new double[number];
     1245    for (i=0;i<number;i++)
     1246      newArray[i]=array[which[i]];
     1247  }
     1248  return newArray;
     1249}
     1250char * whichChar(char * array , int number, const int * which)
     1251{
     1252  char * newArray=NULL;
     1253  if (array&&number) {
     1254    int i ;
     1255    newArray = new char[number];
     1256    for (i=0;i<number;i++)
     1257      newArray[i]=array[which[i]];
     1258  }
     1259  return newArray;
     1260}
     1261unsigned char * whichUnsignedChar(unsigned char * array ,
     1262                                  int number, const int * which)
     1263{
     1264  unsigned char * newArray=NULL;
     1265  if (array&&number) {
     1266    int i ;
     1267    newArray = new unsigned char[number];
     1268    for (i=0;i<number;i++)
     1269      newArray[i]=array[which[i]];
     1270  }
     1271  return newArray;
     1272}
     1273// Replace Clp Matrix (current is not deleted)
     1274void
     1275ClpModel::replaceMatrix( ClpMatrixBase * matrix)
     1276{
     1277  matrix_=matrix;
     1278}
     1279// Subproblem constructor
     1280ClpModel::ClpModel ( const ClpModel * rhs,
     1281                     int numberRows, const int * whichRow,
     1282                     int numberColumns, const int * whichColumn,
     1283                     bool dropNames, bool dropIntegers)
     1284{
     1285#if 0
     1286  // Could be recoded to be faster and take less memory
     1287  // and to allow re-ordering and duplicates etc
     1288  gutsOfCopy(*rhs,true);
     1289  int numberRowsWhole = rhs->numberRows();
     1290  int * delRow = new int[numberRowsWhole];
     1291  memset(delRow,0,numberRowsWhole*sizeof(int));
     1292  int i;
     1293  for (i=0;i<numberRows;i++) {
     1294    int iRow=whichRow[i];
     1295    assert (iRow>=0&&iRow<numberRowsWhole);
     1296    delRow[iRow]=1;
     1297  }
     1298  numberRows=0;
     1299  for (i=0;i<numberRowsWhole;i++) {
     1300    if (delRow[i]==0)
     1301      delRow[numberRows++]=i;
     1302  }
     1303  deleteRows(numberRows,delRow);
     1304  delete [] delRow;
     1305  int numberColumnsWhole = rhs->numberColumns();
     1306  int * delColumn = new int[numberColumnsWhole];
     1307  memset(delColumn,0,numberColumnsWhole*sizeof(int));
     1308  for (i=0;i<numberColumns;i++) {
     1309    int iColumn=whichColumn[i];
     1310    assert (iColumn>=0&&iColumn<numberColumnsWhole);
     1311    delColumn[iColumn]=1;
     1312  }
     1313  numberColumns=0;
     1314  for (i=0;i<numberColumnsWhole;i++) {
     1315    if (delColumn[i]==0)
     1316      delColumn[numberColumns++]=i;
     1317  }
     1318  deleteColumns(numberColumns,delColumn);
     1319  delete [] delColumn;
     1320#else
     1321  defaultHandler_ = rhs->defaultHandler_;
     1322  if (defaultHandler_)
     1323    handler_ = new CoinMessageHandler(*rhs->handler_);
     1324   else
     1325    handler_ = rhs->handler_;
     1326  messages_ = rhs->messages_;
     1327  intParam_[ClpMaxNumIteration] = rhs->intParam_[ClpMaxNumIteration];
     1328  intParam_[ClpMaxNumIterationHotStart] =
     1329    rhs->intParam_[ClpMaxNumIterationHotStart];
     1330
     1331  dblParam_[ClpDualObjectiveLimit] = rhs->dblParam_[ClpDualObjectiveLimit];
     1332  dblParam_[ClpPrimalObjectiveLimit] = rhs->dblParam_[ClpPrimalObjectiveLimit];
     1333  dblParam_[ClpDualTolerance] = rhs->dblParam_[ClpDualTolerance];
     1334  dblParam_[ClpPrimalTolerance] = rhs->dblParam_[ClpPrimalTolerance];
     1335  dblParam_[ClpObjOffset] = rhs->dblParam_[ClpObjOffset];
     1336  dblParam_[ClpMaxSeconds] = rhs->dblParam_[ClpMaxSeconds];
     1337
     1338  strParam_[ClpProbName] = rhs->strParam_[ClpProbName];
     1339
     1340  optimizationDirection_ = rhs->optimizationDirection_;
     1341  objectiveValue_=rhs->objectiveValue_;
     1342  smallElement_ = rhs->smallElement_;
     1343  numberIterations_ = rhs->numberIterations_;
     1344  solveType_ = rhs->solveType_;
     1345  problemStatus_ = rhs->problemStatus_;
     1346  secondaryStatus_ = rhs->secondaryStatus_;
     1347  // check valid lists
     1348  int numberBad=0;
     1349  int i;
     1350  for (i=0;i<numberRows;i++)
     1351    if (whichRow[i]<0||whichRow[i]>=rhs->numberRows_)
     1352      numberBad++;
     1353  if (numberBad)
     1354    throw CoinError("bad row list", "subproblem constructor", "ClpModel");
     1355  numberBad=0;
     1356  for (i=0;i<numberColumns;i++)
     1357    if (whichColumn[i]<0||whichColumn[i]>=rhs->numberColumns_)
     1358      numberBad++;
     1359  if (numberBad)
     1360    throw CoinError("bad column list", "subproblem constructor", "ClpModel");
     1361  numberRows_ = numberRows;
     1362  numberColumns_ = numberColumns;
     1363  if (!dropNames) {
     1364    unsigned int maxLength=0;
     1365    int iRow;
     1366    rowNames_ = std::vector<std::string> ();
     1367    columnNames_ = std::vector<std::string> ();
     1368    rowNames_.reserve(numberRows_);
     1369    for (iRow=0;iRow<numberRows_;iRow++) {
     1370      rowNames_[iRow] = rhs->rowNames_[whichRow[iRow]];
     1371      maxLength = max(maxLength,(unsigned int) strlen(rowNames_[iRow].c_str()));
     1372    }
     1373    int iColumn;
     1374    columnNames_.reserve(numberColumns_);
     1375    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1376      columnNames_[iColumn] = rhs->columnNames_[whichColumn[iColumn]];
     1377      maxLength = max(maxLength,(unsigned int) strlen(columnNames_[iColumn].c_str()));
     1378    }
     1379    lengthNames_=(int) maxLength;
     1380  } else {
     1381    lengthNames_ = 0;
     1382    rowNames_ = std::vector<std::string> ();
     1383    columnNames_ = std::vector<std::string> ();
     1384  }
     1385  if (rhs->integerType_&&!dropIntegers) {
     1386    integerType_ = whichChar(rhs->integerType_,numberColumns,whichColumn);
     1387  } else {
     1388    integerType_ = NULL;
     1389  }
     1390  if (rhs->rowActivity_) {
     1391    rowActivity_=whichDouble(rhs->rowActivity_,numberRows,whichRow);
     1392    dual_=whichDouble(rhs->dual_,numberRows,whichRow);
     1393    columnActivity_=whichDouble(rhs->columnActivity_,numberColumns,
     1394                                whichColumn);
     1395    reducedCost_=whichDouble(rhs->reducedCost_,numberColumns,
     1396                                whichColumn);
     1397  } else {
     1398    rowActivity_=NULL;
     1399    columnActivity_=NULL;
     1400    dual_=NULL;
     1401    reducedCost_=NULL;
     1402  }
     1403  rowLower_=whichDouble(rhs->rowLower_,numberRows,whichRow);
     1404  rowUpper_=whichDouble(rhs->rowUpper_,numberRows,whichRow);
     1405  columnLower_=whichDouble(rhs->columnLower_,numberColumns,whichColumn);
     1406  columnUpper_=whichDouble(rhs->columnUpper_,numberColumns,whichColumn);
     1407  if (rhs->objective_)
     1408    objective_  = rhs->objective_->subsetClone(numberColumns,whichColumn);
     1409  else
     1410    objective_ = NULL;
     1411  rowObjective_=whichDouble(rhs->rowObjective_,numberRows,whichRow);
     1412  // status has to be done in two stages
     1413  status_ = new unsigned char[numberColumns_+numberRows_];
     1414  unsigned char * rowStatus = whichUnsignedChar(rhs->status_+rhs->numberColumns_,
     1415                                                numberRows_,whichRow);
     1416  unsigned char * columnStatus = whichUnsignedChar(rhs->status_,
     1417                                                numberColumns_,whichColumn);
     1418  memcpy(status_+numberColumns_,rowStatus,numberRows_);
     1419  delete [] rowStatus;
     1420  memcpy(status_,columnStatus,numberColumns_);
     1421  delete [] columnStatus;
     1422  ray_ = NULL;
     1423  if (problemStatus_==1&&!secondaryStatus_)
     1424    ray_ = whichDouble (rhs->ray_,numberRows,whichRow);
     1425  else if (problemStatus_==2)
     1426    ray_ = whichDouble (rhs->ray_,numberColumns,whichColumn);
     1427  if (rhs->rowCopy_) {
     1428    rowCopy_ = rhs->rowCopy_->subsetClone(numberRows,whichRow,
     1429                                          numberColumns,whichColumn);
     1430  } else {
     1431    rowCopy_=NULL;
     1432  }
     1433  matrix_=NULL;
     1434  if (rhs->matrix_) {
     1435    matrix_ = rhs->matrix_->subsetClone(numberRows,whichRow,
     1436                                        numberColumns,whichColumn);
     1437  }
     1438#endif
     1439}
     1440// Copies in names
     1441void
     1442ClpModel::copyNames(std::vector<std::string> & rowNames,
     1443                 std::vector<std::string> & columnNames)
     1444{
     1445  unsigned int maxLength=0;
     1446  int iRow;
     1447  rowNames_ = std::vector<std::string> ();
     1448  columnNames_ = std::vector<std::string> ();
     1449  rowNames_.reserve(numberRows_);
     1450  for (iRow=0;iRow<numberRows_;iRow++) {
     1451    rowNames_[iRow] = rowNames[iRow];
     1452    maxLength = max(maxLength,(unsigned int) strlen(rowNames_[iRow].c_str()));
     1453  }
     1454  int iColumn;
     1455  columnNames_.reserve(numberColumns_);
     1456  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1457    columnNames_[iColumn] = columnNames[iColumn];
     1458    maxLength = max(maxLength,(unsigned int) strlen(columnNames_[iColumn].c_str()));
     1459  }
     1460  lengthNames_=(int) maxLength;
    11901461}
    11911462//#############################################################################
  • trunk/ClpNetworkBasis.cpp

    r131 r225  
    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[i];
     685        region2[i]=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  }
  • trunk/ClpNetworkMatrix.cpp

    r124 r225  
    1515#include "ClpPlusMinusOneMatrix.hpp"
    1616#include "ClpMessage.hpp"
     17#include <iostream>
    1718
    1819//#############################################################################
     
    369370  ClpPlusMinusOneMatrix* rowCopy =
    370371    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    371   if (numberInRowArray>0.3*numberRows||!rowCopy) {
     372  bool packed = rowArray->packedMode();
     373  double factor = 0.3;
     374  // We may not want to do by row if there may be cache problems
     375  int numberColumns = model->numberColumns();
     376  // It would be nice to find L2 cache size - for moment 512K
     377  // Be slightly optimistic
     378  if (numberColumns*sizeof(double)>1000000) {
     379    if (numberRows*10<numberColumns)
     380      factor=0.1;
     381    else if (numberRows*4<numberColumns)
     382      factor=0.15;
     383    else if (numberRows*2<numberColumns)
     384      factor=0.2;
     385    //if (model->numberIterations()%50==0)
     386    //printf("%d nonzero\n",numberInRowArray);
     387  }
     388  if (numberInRowArray>factor*numberRows||!rowCopy) {
    372389    // do by column
    373390    int iColumn;
    374     double * markVector = y->denseVector(); // probably empty
     391    assert (!y->getNumElements());
    375392    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         }
     393    if (packed) {
     394      // need to expand pi into y
     395      assert(y->capacity()>=numberRows);
     396      double * piOld = pi;
     397      pi = y->denseVector();
     398      const int * whichRow = rowArray->getIndices();
     399      int i;
     400      // modify pi so can collapse to one loop
     401      for (i=0;i<numberInRowArray;i++) {
     402        int iRow = whichRow[i];
     403        pi[iRow]=scalar*piOld[i];
     404      }
     405      if (trueNetwork_) {
     406        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     407          double value = 0.0;
     408          int iRowM = indices_[j];
     409          int iRowP = indices_[j+1];
     410          value -= pi[iRowM];
     411          value += pi[iRowP];
     412          if (fabs(value)>zeroTolerance) {
     413            array[numberNonZero]=value;
     414            index[numberNonZero++]=iColumn;
     415          }
     416        }
     417      } else {
     418        // skip negative rows
     419        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     420          double value = 0.0;
     421          int iRowM = indices_[j];
     422          int iRowP = indices_[j+1];
     423          if (iRowM>=0)
     424            value -= pi[iRowM];
     425          if (iRowP>=0)
     426            value += pi[iRowP];
     427          if (fabs(value)>zeroTolerance) {
     428            array[numberNonZero]=value;
     429            index[numberNonZero++]=iColumn;
     430          }
     431        }
     432      }
     433      for (i=0;i<numberInRowArray;i++) {
     434        int iRow = whichRow[i];
     435        pi[iRow]=0.0;
    388436      }
    389437    } 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)
     438      if (trueNetwork_) {
     439        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     440          double value = 0.0;
     441          int iRowM = indices_[j];
     442          int iRowP = indices_[j+1];
    397443          value -= scalar*pi[iRowM];
    398         if (iRowP>=0)
    399444          value += scalar*pi[iRowP];
    400         if (fabs(value)>zeroTolerance) {
    401           index[numberNonZero++]=iColumn;
    402           array[iColumn]=value;
     445          if (fabs(value)>zeroTolerance) {
     446            index[numberNonZero++]=iColumn;
     447            array[iColumn]=value;
     448          }
     449        }
     450      } else {
     451        // skip negative rows
     452        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
     453          double value = 0.0;
     454          int iRowM = indices_[j];
     455          int iRowP = indices_[j+1];
     456          if (iRowM>=0)
     457            value -= scalar*pi[iRowM];
     458          if (iRowP>=0)
     459            value += scalar*pi[iRowP];
     460          if (fabs(value)>zeroTolerance) {
     461            index[numberNonZero++]=iColumn;
     462            array[iColumn]=value;
     463          }
    403464        }
    404465      }
    405466    }
    406467    columnArray->setNumElements(numberNonZero);
    407     y->setNumElements(0);
    408468  } else {
    409469    // do by row
     
    430490  int numberToDo = y->getNumElements();
    431491  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       }
     492  bool packed = rowArray->packedMode();
     493  if (packed) {
     494    // Must line up with y
     495    // need to expand pi into y
     496    int numberInRowArray = rowArray->getNumElements();
     497    assert(y->capacity()>=model->numberRows());
     498    double * piOld = pi;
     499    pi = y->denseVector();
     500    const int * whichRow = rowArray->getIndices();
     501    int i;
     502    for (i=0;i<numberInRowArray;i++) {
     503      int iRow = whichRow[i];
     504      pi[iRow]=piOld[i];
     505    }
     506    if (trueNetwork_) {
     507      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     508        int iColumn = which[jColumn];
     509        double value = 0.0;
     510        CoinBigIndex j=iColumn<<1;
     511        int iRowM = indices_[j];
     512        int iRowP = indices_[j+1];
     513        value -= pi[iRowM];
     514        value += pi[iRowP];
     515        array[jColumn]=value;
     516      }
     517    } else {
     518      // skip negative rows
     519      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     520        int iColumn = which[jColumn];
     521        double value = 0.0;
     522        CoinBigIndex j=iColumn<<1;
     523        int iRowM = indices_[j];
     524        int iRowP = indices_[j+1];
     525        if (iRowM>=0)
     526          value -= pi[iRowM];
     527        if (iRowP>=0)
     528          value += pi[iRowP];
     529        array[jColumn]=value;
     530      }
     531    }
     532    for (i=0;i<numberInRowArray;i++) {
     533      int iRow = whichRow[i];
     534      pi[iRow]=0.0;
    445535    }
    446536  } 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)
     537    if (trueNetwork_) {
     538      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     539        int iColumn = which[jColumn];
     540        double value = 0.0;
     541        CoinBigIndex j=iColumn<<1;
     542        int iRowM = indices_[j];
     543        int iRowP = indices_[j+1];
    455544        value -= pi[iRowM];
    456       if (iRowP>=0)
    457545        value += pi[iRowP];
    458       if (fabs(value)>zeroTolerance) {
    459         index[numberNonZero++]=iColumn;
    460         array[iColumn]=value;
     546        if (fabs(value)>zeroTolerance) {
     547          index[numberNonZero++]=iColumn;
     548          array[iColumn]=value;
     549        }
     550      }
     551    } else {
     552      // skip negative rows
     553      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     554        int iColumn = which[jColumn];
     555        double value = 0.0;
     556        CoinBigIndex j=iColumn<<1;
     557        int iRowM = indices_[j];
     558        int iRowP = indices_[j+1];
     559        if (iRowM>=0)
     560          value -= pi[iRowM];
     561        if (iRowP>=0)
     562          value += pi[iRowP];
     563        if (fabs(value)>zeroTolerance) {
     564          index[numberNonZero++]=iColumn;
     565          array[iColumn]=value;
     566        }
    461567      }
    462568    }
     
    542648  return numberElements;
    543649}
     650/* If element NULL returns number of elements in column part of basis,
     651   If not NULL fills in as well */
     652CoinBigIndex
     653ClpNetworkMatrix::fillBasis(const ClpSimplex * model,
     654                                 const int * whichColumn,
     655                                 int numberBasic,
     656                                 int numberColumnBasic,
     657                                 int * indexRowU, int * indexColumnU,
     658                                 double * elementU) const
     659{
     660  int i;
     661  CoinBigIndex numberElements=0;
     662  if (elementU!=NULL) {
     663    if (trueNetwork_) {
     664      for (i=0;i<numberColumnBasic;i++) {
     665        int iColumn = whichColumn[i];
     666        CoinBigIndex j=iColumn<<1;
     667        int iRowM = indices_[j];
     668        int iRowP = indices_[j+1];
     669        indexRowU[numberElements]=iRowM;
     670        indexColumnU[numberElements]=numberBasic;
     671        elementU[numberElements]=-1.0;
     672        indexRowU[numberElements+1]=iRowP;
     673        indexColumnU[numberElements+1]=numberBasic;
     674        elementU[numberElements+1]=1.0;
     675        numberElements+=2;
     676        numberBasic++;
     677      }
     678    } else {
     679      for (i=0;i<numberColumnBasic;i++) {
     680        int iColumn = whichColumn[i];
     681        CoinBigIndex j=iColumn<<1;
     682        int iRowM = indices_[j];
     683        int iRowP = indices_[j+1];
     684        if (iRowM>=0) {
     685          indexRowU[numberElements]=iRowM;
     686          indexColumnU[numberElements]=numberBasic;
     687          elementU[numberElements++]=-1.0;
     688        }
     689        if (iRowP>=0) {
     690          indexRowU[numberElements]=iRowP;
     691          indexColumnU[numberElements]=numberBasic;
     692          elementU[numberElements++]=1.0;
     693        }
     694        numberBasic++;
     695      }
     696    }
     697  } else {
     698    if (trueNetwork_) {
     699      numberElements = 2*numberColumnBasic;
     700    } else {
     701      for (i=0;i<numberColumnBasic;i++) {
     702        int iColumn = whichColumn[i];
     703        CoinBigIndex j=iColumn<<1;
     704        int iRowM = indices_[j];
     705        int iRowP = indices_[j+1];
     706        if (iRowM>=0)
     707          numberElements ++;
     708        if (iRowP>=0)
     709          numberElements ++;
     710      }
     711    }
     712  }
     713  return numberElements;
     714}
    544715/* Unpacks a column into an CoinIndexedvector
    545       Note that model is NOT const.  Bounds and objective could
    546       be modified if doing column generation */
     716 */
    547717void
    548718ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
     
    557727    rowArray->add(iRowP,1.0);
    558728}
     729/* Unpacks a column into an CoinIndexedvector
     730** in packed foramt
     731Note that model is NOT const.  Bounds and objective could
     732be modified if doing column generation (just for this variable) */
     733void
     734ClpNetworkMatrix::unpackPacked(ClpSimplex * model,
     735                            CoinIndexedVector * rowArray,
     736                            int iColumn) const
     737{
     738  int * index = rowArray->getIndices();
     739  double * array = rowArray->denseVector();
     740  int number = 0;
     741  CoinBigIndex j=iColumn<<1;
     742  int iRowM = indices_[j];
     743  int iRowP = indices_[j+1];
     744  if (iRowM>=0) {
     745    array[number]=-1.0;
     746    index[number++]=iRowM;
     747  }
     748  if (iRowP>=0) {
     749    array[number]=1.0;
     750    index[number++]=iRowP;
     751  }
     752  rowArray->setNumElements(number);
     753  rowArray->setPackedMode(true);
     754}
    559755/* Adds multiple of a column into an CoinIndexedvector
    560756      You can use quickAdd to add to vector */
     
    631827void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel)
    632828{
     829  std::cerr<<"deleteCols not implemented in ClpNetworkMatrix"<<std::endl;
    633830  abort();
    634831}
     
    636833void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel)
    637834{
     835  std::cerr<<"deleteRows not implemented in ClpNetworkMatrix"<<std::endl;
    638836  abort();
    639837}
     838/* Given positive integer weights for each row fills in sum of weights
     839   for each column (and slack).
     840   Returns weights vector
     841*/
     842CoinBigIndex *
     843ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
     844{
     845  int numberRows = model->numberRows();
     846  int numberColumns =model->numberColumns();
     847  int number = numberRows+numberColumns;
     848  CoinBigIndex * weights = new CoinBigIndex[number];
     849  int i;
     850  for (i=0;i<numberColumns;i++) {
     851    CoinBigIndex j=i<<1;
     852    CoinBigIndex count=0;
     853    int iRowM = indices_[j];
     854    int iRowP = indices_[j+1];
     855    if (iRowM>=0) {
     856      count += inputWeights[iRowM];
     857    }
     858    if (iRowP>=0) {
     859      count += inputWeights[iRowP];
     860    }
     861    weights[i]=count;
     862  }
     863  for (i=0;i<numberRows;i++) {
     864    weights[i+numberColumns]=inputWeights[i];
     865  }
     866  return weights;
     867}
     868/* Returns largest and smallest elements of both signs.
     869   Largest refers to largest absolute value.
     870*/
     871void
     872ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
     873                       double & smallestPositive, double & largestPositive)
     874{
     875  smallestNegative=-1.0;
     876  largestNegative=-1.0;
     877  smallestPositive=1.0;
     878  largestPositive=1.0;
     879}
  • trunk/ClpNonLinearCost.cpp

    r193 r225  
    4040   later may do dual analysis and even finding duplicate columns
    4141*/
    42 ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model,
    43                                      int numberOriginalColumns)
     42ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model)
    4443{
    4544  model_ = model;
     
    6867  double * cost = model_->costRegion();
    6968
    70   bool forQuadratic = (numberOriginalColumns>0);
    71 
    7269  // For quadratic we need -inf,0,0,+inf
    73   if (!forQuadratic) {
    74     for (iSequence=0;iSequence<numberTotal;iSequence++) {
    75       if (lower[iSequence]>-1.0e20)
    76         put++;
    77       if (upper[iSequence]<1.0e20)
    78         put++;
    79       put += 2;
    80     }
    81   } else {
    82     put = 4*numberTotal;
     70  for (iSequence=0;iSequence<numberTotal;iSequence++) {
     71    if (lower[iSequence]>-1.0e20)
     72      put++;
     73    if (upper[iSequence]<1.0e20)
     74      put++;
     75    put += 2;
    8376  }
    8477
     
    9487  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    9588
    96     if (!forQuadratic||iSequence<numberOriginalColumns) {
    97       if (lower[iSequence]>-COIN_DBL_MAX) {
    98         lower_[put] = -COIN_DBL_MAX;
    99         setInfeasible(put,true);
    100         cost_[put++] = cost[iSequence]-infeasibilityCost;
    101       }
    102       whichRange_[iSequence]=put;
    103       lower_[put] = lower[iSequence];
    104       cost_[put++] = cost[iSequence];
    105       lower_[put] = upper[iSequence];
    106       cost_[put++] = cost[iSequence]+infeasibilityCost;
    107       if (upper[iSequence]<1.0e20) {
    108         lower_[put] = COIN_DBL_MAX;
    109         setInfeasible(put-1,true);
    110         cost_[put++] = 1.0e50;
    111       }
    112       start_[iSequence+1]=put;
    113     } else {
    114       // quadratic  variable
     89    if (lower[iSequence]>-COIN_DBL_MAX) {
    11590      lower_[put] = -COIN_DBL_MAX;
    11691      setInfeasible(put,true);
    117       cost_[put++] = -infeasibilityCost;
    118       whichRange_[iSequence]=put;
    119       lower_[put] = max(-0.5*COIN_DBL_MAX,lower[iSequence]);
    120       cost_[put++] = 0.0;
    121       lower_[put] = min(0.5*COIN_DBL_MAX,upper[iSequence]);
    122       cost_[put++] = infeasibilityCost;
     92      cost_[put++] = cost[iSequence]-infeasibilityCost;
     93    }
     94    whichRange_[iSequence]=put;
     95    lower_[put] = lower[iSequence];
     96    cost_[put++] = cost[iSequence];
     97    lower_[put] = upper[iSequence];
     98    cost_[put++] = cost[iSequence]+infeasibilityCost;
     99    if (upper[iSequence]<1.0e20) {
    123100      lower_[put] = COIN_DBL_MAX;
    124101      setInfeasible(put-1,true);
    125       cost_[put++] = 0.0;
    126       start_[iSequence+1]=put;
    127     }
     102      cost_[put++] = 1.0e50;
     103    }
     104    start_[iSequence+1]=put;
    128105  }
    129106
     
    364341// purpose which will come to me later
    365342void
    366 ClpNonLinearCost::checkInfeasibilities(bool toNearest)
     343ClpNonLinearCost::checkInfeasibilities(double oldTolerance)
    367344{
    368345  numberInfeasibilities_=0;
     
    372349  sumInfeasibilities_ = 0.0;
    373350  double primalTolerance = model_->currentPrimalTolerance();
    374  
    375351  int iSequence;
    376352  double * solution = model_->solutionRegion();
     
    378354  double * lower = model_->lowerRegion();
    379355  double * cost = model_->costRegion();
     356  bool toNearest = oldTolerance<=0.0;
    380357   
    381358  // nonbasic should be at a valid bound
     
    454431      if (!toNearest) {
    455432        // 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);
    458           model_->setStatus(iSequence,ClpSimplex::atLowerBound);
     433        if (fabs(value-upperValue)>oldTolerance*1.0001) {
     434          if (fabs(value-lowerValue)<=oldTolerance*1.0001) {
     435            if  (fabs(value-lowerValue)>primalTolerance)
     436              solution[iSequence]=lowerValue;
     437            model_->setStatus(iSequence,ClpSimplex::atLowerBound);
     438          } else {
     439            model_->setStatus(iSequence,ClpSimplex::superBasic);
     440          }
     441        } else if  (fabs(value-upperValue)>primalTolerance) {
     442          solution[iSequence]=upperValue;
    459443        }
    460444      } else {
    461         if (fabs(value-upperValue)<=fabs(value-lowerValue)) {
    462           solution[iSequence] = upperValue;
    463         } else {
    464           model_->setStatus(iSequence,ClpSimplex::atLowerBound);
    465           solution[iSequence] = lowerValue;
    466         }
     445        // Set to nearest and make at upper bound
     446        int kRange;
     447        iRange=-1;
     448        double nearest = COIN_DBL_MAX;
     449        for (kRange=start; kRange<end;kRange++) {
     450          if (fabs(lower_[kRange]-value)<nearest) {
     451            nearest = fabs(lower_[kRange]-value);
     452            iRange=kRange;
     453          }
     454        }
     455        assert (iRange>=0);
     456        iRange--;
     457        solution[iSequence]=lower_[iRange+1];
    467458      }
    468459      break;
     
    470461      if (!toNearest) {
    471462        // 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);
     463        if (fabs(value-lowerValue)>oldTolerance*1.0001) {
     464          if (fabs(value-upperValue)<=oldTolerance*1.0001) {
     465            if  (fabs(value-upperValue)>primalTolerance)
     466              solution[iSequence]=upperValue;
     467            model_->setStatus(iSequence,ClpSimplex::atLowerBound);
     468          } else {
     469            model_->setStatus(iSequence,ClpSimplex::superBasic);
     470          }
     471        } else if  (fabs(value-lowerValue)>primalTolerance) {
     472          solution[iSequence]=lowerValue;
    475473        }
    476474      } else {
    477         if (fabs(value-lowerValue)<=fabs(value-upperValue)) {
    478           solution[iSequence] = lowerValue;
    479         } else {
    480           model_->setStatus(iSequence,ClpSimplex::atUpperBound);
    481           solution[iSequence] = upperValue;
    482         }
     475        // Set to nearest and make at lower bound
     476        int kRange;
     477        iRange=-1;
     478        double nearest = COIN_DBL_MAX;
     479        for (kRange=start; kRange<end;kRange++) {
     480          if (fabs(lower_[kRange]-value)<nearest) {
     481            nearest = fabs(lower_[kRange]-value);
     482            iRange=kRange;
     483          }
     484        }
     485        assert (iRange>=0);
     486        solution[iSequence]=lower_[iRange];
    483487      }
    484488      break;
    485489    case ClpSimplex::isFixed:
     490      if (toNearest) {
     491        // Set to true fixed
     492        for (iRange=start; iRange<end;iRange++) {
     493          if (lower_[iRange]==lower_[iRange+1])
     494            break;
     495        }
     496        solution[iSequence]=lower_[iRange];
     497      }
    486498      break;
    487499    }
     
    688700  int end = start_[iPivot+1]-1;
    689701  if (!bothWays_) {
    690     for (iRange=start; iRange<end;iRange++) {
    691       if (value<lower_[iRange+1]+primalTolerance) {
    692         // put in better range
    693         if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
    694           iRange++;
    695         break;
     702    // If fixed try and get feasible
     703    if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
     704      iRange =start+1;
     705    } else {
     706      for (iRange=start; iRange<end;iRange++) {
     707        if (value<=lower_[iRange+1]+primalTolerance) {
     708          // put in better range
     709          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
     710            iRange++;
     711          break;
     712        }
    696713      }
    697714    }
     
    755772  return difference;
    756773}
     774/* Sets bounds and cost for outgoing variable
     775   may change value
     776   Returns direction */
     777int
     778ClpNonLinearCost::setOneOutgoing(int iPivot, double & value)
     779{
     780  assert (model_!=NULL);
     781  double primalTolerance = model_->currentPrimalTolerance();
     782  // get where in bound sequence
     783  int iRange;
     784  int currentRange = whichRange_[iPivot];
     785  int start = start_[iPivot];
     786  int end = start_[iPivot+1]-1;
     787  // Set perceived direction out
     788  int direction;
     789  if (value<=lower_[currentRange]+1.001*primalTolerance) {
     790    direction=1;
     791  } else if (value>=lower_[currentRange+1]-1.001*primalTolerance) {
     792    direction=-1;
     793  } else {
     794    // odd
     795    direction=0;
     796  }
     797  // If fixed try and get feasible
     798  if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
     799    iRange =start+1;
     800  } else {
     801    // See if exact
     802    for (iRange=start; iRange<end;iRange++) {
     803      if (value==lower_[iRange+1]) {
     804        // put in better range
     805        if (infeasible(iRange)&&iRange==start)
     806          iRange++;
     807        break;
     808      }
     809    }
     810    if (iRange==end) {
     811      // not exact
     812      for (iRange=start; iRange<end;iRange++) {
     813        if (value<=lower_[iRange+1]+primalTolerance) {
     814          // put in better range
     815          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
     816            iRange++;
     817          break;
     818        }
     819      }
     820    }
     821  }
     822  assert(iRange<end);
     823  whichRange_[iPivot]=iRange;
     824  if (iRange!=currentRange) {
     825    if (infeasible(iRange))
     826      numberInfeasibilities_++;
     827    if (infeasible(currentRange))
     828      numberInfeasibilities_--;
     829  }
     830  double & lower = model_->lowerAddress(iPivot);
     831  double & upper = model_->upperAddress(iPivot);
     832  double & cost = model_->costAddress(iPivot);
     833  lower = lower_[iRange];
     834  upper = lower_[iRange+1];
     835  if (upper==lower) {
     836    value=upper;
     837  } else {
     838    // set correctly
     839    if (fabs(value-lower)<=primalTolerance*1.001){
     840      value = min(value,lower+primalTolerance);
     841    } else if (fabs(value-upper)<=primalTolerance*1.001){
     842      value = max(value,upper-primalTolerance);
     843    } else {
     844      printf("*** variable wandered off bound %g %g %g!\n",
     845             lower,value,upper);
     846      if (value-lower<=upper-value)
     847        value = lower+primalTolerance;
     848      else
     849        value = upper-primalTolerance;
     850    }
     851  }
     852  double difference = cost-cost_[iRange];
     853  cost = cost_[iRange];
     854  changeCost_ += value*difference;
     855  return direction;
     856}
    757857// Returns nearest bound
    758858double
     
    775875  return lower_[jRange];
    776876}
    777 // Sets inside bounds (i.e. non infinite - used in QP
    778 void
    779 ClpNonLinearCost::setBounds(int sequence, double lower, double upper)
    780 {
    781   int start = start_[sequence];
    782   assert(start_[sequence+1]==start+4);
    783   lower_[start+1]=lower;
    784   lower_[start+2]=upper;
    785 }
    786 
     877
  • trunk/ClpObjective.cpp

    r124 r225  
    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
  • trunk/ClpPackedMatrix.cpp

    r145 r225  
    1111#include "ClpSimplex.hpp"
    1212#include "ClpFactorization.hpp"
     13#include "ClpQuadraticObjective.hpp"
    1314// at end to get min/max!
    1415#include "ClpPackedMatrix.hpp"
     
    9091{
    9192  return new ClpPackedMatrix(*this);
     93}
     94/* Subset clone (without gaps).  Duplicates are allowed
     95   and order is as given */
     96ClpMatrixBase *
     97ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
     98                              int numberColumns,
     99                              const int * whichColumns) const
     100{
     101  return new ClpPackedMatrix(*this, numberRows, whichRows,
     102                                   numberColumns, whichColumns);
     103}
     104/* Subset constructor (without gaps).  Duplicates are allowed
     105   and order is as given */
     106ClpPackedMatrix::ClpPackedMatrix (
     107                       const ClpPackedMatrix & rhs,
     108                       int numberRows, const int * whichRows,
     109                       int numberColumns, const int * whichColumns)
     110: ClpMatrixBase(rhs)
     111{
     112  matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
     113                                 numberColumns,whichColumns);
     114  zeroElements_ = rhs.zeroElements_;
     115}
     116ClpPackedMatrix::ClpPackedMatrix (
     117                       const CoinPackedMatrix & rhs,
     118                       int numberRows, const int * whichRows,
     119                       int numberColumns, const int * whichColumns)
     120: ClpMatrixBase()
     121{
     122  matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
     123                                 numberColumns,whichColumns);
     124  zeroElements_ = false;
     125  setType(1);
    92126}
    93127
     
    156190                       const double * columnScale) const
    157191{
    158   int iRow,iColumn;
    159   // get matrix data pointers
    160   const int * row = matrix_->getIndices();
    161   const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    162   const int * columnLength = matrix_->getVectorLengths();
    163   const double * elementByColumn = matrix_->getElements();
    164   int numberColumns = matrix_->getNumCols();
    165   //memset(y,0,matrix_->getNumRows()*sizeof(double));
    166   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    167     CoinBigIndex j;
    168     double value = x[iColumn];
    169     if (value) {
    170       // scaled
    171       value *= scalar*columnScale[iColumn];
    172       for (j=columnStart[iColumn];
    173            j<columnStart[iColumn]+columnLength[iColumn];j++) {
    174         iRow=row[j];
    175         y[iRow] += value*elementByColumn[j];
    176       }
    177     }
    178   }
    179   int numberRows = getNumRows();
    180   for (iRow=0;iRow<numberRows;iRow++) {
    181     y[iRow] *= rowScale[iRow];
     192  if (rowScale) {
     193    int iRow,iColumn;
     194    // get matrix data pointers
     195    const int * row = matrix_->getIndices();
     196    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     197    const int * columnLength = matrix_->getVectorLengths();
     198    const double * elementByColumn = matrix_->getElements();
     199    int numberColumns = matrix_->getNumCols();
     200    //memset(y,0,matrix_->getNumRows()*sizeof(double));
     201    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     202      CoinBigIndex j;
     203      double value = x[iColumn];
     204      if (value) {
     205        // scaled
     206        value *= scalar*columnScale[iColumn];
     207        for (j=columnStart[iColumn];
     208             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     209          iRow=row[j];
     210          y[iRow] += value*elementByColumn[j];
     211        }
     212      }
     213    }
     214    int numberRows = getNumRows();
     215    for (iRow=0;iRow<numberRows;iRow++) {
     216      y[iRow] *= rowScale[iRow];
     217    }
     218  } else {
     219    times(scalar,x,y);
    182220  }
    183221}
     
    188226                                 const double * columnScale) const
    189227{
    190   int iColumn;
    191   // get matrix data pointers
    192   const int * row = matrix_->getIndices();
    193   const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    194   const int * columnLength = matrix_->getVectorLengths();
    195   const double * elementByColumn = matrix_->getElements();
    196   int numberColumns = matrix_->getNumCols();
    197   //memset(y,0,numberColumns*sizeof(double));
    198   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    199     CoinBigIndex j;
    200     double value=0.0;
    201     // scaled
    202     for (j=columnStart[iColumn];
    203          j<columnStart[iColumn]+columnLength[iColumn];j++) {
    204       int jRow=row[j];
     228  if (rowScale) {
     229    int iColumn;
     230    // get matrix data pointers
     231    const int * row = matrix_->getIndices();
     232    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     233    const int * columnLength = matrix_->getVectorLengths();
     234    const double * elementByColumn = matrix_->getElements();
     235    int numberColumns = matrix_->getNumCols();
     236    //memset(y,0,numberColumns*sizeof(double));
     237    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     238      CoinBigIndex j;
     239      double value=0.0;
     240      // scaled
     241      for (j=columnStart[iColumn];
     242           j<columnStart[iColumn]+columnLength[iColumn];j++) {
     243        int jRow=row[j];
    205244      value += x[jRow]*elementByColumn[j]*rowScale[jRow];
    206     }
    207     y[iColumn] += value*scalar*columnScale[iColumn];
     245      }
     246      y[iColumn] += value*scalar*columnScale[iColumn];
     247    }
     248  } else {
     249    transposeTimes(scalar,x,y);
    208250  }
    209251}
     
    227269  ClpPackedMatrix* rowCopy =
    228270    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
    229   if (numberInRowArray>0.333*numberRows||!rowCopy) {
     271  bool packed = rowArray->packedMode();
     272  double factor = 0.3;
     273  // We may not want to do by row if there may be cache problems
     274  int numberColumns = model->numberColumns();
     275  // It would be nice to find L2 cache size - for moment 512K
     276  // Be slightly optimistic
     277  if (numberColumns*sizeof(double)>1000000) {
     278    if (numberRows*10<numberColumns)
     279      factor=0.1;
     280    else if (numberRows*4<numberColumns)
     281      factor=0.15;
     282    else if (numberRows*2<numberColumns)
     283      factor=0.2;
     284    //if (model->numberIterations()%50==0)
     285    //printf("%d nonzero\n",numberInRowArray);
     286  }
     287  if (numberInRowArray>factor*numberRows||!rowCopy) {
    230288    // do by column
    231289    int iColumn;
     
    238296    int numberColumns = model->numberColumns();
    239297    if (!y->getNumElements()) {
    240       if (!rowScale) {
    241         if (scalar==-1.0) {
     298      if (packed) {
     299        // need to expand pi into y
     300        assert(y->capacity()>=numberRows);
     301        double * piOld = pi;
     302        pi = y->denseVector();
     303        const int * whichRow = rowArray->getIndices();
     304        int i;
     305        if (!rowScale) {
     306          // modify pi so can collapse to one loop
     307          for (i=0;i<numberInRowArray;i++) {
     308            int iRow = whichRow[i];
     309            pi[iRow]=scalar*piOld[i];
     310          }
    242311          for (iColumn=0;iColumn<numberColumns;iColumn++) {
    243312            double value = 0.0;
     
    249318            }
    250319            if (fabs(value)>zeroTolerance) {
     320              array[numberNonZero]=value;
    251321              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;
    267322            }
    268323          }
    269324        } 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) {
     325          // scaled
     326          // modify pi so can collapse to one loop
     327          for (i=0;i<numberInRowArray;i++) {
     328            int iRow = whichRow[i];
     329            pi[iRow]=scalar*piOld[i]*rowScale[iRow];
     330          }
    288331          for (iColumn=0;iColumn<numberColumns;iColumn++) {
    289332            double value = 0.0;
     
    293336                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
    294337              int iRow = row[j];
    295               value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     338              value += pi[iRow]*elementByColumn[j];
    296339            }
    297340            value *= columnScale[iColumn];
    298341            if (fabs(value)>zeroTolerance) {
     342              array[numberNonZero]=value;
    299343              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;
     344            }
     345          }
     346        }
     347        // zero out
     348        for (i=0;i<numberInRowArray;i++) {
     349          int iRow = whichRow[i];
     350          pi[iRow]=0.0;
     351        }
     352      } else {
     353        if (!rowScale) {
     354          if (scalar==-1.0) {
     355            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     356              double value = 0.0;
     357              CoinBigIndex j;
     358              for (j=columnStart[iColumn];
     359                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     360                int iRow = row[j];
     361                value += pi[iRow]*elementByColumn[j];
     362              }
     363              if (fabs(value)>zeroTolerance) {
     364                index[numberNonZero++]=iColumn;
     365                array[iColumn]=-value;
     366              }
     367            }
     368          } else if (scalar==1.0) {
     369            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     370              double value = 0.0;
     371              CoinBigIndex j;
     372              for (j=columnStart[iColumn];
     373                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     374                int iRow = row[j];
     375                value += pi[iRow]*elementByColumn[j];
     376              }
     377              if (fabs(value)>zeroTolerance) {
     378                index[numberNonZero++]=iColumn;
     379                array[iColumn]=value;
     380              }
     381            }
     382          } else {
     383            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     384              double value = 0.0;
     385              CoinBigIndex j;
     386              for (j=columnStart[iColumn];
     387                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     388                int iRow = row[j];
     389                value += pi[iRow]*elementByColumn[j];
     390              }
     391              value *= scalar;
     392              if (fabs(value)>zeroTolerance) {
     393                index[numberNonZero++]=iColumn;
     394                array[iColumn]=value;
     395              }
    317396            }
    318397          }
    319398        } 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;
     399          // scaled
     400          if (scalar==-1.0) {
     401            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     402              double value = 0.0;
     403              CoinBigIndex j;
     404              const double * columnScale = model->columnScale();
     405              for (j=columnStart[iColumn];
     406                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     407                int iRow = row[j];
     408                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     409              }
     410              value *= columnScale[iColumn];
     411              if (fabs(value)>zeroTolerance) {
     412                index[numberNonZero++]=iColumn;
     413                array[iColumn]=-value;
     414              }
     415            }
     416          } else if (scalar==1.0) {
     417            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     418              double value = 0.0;
     419              CoinBigIndex j;
     420              const double * columnScale = model->columnScale();
     421              for (j=columnStart[iColumn];
     422                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     423                int iRow = row[j];
     424                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     425              }
     426              value *= columnScale[iColumn];
     427              if (fabs(value)>zeroTolerance) {
     428                index[numberNonZero++]=iColumn;
     429                array[iColumn]=value;
     430              }
     431            }
     432          } else {
     433            for (iColumn=0;iColumn<numberColumns;iColumn++) {
     434              double value = 0.0;
     435              CoinBigIndex j;
     436              const double * columnScale = model->columnScale();
     437              for (j=columnStart[iColumn];
     438                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
     439                int iRow = row[j];
     440                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     441              }
     442              value *= scalar*columnScale[iColumn];
     443              if (fabs(value)>zeroTolerance) {
     444                index[numberNonZero++]=iColumn;
     445                array[iColumn]=value;
     446              }
    333447            }
    334448          }
     
    336450      }
    337451    } else {
     452      assert(!packed);
    338453      double * markVector = y->denseVector(); // not empty
    339454      if (!rowScale) {
     
    381496    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
    382497  }
     498  if (packed)
     499    columnArray->setPackedMode(true);
     500  if (0) {
     501    columnArray->checkClean();
     502    int numberNonZero=columnArray->getNumElements();;
     503    int * index = columnArray->getIndices();
     504    double * array = columnArray->denseVector();
     505    int i;
     506    for (i=0;i<numberNonZero;i++) {
     507      int j=index[i];
     508      double value;
     509      if (packed)
     510        value=array[i];
     511      else
     512        value=array[j];
     513      printf("Ti %d %d %g\n",i,j,value);
     514    }
     515  }
    383516}
    384517/* Return <code>x * A + y</code> in <code>z</code>.
     
    402535  const double * element = getElements();
    403536  const int * whichRow = rowArray->getIndices();
     537  bool packed = rowArray->packedMode();
    404538  if (numberInRowArray>2||y->getNumElements()) {
    405539    // do by rows
    406540    // ** Row copy is already scaled
    407541    int iRow;
    408     double * markVector = y->denseVector(); // probably empty .. but
    409542    int * mark = y->getIndices();
    410543    int numberOriginal=y->getNumElements();
    411544    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;
    430       for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    431         int iColumn = column[j];
    432         if (!marked[iColumn]) {
    433           marked[iColumn]=1;
     545    if (packed) {
     546      assert(!numberOriginal);
     547      numberNonZero=0;
     548      // and set up mark as char array
     549      char * marked = (char *) (index+columnArray->capacity());
     550      double * array2 = y->denseVector();
     551#ifdef CLP_DEBUG
     552      int numberColumns = model->numberColumns();
     553      for (i=0;i<numberColumns;i++) {
     554        assert(!marked[i]);
     555        assert(!array2[i]);
     556      }
     557#endif     
     558      for (i=0;i<numberInRowArray;i++) {
     559        iRow = whichRow[i];
     560        double value = pi[i]*scalar;
     561        CoinBigIndex j;
     562        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     563          int iColumn = column[j];
     564          if (!marked[iColumn]) {
     565            marked[iColumn]=1;
     566            index[numberNonZero++]=iColumn;
     567          }
     568          array2[iColumn] += value*element[j];
     569        }
     570      }
     571      // get rid of tiny values and zero out marked
     572      numberOriginal=numberNonZero;
     573      numberNonZero=0;
     574      for (i=0;i<numberOriginal;i++) {
     575        int iColumn = index[i];
     576        if (marked[iColumn]) {
     577          double value = array2[iColumn];
     578          array2[iColumn]=0.0;
     579          marked[iColumn]=0;
     580          if (fabs(value)>zeroTolerance) {
     581            array[numberNonZero]=value;
     582            index[numberNonZero++]=iColumn;
     583          }
     584        }
     585      }
     586    } else {
     587      double * markVector = y->denseVector(); // probably empty .. but
     588      for (i=0;i<numberOriginal;i++) {
     589        int iColumn = mark[i];
     590        index[i]=iColumn;
     591        array[iColumn]=markVector[iColumn];
     592        markVector[iColumn]=0.0;
     593      }
     594      numberNonZero=numberOriginal;
     595      // and set up mark as char array
     596      char * marked = (char *) markVector;
     597      for (i=0;i<numberOriginal;i++) {
     598        int iColumn = index[i];
     599        marked[iColumn]=0;
     600      }
     601     
     602      for (i=0;i<numberInRowArray;i++) {
     603        iRow = whichRow[i];
     604        double value = pi[iRow]*scalar;
     605        CoinBigIndex j;
     606        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     607          int iColumn = column[j];
     608          if (!marked[iColumn]) {
     609            marked[iColumn]=1;
     610            index[numberNonZero++]=iColumn;
     611          }
     612          array[iColumn] += value*element[j];
     613        }
     614      }
     615      // get rid of tiny values and zero out marked
     616      numberOriginal=numberNonZero;
     617      numberNonZero=0;
     618      for (i=0;i<numberOriginal;i++) {
     619        int iColumn = index[i];
     620        marked[iColumn]=0;
     621        if (fabs(array[iColumn])>zeroTolerance) {
    434622          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;
     623        } else {
     624          array[iColumn]=0.0;
     625        }
    449626      }
    450627    }
    451628  } else if (numberInRowArray==2) {
    452629    // do by rows when two rows
    453     int iRow;
    454630    int numberOriginal;
    455631    int i;
     632    CoinBigIndex j;
    456633    numberNonZero=0;
    457634
    458635    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
     636    if (packed) {
     637      int iRow0 = whichRow[0];
     638      int iRow1 = whichRow[1];
     639      double pi0 = pi[0];
     640      double pi1 = pi[1];
     641      if (rowStart[iRow0+1]-rowStart[iRow0]>
     642          rowStart[iRow1+1]-rowStart[iRow1]) {
     643        // do one with fewer first
     644        iRow0=iRow1;
     645        iRow1=whichRow[0];
     646        pi0=pi1;
     647        pi1=pi[0];
     648      }
     649      // and set up mark as char array
     650      char * marked = (char *) (index+columnArray->capacity());
     651      int * lookup = y->getIndices();
     652      value = pi0*scalar;
     653      for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
     654        int iColumn = column[j];
     655        double value2 = value*element[j];
     656        array[numberNonZero] = value2;
     657        marked[iColumn]=1;
     658        lookup[iColumn]=numberNonZero;
    477659        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) {
     660      }
     661      numberOriginal = numberNonZero;
     662      value = pi1*scalar;
     663      for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
     664        int iColumn = column[j];
     665        double value2 = value*element[j];
     666        // I am assuming no zeros in matrix
     667        if (marked[iColumn]) {
     668          int iLookup = lookup[iColumn];
     669          array[iLookup] += value2;
     670        } else {
     671          if (fabs(value2)>zeroTolerance) {
     672            array[numberNonZero] = value2;
     673            index[numberNonZero++]=iColumn;
     674          }
     675        }
     676      }
     677      // get rid of tiny values and zero out marked
     678      int nDelete=0;
     679      for (i=0;i<numberOriginal;i++) {
     680        int iColumn = index[i];
     681        marked[iColumn]=0;
     682        if (fabs(array[i])<=zeroTolerance)
     683          nDelete++;
     684      }
     685      if (nDelete) {
     686        numberOriginal=numberNonZero;
     687        numberNonZero=0;
     688        for (i=0;i<numberOriginal;i++) {
     689          int iColumn = index[i];
     690          double value = array[i];
     691          array[i]=0.0;
     692          if (fabs(value)>zeroTolerance) {
     693            array[numberNonZero]=value;
     694            index[numberNonZero++]=iColumn;
     695          }
     696        }
     697      }
     698    } else {
     699      int iRow = whichRow[0];
     700      value = pi[iRow]*scalar;
     701      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     702        int iColumn = column[j];
     703        double value2 = value*element[j];
    486704        index[numberNonZero++]=iColumn;
    487       } else {
    488         array[iColumn]=0.0;
     705        array[iColumn] = value2;
     706      }
     707      iRow = whichRow[1];
     708      value = pi[iRow]*scalar;
     709      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     710        int iColumn = column[j];
     711        double value2 = value*element[j];
     712        // I am assuming no zeros in matrix
     713        if (array[iColumn])
     714          value2 += array[iColumn];
     715        else
     716          index[numberNonZero++]=iColumn;
     717        array[iColumn] = value2;
     718      }
     719      // get rid of tiny values and zero out marked
     720      numberOriginal=numberNonZero;
     721      numberNonZero=0;
     722      for (i=0;i<numberOriginal;i++) {
     723        int iColumn = index[i];
     724        if (fabs(array[iColumn])>zeroTolerance) {
     725          index[numberNonZero++]=iColumn;
     726        } else {
     727          array[iColumn]=0.0;
     728        }
    489729      }
    490730    }
     
    493733    int iRow=rowArray->getIndices()[0];
    494734    numberNonZero=0;
    495     double value = pi[iRow]*scalar;
    496735    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;
     736    if (packed) {
     737      double value = pi[0]*scalar;
     738      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     739        int iColumn = column[j];
     740        double value2 = value*element[j];
     741        if (fabs(value2)>zeroTolerance) {
     742          array[numberNonZero] = value2;
     743          index[numberNonZero++]=iColumn;
     744        }
     745      }
     746    } else {
     747      double value = pi[iRow]*scalar;
     748      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     749        int iColumn = column[j];
     750        double value2 = value*element[j];
     751        if (fabs(value2)>zeroTolerance) {
     752          index[numberNonZero++]=iColumn;
     753          array[iColumn] = value2;
     754        }
    503755      }
    504756    }
     
    532784  int numberToDo = y->getNumElements();
    533785  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       }
     786  bool packed = rowArray->packedMode();
     787  if (packed) {
     788    // need to expand pi into y
     789    int numberInRowArray = rowArray->getNumElements();
     790    assert(y->capacity()>=model->numberRows());
     791    double * piOld = pi;
     792    pi = y->denseVector();
     793    const int * whichRow = rowArray->getIndices();
     794    int i;
     795    // Do NOT squash small elements - must line up with  y
     796    if (!rowScale) {
     797      for (i=0;i<numberInRowArray;i++) {
     798        int iRow = whichRow[i];
     799        pi[iRow]=piOld[i];
     800      }
     801      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     802        int iColumn = which[jColumn];
     803        double value = 0.0;
     804        CoinBigIndex j;
     805        for (j=columnStart[iColumn];
     806             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     807          int iRow = row[j];
     808          value += pi[iRow]*elementByColumn[j];
     809        }
     810        array[jColumn]=value;
     811      }
     812    } else {
     813      // scaled
     814      for (i=0;i<numberInRowArray;i++) {
     815        int iRow = whichRow[i];
     816        pi[iRow]=rowScale[iRow]*piOld[i];
     817      }
     818      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     819        int iColumn = which[jColumn];
     820        double value = 0.0;
     821        CoinBigIndex j;
     822        const double * columnScale = model->columnScale();
     823        for (j=columnStart[iColumn];
     824             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     825          int iRow = row[j];
     826          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     827        }
     828        value *= columnScale[iColumn];
     829        array[jColumn]=value;
     830      }
     831    }
     832    for (i=0;i<numberInRowArray;i++) {
     833      int iRow = whichRow[i];
     834      pi[iRow]=0.0;
    548835    }
    549836  } 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;
     837    if (!rowScale) {
     838      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     839        int iColumn = which[jColumn];
     840        double value = 0.0;
     841        CoinBigIndex j;
     842        for (j=columnStart[iColumn];
     843             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     844          int iRow = row[j];
     845          value += pi[iRow]*elementByColumn[j];
     846        }
     847        if (fabs(value)>zeroTolerance) {
     848          index[numberNonZero++]=iColumn;
     849          array[iColumn]=value;
     850        }
     851      }
     852    } else {
     853      // scaled
     854      for (jColumn=0;jColumn<numberToDo;jColumn++) {
     855        int iColumn = which[jColumn];
     856        double value = 0.0;
     857        CoinBigIndex j;
     858        const double * columnScale = model->columnScale();
     859        for (j=columnStart[iColumn];
     860             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     861          int iRow = row[j];
     862          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
     863        }
     864        value *= columnScale[iColumn];
     865        if (fabs(value)>zeroTolerance) {
     866          index[numberNonZero++]=iColumn;
     867          array[iColumn]=value;
     868        }
    565869      }
    566870    }
     
    681985        }
    682986      }
     987    }
     988  }
     989  return numberElements;
     990}
     991/* If element NULL returns number of elements in column part of basis,
     992   If not NULL fills in as well */
     993CoinBigIndex
     994ClpPackedMatrix::fillBasis(const ClpSimplex * model,
     995                           const int * whichColumn,
     996                           int numberBasic,
     997                           int numberColumnBasic,
     998                           int * indexRowU, int * indexColumnU,
     999                           double * elementU) const
     1000{
     1001  const int * columnLength = matrix_->getVectorLengths();
     1002  int i;
     1003  CoinBigIndex numberElements=0;
     1004  if (elementU!=NULL) {
     1005    // fill
     1006    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1007    const double * rowScale = model->rowScale();
     1008    const int * row = matrix_->getIndices();
     1009    const double * elementByColumn = matrix_->getElements();
     1010    if (!zeroElements_) {
     1011      if (!rowScale) {
     1012        // no scaling
     1013        for (i=0;i<numberColumnBasic;i++) {
     1014          int iColumn = whichColumn[i];
     1015          CoinBigIndex j;
     1016          for (j=columnStart[iColumn];
     1017               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1018            indexRowU[numberElements]=row[j];
     1019            indexColumnU[numberElements]=numberBasic;
     1020            elementU[numberElements++]=elementByColumn[j];
     1021          }
     1022          numberBasic++;
     1023        }
     1024      } else {
     1025        // scaling
     1026        const double * columnScale = model->columnScale();
     1027        for (i=0;i<numberColumnBasic;i++) {
     1028          int iColumn = whichColumn[i];
     1029          CoinBigIndex j;
     1030          double scale = columnScale[iColumn];
     1031          for (j=columnStart[iColumn];
     1032               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1033            int iRow = row[j];
     1034            indexRowU[numberElements]=iRow;
     1035            indexColumnU[numberElements]=numberBasic;
     1036            elementU[numberElements++]=
     1037              elementByColumn[j]*scale*rowScale[iRow];
     1038          }
     1039          numberBasic++;
     1040        }
     1041      }
     1042    } else {
     1043      // there are zero elements so need to look more closely
     1044      if (!rowScale) {
     1045        // no scaling
     1046        for (i=0;i<numberColumnBasic;i++) {
     1047          int iColumn = whichColumn[i];
     1048          CoinBigIndex j;
     1049          for (j=columnStart[iColumn];
     1050               j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1051            double value = elementByColumn[j];
     1052            if (value) {
     1053              indexRowU[numberElements]=row[j];
     1054              indexColumnU[numberElements]=numberBasic;
     1055              elementU[numberElements++]=value;
     1056            }
     1057          }
     1058          numberBasic++;
     1059        }
     1060      } else {
     1061        // scaling
     1062        const double * columnScale = model->columnScale();
     1063        for (i=0;i<numberColumnBasic;i++) {
     1064          int iColumn = whichColumn[i];
     1065          CoinBigIndex j;
     1066          double scale = columnScale[iColumn];
     1067          for (j=columnStart[iColumn];
     1068               j<columnStart[iColumn]+columnLength[i];j++) {
     1069            double value = elementByColumn[j];
     1070            if (value) {
     1071              int iRow = row[j];
     1072              indexRowU[numberElements]=iRow;
     1073              indexColumnU[numberElements]=numberBasic;
     1074              elementU[numberElements++]=value*scale*rowScale[iRow];
     1075            }
     1076          }
     1077          numberBasic++;
     1078        }
     1079      }
     1080    }
     1081  } else {
     1082    // just count - can be over so ignore zero problem
     1083    for (i=0;i<numberColumnBasic;i++) {
     1084      int iColumn = whichColumn[i];
     1085      numberElements += columnLength[iColumn];
    6831086    }
    6841087  }
     
    7231126  // mark empty and fixed columns
    7241127  // also see if worth scaling
    725   assert (model->scalingFlag()==1); // dynamic not implemented
     1128  assert (model->scalingFlag()<4); // dynamic not implemented
    7261129  double largest=0.0;
    7271130  double smallest=1.0e50;
     
    7611164    return 1;
    7621165  } else {
     1166      // need to scale
     1167    int scalingMethod = model->scalingFlag();
     1168    if (scalingMethod==3) {
     1169      // Choose between 1 and 2
     1170      if (smallest<1.0e-5||smallest*largest<1.0)
     1171        scalingMethod=1;
     1172      else
     1173        scalingMethod=2;
     1174    }
    7631175    // and see if there any empty rows
    7641176    for (iRow=0;iRow<numberRows;iRow++) {
     
    7781190    ClpFillN ( rowScale, numberRows,1.0);
    7791191    ClpFillN ( columnScale, numberColumns,1.0);
    780     int numberPass=3;
    781     double overallLargest;
    782     double overallSmallest;
    783     while (numberPass) {
    784       overallLargest=0.0;
    785       overallSmallest=1.0e50;
    786       numberPass--;
    787       // Geometric mean on row scales
     1192    double overallLargest=-1.0e-30;
     1193    double overallSmallest=1.0e30;
     1194    if (scalingMethod==1) {
     1195      // Maximum in each row
    7881196      for (iRow=0;iRow<numberRows;iRow++) {
    7891197        if (usefulRow[iRow]) {
    7901198          CoinBigIndex j;
    7911199          largest=1.0e-20;
    792           smallest=1.0e50;
    7931200          for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    7941201            int iColumn = column[j];
    7951202            if (usefulColumn[iColumn]) {
    7961203              double value = fabs(element[j]);
     1204              largest = max(largest,value);
     1205            }
     1206          }
     1207          rowScale[iRow]=1.0/largest;
     1208          overallLargest = max(overallLargest,largest);
     1209          overallSmallest = min(overallSmallest,largest);
     1210        }
     1211      }
     1212    } else {
     1213      assert(scalingMethod==2);
     1214      int numberPass=3;
     1215#ifdef USE_OBJECTIVE
     1216      // This will be used to help get scale factors
     1217      double * objective = new double[numberColumns];
     1218      memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
     1219      double objScale=1.0;
     1220#endif
     1221      while (numberPass) {
     1222        overallLargest=0.0;
     1223        overallSmallest=1.0e50;
     1224        numberPass--;
     1225        // Geometric mean on row scales
     1226        for (iRow=0;iRow<numberRows;iRow++) {
     1227          if (usefulRow[iRow]) {
     1228            CoinBigIndex j;
     1229            largest=1.0e-20;
     1230            smallest=1.0e50;
     1231            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     1232              int iColumn = column[j];
     1233              if (usefulColumn[iColumn]) {
     1234                double value = fabs(element[j]);
     1235                // Don't bother with tiny elements
     1236                if (value>1.0e-30) {
     1237                  value *= columnScale[iColumn];
     1238                  largest = max(largest,value);
     1239                  smallest = min(smallest,value);
     1240                }
     1241              }
     1242            }
     1243            rowScale[iRow]=1.0/sqrt(smallest*largest);
     1244            overallLargest = max(largest*rowScale[iRow],overallLargest);
     1245            overallSmallest = min(smallest*rowScale[iRow],overallSmallest);
     1246          }
     1247        }
     1248#ifdef USE_OBJECTIVE
     1249        largest=1.0e-20;
     1250        smallest=1.0e50;
     1251        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1252          if (usefulColumn[iColumn]) {
     1253            double value = fabs(objective[iColumn]);
     1254            // Don't bother with tiny elements
     1255            if (value>1.0e-30) {
     1256              value *= columnScale[iColumn];
     1257              largest = max(largest,value);
     1258              smallest = min(smallest,value);
     1259            }
     1260          }
     1261        }
     1262        objScale=1.0/sqrt(smallest*largest);
     1263#endif
     1264        model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
     1265          <<overallSmallest
     1266          <<overallLargest
     1267          <<CoinMessageEol;
     1268        // skip last column round
     1269        if (numberPass==1)
     1270          break;
     1271        // Geometric mean on column scales
     1272        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1273          if (usefulColumn[iColumn]) {
     1274            CoinBigIndex j;
     1275            largest=1.0e-20;
     1276            smallest=1.0e50;
     1277            for (j=columnStart[iColumn];
     1278                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1279              iRow=row[j];
     1280              double value = fabs(elementByColumn[j]);
    7971281              // Don't bother with tiny elements
    798               if (value>1.0e-30) {
    799                 value *= columnScale[iColumn];
     1282              if (value>1.0e-30&&usefulRow[iRow]) {
     1283                value *= rowScale[iRow];
    8001284                largest = max(largest,value);
    8011285                smallest = min(smallest,value);
    8021286              }
    8031287            }
    804           }
    805           rowScale[iRow]=1.0/sqrt(smallest*largest);
    806           overallLargest = max(largest*rowScale[iRow],overallLargest);
    807           overallSmallest = min(smallest*rowScale[iRow],overallSmallest);
    808         }
    809       }
    810       model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
    811         <<overallSmallest
    812         <<overallLargest
    813         <<CoinMessageEol;
    814       // skip last column round
    815       if (numberPass==1)
    816         break;
    817       // Geometric mean on column scales
    818       for (iColumn=0;iColumn<numberColumns;iColumn++) {
    819         if (usefulColumn[iColumn]) {
    820           CoinBigIndex j;
    821           largest=1.0e-20;
    822           smallest=1.0e50;
    823           for (j=columnStart[iColumn];
    824                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    825             iRow=row[j];
    826             double value = fabs(elementByColumn[j]);
    827             // Don't bother with tiny elements
    828             if (value>1.0e-30&&usefulRow[iRow]) {
    829               value *= rowScale[iRow];
     1288#ifdef USE_OBJECTIVE
     1289            if (fabs(objective[iColumn])>1.0e-30) {
     1290              double value = fabs(objective[iColumn])*objScale;
    8301291              largest = max(largest,value);
    8311292              smallest = min(smallest,value);
    8321293            }
    833           }
    834           columnScale[iColumn]=1.0/sqrt(smallest*largest);
    835         }
    836       }
    837     }
    838     // final pass to scale columns so largest is 1.0
    839     overallLargest=0.0;
     1294#endif
     1295            columnScale[iColumn]=1.0/sqrt(smallest*largest);
     1296          }
     1297        }
     1298      }
     1299#ifdef USE_OBJECTIVE
     1300      delete [] objective;
     1301      printf("obj scale %g - use it if you want\n",objScale);
     1302#endif
     1303    }
     1304    // If ranges will make horrid then scale
     1305    double tolerance = 5.0*model->primalTolerance();
     1306    for (iRow=0;iRow<numberRows;iRow++) {
     1307      if (usefulRow[iRow]) {
     1308        double difference = rowUpper[iRow]-rowLower[iRow];
     1309        double scaledDifference = difference*rowScale[iRow];
     1310        if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
     1311          // make gap larger
     1312          rowScale[iRow] *= 1.0e-4/scaledDifference;
     1313          //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
     1314          // scaledDifference,difference*rowScale[iRow]);
     1315        }
     1316      }
     1317    }
     1318    // final pass to scale columns so largest is reasonable
     1319    // See what smallest will be if largest is 1.0
    8401320    overallSmallest=1.0e50;
    8411321    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     
    8531333          }
    8541334        }
    855         columnScale[iColumn]=1.0/largest;
    856         overallLargest = max(overallLargest,largest*columnScale[iColumn]);
     1335        if (overallSmallest*largest>smallest)
     1336          overallSmallest = smallest/largest;
     1337      }
     1338    }
     1339#define RANDOMIZE
     1340#ifdef RANDOMIZE
     1341    // randomize by up to 10%
     1342    for (iRow=0;iRow<numberRows;iRow++) {
     1343      double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
     1344      rowScale[iRow] *= (1.0+0.1*value);
     1345    }
     1346#endif
     1347    overallLargest=1.0;
     1348    if (overallSmallest<1.0e-1)
     1349      overallLargest = 1.0/sqrt(overallSmallest);
     1350    overallLargest = min(1000.0,overallLargest);
     1351    overallSmallest=1.0e50;
     1352    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1353      if (usefulColumn[iColumn]) {
     1354        CoinBigIndex j;
     1355        largest=1.0e-20;
     1356        smallest=1.0e50;
     1357        for (j=columnStart[iColumn];
     1358             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1359          iRow=row[j];
     1360          if(elementByColumn[j]&&usefulRow[iRow]) {
     1361            double value = fabs(elementByColumn[j]*rowScale[iRow]);
     1362            largest = max(largest,value);
     1363            smallest = min(smallest,value);
     1364          }
     1365        }
     1366        columnScale[iColumn]=overallLargest/largest;
     1367#ifdef RANDOMIZE
     1368        double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
     1369        columnScale[iColumn] *= (1.0+0.1*value);
     1370#endif
     1371        double difference = columnUpper[iColumn]-columnLower[iColumn];
     1372        double scaledDifference = difference*columnScale[iColumn];
     1373        if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
     1374          // make gap larger
     1375          columnScale[iColumn] *= 1.0e-4/scaledDifference;
     1376          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
     1377          // scaledDifference,difference*columnScale[iColumn]);
     1378        }
    8571379        overallSmallest = min(overallSmallest,smallest*columnScale[iColumn]);
    8581380      }
     
    8621384      <<overallLargest
    8631385      <<CoinMessageEol;
    864    
    8651386    delete [] usefulRow;
    8661387    delete [] usefulColumn;
     1388    // If quadratic then make symmetric
     1389    ClpObjective * obj = model->objectiveAsObject();
     1390    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
     1391    if (quadraticObj) {
     1392      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     1393      int numberXColumns = quadratic->getNumCols();
     1394      if (numberXColumns<numberColumns) {
     1395        // we assume symmetric
     1396        int numberQuadraticColumns=0;
     1397        int i;
     1398        //const int * columnQuadratic = quadratic->getIndices();
     1399        //const int * columnQuadraticStart = quadratic->getVectorStarts();
     1400        const int * columnQuadraticLength = quadratic->getVectorLengths();
     1401        for (i=0;i<numberXColumns;i++) {
     1402          int length=columnQuadraticLength[i];
     1403#ifndef CORRECT_COLUMN_COUNTS
     1404          length=1;
     1405#endif
     1406          if (length)
     1407            numberQuadraticColumns++;
     1408        }
     1409        int numberXRows = numberRows-numberQuadraticColumns;
     1410        numberQuadraticColumns=0;
     1411        for (i=0;i<numberXColumns;i++) {
     1412          int length=columnQuadraticLength[i];
     1413#ifndef CORRECT_COLUMN_COUNTS
     1414          length=1;
     1415#endif
     1416          if (length) {
     1417            rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
     1418            numberQuadraticColumns++;
     1419          }
     1420        }   
     1421        int numberQuadraticRows=0;
     1422        for (i=0;i<numberXRows;i++) {
     1423          // See if any in row quadratic
     1424          int j;
     1425          int numberQ=0;
     1426          for (j=rowStart[i];j<rowStart[i+1];j++) {
     1427            int iColumn = column[j];
     1428            if (columnQuadraticLength[iColumn])
     1429              numberQ++;
     1430          }
     1431#ifndef CORRECT_ROW_COUNTS
     1432          numberQ=1;
     1433#endif
     1434          if (numberQ) {
     1435            columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
     1436            numberQuadraticRows++;
     1437          }
     1438        }
     1439        // and make sure Sj okay
     1440        for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
     1441          CoinBigIndex j=columnStart[iColumn];
     1442          assert(columnLength[iColumn]==1);
     1443          int iRow=row[j];
     1444          double value = fabs(elementByColumn[j]*rowScale[iRow]);
     1445          columnScale[iColumn]=1.0/value;
     1446        }
     1447      }
     1448    }
    8671449    model->setRowScale(rowScale);
    8681450    model->setColumnScale(columnScale);
     
    8931475}
    8941476/* Unpacks a column into an CoinIndexedvector
    895       Note that model is NOT const.  Bounds and objective could
    896       be modified if doing column generation */
     1477 */
    8971478void
    8981479ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
     
    9181499      rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
    9191500    }
     1501  }
     1502}
     1503/* Unpacks a column into a CoinIndexedVector
     1504** in packed format
     1505Note that model is NOT const.  Bounds and objective could
     1506be modified if doing column generation (just for this variable) */
     1507void
     1508ClpPackedMatrix::unpackPacked(ClpSimplex * model,
     1509                            CoinIndexedVector * rowArray,
     1510                            int iColumn) const
     1511{
     1512  const double * rowScale = model->rowScale();
     1513  const int * row = matrix_->getIndices();
     1514  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1515  const int * columnLength = matrix_->getVectorLengths();
     1516  const double * elementByColumn = matrix_->getElements();
     1517  CoinBigIndex i;
     1518  if (!rowScale) {
     1519    int j=columnStart[iColumn];
     1520    rowArray->createPacked(columnLength[iColumn],
     1521                           row+j,elementByColumn+j);
     1522  } else {
     1523    // apply scaling
     1524    double scale = model->columnScale()[iColumn];
     1525    int * index = rowArray->getIndices();
     1526    double * array = rowArray->denseVector();
     1527    int number = 0;
     1528    for (i=columnStart[iColumn];
     1529         i<columnStart[iColumn]+columnLength[iColumn];i++) {
     1530      int iRow = row[i];
     1531      array[number]=elementByColumn[i]*scale*rowScale[iRow];
     1532      index[number++]=iRow;
     1533    }
     1534    rowArray->setNumElements(number);
     1535    rowArray->setPackedMode(true);
    9201536  }
    9211537}
     
    9541570*/
    9551571bool
    956 ClpPackedMatrix::allElementsInRange(ClpSimplex * model,
     1572ClpPackedMatrix::allElementsInRange(ClpModel * model,
    9571573                                    double smallest, double largest)
    9581574{
     
    10371653  return true;
    10381654}
     1655/* Given positive integer weights for each row fills in sum of weights
     1656   for each column (and slack).
     1657   Returns weights vector
     1658*/
     1659CoinBigIndex *
     1660ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
     1661{
     1662  int numberRows = model->numberRows();
     1663  int numberColumns =model->numberColumns();
     1664  int number = numberRows+numberColumns;
     1665  CoinBigIndex * weights = new CoinBigIndex[number];
     1666  // get matrix data pointers
     1667  const int * row = matrix_->getIndices();
     1668  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1669  const int * columnLength = matrix_->getVectorLengths();
     1670  int i;
     1671  for (i=0;i<numberColumns;i++) {
     1672    CoinBigIndex j;
     1673    CoinBigIndex count=0;
     1674    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     1675      int iRow=row[j];
     1676      count += inputWeights[iRow];
     1677    }
     1678    weights[i]=count;
     1679  }
     1680  for (i=0;i<numberRows;i++) {
     1681    weights[i+numberColumns]=inputWeights[i];
     1682  }
     1683  return weights;
     1684}
     1685/* Returns largest and smallest elements of both signs.
     1686   Largest refers to largest absolute value.
     1687*/
     1688void
     1689ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
     1690                       double & smallestPositive, double & largestPositive)
     1691{
     1692  smallestNegative=-COIN_DBL_MAX;
     1693  largestNegative=0.0;
     1694  smallestPositive=COIN_DBL_MAX;
     1695  largestPositive=0.0;
     1696  int numberColumns = matrix_->getNumCols();
     1697  // get matrix data pointers
     1698  const double * elementByColumn = matrix_->getElements();
     1699  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1700  const int * columnLength = matrix_->getVectorLengths();
     1701  int i;
     1702  for (i=0;i<numberColumns;i++) {
     1703    CoinBigIndex j;
     1704    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     1705      double value = elementByColumn[j];
     1706      if (value>0.0) {
     1707        smallestPositive = min(smallestPositive,value);
     1708        largestPositive = max(largestPositive,value);
     1709      } else if (value<0.0) {
     1710        smallestNegative = max(smallestNegative,value);
     1711        largestNegative = min(largestNegative,value);
     1712      }
     1713    }
     1714  }
     1715}
    10391716
    10401717
  • trunk/ClpPlusMinusOneMatrix.cpp

    r124 r225  
    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));
     
    118119    delete [] indices_;
    119120    // put in message
    120     printf("Not all +-1 - can test if indices_ null\n");
    121121    indices_=NULL;
    122122    numberRows_=0;
     
    165165    columnOrdered_=rhs.columnOrdered_;
    166166    if (numberColumns_) {
    167       indices_ = new int [ 2*numberColumns_];
    168       memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
     167      int numberElements = rhs.startPositive_[numberColumns_];
     168      indices_ = new int [ numberElements];
     169      memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
    169170      startPositive_ = new int [ numberColumns_+1];
    170171      memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
     
    182183  return new ClpPlusMinusOneMatrix(*this);
    183184}
     185/* Subset clone (without gaps).  Duplicates are allowed
     186   and order is as given */
     187ClpMatrixBase *
     188ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
     189                              int numberColumns,
     190                              const int * whichColumns) const
     191{
     192  return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
     193                                   numberColumns, whichColumns);
     194}
     195/* Subset constructor (without gaps).  Duplicates are allowed
     196   and order is as given */
     197ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
     198                       const ClpPlusMinusOneMatrix & rhs,
     199                       int numberRows, const int * whichRow,
     200                       int numberColumns, const int * whichColumn)
     201: ClpMatrixBase(rhs)
     202
     203  elements_ = NULL;
     204  startPositive_ = NULL;
     205  startNegative_ = NULL;
     206  lengths_=NULL;
     207  indices_=NULL;
     208  numberRows_=0;
     209  numberColumns_=0;
     210  columnOrdered_=rhs.columnOrdered_;
     211  if (numberRows<=0||numberColumns<=0) {
     212    startPositive_ = new int[1];
     213    startPositive_[0] = 0;
     214  } else {
     215    numberColumns_ = numberColumns;
     216    numberRows_ = numberRows;
     217    const int * index1 = rhs.indices_;
     218    int * startPositive1 = rhs.startPositive_;
     219
     220    int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
     221    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
     222    int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
     223    int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
     224    // Also swap incoming if not column ordered
     225    if (!columnOrdered_) {
     226      int temp1 = numberRows;
     227      numberRows = numberColumns;
     228      numberColumns = temp1;
     229      const int * temp2;
     230      temp2 = whichRow;
     231      whichRow = whichColumn;
     232      whichColumn = temp2;
     233    }
     234    // Throw exception if rhs empty
     235    if (numberMajor1 <= 0 || numberMinor1 <= 0)
     236      throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
     237    // Array to say if an old row is in new copy
     238    int * newRow = new int [numberMinor1];
     239    int iRow;
     240    for (iRow=0;iRow<numberMinor1;iRow++)
     241      newRow[iRow] = -1;
     242    // and array for duplicating rows
     243    int * duplicateRow = new int [numberMinor];
     244    int numberBad=0;
     245    for (iRow=0;iRow<numberMinor;iRow++) {
     246      duplicateRow[iRow] = -1;
     247      int kRow = whichRow[iRow];
     248      if (kRow>=0  && kRow < numberMinor1) {
     249        if (newRow[kRow]<0) {
     250          // first time
     251          newRow[kRow]=iRow;
     252        } else {
     253          // duplicate
     254          int lastRow = newRow[kRow];
     255          newRow[kRow]=iRow;
     256          duplicateRow[iRow] = lastRow;
     257        }
     258      } else {
     259        // bad row
     260        numberBad++;
     261      }
     262    }
     263
     264    if (numberBad)
     265      throw CoinError("bad minor entries",
     266                      "subset constructor", "ClpPlusMinusOneMatrix");
     267    // now get size and check columns
     268    int size = 0;
     269    int iColumn;
     270    numberBad=0;
     271    for (iColumn=0;iColumn<numberMajor;iColumn++) {
     272      int kColumn = whichColumn[iColumn];
     273      if (kColumn>=0  && kColumn <numberMajor1) {
     274        int i;
     275        for (i=startPositive1[kColumn];i<startPositive1[kColumn+1];i++) {
     276          int kRow = index1[i];
     277          kRow = newRow[kRow];
     278          while (kRow>=0) {
     279            size++;
     280            kRow = duplicateRow[kRow];
     281          }
     282        }
     283      } else {
     284        // bad column
     285        numberBad++;
     286        printf("%d %d %d %d\n",iColumn,numberMajor,numberMajor1,kColumn);
     287      }
     288    }
     289    if (numberBad)
     290      throw CoinError("bad major entries",
     291                      "subset constructor", "ClpPlusMinusOneMatrix");
     292    // now create arrays
     293    startPositive_ = new int [numberMajor+1];
     294    startNegative_ = new int [numberMajor];
     295    indices_ = new int[size];
     296    // and fill them
     297    size = 0;
     298    startPositive_[0]=0;
     299    int * startNegative1 = rhs.startNegative_;
     300    for (iColumn=0;iColumn<numberMajor;iColumn++) {
     301      int kColumn = whichColumn[iColumn];
     302      int i;
     303      for (i=startPositive1[kColumn];i<startNegative1[kColumn];i++) {
     304        int kRow = index1[i];
     305        kRow = newRow[kRow];
     306        while (kRow>=0) {
     307          indices_[size++] = kRow;
     308          kRow = duplicateRow[kRow];
     309        }
     310      }
     311      startNegative_[iColumn] = size;
     312      for (;i<startPositive1[kColumn+1];i++) {
     313        int kRow = index1[i];
     314        kRow = newRow[kRow];
     315        while (kRow>=0) {
     316          indices_[size++] = kRow;
     317          kRow = duplicateRow[kRow];
     318        }
     319      }
     320      startPositive_[iColumn+1] = size;
     321    }
     322    delete [] newRow;
     323    delete [] duplicateRow;
     324  }
     325}
     326
    184327
    185328/* Returns a new matrix in reverse order without gaps */
     
    323466  double zeroTolerance = model->factorization()->zeroTolerance();
    324467  int numberRows = model->numberRows();
     468  bool packed = rowArray->packedMode();
    325469  ClpPlusMinusOneMatrix* rowCopy =
    326470    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    327   if (numberInRowArray>0.3*numberRows||!rowCopy) {
     471  double factor = 0.3;
     472  // We may not want to do by row if there may be cache problems
     473  int numberColumns = model->numberColumns();
     474  // It would be nice to find L2 cache size - for moment 512K
     475  // Be slightly optimistic
     476  if (numberColumns*sizeof(double)>1000000) {
     477    if (numberRows*10<numberColumns)
     478      factor=0.1;
     479    else if (numberRows*4<numberColumns)
     480      factor=0.15;
     481    else if (numberRows*2<numberColumns)
     482      factor=0.2;
     483  }
     484  if (numberInRowArray>factor*numberRows||!rowCopy) {
     485    assert (!y->getNumElements());
    328486    // do by column
     487    // Need to expand if packed mode
    329488    int iColumn;
    330     double * markVector = y->denseVector(); // probably empty
    331489    CoinBigIndex j=0;
    332490    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;
     491    if (packed) {
     492      // need to expand pi into y
     493      assert(y->capacity()>=numberRows);
     494      double * piOld = pi;
     495      pi = y->denseVector();
     496      const int * whichRow = rowArray->getIndices();
     497      int i;
     498      // modify pi so can collapse to one loop
     499      for (i=0;i<numberInRowArray;i++) {
     500        int iRow = whichRow[i];
     501        pi[iRow]=scalar*piOld[i];
     502      }
     503      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     504        double value = 0.0;
     505        for (;j<startNegative_[iColumn];j++) {
     506          int iRow = indices_[j];
     507          value += pi[iRow];
     508        }
     509        for (;j<startPositive_[iColumn+1];j++) {
     510          int iRow = indices_[j];
     511          value -= pi[iRow];
     512        }
     513        if (fabs(value)>zeroTolerance) {
     514          array[numberNonZero]=value;
     515          index[numberNonZero++]=iColumn;
     516        }
     517      }
     518      for (i=0;i<numberInRowArray;i++) {
     519        int iRow = whichRow[i];
     520        pi[iRow]=0.0;
     521      }
     522    } else {
     523      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     524        double value = 0.0;
     525        for (;j<startNegative_[iColumn];j++) {
     526          int iRow = indices_[j];
     527          value += pi[iRow];
     528        }
     529        for (;j<startPositive_[iColumn+1];j++) {
     530          int iRow = indices_[j];
     531          value -= pi[iRow];
     532        }
     533        value *= scalar;
     534        if (fabs(value)>zeroTolerance) {
     535          index[numberNonZero++]=iColumn;
     536          array[iColumn]=value;
     537        }
    349538      }
    350539    }
    351540    columnArray->setNumElements(numberNonZero);
    352     y->setNumElements(0);
    353541  } else {
    354542    // do by row
     
    376564  const CoinBigIndex * startNegative = startNegative_;
    377565  const int * whichRow = rowArray->getIndices();
     566  bool packed = rowArray->packedMode();
    378567  if (numberInRowArray>2||y->getNumElements()) {
    379568    // do by rows
     
    383572    int numberOriginal=y->getNumElements();
    384573    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;
     574    if (packed) {
     575      assert(!numberOriginal);
     576      numberNonZero=0;
     577      // and set up mark as char array
     578      char * marked = (char *) (index+columnArray->capacity());
     579      double * array2 = y->denseVector();
     580#ifdef CLP_DEBUG
     581      int numberColumns = model->numberColumns();
     582      for (i=0;i<numberColumns;i++) {
     583        assert(!marked[i]);
     584        assert(!array2[i]);
     585      }
     586#endif
     587      for (i=0;i<numberInRowArray;i++) {
     588        iRow = whichRow[i];
     589        double value = pi[i]*scalar;
     590        CoinBigIndex j;
     591        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     592          int iColumn = column[j];
     593          if (!marked[iColumn]) {
     594            marked[iColumn]=1;
     595            index[numberNonZero++]=iColumn;
     596          }
     597          array2[iColumn] += value;
     598        }
     599        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     600          int iColumn = column[j];
     601          if (!marked[iColumn]) {
     602            marked[iColumn]=1;
     603            index[numberNonZero++]=iColumn;
     604          }
     605          array2[iColumn] -= value;
     606        }
     607      }
     608      // get rid of tiny values and zero out marked
     609      numberOriginal=numberNonZero;
     610      numberNonZero=0;
     611      for (i=0;i<numberOriginal;i++) {
     612        int iColumn = index[i];
     613        if (marked[iColumn]) {
     614          double value = array2[iColumn];
     615          array2[iColumn]=0.0;
     616          marked[iColumn]=0;
     617          if (fabs(value)>zeroTolerance) {
     618            array[numberNonZero]=value;
     619            index[numberNonZero++]=iColumn;
     620          }
     621        }
     622      }
     623    } else {     
     624      for (i=0;i<numberOriginal;i++) {
     625        int iColumn = mark[i];
     626        index[i]=iColumn;
     627        array[iColumn]=markVector[iColumn];
     628        markVector[iColumn]=0.0;
     629      }
     630      numberNonZero=numberOriginal;
     631      // and set up mark as char array
     632      char * marked = (char *) markVector;
     633      for (i=0;i<numberOriginal;i++) {
     634        int iColumn = index[i];
     635        marked[iColumn]=0;
     636      }
     637      for (i=0;i<numberInRowArray;i++) {
     638        iRow = whichRow[i];
     639        double value = pi[iRow]*scalar;
     640        CoinBigIndex j;
     641        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     642          int iColumn = column[j];
     643          if (!marked[iColumn]) {
     644            marked[iColumn]=1;
     645            index[numberNonZero++]=iColumn;
     646          }
     647          array[iColumn] += value;
     648        }
     649        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     650          int iColumn = column[j];
     651          if (!marked[iColumn]) {
     652            marked[iColumn]=1;
     653            index[numberNonZero++]=iColumn;
     654          }
     655          array[iColumn] -= value;
     656        }
     657      }
     658      // get rid of tiny values and zero out marked
     659      numberOriginal=numberNonZero;
     660      numberNonZero=0;
     661      for (i=0;i<numberOriginal;i++) {
     662        int iColumn = index[i];
     663        marked[iColumn]=0;
     664        if (fabs(array[iColumn])>zeroTolerance) {
    406665          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;
     666        } else {
     667          array[iColumn]=0.0;
     668        }
    429669      }
    430670    }
    431671  } else if (numberInRowArray==2) {
    432     // do by rows when two rows (do longer first)
     672    /* do by rows when two rows (do longer first when not packed
     673       and shorter first if packed */
    433674    int iRow0 = whichRow[0];
    434675    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;
     676    int j;
     677    if (packed) {
     678      double pi0 = pi[0];
     679      double pi1 = pi[1];
     680      if (startPositive[iRow0+1]-startPositive[iRow0]>
     681          startPositive[iRow1+1]-startPositive[iRow1]) {
     682        int temp = iRow0;
     683        iRow0 = iRow1;
     684        iRow1 = temp;
     685        pi0=pi1;
     686        pi1=pi[0];
     687      }
     688      // and set up mark as char array
     689      char * marked = (char *) (index+columnArray->capacity());
     690      int * lookup = y->getIndices();
     691      double value = pi0*scalar;
     692      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
     693        int iColumn = column[j];
     694        array[numberNonZero] = value;
     695        marked[iColumn]=1;
     696        lookup[iColumn]=numberNonZero;
    465697        index[numberNonZero++]=iColumn;
    466698      }
    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;
     699      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
     700        int iColumn = column[j];
     701        array[numberNonZero] = -value;
     702        marked[iColumn]=1;
     703        lookup[iColumn]=numberNonZero;
    476704        index[numberNonZero++]=iColumn;
    477705      }
    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) {
     706      int numberOriginal = numberNonZero;
     707      value = pi1*scalar;
     708      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
     709        int iColumn = column[j];
     710        if (marked[iColumn]) {
     711          int iLookup = lookup[iColumn];
     712          array[iLookup] += value;
     713        } else {
     714          if (fabs(value)>zeroTolerance) {
     715            array[numberNonZero] = value;
     716            index[numberNonZero++]=iColumn;
     717          }
     718        }
     719      }
     720      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
     721        int iColumn = column[j];
     722        if (marked[iColumn]) {
     723          int iLookup = lookup[iColumn];
     724          array[iLookup] -= value;
     725        } else {
     726          if (fabs(value)>zeroTolerance) {
     727            array[numberNonZero] = -value;
     728            index[numberNonZero++]=iColumn;
     729          }
     730        }
     731      }
     732      // get rid of tiny values and zero out marked
     733      int nDelete=0;
     734      for (j=0;j<numberOriginal;j++) {
     735        int iColumn = index[j];
     736        marked[iColumn]=0;
     737        if (fabs(array[j])<=zeroTolerance)
     738          nDelete++;
     739      }
     740      if (nDelete) {
     741        numberOriginal=numberNonZero;
     742        numberNonZero=0;
     743        for (j=0;j<numberOriginal;j++) {
     744          int iColumn = index[j];
     745          double value = array[j];
     746          array[j]=0.0;
     747          if (fabs(value)>zeroTolerance) {
     748            array[numberNonZero]=value;
     749            index[numberNonZero++]=iColumn;
     750          }
     751        }
     752      }
     753    } else {
     754      if (startPositive[iRow0+1]-startPositive[iRow0]<
     755          startPositive[iRow1+1]-startPositive[iRow1]) {
     756        int temp = iRow0;
     757        iRow0 = iRow1;
     758        iRow1 = temp;
     759      }
     760      int numberOriginal;
     761      int i;
     762      numberNonZero=0;
     763      double value;
     764      value = pi[iRow0]*scalar;
     765      CoinBigIndex j;
     766      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
     767        int iColumn = column[j];
    486768        index[numberNonZero++]=iColumn;
    487       } else {
    488         array[iColumn]=0.0;
     769        array[iColumn] = value;
     770      }
     771      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
     772        int iColumn = column[j];
     773        index[numberNonZero++]=iColumn;
     774        array[iColumn] = -value;
     775      }
     776      value = pi[iRow1]*scalar;
     777      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
     778        int iColumn = column[j];
     779        double value2= array[iColumn];
     780        if (value2) {
     781          value2 += value;
     782        } else {
     783          value2 = value;
     784          index[numberNonZero++]=iColumn;
     785        }
     786        array[iColumn] = value2;
     787      }
     788      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
     789        int iColumn = column[j];
     790        double value2= array[iColumn];
     791        if (value2) {
     792          value2 -= value;
     793        } else {
     794          value2 = -value;
     795          index[numberNonZero++]=iColumn;
     796        }
     797        array[iColumn] = value2;
     798      }
     799      // get rid of tiny values and zero out marked
     800      numberOriginal=numberNonZero;
     801      numberNonZero=0;
     802      for (i=0;i<numberOriginal;i++) {
     803        int iColumn = index[i];
     804        if (fabs(array[iColumn])>zeroTolerance) {
     805          index[numberNonZero++]=iColumn;
     806        } else {
     807          array[iColumn]=0.0;
     808        }
    489809      }
    490810    }
     
    495815    double value;
    496816    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;
     817    CoinBigIndex j;
     818    if (packed) {
     819      value = pi[0]*scalar;
     820      if (fabs(value)>zeroTolerance) {
     821        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     822          int iColumn = column[j];
     823          array[numberNonZero] = value;
     824          index[numberNonZero++]=iColumn;
     825        }
     826        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     827          int iColumn = column[j];
     828          array[numberNonZero] = -value;
     829          index[numberNonZero++]=iColumn;
     830        }
     831      }
     832    } else {
     833      value = pi[iRow]*scalar;
     834      if (fabs(value)>zeroTolerance) {
     835        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     836          int iColumn = column[j];
     837          array[iColumn] = value;
     838          index[numberNonZero++]=iColumn;
     839        }
     840        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     841          int iColumn = column[j];
     842          array[iColumn] = -value;
     843          index[numberNonZero++]=iColumn;
     844        }
    509845      }
    510846    }
     
    532868  int numberToDo = y->getNumElements();
    533869  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;
     870  bool packed = rowArray->packedMode();
     871  if (packed) {
     872    // need to expand pi into y
     873    int numberInRowArray = rowArray->getNumElements();
     874    assert(y->capacity()>=model->numberRows());
     875    double * piOld = pi;
     876    pi = y->denseVector();
     877    const int * whichRow = rowArray->getIndices();
     878    int i;
     879    for (i=0;i<numberInRowArray;i++) {
     880      int iRow = whichRow[i];
     881      pi[iRow]=piOld[i];
     882    }
     883    // Must line up with y
     884    for (jColumn=0;jColumn<numberToDo;jColumn++) {
     885      int iColumn = which[jColumn];
     886      double value = 0.0;
     887      CoinBigIndex j=startPositive_[iColumn];
     888      for (;j<startNegative_[iColumn];j++) {
     889        int iRow = indices_[j];
     890        value += pi[iRow];
     891      }
     892      for (;j<startPositive_[iColumn+1];j++) {
     893        int iRow = indices_[j];
     894        value -= pi[iRow];
     895      }
     896      array[jColumn]=value;
     897    }
     898    for (i=0;i<numberInRowArray;i++) {
     899      int iRow = whichRow[i];
     900      pi[iRow]=0.0;
     901    }
     902  } else {
     903    for (jColumn=0;jColumn<numberToDo;jColumn++) {
     904      int iColumn = which[jColumn];
     905      double value = 0.0;
     906      CoinBigIndex j=startPositive_[iColumn];
     907      for (;j<startNegative_[iColumn];j++) {
     908        int iRow = indices_[j];
     909        value += pi[iRow];
     910      }
     911      for (;j<startPositive_[iColumn+1];j++) {
     912        int iRow = indices_[j];
     913        value -= pi[iRow];
     914      }
     915      if (fabs(value)>zeroTolerance) {
     916        index[numberNonZero++]=iColumn;
     917        array[iColumn]=value;
     918      }
    549919    }
    550920  }
     
    598968  return numberElements;
    599969}
     970/* If element NULL returns number of elements in column part of basis,
     971   If not NULL fills in as well */
     972CoinBigIndex
     973ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
     974                                 const int * whichColumn,
     975                                 int numberBasic,
     976                                 int numberColumnBasic,
     977                                 int * indexRowU, int * indexColumnU,
     978                                 double * elementU) const
     979{
     980  int i;
     981  CoinBigIndex numberElements=0;
     982  if (elementU!=NULL) {
     983    assert (columnOrdered_);
     984    for (i=0;i<numberColumnBasic;i++) {
     985      int iColumn = whichColumn[i];
     986      CoinBigIndex j=startPositive_[iColumn];
     987      for (;j<startNegative_[iColumn];j++) {
     988        int iRow = indices_[j];
     989        indexRowU[numberElements]=iRow;
     990        indexColumnU[numberElements]=numberBasic;
     991        elementU[numberElements++]=1.0;
     992      }
     993      for (;j<startPositive_[iColumn+1];j++) {
     994        int iRow = indices_[j];
     995        indexRowU[numberElements]=iRow;
     996        indexColumnU[numberElements]=numberBasic;
     997        elementU[numberElements++]=-1.0;
     998      }
     999      numberBasic++;
     1000    }
     1001  } else {
     1002    for (i=0;i<numberColumnBasic;i++) {
     1003      int iColumn = whichColumn[i];
     1004      numberElements += startPositive_[iColumn+1]-startPositive_[iColumn];
     1005    }
     1006  }
     1007  return numberElements;
     1008}
    6001009/* Unpacks a column into an CoinIndexedvector
    601       Note that model is NOT const.  Bounds and objective could
    602       be modified if doing column generation */
     1010 */
    6031011void
    6041012ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
     
    6151023    rowArray->add(iRow,-1.0);
    6161024  }
     1025}
     1026/* Unpacks a column into an CoinIndexedvector
     1027** in packed foramt
     1028Note that model is NOT const.  Bounds and objective could
     1029be modified if doing column generation (just for this variable) */
     1030void
     1031ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * model,
     1032                            CoinIndexedVector * rowArray,
     1033                            int iColumn) const
     1034{
     1035  int * index = rowArray->getIndices();
     1036  double * array = rowArray->denseVector();
     1037  int number = 0;
     1038  CoinBigIndex j=startPositive_[iColumn];
     1039  for (;j<startNegative_[iColumn];j++) {
     1040    int iRow = indices_[j];
     1041    array[number]=1.0;
     1042    index[number++]=iRow;
     1043  }
     1044  for (;j<startPositive_[iColumn+1];j++) {
     1045    int iRow = indices_[j];
     1046    array[number]=-1.0;
     1047    index[number++]=iRow;
     1048  }
     1049  rowArray->setNumElements(number);
     1050  rowArray->setPackedMode(true);
    6171051}
    6181052/* Adds multiple of a column into an CoinIndexedvector
     
    6931127ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel)
    6941128{
    695   abort();
     1129  int iColumn;
     1130  int newSize=startPositive_[numberColumns_];;
     1131  int numberBad=0;
     1132  // Use array to make sure we can have duplicates
     1133  int * which = new int[numberColumns_];
     1134  memset(which,0,numberColumns_*sizeof(int));
     1135  int nDuplicate=0;
     1136  for (iColumn=0;iColumn<numDel;iColumn++) {
     1137    int jColumn = indDel[iColumn];
     1138    if (jColumn<0||jColumn>=numberColumns_) {
     1139      numberBad++;
     1140    } else {
     1141      newSize -= startPositive_[jColumn+1]-startPositive_[jColumn];
     1142      if (which[jColumn])
     1143        nDuplicate++;
     1144      else
     1145        which[jColumn]=1;
     1146    }
     1147  }
     1148  if (numberBad)
     1149    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
     1150  int newNumber = numberColumns_-numDel+nDuplicate;
     1151  // Get rid of temporary arrays
     1152  delete [] lengths_;
     1153  lengths_=NULL;
     1154  delete [] elements_;
     1155  elements_= NULL;
     1156  int * newPositive = new int [newNumber+1];
     1157  int * newNegative = new int [newNumber];
     1158  int * newIndices = new int [newSize];
     1159  newNumber=0;
     1160  newSize=0;
     1161  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1162    if (!which[iColumn]) {
     1163      int start,end;
     1164      int i;
     1165      start = startPositive_[iColumn];
     1166      end=startNegative_[iColumn];
     1167      newPositive[newNumber]=newSize;
     1168      for (i=start;i<end;i++)
     1169        newIndices[newSize++]=indices_[i];
     1170      start = startNegative_[iColumn];
     1171      end=startPositive_[iColumn+1];
     1172      newNegative[newNumber++]=newSize;
     1173      for (i=start;i<end;i++)
     1174        newIndices[newSize++]=indices_[i];
     1175    }
     1176  }
     1177  newPositive[newNumber]=newSize;
     1178  delete [] which;
     1179  delete [] startPositive_;
     1180  startPositive_= newPositive;
     1181  delete [] startNegative_;
     1182  startNegative_= newNegative;
     1183  delete [] indices_;
     1184  indices_= newIndices;
     1185  numberColumns_ = newNumber;
    6961186}
    6971187/* Delete the rows whose indices are listed in <code>indDel</code>. */
     
    6991189ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel)
    7001190{
    701   abort();
     1191  int iRow;
     1192  int numberBad=0;
     1193  // Use array to make sure we can have duplicates
     1194  int * which = new int[numberRows_];
     1195  memset(which,0,numberRows_*sizeof(int));
     1196  int nDuplicate=0;
     1197  for (iRow=0;iRow<numDel;iRow++) {
     1198    int jRow = indDel[iRow];
     1199    if (jRow<0||jRow>=numberRows_) {
     1200      numberBad++;
     1201    } else {
     1202      if (which[jRow])
     1203        nDuplicate++;
     1204      else
     1205        which[jRow]=1;
     1206    }
     1207  }
     1208  if (numberBad)
     1209    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
     1210  int iElement;
     1211  int numberElements=startPositive_[numberColumns_];
     1212  int newSize=0;
     1213  for (iElement=0;iElement<numberElements;iElement++) {
     1214    iRow = indices_[iElement];
     1215    if (!which[iRow])
     1216      newSize++;
     1217  }
     1218  int newNumber = numberRows_-numDel+nDuplicate;
     1219  // Get rid of temporary arrays
     1220  delete [] lengths_;
     1221  lengths_=NULL;
     1222  delete [] elements_;
     1223  elements_= NULL;
     1224  int * newIndices = new int [newSize];
     1225  newSize=0;
     1226  int iColumn;