Changeset 1266 for trunk


Ignore:
Timestamp:
Sep 4, 2008 4:46:53 PM (12 years ago)
Author:
forrest
Message:

trying to make faster

Location:
trunk/Clp/src
Files:
15 edited

Legend:

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

    r1242 r1266  
    15041504  parameters[numberParameters-1].setDoubleValue(0.0);
    15051505  parameters[numberParameters++]=
     1506    CbcOrClpParam("dextra5","Extra double parameter 5",
     1507                  -COIN_DBL_MAX,COIN_DBL_MAX,DEXTRA5,false);
     1508  parameters[numberParameters-1].setDoubleValue(0.0);
     1509  parameters[numberParameters++]=
    15061510      CbcOrClpParam("Dins","Whether to try Distance Induced Neighborhood Search",
    15071511                    "off",DINS);
     
    15781582                    "off",DIVINGS);
    15791583  parameters[numberParameters-1].append("on");
     1584  parameters[numberParameters-1].append("do");
    15801585  parameters[numberParameters-1].setLonghelp
    15811586    (
     
    15881593                    "off",DIVINGC);
    15891594  parameters[numberParameters-1].append("on");
     1595  parameters[numberParameters-1].append("do");
    15901596  parameters[numberParameters++]=
    15911597      CbcOrClpParam("DivingF!ractional","Whether to try DiveFractional",
    15921598                    "off",DIVINGF);
    15931599  parameters[numberParameters-1].append("on");
     1600  parameters[numberParameters-1].append("do");
    15941601  parameters[numberParameters++]=
    15951602      CbcOrClpParam("DivingG!uided","Whether to try DiveGuided",
    15961603                    "off",DIVINGG);
    15971604  parameters[numberParameters-1].append("on");
     1605  parameters[numberParameters-1].append("do");
    15981606  parameters[numberParameters++]=
    15991607      CbcOrClpParam("DivingL!ineSearch","Whether to try DiveLineSearch",
    16001608                    "off",DIVINGL);
    16011609  parameters[numberParameters-1].append("on");
     1610  parameters[numberParameters-1].append("do");
    16021611  parameters[numberParameters++]=
    16031612      CbcOrClpParam("DivingP!seudoCost","Whether to try DivePseudoCost",
    16041613                    "off",DIVINGP);
    16051614  parameters[numberParameters-1].append("on");
     1615  parameters[numberParameters-1].append("do");
    16061616  parameters[numberParameters++]=
    16071617      CbcOrClpParam("DivingV!ectorLength","Whether to try DiveVectorLength",
    16081618                    "off",DIVINGV);
    16091619  parameters[numberParameters-1].append("on");
     1620  parameters[numberParameters-1].append("do");
    16101621  parameters[numberParameters++]=
    16111622    CbcOrClpParam("doH!euristic","Do heuristics before any preprocessing",
     
    17551766                  1.0,1.0e15,FAKEBOUND,false);
    17561767#ifdef COIN_HAS_CBC
    1757     parameters[numberParameters++]=
    1758       CbcOrClpParam("feas!ibilityPump","Whether to try Feasibility Pump",
    1759                     "off",FPUMP);
    1760     parameters[numberParameters-1].append("on");
    1761     parameters[numberParameters-1].append("do");
     1768  parameters[numberParameters++]=
     1769    CbcOrClpParam("feas!ibilityPump","Whether to try Feasibility Pump",
     1770                  "off",FPUMP);
     1771  parameters[numberParameters-1].append("on");
     1772  parameters[numberParameters-1].append("do");
    17621773  parameters[numberParameters-1].setLonghelp
    17631774    (
     
    22272238  parameters[numberParameters++]=
    22282239    CbcOrClpParam("passT!reeCuts","Number of cut passes in tree",
    2229                   -999999,999999,CUTPASSINTREE);
     2240                  -9999999,9999999,CUTPASSINTREE);
    22302241  parameters[numberParameters-1].setLonghelp
    22312242    (
     
    24442455\t>=1000 use index+1 as number of large loops\n\
    24452456\t>=100 use 0.05 objvalue as increment\n\
    2446 \t>=10 use +0.1 objvalue for cutoff (add)\n\
     2457\t%100 == 10,20 etc for experimentation\n\
    24472458\t1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds"
    24482459     );
     
    27262737     );
    27272738#ifdef COIN_HAS_CBC
     2739  parameters[numberParameters++]=
     2740    CbcOrClpParam("strat!egy","Switches on groups of features",
     2741                  0,2,STRATEGY);
     2742  parameters[numberParameters-1].setLonghelp
     2743    (
     2744     "This turns on newer features. \
     2745Use 0 for easy problems, 1 is default, 2 is aggressive. \
     27461 uses Gomory cuts using tolerance of 0.01 at root, \
     2747does a possible restart after 100 nodes if can fix many \
     2748and activates a diving and RINS heuristic and makes feasibility pump \
     2749more aggressive. \
     2750This does not apply to unit tests (where 'experiment' may have similar effects)."
     2751     );
     2752  parameters[numberParameters-1].setIntValue(1);
    27282753  parameters[numberParameters++]=
    27292754    CbcOrClpParam("strengthen","Create strengthened problem",
  • trunk/Clp/src/CbcOrClpParam.hpp

    r1242 r1266  
    5353   
    5454    DJFIX = 81, GAPRATIO,TIGHTENFACTOR,PRESOLVETOLERANCE,OBJSCALE2,
    55     DEXTRA1, DEXTRA2, DEXTRA3, DEXTRA4,
     55    DEXTRA1, DEXTRA2, DEXTRA3, DEXTRA4, DEXTRA5,
    5656
    5757    SOLVERLOGLEVEL=101,
     
    6666    NUMBERMINI,MIPOPTIONS,MOREMIPOPTIONS,MAXHOTITS,FPUMPITS,MAXSOLS,
    6767    FPUMPTUNE,TESTOSI,EXTRA1,EXTRA2,EXTRA3,EXTRA4,CUTPASSINTREE,
    68     THREADS,CUTPASS,VUBTRY,DENSE,EXPERIMENT,DIVEOPT,
     68    THREADS,CUTPASS,VUBTRY,DENSE,EXPERIMENT,DIVEOPT,STRATEGY,
    6969#ifdef COIN_HAS_CBC
    7070    LOGLEVEL ,
  • trunk/Clp/src/ClpFactorization.hpp

    r1197 r1266  
    176176  inline int numberDense() const
    177177  { if (coinFactorizationA_) return coinFactorizationA_->numberDense(); else return 0 ;}
    178 #if 0
     178#if 1
    179179  /// Returns number in U area
    180180  inline CoinBigIndex numberElementsU (  ) const {
    181     if (coinFactorizationA_) return coinFactorizationA_->numberElementsU(); else return 0 ;
     181    if (coinFactorizationA_) return coinFactorizationA_->numberElementsU(); else return -1 ;
    182182  }
    183183  /// Returns number in L area
    184184  inline CoinBigIndex numberElementsL (  ) const {
    185     if (coinFactorizationA_) return coinFactorizationA_->numberElementsL(); else return 0 ;
     185    if (coinFactorizationA_) return coinFactorizationA_->numberElementsL(); else return -1 ;
    186186  }
    187187  /// Returns number in R area
  • trunk/Clp/src/ClpModel.cpp

    r1264 r1266  
    306306    // later may want to keep as unknown class
    307307    CoinPackedMatrix matrix2;
     308    matrix2.setExtraGap(0.0);
     309    matrix2.setExtraMajor(0.0);
    308310    matrix2.reverseOrderedCopyOf(*matrix.getPackedMatrix());
    309311    matrix.releasePackedMatrix();
     
    319321                                const double * rowObjective)
    320322{
     323  ClpPackedMatrix* clpMatrix =
     324    dynamic_cast< ClpPackedMatrix*>(matrix_);
     325  bool special = (clpMatrix) ? clpMatrix->wantsSpecialColumnCopy() : false;
    321326  gutsOfLoadModel(matrix.getNumRows(),matrix.getNumCols(),
    322327                  collb, colub, obj, rowlb, rowub, rowObjective);
    323328  if (matrix.isColOrdered()) {
    324329    matrix_=new ClpPackedMatrix(matrix);
     330    if (special) {
     331      clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix_);
     332      clpMatrix->makeSpecialColumnCopy();
     333    }
    325334  } else {
    326335    CoinPackedMatrix matrix2;
     336    matrix2.setExtraGap(0.0);
     337    matrix2.setExtraMajor(0.0);
    327338    matrix2.reverseOrderedCopyOf(matrix);
    328339    matrix_=new ClpPackedMatrix(matrix2);
     
    37263737  if (!scaledMatrix_||!rowScale_) {
    37273738    if (rowScale_)
    3728       matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_);
     3739      matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_,NULL);
    37293740    else
    37303741      matrix_->transposeTimes(scalar,x,y);
     
    38063817  CoinModel * coinModel = new CoinModel();
    38073818  CoinPackedMatrix matrixByRow;
     3819  matrixByRow.setExtraGap(0.0);
     3820  matrixByRow.setExtraMajor(0.0);
    38083821  matrixByRow.reverseOrderedCopyOf(*matrix());
    38093822  coinModel->setObjectiveOffset(objectiveOffset());
     
    39663979  // Get a row copy in standard format
    39673980  CoinPackedMatrix * copy = new CoinPackedMatrix();
     3981  copy->setExtraGap(0.0);
     3982  copy->setExtraMajor(0.0);
    39683983  copy->reverseOrderedCopyOf(*columnCopy);
    39693984  // make sure ordered and no gaps
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1264 r1266  
    6060  numberActiveColumns_ = rhs.numberActiveColumns_;
    6161  flags_ = rhs.flags_&(~2);
    62   int numberRows = getNumRows();
     62  int numberRows = matrix_->getNumRows();
    6363  if (rhs.rhsOffset_&&numberRows) {
    6464    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
     
    7373  }
    7474  if (rhs.columnCopy_) {
    75     assert ((flags_&8)!=0);
     75    assert ((flags_&(8+16))==8+16);
    7676    columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
    7777  } else {
     
    137137    }
    138138    if (rhs.columnCopy_) {
    139       assert ((flags_&8)!=0);
     139      assert ((flags_&(8+16))==8+16);
    140140      columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
    141141    } else {
     
    224224  const int * row = matrix_->getIndices();
    225225  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    226   const int * columnLength = matrix_->getVectorLengths();
    227226  const double * elementByColumn = matrix_->getElements();
    228227  //memset(y,0,matrix_->getNumRows()*sizeof(double));
    229   for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    230     CoinBigIndex j;
    231     double value = scalar*x[iColumn];
    232     if (value) {
    233       for (j=columnStart[iColumn];
    234            j<columnStart[iColumn]+columnLength[iColumn];j++) {
    235         iRow=row[j];
    236         y[iRow] += value*elementByColumn[j];
     228  assert (((flags_&2)!=0)==matrix_->hasGaps());
     229  if (!(flags_&2)) {
     230    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     231      CoinBigIndex j;
     232      double value = x[iColumn];
     233      if (value) {
     234        CoinBigIndex start = columnStart[iColumn];
     235        CoinBigIndex end = columnStart[iColumn+1];
     236        value *= scalar;
     237        for (j=start; j<end; j++) {
     238          iRow=row[j];
     239          y[iRow] += value*elementByColumn[j];
     240        }
     241      }
     242    }
     243  } else {
     244    const int * columnLength = matrix_->getVectorLengths();
     245    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     246      CoinBigIndex j;
     247      double value = x[iColumn];
     248      if (value) {
     249        CoinBigIndex start = columnStart[iColumn];
     250        CoinBigIndex end = start + columnLength[iColumn];
     251        value *= scalar;
     252        for (j=start; j<end; j++) {
     253          iRow=row[j];
     254          y[iRow] += value*elementByColumn[j];
     255        }
    237256      }
    238257    }
     
    247266  const int * row = matrix_->getIndices();
    248267  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    249   const int * columnLength = matrix_->getVectorLengths();
    250268  const double * elementByColumn = matrix_->getElements();
    251269  if (!(flags_&2)) {
    252     if (scalar==1.0) {
    253       CoinBigIndex start=columnStart[0];
    254       for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    255         CoinBigIndex j;
    256         CoinBigIndex next=columnStart[iColumn+1];
    257         double value=y[iColumn];
    258         for (j=start;j<next;j++) {
    259           int jRow=row[j];
    260           value += x[jRow]*elementByColumn[j];
    261         }
    262         start=next;
    263         y[iColumn] = value;
    264       }
    265     } else if (scalar==-1.0) {
     270    if (scalar==-1.0) {
    266271      CoinBigIndex start=columnStart[0];
    267272      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     
    291296    }
    292297  } else {
     298    const int * columnLength = matrix_->getVectorLengths();
    293299    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    294300      CoinBigIndex j;
    295301      double value=0.0;
    296       for (j=columnStart[iColumn];
    297            j<columnStart[iColumn]+columnLength[iColumn];j++) {
     302      CoinBigIndex start = columnStart[iColumn];
     303      CoinBigIndex end = start + columnLength[iColumn];
     304      for (j=start; j<end; j++) {
    298305        int jRow=row[j];
    299306        value += x[jRow]*elementByColumn[j];
     
    305312void
    306313ClpPackedMatrix::times(double scalar,
    307                        const double * x, double * y,
    308                        const double * rowScale,
    309                        const double * columnScale) const
     314                       const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
     315                       const double * COIN_RESTRICT rowScale,
     316                       const double * COIN_RESTRICT columnScale) const
    310317{
    311318  if (rowScale) {
    312319    int iRow,iColumn;
    313320    // get matrix data pointers
    314     const int * row = matrix_->getIndices();
    315     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    316     const int * columnLength = matrix_->getVectorLengths();
    317     const double * elementByColumn = matrix_->getElements();
    318     for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    319       CoinBigIndex j;
    320       double value = x[iColumn];
    321       if (value) {
    322         // scaled
    323         value *= scalar*columnScale[iColumn];
    324         for (j=columnStart[iColumn];
    325              j<columnStart[iColumn]+columnLength[iColumn];j++) {
    326           iRow=row[j];
    327           y[iRow] += value*elementByColumn[j]*rowScale[iRow];
     321    const int * COIN_RESTRICT row = matrix_->getIndices();
     322    const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     323    const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     324    if (!(flags_&2)) {
     325      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     326        double value = x[iColumn];
     327        if (value) {
     328          // scaled
     329          value *= scalar*columnScale[iColumn];
     330          CoinBigIndex start = columnStart[iColumn];
     331          CoinBigIndex end = columnStart[iColumn+1];
     332          CoinBigIndex j;
     333          for (j=start; j<end; j++) {
     334            iRow=row[j];
     335            y[iRow] += value*elementByColumn[j]*rowScale[iRow];
     336          }
     337        }
     338      }
     339    } else {
     340      const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     341      for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     342        double value = x[iColumn];
     343        if (value) {
     344          // scaled
     345          value *= scalar*columnScale[iColumn];
     346          CoinBigIndex start = columnStart[iColumn];
     347          CoinBigIndex end = start + columnLength[iColumn];
     348          CoinBigIndex j;
     349          for (j=start; j<end; j++) {
     350            iRow=row[j];
     351            y[iRow] += value*elementByColumn[j]*rowScale[iRow];
     352          }
    328353        }
    329354      }
     
    335360void
    336361ClpPackedMatrix::transposeTimes( double scalar,
    337                                  const double * x, double * y,
    338                                  const double * rowScale,
    339                                  const double * columnScale,
    340                                  double * spare) const
     362                                 const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
     363                                 const double * COIN_RESTRICT rowScale,
     364                                 const double * COIN_RESTRICT columnScale,
     365                                 double * COIN_RESTRICT spare) const
    341366{
    342367  if (rowScale) {
    343368    int iColumn;
    344369    // get matrix data pointers
    345     const int * row = matrix_->getIndices();
    346     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    347     const int * columnLength = matrix_->getVectorLengths();
    348     const double * elementByColumn = matrix_->getElements();
     370    const int * COIN_RESTRICT row = matrix_->getIndices();
     371    const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     372    const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     373    const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    349374    if (!spare) {
    350375      if (!(flags_&2)) {
    351376        CoinBigIndex start=columnStart[0];
    352         for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    353           CoinBigIndex j;
    354           CoinBigIndex next=columnStart[iColumn+1];
    355           double value=0.0;
    356           // scaled
    357           for (j=start;j<next;j++) {
    358             int jRow=row[j];
    359             value += x[jRow]*elementByColumn[j]*rowScale[jRow];
    360           }
    361           start=next;
    362           y[iColumn] += value*scalar*columnScale[iColumn];
     377        if (scalar==-1.0) {
     378          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     379            CoinBigIndex j;
     380            CoinBigIndex next=columnStart[iColumn+1];
     381            double value=0.0;
     382            // scaled
     383            for (j=start;j<next;j++) {
     384              int jRow=row[j];
     385              value += x[jRow]*elementByColumn[j]*rowScale[jRow];
     386            }
     387            start=next;
     388            y[iColumn] -= value*columnScale[iColumn];
     389          }
     390        } else {       
     391          for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     392            CoinBigIndex j;
     393            CoinBigIndex next=columnStart[iColumn+1];
     394            double value=0.0;
     395            // scaled
     396            for (j=start;j<next;j++) {
     397              int jRow=row[j];
     398              value += x[jRow]*elementByColumn[j]*rowScale[jRow];
     399            }
     400            start=next;
     401            y[iColumn] += value*scalar*columnScale[iColumn];
     402          }
    363403        }
    364404      } else {
     
    378418      // can use spare region
    379419      int iRow;
    380       int numberRows = getNumRows();
    381       for (iRow=0;iRow<numberRows;iRow++)
    382         spare[iRow] = x[iRow]*rowScale[iRow];
     420      int numberRows = matrix_->getNumRows();
     421      for (iRow=0;iRow<numberRows;iRow++) {
     422        double value = x[iRow];
     423        if (value)
     424          spare[iRow] = value*rowScale[iRow];
     425        else
     426          spare[iRow]=0.0;
     427      }
    383428      if (!(flags_&2)) {
    384429        CoinBigIndex start=columnStart[0];
     
    414459  } else {
    415460    transposeTimes(scalar,x,y);
     461  }
     462}
     463void
     464ClpPackedMatrix::transposeTimesSubset( int number,
     465                                       const int * which,
     466                                       const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
     467                                       const double * COIN_RESTRICT rowScale,
     468                                       const double * COIN_RESTRICT columnScale,
     469                                       double * COIN_RESTRICT spare) const
     470{
     471  assert (rowScale);
     472  // get matrix data pointers
     473  const int * COIN_RESTRICT row = matrix_->getIndices();
     474  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     475  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     476  if (!spare) {
     477    for (int jColumn=0;jColumn<number;jColumn++) {
     478      int iColumn = which[jColumn];
     479      CoinBigIndex j;
     480      CoinBigIndex start=columnStart[iColumn];
     481      CoinBigIndex next=columnStart[iColumn+1];
     482      double value=0.0;
     483      for (j=start;j<next;j++) {
     484        int jRow=row[j];
     485        value += x[jRow]*elementByColumn[j]*rowScale[jRow];
     486      }
     487      y[iColumn] -= value*columnScale[iColumn];
     488    }
     489  } else {
     490    // can use spare region
     491    int iRow;
     492    int numberRows = matrix_->getNumRows();
     493    for (iRow=0;iRow<numberRows;iRow++) {
     494      double value = x[iRow];
     495      if (value)
     496        spare[iRow] = value*rowScale[iRow];
     497      else
     498        spare[iRow]=0.0;
     499    }
     500    for (int jColumn=0;jColumn<number;jColumn++) {
     501      int iColumn = which[jColumn];
     502      CoinBigIndex j;
     503      CoinBigIndex start=columnStart[iColumn];
     504      CoinBigIndex next=columnStart[iColumn+1];
     505      double value=0.0;
     506      for (j=start;j<next;j++) {
     507        int jRow=row[j];
     508        value += spare[jRow]*elementByColumn[j];
     509      }
     510      y[iColumn] -= value*columnScale[iColumn];
     511    }
    416512  }
    417513}
     
    439535  ClpPackedMatrix* rowCopy =
    440536    static_cast< ClpPackedMatrix*>(model->rowCopy());
    441 #endif
     537#endif 
    442538  bool packed = rowArray->packedMode();
    443539  double factor = (numberRows<100) ? 0.25 : 0.35;
     540  factor=0.5;
    444541  // We may not want to do by row if there may be cache problems
    445542  // It would be nice to find L2 cache size - for moment 512K
     
    463560  if (packed&&rowCopy_&&numberInRowArray>2&&numberInRowArray>factor2*numberRows&&
    464561      numberInRowArray<0.9*numberRows&&0) {
    465     rowCopy_->transposeTimes(model,rowCopy->getPackedMatrix(),rowArray,y,columnArray);
     562    rowCopy_->transposeTimes(model,rowCopy->matrix_,rowArray,y,columnArray);
    466563    return;
    467564  }
     
    643740  }
    644741}
     742//static int xA=0;
     743//static int xB=0;
     744//static int xC=0;
     745//static int xD=0;
     746//static double yA=0.0;
     747//static double yC=0.0;
    645748/* Return <code>x * scalar * A + y</code> in <code>z</code>.
    646749   Note - If x packed mode - then z packed mode
     
    653756                                        CoinIndexedVector * columnArray) const
    654757{
    655   double * pi = rowArray->denseVector();
     758  double * COIN_RESTRICT pi = rowArray->denseVector();
    656759  int numberNonZero=0;
    657   int * index = columnArray->getIndices();
    658   double * array = columnArray->denseVector();
     760  int * COIN_RESTRICT index = columnArray->getIndices();
     761  double * COIN_RESTRICT array = columnArray->denseVector();
    659762  int numberInRowArray = rowArray->getNumElements();
    660763  // maybe I need one in OsiSimplex
     
    664767  int iColumn;
    665768  // get matrix data pointers
    666   const int * row = matrix_->getIndices();
    667   const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    668   const double * elementByColumn = matrix_->getElements();
    669   const double * rowScale = model->rowScale();
     769  const int * COIN_RESTRICT row = matrix_->getIndices();
     770  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     771  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     772  const double * COIN_RESTRICT rowScale = model->rowScale();
    670773  assert (!y->getNumElements());
    671774  assert (numberActiveColumns_>0);
     
    675778    double * piOld = pi;
    676779    pi = y->denseVector();
    677     const int * whichRow = rowArray->getIndices();
     780    const int * COIN_RESTRICT whichRow = rowArray->getIndices();
    678781    int i;
    679782    if (!rowScale) {
    680783      // modify pi so can collapse to one loop
    681784      if (scalar==-1.0) {
     785        //yA += numberInRowArray;
    682786        for (i=0;i<numberInRowArray;i++) {
    683787          int iRow = whichRow[i];
     
    695799                                           zeroTolerance);
    696800        columnArray->setNumElements(numberNonZero);
     801        //xA++;
    697802      } else {
    698803        columnCopy_->transposeTimes(model,pi,columnArray);
    699804        numberNonZero = columnArray->getNumElements();
     805        //xB++;
    700806      }
    701807    } else {
     
    703809      // modify pi so can collapse to one loop
    704810      if (scalar==-1.0) {
     811        //yC += numberInRowArray;
    705812        for (i=0;i<numberInRowArray;i++) {
    706813          int iRow = whichRow[i];
     
    720827                                                 zeroTolerance);
    721828        columnArray->setNumElements(numberNonZero);
     829        //xC++;
    722830      } else {
    723831        columnCopy_->transposeTimes(model,pi,columnArray);
    724832        numberNonZero = columnArray->getNumElements();
     833        //xD++;
    725834      }
    726835    }
     
    735844      CoinZeroN(pi,numberRows);
    736845    }
     846    //int kA=xA+xB;
     847    //int kC=xC+xD;
     848    //if ((kA+kC)%10000==0)
     849    //printf("AA %d %d %g, CC %d %d %g\n",
     850    //     xA,xB,kA ? yA/(double)(kA): 0.0,xC,xD,kC ? yC/(double) (kC) :0.0);
    737851  } else {
    738852    if (!rowScale) {
     
    876990  // maybe I need one in OsiSimplex
    877991  double zeroTolerance = model->factorization()->zeroTolerance();
    878   const int * column = getIndices();
     992  const int * column = matrix_->getIndices();
    879993  const CoinBigIndex * rowStart = getVectorStarts();
    880994  const double * element = getElements();
     
    8921006      assert (!y->getNumElements());
    8931007      // and set up mark as char array
    894       char * marked = (char *) (index+columnArray->capacity());
    895       int * lookup = y->getIndices();
     1008      //char * marked = (char *) (index+columnArray->capacity());
     1009      //int * lookup = y->getIndices();
     1010      int numberColumns = matrix_->getNumCols();
     1011      //int numberRows = matrix_->getNumRows();
    8961012#ifndef NDEBUG
    897       int numberColumns = matrix_->getNumCols();
    898       for (int i=0;i<numberColumns;i++)
    899         assert(!marked[i]);
    900 #endif
    901       numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
    902                                                  lookup,marked,zeroTolerance,scalar);
     1013      //for (int i=0;i<numberColumns;i++)
     1014      //assert(!marked[i]);
     1015#endif
     1016      //if (numberInRowArray>0)
     1017        numberNonZero=gutsOfTransposeTimesByRowGEK(rowArray,index,array,
     1018                                                   numberColumns,zeroTolerance,scalar);
     1019        //else
     1020        //numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
     1021        //                                 lookup,marked,zeroTolerance,scalar);
    9031022      columnArray->setNumElements(numberNonZero);
    9041023    } else {
     
    10141133  int numberNonZero=0;
    10151134  // get matrix data pointers
    1016   const int * row = matrix_->getIndices();
    1017   const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1018   const double * elementByColumn = matrix_->getElements();
     1135  const int * COIN_RESTRICT row = matrix_->getIndices();
     1136  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1137  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    10191138  double value = 0.0;
    10201139  CoinBigIndex j;
     
    10541173  int numberNonZero=0;
    10551174  // get matrix data pointers
    1056   const int * row = matrix_->getIndices();
    1057   const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1058   const double * elementByColumn = matrix_->getElements();
     1175  const int * COIN_RESTRICT row = matrix_->getIndices();
     1176  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1177  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    10591178  double value = 0.0;
    10601179  double scale=columnScale[0];
     
    10881207  return numberNonZero;
    10891208}
    1090 // Meat of transposeTimes by row n > 2 if packed - returns number nonzero
     1209// Meat of transposeTimes by row n > K if packed - returns number nonzero
    10911210int
    1092 ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector,
     1211ClpPackedMatrix::gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector,
    10931212                                              int * COIN_RESTRICT index,
    10941213                                              double * COIN_RESTRICT output,
    1095                                               int * COIN_RESTRICT lookup,
    1096                                               char * COIN_RESTRICT marked,
     1214                                              int numberColumns,
    10971215                                              const double tolerance,
    10981216                                              const double scalar) const
    10991217{
    1100   const double * pi = piVector->denseVector();
    1101   int numberNonZero=0;
     1218  const double * COIN_RESTRICT pi = piVector->denseVector();
    11021219  int numberInRowArray = piVector->getNumElements();
    1103   const int * column = getIndices();
    1104   const CoinBigIndex * rowStart = getVectorStarts();
    1105   const double * element = getElements();
    1106   const int * whichRow = piVector->getIndices();
    1107 #ifndef NDEBUG
    1108   int maxColumn=0;
    1109 #endif
     1220  const int * COIN_RESTRICT column = matrix_->getIndices();
     1221  const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
     1222  const double * COIN_RESTRICT element = matrix_->getElements();
     1223  const int * COIN_RESTRICT whichRow = piVector->getIndices();
    11101224  // ** Row copy is already scaled
    1111   int iRow;
    1112   int i;
    1113   for (i=0;i<numberInRowArray;i++) {
    1114     iRow = whichRow[i];
     1225  for (int i=0;i<numberInRowArray;i++) {
     1226    int iRow = whichRow[i];
    11151227    double value = pi[i]*scalar;
    11161228    CoinBigIndex j;
     1229    // could do by twos
    11171230    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    11181231      int iColumn = column[j];
    1119 #ifndef NDEBUG
    1120       maxColumn=CoinMax(maxColumn,iColumn);
    1121 #endif
    11221232      double elValue = element[j];
    11231233      elValue *= value;
    1124       if (marked[iColumn]) {
    1125         int k = lookup[iColumn];
    1126         output[k] += elValue;
    1127       } else {
    1128         output[numberNonZero] = elValue;
    1129         marked[iColumn]=1;
    1130         lookup[iColumn]=numberNonZero;
    1131         index[numberNonZero++]=iColumn;
     1234      output[iColumn] += elValue;
     1235    }
     1236  }
     1237  // get rid of tiny values and count
     1238  int numberNonZero=0;
     1239  for (int i=0;i<numberColumns;i++) {
     1240    double value = output[i];
     1241    if (value) {
     1242      output[i]=0.0;
     1243      if (fabs(value)>tolerance) {
     1244        output[numberNonZero]=value;
     1245        index[numberNonZero++] = i;
    11321246      }
    11331247    }
    11341248  }
    11351249#ifndef NDEBUG
    1136   int saveN = numberNonZero;
    1137 #endif
    1138   // get rid of tiny values and zero out lookup
    1139   for (i=0;i<numberNonZero;i++) {
    1140     int iColumn = index[i];
    1141     marked[iColumn]=0;
    1142     double value = output[i];
    1143     if (fabs(value)<=tolerance) {
    1144       while (fabs(value)<=tolerance) {
    1145         numberNonZero--;
    1146         value = output[numberNonZero];
    1147         iColumn = index[numberNonZero];
    1148         marked[iColumn]=0;
    1149         if (i<numberNonZero) {
    1150           output[numberNonZero]=0.0;
    1151           output[i] = value;
    1152           index[i] = iColumn;
    1153         } else {
    1154           output[i]=0.0;
    1155           value =1.0; // to force end of while
    1156         }
    1157       }
    1158     }
    1159   }
    1160 #ifndef NDEBUG
    1161   for (i=numberNonZero;i<saveN;i++)
     1250  for (int i=numberNonZero;i<numberColumns;i++)
    11621251    assert(!output[i]);
    1163   for (i=0;i<=maxColumn;i++)
    1164     assert (!marked[i]);
    11651252#endif
    11661253  return numberNonZero;
     
    11751262  int * index = output->getIndices();
    11761263  double * array = output->denseVector();
    1177   const int * column = getIndices();
    1178   const CoinBigIndex * rowStart = getVectorStarts();
    1179   const double * element = getElements();
     1264  const int * column = matrix_->getIndices();
     1265  const CoinBigIndex * rowStart = matrix_->getVectorStarts();
     1266  const double * element = matrix_->getElements();
    11801267  const int * whichRow = piVector->getIndices();
    11811268  int iRow0 = whichRow[0];
     
    12731360  int * index = output->getIndices();
    12741361  double * array = output->denseVector();
    1275   const int * column = getIndices();
    1276   const CoinBigIndex * rowStart = getVectorStarts();
    1277   const double * element = getElements();
     1362  const int * column = matrix_->getIndices();
     1363  const CoinBigIndex * rowStart = matrix_->getVectorStarts();
     1364  const double * element = matrix_->getElements();
    12781365  int iRow=piVector->getIndices()[0];
    12791366  numberNonZero=0;
     
    13011388{
    13021389  columnArray->clear();
    1303   double * pi = rowArray->denseVector();
    1304   double * array = columnArray->denseVector();
     1390  double * COIN_RESTRICT pi = rowArray->denseVector();
     1391  double * COIN_RESTRICT array = columnArray->denseVector();
    13051392  int jColumn;
    13061393  // get matrix data pointers
    1307   const int * row = matrix_->getIndices();
    1308   const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1309   const int * columnLength = matrix_->getVectorLengths();
    1310   const double * elementByColumn = matrix_->getElements();
    1311   const double * rowScale = model->rowScale();
     1394  const int * COIN_RESTRICT row = matrix_->getIndices();
     1395  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1396  const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
     1397  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1398  const double * COIN_RESTRICT rowScale = model->rowScale();
    13121399  int numberToDo = y->getNumElements();
    1313   const int * which = y->getIndices();
     1400  const int * COIN_RESTRICT which = y->getIndices();
    13141401  assert (!rowArray->packedMode());
    13151402  columnArray->setPacked();
     
    18231910void
    18241911ClpPackedMatrix::fillBasis(ClpSimplex * model,
    1825                          const int * whichColumn,
    1826                          int & numberColumnBasic,
    1827                          int * indexRowU, int * start,
    1828                          int * rowCount, int * columnCount,
    1829                          double * elementU)
    1830 {
    1831   const int * columnLength = matrix_->getVectorLengths();
     1912                           const int * COIN_RESTRICT whichColumn,
     1913                           int & numberColumnBasic,
     1914                           int * COIN_RESTRICT indexRowU,
     1915                           int * COIN_RESTRICT start,
     1916                           int * COIN_RESTRICT rowCount,
     1917                           int * COIN_RESTRICT columnCount,
     1918                           double * COIN_RESTRICT elementU)
     1919{
     1920  const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
    18321921  int i;
    18331922  CoinBigIndex numberElements=start[0];
    18341923  // fill
    1835   const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1836   const double * rowScale = model->rowScale();
    1837   const int * row = matrix_->getIndices();
    1838   const double * elementByColumn = matrix_->getElements();
     1924  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1925  const double * COIN_RESTRICT rowScale = model->rowScale();
     1926  const int * COIN_RESTRICT row = matrix_->getIndices();
     1927  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
    18391928  ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
    18401929  if (scaledMatrix&&true) {
    1841     columnLength = scaledMatrix->getVectorLengths();
    1842     columnStart = scaledMatrix->getVectorStarts();
     1930    columnLength = scaledMatrix->matrix_->getVectorLengths();
     1931    columnStart = scaledMatrix->matrix_->getVectorStarts();
    18431932    rowScale = NULL;
    1844     row = scaledMatrix->getIndices();
    1845     elementByColumn = scaledMatrix->getElements();
     1933    row = scaledMatrix->matrix_->getIndices();
     1934    elementByColumn = scaledMatrix->matrix_->getElements();
    18461935  } 
    18471936  if ((flags_&1)==0) {
     
    18501939      for (i=0;i<numberColumnBasic;i++) {
    18511940        int iColumn = whichColumn[i];
    1852         CoinBigIndex j;
    1853         for (j=columnStart[iColumn];
    1854              j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1941        int length = columnLength[iColumn];
     1942        CoinBigIndex startThis = columnStart[iColumn];
     1943        columnCount[i]=length;
     1944        CoinBigIndex endThis = startThis+length;
     1945        for (CoinBigIndex j=startThis;j<endThis;j++) {
    18551946          int iRow=row[j];
    18561947          indexRowU[numberElements]=iRow;
     
    18591950        }
    18601951        start[i+1]=numberElements;
    1861         columnCount[i]=columnLength[iColumn];
    18621952      }
    18631953    } else {
    18641954      // scaling
    1865       const double * columnScale = model->columnScale();
     1955      const double * COIN_RESTRICT columnScale = model->columnScale();
    18661956      for (i=0;i<numberColumnBasic;i++) {
    18671957        int iColumn = whichColumn[i];
    1868         CoinBigIndex j;
    18691958        double scale = columnScale[iColumn];
    1870         for (j=columnStart[iColumn];
    1871              j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1959        int length = columnLength[iColumn];
     1960        CoinBigIndex startThis = columnStart[iColumn];
     1961        columnCount[i]=length;
     1962        CoinBigIndex endThis = startThis+length;
     1963        for (CoinBigIndex j=startThis;j<endThis;j++) {
    18721964          int iRow = row[j];
    18731965          indexRowU[numberElements]=iRow;
     
    18771969        }
    18781970        start[i+1]=numberElements;
    1879         columnCount[i]=columnLength[iColumn];
    18801971      }
    18811972    }
     
    28632954  return true;
    28642955}
     2956int
     2957ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector,
     2958                                              int * COIN_RESTRICT index,
     2959                                              double * COIN_RESTRICT output,
     2960                                              int * COIN_RESTRICT lookup,
     2961                                              char * COIN_RESTRICT marked,
     2962                                              const double tolerance,
     2963                                              const double scalar) const
     2964{
     2965  const double * COIN_RESTRICT pi = piVector->denseVector();
     2966  int numberNonZero=0;
     2967  int numberInRowArray = piVector->getNumElements();
     2968  const int * COIN_RESTRICT column = matrix_->getIndices();
     2969  const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
     2970  const double * COIN_RESTRICT element = matrix_->getElements();
     2971  const int * COIN_RESTRICT whichRow = piVector->getIndices();
     2972#ifndef NDEBUG
     2973  int maxColumn=0;
     2974#endif
     2975  // ** Row copy is already scaled
     2976  int iRow;
     2977  int i;
     2978  for (i=0;i<numberInRowArray;i++) {
     2979    iRow = whichRow[i];
     2980    double value = pi[i]*scalar;
     2981    CoinBigIndex j;
     2982    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     2983      int iColumn = column[j];
     2984#ifndef NDEBUG
     2985      maxColumn=CoinMax(maxColumn,iColumn);
     2986#endif
     2987      double elValue = element[j];
     2988      elValue *= value;
     2989      if (marked[iColumn]) {
     2990        int k = lookup[iColumn];
     2991        output[k] += elValue;
     2992      } else {
     2993        output[numberNonZero] = elValue;
     2994        marked[iColumn]=1;
     2995        lookup[iColumn]=numberNonZero;
     2996        index[numberNonZero++]=iColumn;
     2997      }
     2998    }
     2999  }
     3000#ifndef NDEBUG
     3001  int saveN = numberNonZero;
     3002#endif
     3003  // get rid of tiny values and zero out lookup
     3004  for (i=0;i<numberNonZero;i++) {
     3005    int iColumn = index[i];
     3006    marked[iColumn]=0;
     3007    double value = output[i];
     3008    if (fabs(value)<=tolerance) {
     3009      while (fabs(value)<=tolerance) {
     3010        numberNonZero--;
     3011        value = output[numberNonZero];
     3012        iColumn = index[numberNonZero];
     3013        marked[iColumn]=0;
     3014        if (i<numberNonZero) {
     3015          output[numberNonZero]=0.0;
     3016          output[i] = value;
     3017          index[i] = iColumn;
     3018        } else {
     3019          output[i]=0.0;
     3020          value =1.0; // to force end of while
     3021        }
     3022      }
     3023    }
     3024  }
     3025#ifndef NDEBUG
     3026  for (i=numberNonZero;i<saveN;i++)
     3027    assert(!output[i]);
     3028  for (i=0;i<=maxColumn;i++)
     3029    assert (!marked[i]);
     3030#endif
     3031  return numberNonZero;
     3032}
    28653033/* Given positive integer weights for each row fills in sum of weights
    28663034   for each column (and slack).
     
    32013369  assert (rowCopy!=NULL);
    32023370
    3203   const int * column = rowCopy->getIndices();
    3204   const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    3205   const int * rowLength = rowCopy->getVectorLengths();
    3206   const double * element = rowCopy->getElements();
    3207   for (int i=0;i<rowCopy->getNumRows();i++) {
     3371  const int * column = rowCopy->matrix_->getIndices();
     3372  const CoinBigIndex * rowStart = rowCopy->matrix_->getVectorStarts();
     3373  const int * rowLength = rowCopy->matrix_->getVectorLengths();
     3374  const double * element = rowCopy->matrix_->getElements();
     3375  int numberRows = rowCopy->matrix_->getNumRows();
     3376  for (int i=0;i<numberRows;i++) {
    32083377    if (!rowLength[i])
    32093378      printf("zero row %d\n",i);
     
    32223391    // need to replace row by row
    32233392    int numberRows = model->numberRows();
     3393#ifndef NDEBUG
    32243394    int numberColumns = matrix_->getNumCols();
    3225     double * newElement = new double[numberColumns];
     3395#endif
    32263396    ClpMatrixBase * rowCopyBase=model->rowCopy();
    32273397#ifndef NO_RTTI
     
    32373407    const int * column = rowCopy->getIndices();
    32383408    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    3239     const double * element = rowCopy->getElements();
     3409    double * element = rowCopy->getMutableElements();
    32403410    const double * rowScale = model->rowScale();
    32413411    const double * columnScale = model->columnScale();
     
    32443414      CoinBigIndex j;
    32453415      double scale = rowScale[iRow];
    3246       const double * elementsInThisRow = element + rowStart[iRow];
     3416      double * elementsInThisRow = element + rowStart[iRow];
    32473417      const int * columnsInThisRow = column + rowStart[iRow];
    32483418      int number = rowStart[iRow+1]-rowStart[iRow];
     
    32503420      for (j=0;j<number;j++) {
    32513421        int iColumn = columnsInThisRow[j];
    3252         newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
    3253       }
    3254       rowCopy->replaceVector(iRow,number,newElement);
    3255     }
    3256     delete [] newElement;
     3422        elementsInThisRow[j] *= scale*columnScale[iColumn];
     3423      }
     3424    }
    32573425  }
    32583426}
     
    32643432{
    32653433  // need to replace column by column
     3434#ifndef NDEBUG
    32663435  int numberRows = model->numberRows();
     3436#endif
    32673437  int numberColumns = matrix_->getNumCols();
    3268   double * newElement = new double[numberRows];
    32693438  ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
    32703439  const int * row = copy->getIndices();
    32713440  const CoinBigIndex * columnStart = copy->getVectorStarts();
    32723441  const int * length = copy->getVectorLengths();
    3273   const double * element = copy->getElements();
     3442  double * element = copy->getMutableElements();
    32743443  const double * rowScale = model->rowScale();
    32753444  const double * columnScale = model->columnScale();
     
    32783447    CoinBigIndex j;
    32793448    double scale = columnScale[iColumn];
    3280     const double * elementsInThisColumn = element + columnStart[iColumn];
     3449    double * elementsInThisColumn = element + columnStart[iColumn];
    32813450    const int * rowsInThisColumn = row + columnStart[iColumn];
    32823451    int number = length[iColumn];
     
    32843453    for (j=0;j<number;j++) {
    32853454      int iRow = rowsInThisColumn[j];
    3286       newElement[j] = elementsInThisColumn[j]*scale*rowScale[iRow];
    3287     }
    3288     copy->replaceVector(iColumn,number,newElement);
    3289   }
    3290   delete [] newElement;
     3455      elementsInThisColumn[j] *= scale*rowScale[iRow];
     3456    }
     3457  }
    32913458  return copy;
    32923459}
     
    33803547    if (matrix_->isColOrdered()&&numberOther>matrix_->getNumCols())
    33813548      matrix_->setDimensions(-1,numberOther);
    3382     numberErrors=matrix_->appendRows(number,starts,index,element,numberOther);
     3549    if (!matrix_->isColOrdered()||numberOther>=0||matrix_->getExtraGap()) {
     3550      numberErrors=matrix_->appendRows(number,starts,index,element,numberOther);
     3551    } else {
     3552      //CoinPackedMatrix mm(*matrix_);
     3553      matrix_->appendMinorFast(number,starts,index,element);
     3554      //mm.appendRows(number,starts,index,element,numberOther);
     3555      //if (!mm.isEquivalent(*matrix_)) {
     3556      //printf("bad append\n");
     3557      //abort();
     3558      //}
     3559    } 
    33833560  } else {
    33843561    // columns
     
    34093586{
    34103587  delete columnCopy_;
    3411   if ((flags_&8)!=0)
     3588  if ((flags_&16)!=0) {
    34123589    columnCopy_ = new ClpPackedMatrix3(model,matrix_);
    3413   else
     3590    flags_ |= 8;
     3591  } else {
    34143592    columnCopy_=NULL;
     3593  }
    34153594}
    34163595// Say we don't want special column copy
     
    34183597ClpPackedMatrix::releaseSpecialColumnCopy()
    34193598{
    3420   flags_ &= ~8;
     3599  flags_ &= ~(8+16);
    34213600  delete columnCopy_;
    34223601  columnCopy_=NULL;
     
    43434522#define MAXBLOCK 100
    43444523#define MAXUNROLL 10
    4345   numberColumns_ = columnCopy->getNumCols();
     4524  numberColumns_ = model->getNumCols();
     4525  int numberColumns = columnCopy->getNumCols();
     4526  assert (numberColumns_ >= numberColumns);
    43464527  int numberRows = columnCopy->getNumRows();
    43474528  int * counts = new int[numberRows+1];
     
    43554536  const int * columnLength = columnCopy->getVectorLengths();
    43564537  const double * elementByColumn = columnCopy->getElements();
    4357   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     4538  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    43584539    CoinBigIndex start = columnStart[iColumn];
    43594540    int n = columnLength[iColumn];
     
    43694550    counts[n]++;
    43704551  }
     4552  counts[0] += numberColumns_-numberColumns;
    43714553  int nZeroColumns = counts[0];
    43724554  counts[0]=-1;
    43734555  numberColumns_ -= nZeroColumns;
    4374   column_ = new int [2*numberColumns_];
     4556  column_ = new int [2*numberColumns_+nZeroColumns];
    43754557  int * lookup = column_ + numberColumns_;
    43764558  row_ = new int[nels];
     
    44184600    }
    44194601  }
     4602  for (iColumn=numberColumns;iColumn<numberColumns_;iColumn++)
     4603    lookup[iColumn]=-1;
    44204604  // fill
    44214605  start_[0]=0;
     
    44234607  nInOdd=0;
    44244608  const double * columnScale = model->columnScale();
    4425   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     4609  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    44264610    CoinBigIndex start = columnStart[iColumn];
    44274611    int n = columnLength[iColumn];
     
    44674651        start_[nOdd]=nInOdd;
    44684652      }
     4653    } else {
     4654      // zero column
     4655      lookup[iColumn]=-1;
    44694656    }
    44704657  }
     
    46544841  int kB;
    46554842  if (moveUp) {
     4843    // May already be in correct place (e.g. fixed basic leaving basis)
     4844    if (kA>=lastPrice)
     4845      return;
    46564846    kB=lastPrice-1;
    46574847    block->numberPrice_--;
    46584848  } else {
     4849    assert (kA>=lastPrice);
    46594850    kB=lastPrice;
    46604851    block->numberPrice_++;
     
    46784869    elementB[i]=tempE;
    46794870  }
    4680 #if 1
    46814871#ifndef NDEBUG
    46824872  // check
     
    46964886    assert (lookup[iColumn]==i);
    46974887  }
    4698 #endif
    46994888#endif
    47004889}
  • trunk/Clp/src/ClpPackedMatrix.hpp

    r1200 r1266  
    2323   /// Return a complete CoinPackedMatrix
    2424  virtual CoinPackedMatrix * getPackedMatrix() const { return matrix_;}
    25     /** Whether the packed matrix is column major ordered or not. */
    26     virtual bool isColOrdered() const { return matrix_->isColOrdered(); }
     25  /** Whether the packed matrix is column major ordered or not. */
     26  virtual bool isColOrdered() const { return matrix_->isColOrdered(); }
    2727   /** Number of entries in the packed matrix. */
    2828  virtual  CoinBigIndex getNumElements() const
    2929  { return matrix_->getNumElements(); }
    30    /** Number of columns. */
    31    virtual int getNumCols() const { return matrix_->getNumCols(); }
    32    /** Number of rows. */
     30  /** Number of columns. */
     31  virtual int getNumCols() const { return matrix_->getNumCols(); }
     32  /** Number of rows. */
    3333  virtual int getNumRows() const { return matrix_->getNumRows(); }
    34 
    35    /** A vector containing the elements in the packed matrix. Note that there
    36         might be gaps in this list, entries that do not belong to any
    37         major-dimension vector. To get the actual elements one should look at
    38         this vector together with vectorStarts and vectorLengths. */
    39    virtual const double * getElements() const
     34 
     35  /** A vector containing the elements in the packed matrix. Note that there
     36      might be gaps in this list, entries that do not belong to any
     37      major-dimension vector. To get the actual elements one should look at
     38      this vector together with vectorStarts and vectorLengths. */
     39  virtual const double * getElements() const
    4040  { return matrix_->getElements();}
    4141  /// Mutable elements
     
    199199                                const double * columnScale,
    200200                                double * spare=NULL) const;
     201    /** Return <code>y - pi * A</code> in <code>y</code>.
     202        @pre <code>pi</code> must be of size <code>numRows()</code>
     203        @pre <code>y</code> must be of size <code>numColumns()</code>
     204        This just does subset (but puts in correct place in y) */
     205  void transposeTimesSubset( int number,
     206                             const int * which,
     207                             const double * pi, double * y,
     208                             const double * rowScale,
     209                             const double * columnScale,
     210                             double * spare=NULL) const;
    201211    /** Return <code>x * scalar * A + y</code> in <code>z</code>.
    202212        Can use y as temporary array (will be empty at end)
     
    267277  /// Say we want special column copy
    268278  inline void makeSpecialColumnCopy()
    269   { flags_ |= 8;}
     279  { flags_ |= 16;}
    270280  /// Say we don't want special column copy
    271281  void releaseSpecialColumnCopy();
     
    275285  /// Do we want special column copy
    276286  inline bool wantsSpecialColumnCopy() const
    277   { return ((flags_&8)!=0);}
     287  { return ((flags_&16)!=0);}
    278288  /// Flags
    279289  inline int flags() const
     
    337347                                 double * COIN_RESTRICT array,
    338348                                 const double tolerance) const;
     349  /// Meat of transposeTimes by row n > K if packed - returns number nonzero
     350  int gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector,
     351                                   int * COIN_RESTRICT index,
     352                                   double * COIN_RESTRICT output,
     353                                   int numberColumns,
     354                                   const double tolerance,
     355                                   const double scalar) const;
    339356  /// Meat of transposeTimes by row n > 2 if packed - returns number nonzero
    340357  int gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector,
     
    370387      4 - has special row copy
    371388      8 - has special column copy
     389      16 - wants special column copy
    372390  */
    373391  int flags_;
  • trunk/Clp/src/ClpPresolve.cpp

    r1197 r1266  
    441441    const bool doubleton = doDoubleton();
    442442    const bool tripleton = doTripleton();
     443#define NO_FORCING
     444#ifndef NO_FORCING
    443445    const bool forcing = doForcing();
     446#endif
    444447    const bool ifree = doImpliedFree();
    445448    const bool zerocost = doTighten();
     
    477480#endif
    478481    if (dupcol) {
    479       //paction_ = dupcol_action::presolve(prob, paction_);
     482      // maybe allow integer columns to be checked
     483      if ((presolveActions_&512)!=0)
     484        prob->setPresolveOptions(prob->presolveOptions()|1);
     485      paction_ = dupcol_action::presolve(prob, paction_);
     486    }
     487    if (duprow) {
     488      paction_ = duprow_action::presolve(prob, paction_);
     489    }
     490    if (doGubrow()) {
     491      paction_ = gubrow_action::presolve(prob, paction_);
    480492    }
    481493
     
    494506      const CoinPresolveAction * const paction0 = paction_;
    495507      // look for substitutions with no fill
    496       int fill_level=2;
    497       //fill_level=10;
    498       //printf("** fill_level == 10 !\n");
     508      int fill_level=3;
     509#define IMPLIED 3
     510#define IMPLIED2 99
     511#if IMPLIED!=3
     512#if IMPLIED>2&&IMPLIED<11
     513      fill_level=IMPLIED;
     514      printf("** fill_level == %d !\n",fill_level);
     515#endif
     516#if IMPLIED>11&&IMPLIED<21
     517      fill_level=-(IMPLIED-10);
     518      printf("** fill_level == %d !\n",fill_level);
     519#endif
     520#endif
    499521      int whichPass=0;
    500522      while (1) {
     
    534556            break;
    535557        }
    536 
     558#ifndef NO_FORCING
    537559        if (forcing) {
    538560          paction_ = forcing_constraint_action::presolve(prob, paction_);
     
    540562            break;
    541563        }
     564#endif
    542565
    543566        if (ifree&&(whichPass%5)==1) {
    544         paction_ = implied_free_action::presolve(prob, paction_,fill_level);
     567          paction_ = implied_free_action::presolve(prob, paction_,fill_level);
    545568        if (prob->status_)
    546569          break;
     
    642665          const CoinPresolveAction * const paction2 = paction_;
    643666          if (ifree) {
    644             //int fill_level=0; // switches off substitution
    645             paction_ = implied_free_action::presolve(prob, paction_,fill_level);
     667#if IMPLIED2 ==0
     668            int fill_level=0; // switches off substitution
     669#elif IMPLIED2!=99
     670            int fill_level=IMPLIED2;
     671#endif
     672            if ((itry&1)==0)
     673              paction_ = implied_free_action::presolve(prob, paction_,fill_level);
    646674            if (prob->status_)
    647675              break;
     
    652680      } else if (ifree) {
    653681        // just free
     682#if IMPLIED2 ==0
    654683        int fill_level=0; // switches off substitution
     684#elif IMPLIED2!=99
     685        int fill_level=IMPLIED2;
     686#endif
    655687        paction_ = implied_free_action::presolve(prob, paction_,fill_level);
    656688        if (prob->status_)
     
    10551087  numberNextRowsToDo_(0),
    10561088  presolveOptions_(0)
    1057 
    10581089{
    10591090  const int bufsize = bulk0_;
     
    10901121  // same thing for row rep
    10911122  CoinPackedMatrix * mRow = new CoinPackedMatrix();
     1123  mRow->setExtraGap(0.0);
     1124  mRow->setExtraMajor(0.0);
    10921125  mRow->reverseOrderedCopyOf(*m);
    1093   mRow->removeGaps();
    1094   mRow->setExtraGap(0.0);
     1126  //mRow->removeGaps();
     1127  //mRow->setExtraGap(0.0);
    10951128
    10961129  // Now get rid of matrix
     
    12071240  mcstrt_[ncols_] = bufsize-1;
    12081241  mrstrt_[nrows_] = bufsize-1;
     1242  // Allocate useful arrays
     1243  initializeStuff();
    12091244
    12101245#if     PRESOLVE_CONSISTENCY
     
    15071542    // Do presolve
    15081543    paction_ = presolve(&prob);
     1544    // Get rid of useful arrays
     1545    prob.deleteStuff();
    15091546
    15101547    result =0;
  • trunk/Clp/src/ClpPresolve.hpp

    r1055 r1266  
    123123  inline void setDoSingletonColumn(bool doSingleton)
    124124  { if (doSingleton) presolveActions_  &= ~512; else presolveActions_ |= 512;}
     125  /// Whether we want to do gubrow part of presolve
     126  inline bool doGubrow() const
     127  { return (presolveActions_&1024)==0;}
     128  inline void setDoGubrow(bool doGubrow)
     129  { if (doGubrow) presolveActions_  &= ~1024; else presolveActions_ |= 1024;}
    125130  /// Set whole group
    126131  inline int presolveActions() const
  • trunk/Clp/src/ClpSimplex.cpp

    r1246 r1266  
    10601060      rowReducedCost_[iRow]=value;
    10611061    }
    1062     ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    10631062    // can use work if problem scaled (for better cache)
    1064     if (numberRows_>10000)
    1065       matrix_->transposeTimes(-1.0,dual_,reducedCostWork_,
    1066                               rowScale_,columnScale_,work);
    1067     else
    1068       matrix_->transposeTimes(-1.0,dual_,reducedCostWork_,
    1069                               rowScale_,columnScale_,NULL);
     1063    ClpPackedMatrix* clpMatrix =
     1064      dynamic_cast< ClpPackedMatrix*>(matrix_);
     1065    if (clpMatrix&&rowScale_&&(clpMatrix->flags()&2)==0) {
     1066      CoinIndexedVector * cVector = columnArray_[0];
     1067      int * whichColumn = cVector->getIndices();
     1068      assert (!cVector->getNumElements());
     1069      int n=0;
     1070      for (int i=0;i<numberColumns_;i++) {
     1071        if (getColumnStatus(i)!=basic) {
     1072          whichColumn[n++]=i;
     1073          reducedCostWork_[i]=objectiveWork_[i];
     1074        } else {
     1075          reducedCostWork_[i]=0.0;
     1076        }
     1077      }
     1078      if (numberRows_>4000)
     1079        clpMatrix->transposeTimesSubset(n,whichColumn,dual_,reducedCostWork_,
     1080                                rowScale_,columnScale_,work);
     1081      else
     1082        clpMatrix->transposeTimesSubset(n,whichColumn,dual_,reducedCostWork_,
     1083                                rowScale_,columnScale_,NULL);
     1084    } else {
     1085      ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
     1086      if (numberRows_>4000)
     1087        matrix_->transposeTimes(-1.0,dual_,reducedCostWork_,
     1088                                rowScale_,columnScale_,work);
     1089      else
     1090        matrix_->transposeTimes(-1.0,dual_,reducedCostWork_,
     1091                                rowScale_,columnScale_,NULL);
     1092    }
    10701093    ClpFillN(work,numberRows_,0.0);
    10711094    // Extended duals and check dual infeasibility
     
    16731696    handler_->message()<<CoinMessageEol;
    16741697  }
    1675   //#define COMPUTE_INT_INFEAS
     1698#ifdef COIN_FACTORIZATION_INFO
     1699#define COMPUTE_INT_INFEAS
     1700#endif
    16761701#ifdef COMPUTE_INT_INFEAS
    16771702  if (userPointer_) {
     
    27182743  }
    27192744}
     2745//static int x_gaps[4]={0,0,0,0};
    27202746//static int scale_times[]={0,0,0,0};
    27212747bool
     
    28652891    if (oldMatrix)
    28662892      checkType = 14;
    2867     if ((specialOptions_&COIN_CBC_USING_CLP)!=0)
     2893    bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
     2894    if (inCbcOrOther)
    28682895      checkType -= 4; // don't check for duplicates
    28692896    if (!matrix_->allElementsInRange(this,smallElement_,1.0e20,checkType)) {
     
    30613088  }
    30623089  if (what==63) {
     3090#if 0
     3091    {
     3092      x_gaps[0]++;
     3093      ClpPackedMatrix* clpMatrix =
     3094        dynamic_cast< ClpPackedMatrix*>(matrix_);
     3095      if (clpMatrix) {
     3096        if (!clpMatrix->getPackedMatrix()->hasGaps())
     3097          x_gaps[1]++;
     3098        if ((clpMatrix->flags()&2)==0)
     3099          x_gaps[3]++;
     3100      } else {
     3101        x_gaps[2]++;
     3102      }
     3103      if ((x_gaps[0]%1000)==0)
     3104        printf("create %d times, no gaps %d times - not clp %d times - flagged %d\n",
     3105               x_gaps[0],x_gaps[1],x_gaps[2],x_gaps[3]);
     3106    }
     3107#endif
    30633108    if (newArrays&&(specialOptions_&65536)==0) {
    30643109      delete [] cost_;
     
    49114956  // Get a row copy in standard format
    49124957  CoinPackedMatrix copy;
     4958  copy.setExtraGap(0.0);
     4959  copy.setExtraMajor(0.0);
    49134960  copy.reverseOrderedCopyOf(*matrix());
    49144961  // Matrix may have been created so get rid of it
     
    73347381        // Get a row copy in standard format
    73357382        CoinPackedMatrix copy;
     7383        copy.setExtraGap(0.0);
     7384        copy.setExtraMajor(0.0);
    73367385        copy.reverseOrderedCopyOf(*columnCopy);
    73377386        // get matrix data pointers
     
    81928241ClpSimplex::restoreData(ClpDataSave saved)
    81938242{
    8194   factorization_->sparseThreshold(saved.sparseThreshold_);
     8243  //factorization_->sparseThreshold(saved.sparseThreshold_);
    81958244  factorization_->pivotTolerance(saved.pivotTolerance_);
    81968245  perturbation_ = saved.perturbation_;
  • trunk/Clp/src/ClpSimplex.hpp

    r1264 r1266  
    925925      4 bit - keep arrays like upper_ around
    926926      8 bit - if factorization kept can still declare optimal at once
     927      16 bit - if checking replaceColumn accuracy before updating
    927928  */
    928929  inline int moreSpecialOptions() const
  • trunk/Clp/src/ClpSimplexDual.cpp

    r1243 r1266  
    295295    numberFake_ =0; // Number of variables at fake bounds
    296296    numberChanged_ =0; // Number of variables with changed costs
    297     changeBounds(true,NULL,objectiveChange);
     297    changeBounds(1,NULL,objectiveChange);
    298298   
    299299    if (!ifValuesPass) {
     
    357357  finish(startFinishOptions);
    358358}
     359#ifdef CLP_INVESTIGATE
     360static int z_reason[7]={0,0,0,0,0,0,0};
     361static int z_thinks=-1;
     362#endif
    359363void
    360364ClpSimplexDual::gutsOfDual(int ifValuesPass,double * & saveDuals,int initialStatus,
    361365                           ClpDataSave & data)
    362366{
     367#ifdef CLP_INVESTIGATE
     368  z_reason[0]++;
     369  z_thinks=-1;
     370  int nPivots=9999;
     371#endif
    363372  int lastCleaned=0; // last time objective or bounds cleaned up
    364373 
     
    487496    // Do iterations
    488497    int returnCode=whileIterating(saveDuals,ifValuesPass);
     498#ifdef CLP_INVESTIGATE
     499    nPivots=factorization_->pivots();
     500#endif
    489501    if (returnCode==-2)
    490502      factorType=3;
    491503  }
     504#ifdef CLP_INVESTIGATE
     505  if (z_thinks!=-1) {
     506    assert (z_thinks<4);
     507    if ((!factorization_->pivots()&&nPivots<20)&&z_thinks>=0&&z_thinks<2)
     508      z_thinks += 4;
     509    z_reason[1+z_thinks]++;
     510  }
     511  if ((z_reason[0]%1000)==0) {
     512    printf("Reason");
     513    for (int i=0;i<7;i++)
     514      printf(" %d",z_reason[i]);
     515    printf("\n");
     516  }
     517#endif
    492518}
    493519int
     
    495521{
    496522  algorithm_ = -1;
    497  
     523  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
    498524  // save data
    499525  ClpDataSave data = saveData();
     
    618644    numberFake_ =0; // Number of variables at fake bounds
    619645    numberChanged_ =0; // Number of variables with changed costs
    620     changeBounds(true,NULL,objectiveChange);
     646    changeBounds(1,NULL,objectiveChange);
    621647   
    622648    int lastCleaned=0; // last time objective or bounds cleaned up
     
    772798ClpSimplexDual::whileIterating(double * & givenDuals,int ifValuesPass)
    773799{
     800#ifdef CLP_INVESTIGATE
     801  z_thinks=-1;
     802#endif
    774803#if 0
    775804  if (!numberIterations_&&auxiliaryModel_) {
     
    13471376                                                         rowArray_[1],
    13481377                                                         pivotRow_,
    1349                                                          alpha_);
     1378                                                         alpha_,
     1379                                                     (moreSpecialOptions_&16)!=0);
    13501380        // if no pivots, bad update but reasonable alpha - take and invert
    13511381        if (updateStatus==2&&
     
    13621392          dualRowPivot_->unrollWeights();
    13631393          // later we may need to unwind more e.g. fake bounds
    1364           if (factorization_->pivots()) {
     1394          if (factorization_->pivots()&&
     1395              ((moreSpecialOptions_&16)==0||factorization_->pivots()>4)) {
    13651396            problemStatus_=-2; // factorize now
    13661397            returnCode=-2;
     1398            moreSpecialOptions_ |= 16;
    13671399            break;
    13681400          } else {
     
    15471579        }
    15481580      } else {
     1581#ifdef CLP_INVESTIGATE
     1582        z_thinks=1;
     1583#endif
    15491584        // no incoming column is valid
    15501585        pivotRow_=-1;
     
    15531588          printf("** no column pivot\n");
    15541589#endif
    1555         if (factorization_->pivots()<=CoinMax(dontFactorizePivots_,5)&&acceptablePivot_<=1.0e-8) {
     1590        if (factorization_->pivots()<2&&acceptablePivot_<=1.0e-8) {
    15561591          //&&goodAccuracy()) {
    15571592          // If not in branch and bound etc save ray
     
    15651600          // If we have just factorized and infeasibility reasonable say infeas
    15661601          if (((specialOptions_&4096)!=0||bestPossiblePivot<1.0e-11)&&dualBound_>1.0e8) {
    1567             if (valueOut_>upperOut_+1.0e-3||valueOut_<lowerOut_-1.0e-3
     1602            if (valueOut_>upperOut_+1.0e-4||valueOut_<lowerOut_-1.0e-4
    15681603                || (specialOptions_&64)==0) {
    15691604              // say infeasible
     
    16611696      }
    16621697    } else {
     1698#ifdef CLP_INVESTIGATE
     1699      z_thinks=0;
     1700#endif
    16631701      // no pivot row
    16641702#ifdef CLP_DEBUG
     
    16971735      returnCode=0;
    16981736      if (numberPivots<=CoinMax(dontFactorizePivots_,20)&&
    1699           (specialOptions_&2048)!=0&&!numberChanged_
    1700           &&dualBound_>1.0e8) {
     1737          (specialOptions_&2048)!=0&&(true||!numberChanged_||perturbation_==101)
     1738          &&dualBound_>=1.0e8) {
    17011739        specialCase=true;
    17021740        // as dual bound high - should be okay
     
    17391777            problemStatus_=-5;
    17401778          } else {
    1741             if (numberPivots) {
     1779            problemStatus_=0;
     1780            // May be perturbed
     1781            if (perturbation_==101||numberChanged_) {
     1782              perturbation_=102; // stop any perturbations
     1783              //double changeCost;
     1784              //changeBounds(1,NULL,changeCost);
     1785              createRim4(false);
     1786              // make sure duals are current
     1787              computeDuals(givenDuals);
     1788              checkDualSolution();
     1789              if (numberDualInfeasibilities_) {
     1790                problemStatus_=-3;
     1791              } else {
     1792                computeObjectiveValue(true);
     1793              }
     1794            } else if (numberPivots) {
     1795              computeObjectiveValue(true);
     1796            }
     1797            if (numberPivots<-1000) {
    17421798              // objective may be wrong
    17431799              objectiveValue_ = innerProduct(cost_,numberColumns_+numberRows_,solution_);
     
    17751831              }
    17761832            }
    1777             problemStatus_=0;
    17781833            sumPrimalInfeasibilities_=0.0;
    1779             if ((specialOptions_&(1024+16384))!=0) {
     1834            if ((specialOptions_&(1024+16384))!=0&&!problemStatus_) {
    17801835              CoinIndexedVector * arrayVector = rowArray_[1];
    17811836              arrayVector->clear();
     
    24112466// and cost of change vector
    24122467int
    2413 ClpSimplexDual::changeBounds(bool initialize,
     2468ClpSimplexDual::changeBounds(int initialize,
    24142469                                 CoinIndexedVector * outputArray,
    24152470                                 double & changeCost)
     
    25082563    }
    25092564    return numberInfeasibilities;
    2510   } else {
     2565  } else if (initialize==1) {
    25112566    int iSequence;
    25122567     
     
    25592614
    25602615    return 1;
     2616  } else {
     2617    // just reset changed ones
     2618    if (columnScale_) {
     2619      int iSequence;
     2620      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
     2621        FakeBound fakeStatus = getFakeBound(iSequence);
     2622        if (fakeStatus!=noFake) {
     2623          if (((int) fakeStatus&1)!=0) {
     2624            // lower
     2625            double value = columnLower_[iSequence];
     2626            if (value>-1.0e30) {
     2627              double multiplier = rhsScale_/columnScale_[iSequence];
     2628              value *= multiplier;
     2629            }
     2630            columnLowerWork_[iSequence]=value;
     2631          }
     2632          if (((int) fakeStatus&2)!=0) {
     2633            // upper
     2634            double value = columnUpper_[iSequence];
     2635            if (value<1.0e30) {
     2636              double multiplier = rhsScale_/columnScale_[iSequence];
     2637              value *= multiplier;
     2638            }
     2639            columnUpperWork_[iSequence]=value;
     2640          }
     2641        }
     2642      }
     2643      for (iSequence=0;iSequence<numberRows_;iSequence++) {
     2644        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
     2645        if (fakeStatus!=noFake) {
     2646          if (((int) fakeStatus&1)!=0) {
     2647            // lower
     2648            double value = rowLower_[iSequence];
     2649            if (value>-1.0e30) {
     2650              double multiplier = rhsScale_*rowScale_[iSequence];
     2651              value *= multiplier;
     2652            }
     2653            rowLowerWork_[iSequence]=value;
     2654          }
     2655          if (((int) fakeStatus&2)!=0) {
     2656            // upper
     2657            double value = rowUpper_[iSequence];
     2658            if (value<1.0e30) {
     2659              double multiplier = rhsScale_*rowScale_[iSequence];
     2660              value *= multiplier;
     2661            }
     2662            rowUpperWork_[iSequence]=value;
     2663          }
     2664        }
     2665      }
     2666    } else {
     2667      int iSequence;
     2668      for (iSequence=0;iSequence<numberColumns_;iSequence++) {
     2669        FakeBound fakeStatus = getFakeBound(iSequence);
     2670        if (((int) fakeStatus&1)!=0) {
     2671          // lower
     2672          columnLowerWork_[iSequence]=columnLower_[iSequence];
     2673        }
     2674        if (((int) fakeStatus&2)!=0) {
     2675          // upper
     2676          columnUpperWork_[iSequence]=columnUpper_[iSequence];
     2677        }
     2678      }
     2679      for (iSequence=0;iSequence<numberRows_;iSequence++) {
     2680        FakeBound fakeStatus = getFakeBound(iSequence+numberColumns_);
     2681        if (((int) fakeStatus&1)!=0) {
     2682          // lower
     2683          rowLowerWork_[iSequence]=rowLower_[iSequence];
     2684        }
     2685        if (((int) fakeStatus&2)!=0) {
     2686          // upper
     2687          rowUpperWork_[iSequence]=rowUpper_[iSequence];
     2688        }
     2689      }
     2690    }
     2691    return 0;
    25612692  }
    25622693}
     
    34323563                                      int ifValuesPass)
    34333564{
     3565#ifdef CLP_INVESTIGATE
     3566  if (z_thinks>0&&z_thinks<2)
     3567    z_thinks+=2;
     3568#endif
    34343569  // If lots of iterations then adjust costs if large ones
    34353570  if (numberIterations_>4*(numberRows_+numberColumns_)&&objectiveScale_==1.0) {
     
    38674002      progress_.modifyObjective(objectiveValue_
    38684003                               -sumDualInfeasibilities_*dualBound_);
     4004      // see if cutoff reached
     4005      double limit = 0.0;
     4006      getDblParam(ClpDualObjectiveLimit, limit);
     4007      if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
     4008         limit+1.0e-7+1.0e-8*fabs(limit)&&!numberAtFakeBound()) {
     4009        //looks infeasible on objective
     4010        if (perturbation_==101) {
     4011          perturbation_=102; // stop any perturbations
     4012          cleanDuals=1;
     4013          // make sure fake bounds are back
     4014          changeBounds(1,NULL,changeCost);
     4015          createRim4(false);
     4016          // make sure duals are current
     4017          computeDuals(givenDuals);
     4018          checkDualSolution();
     4019        }
     4020      }
    38694021      if (primalFeasible()&&!givenDuals) {
    38704022        // may be optimal - or may be bounds are wrong
     
    38754027        double * saveColumnSolution = NULL;
    38764028        double * saveRowSolution = NULL;
    3877         bool inCbc = (specialOptions_&0x01000000)!=0;
     4029        bool inCbc = (specialOptions_&(0x01000000|16384))!=0;
    38784030        if (!inCbc) {
    38794031          saveColumnSolution = CoinCopyOfArray(columnActivityWork_,numberColumns_);
    38804032          saveRowSolution = CoinCopyOfArray(rowActivityWork_,numberRows_);
    38814033        }
    3882         numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
     4034        numberChangedBounds=changeBounds(0,rowArray_[3],changeCost);
    38834035        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
    38844036          //looks optimal - do we need to reset tolerance
     
    38884040            // make sure fake bounds are back
    38894041            //computeObjectiveValue();
    3890             changeBounds(true,NULL,changeCost);
     4042            changeBounds(1,NULL,changeCost);
    38914043            //computeObjectiveValue();
    38924044            createRim4(false);
     
    38944046            computeDuals(givenDuals);
    38954047            checkDualSolution();
    3896             //if (numberDualInfeasibilities_)
     4048#define DUAL_TRY_FASTER
     4049#ifdef DUAL_TRY_FASTER
     4050            if (numberDualInfeasibilities_) {
     4051#endif
    38974052              numberChanged_=1; // force something to happen
    3898             //else
    3899             //computeObjectiveValue();
     4053#ifdef DUAL_TRY_FASTER
     4054            } else {
     4055              //double value = objectiveValue_;
     4056              computeObjectiveValue(true);
     4057              //printf("old %g new %g\n",value,objectiveValue_);
     4058              //numberChanged_=1;
     4059            }
     4060#endif
    39004061          }
    39014062          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4&&
     
    40004161      if (problemStatus_==-4||problemStatus_==-5) {
    40014162        // may be infeasible - or may be bounds are wrong
    4002         numberChangedBounds=changeBounds(false,NULL,changeCost);
     4163        numberChangedBounds=changeBounds(0,NULL,changeCost);
    40034164        /* Should this be here as makes no difference to being feasible.
    40044165           But seems to make a difference to run times. */
     
    40084169          numberChangedBounds=1;
    40094170          // make sure fake bounds are back
    4010           changeBounds(true,NULL,changeCost);
     4171          changeBounds(1,NULL,changeCost);
    40114172          createRim4(false);
    40124173        }
     
    40924253          */
    40934254          int saveNumberFake = numberFake_;
    4094           changeBounds(true,NULL,changeCost);
     4255          changeBounds(1,NULL,changeCost);
    40954256          numberFake_ += saveNumberFake;
    40964257          cleanDuals=2;
     
    42054366    }
    42064367  }
    4207   if (type==0||type==1) {
    4208     if (!type) {
    4209       // create save arrays
    4210       delete [] saveStatus_;
    4211       delete [] savedSolution_;
    4212       saveStatus_ = new unsigned char [numberRows_+numberColumns_];
    4213       savedSolution_ = new double [numberRows_+numberColumns_];
    4214     }
    4215     // save arrays
    4216     CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
    4217     CoinMemcpyN(rowActivityWork_,
    4218            numberRows_,savedSolution_+numberColumns_);
    4219     CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
    4220     // save extra stuff
    4221     int dummy;
    4222     matrix_->generalExpanded(this,5,dummy);
    4223   }
    4224   if (weightsSaved) {
    4225     // restore weights (if saved) - also recompute infeasibility list
    4226     if (tentativeStatus>-3)
    4227       dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
    4228     else
    4229       dualRowPivot_->saveWeights(this,3);
    4230   }
    42314368  // unflag all variables (we may want to wait a bit?)
    42324369  if ((tentativeStatus!=-2&&tentativeStatus!=-1)&&unflagVariables) {
     
    42664403#endif
    42674404      }
     4405    }
     4406  }
     4407  if (problemStatus_<0) {
     4408    if (type==0||type==1) {
     4409      if (!type) {
     4410        // create save arrays
     4411        delete [] saveStatus_;
     4412        delete [] savedSolution_;
     4413        saveStatus_ = new unsigned char [numberRows_+numberColumns_];
     4414        savedSolution_ = new double [numberRows_+numberColumns_];
     4415      }
     4416      // save arrays
     4417      CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
     4418      CoinMemcpyN(rowActivityWork_,
     4419                  numberRows_,savedSolution_+numberColumns_);
     4420      CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
     4421      // save extra stuff
     4422      int dummy;
     4423      matrix_->generalExpanded(this,5,dummy);
     4424    }
     4425    if (weightsSaved) {
     4426      // restore weights (if saved) - also recompute infeasibility list
     4427      if (tentativeStatus>-3)
     4428        dualRowPivot_->saveWeights(this,(type <2) ? 2 : 4);
     4429      else
     4430        dualRowPivot_->saveWeights(this,3);
    42684431    }
    42694432  }
     
    43554518  }
    43564519  lastGoodIteration_ = numberIterations_;
     4520  if (numberIterations_>lastBadIteration_+100)
     4521    moreSpecialOptions_ &= ~16; // clear check accuracy flag
    43574522  if (problemStatus_<0) {
    43584523    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
     
    46014766  }
    46024767  int iColumn;
     4768  const int * columnLength = matrix_->getVectorLengths();
    46034769  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    46044770    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
    4605       int length = matrix_->getVectorLength(iColumn);
     4771      int length = columnLength[iColumn];
    46064772      if (length>2) {
    46074773        maxLength = CoinMax(maxLength,length);
     
    46174783  }
    46184784  bool uniformChange=false;
     4785  bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
    46194786  if (perturbation_>50) {
    46204787    // Experiment
    46214788    // maximumFraction could be 1.0e-10 to 1.0
    46224789    double m[]={1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0};
    4623     maximumFraction = m[CoinMin(perturbation_-51,10)];
     4790    int whichOne = perturbation_-51;
     4791    //if (inCbcOrOther&&whichOne>0)
     4792    //whichOne--;
     4793    maximumFraction = m[CoinMin(whichOne,10)];
     4794  } else if (inCbcOrOther) {
     4795    //maximumFraction = 1.0e-6;
    46244796  }
    46254797  int iRow;
     
    48114983      }
    48124984      if (value) {
    4813         int length = matrix_->getVectorLength(iColumn);
     4985        int length = columnLength[iColumn];
    48144986        if (length>3) {
    48154987          length = (int) ((double) length * factor);
     
    52055377  numberFake_ =0; // Number of variables at fake bounds
    52065378  numberChanged_ =0; // Number of variables with changed costs
    5207   changeBounds(true,NULL,objectiveChange);
     5379  changeBounds(1,NULL,objectiveChange);
    52085380
    52095381  problemStatus_ = -1;
     
    53305502  // Say not in fast dual
    53315503  specialOptions_ &= ~16384;
    5332   assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>1.0e8)
     5504  assert(!numberFake_||((specialOptions_&(2048|4096))!=0&&dualBound_>=1.0e8)
    53335505         ||returnCode||problemStatus_); // all bounds should be okay
     5506  if (numberFake_>0) {
     5507    // Set back
     5508    double dummy;
     5509    changeBounds(2,NULL,dummy);
     5510  }
    53345511  // Restore any saved stuff
    53355512  restoreData(data);
     
    59026079  // Get a row copy in standard format
    59036080  CoinPackedMatrix copy;
     6081  copy.setExtraGap(0.0);
     6082  copy.setExtraMajor(0.0);
    59046083  copy.reverseOrderedCopyOf(*columnCopy);
    59056084  // get matrix data pointers
     
    60636242  createRim1(false);
    60646243  double dummyChangeCost=0.0;
    6065   changeBounds(true,rowArray_[2],dummyChangeCost);
     6244  changeBounds(1,rowArray_[2],dummyChangeCost);
    60666245  // throw away change
    60676246  for (int i=0;i<4;i++)
  • trunk/Clp/src/ClpSimplexDual.hpp

    r1243 r1266  
    242242      Fills in changeVector which can be used to see if unbounded
    243243      and cost of change vector
    244   */
    245   int changeBounds(bool initialize,CoinIndexedVector * outputArray,
     244      If 2 sets to original (just changed)
     245  */
     246  int changeBounds(int initialize,CoinIndexedVector * outputArray,
    246247                   double & changeCost);
    247248  /** As changeBounds but just changes new bounds for a single variable.
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r1218 r1266  
    180180
    181181  algorithm_ = +1;
    182   //specialOptions_ |= 4;
     182  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
    183183
    184184  // save data
     
    656656        problemStatus_=-2; //
    657657      } else if (returnCode==-5) {
    658         // something flagged - continue;
     658        if ((moreSpecialOptions_&16)==0&&factorization_->pivots()) {
     659          moreSpecialOptions_ |= 16;
     660          problemStatus_=-2;
     661        }
     662        // otherwise something flagged - continue;
    659663      } else if (returnCode==2) {
    660664        problemStatus_=-5; // looks unbounded
     
    13871391  }
    13881392  lastGoodIteration_ = numberIterations_;
     1393  if (numberIterations_>lastBadIteration_+100)
     1394    moreSpecialOptions_ &= ~16; // clear check accuracy flag
    13891395  if (goToDual)
    13901396    problemStatus_=10; // try dual
     
    27182724                                                     rowArray_[1],
    27192725                                                     pivotRow_,
    2720                                                      alpha_);
     2726                                                     alpha_,
     2727                                                     (moreSpecialOptions_&16)!=0);
    27212728
    27222729      // if no pivots, bad update but reasonable alpha - take and invert
  • trunk/Clp/src/ClpSolve.cpp

    r1264 r1266  
    465465      // switch on singletons to slacks
    466466      pinfo.setDoSingletonColumn(true);
     467      // gub stuff for testing
     468      pinfo.setDoGubrow(true);
    467469    }
    468470#ifndef CLP_NO_STD
     
    517519  int doSlp=0;
    518520  int primalStartup=1;
     521  bool tryItSave = false;
    519522  // switch to primal from automatic if just one cost entry
    520523  if (method==ClpSolve::automatic&&model2->numberColumns()>5000&&
    521524      (specialOptions_&1024)!=0) {
    522525    int numberColumns = model2->numberColumns();
     526    int numberRows = model2->numberRows();
    523527    const double * obj = model2->objective();
    524528    int nNon=0;
     
    532536#endif
    533537      method=ClpSolve::usePrimal;
     538      tryItSave= numberRows>200&&numberColumns>2000&&
     539        (numberColumns>2*numberRows||(specialOptions_&1024)!=0);
    534540    }
    535541  }
     
    742748    model2->setFactorizationFrequency(this->factorizationFrequency());
    743749  }
    744   bool tryItSave = false;
    745750  if (method==ClpSolve::automatic) {
    746751    if (doSprint==0&&doIdiot==0) {
     
    811816        }
    812817        bool tryIt= numberRows>200&&numberColumns>2000&&
    813           (numberColumns>2*numberRows||(method==ClpSolve::automatic&&(specialOptions_&1024)!=0));
     818          (numberColumns>2*numberRows||(method!=ClpSolve::useDual&&(specialOptions_&1024)!=0));
    814819        tryItSave = tryIt;
    815820        if (numberRows<1000&&numberColumns<3000)
  • trunk/Clp/src/Idiot.cpp

    r1264 r1266  
    11601160  }
    11611161  double fixTolerance=IDIOT_FIX_TOLERANCE;
    1162   // double startTime = CoinCpuTime();          unused?
     1162#ifdef COIN_DEVELOP
     1163  double startTime = CoinCpuTime();             
     1164#endif
    11631165  ClpSimplex * saveModel=NULL;
    11641166  ClpMatrixBase * matrix = model_->clpMatrix();
Note: See TracChangeset for help on using the changeset viewer.