Changeset 200 for branches/pre


Ignore:
Timestamp:
Aug 22, 2003 6:31:48 AM (18 years ago)
Author:
forrest
Message:

quadratic stuff

Location:
branches/pre
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • branches/pre/ClpLinearObjective.cpp

    r180 r200  
    110110// Returns gradient
    111111double * 
    112 ClpLinearObjective::gradient(const double * solution, double & offset)
     112ClpLinearObjective::gradient(const double * solution, double & offset,bool refresh)
    113113{
    114114  offset=0.0;

  • branches/pre/ClpModel.cpp

    r199 r200  
    11941194  assert ((dynamic_cast< ClpLinearObjective*>(objective_)));
    11951195  double offset;
    1196   ClpObjective * obj = new ClpQuadraticObjective(objective_->gradient(NULL,offset),numberColumns,
     1196  ClpObjective * obj = new ClpQuadraticObjective(objective_->gradient(NULL,offset,false),numberColumns,
    11971197                                             start,column,element);
    11981198  delete objective_;
    … …  
    12071207  double offset;
    12081208  ClpQuadraticObjective * obj =
    1209     new ClpQuadraticObjective(objective_->gradient(NULL,offset),numberColumns_,
     1209    new ClpQuadraticObjective(objective_->gradient(NULL,offset,false),numberColumns_,
    12101210                                                 NULL,NULL,NULL);
    12111211  delete objective_;

  • branches/pre/ClpPackedMatrix.cpp

    r199 r200  
    12591259      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    12601260      int numberXColumns = quadratic->getNumCols();
    1261       int numberXRows = numberRows-numberXColumns;
    12621261      if (numberXColumns<numberColumns) {
     1262        // we assume symmetric
     1263        int numberQuadraticColumns=0;
    12631264        int i;
    1264         for (i=0;i<numberXColumns;i++)
    1265           rowScale[i+numberXRows] = columnScale[i];
    1266         for (i=0;i<numberXRows;i++)
    1267           columnScale[i+numberXColumns] = rowScale[i];
     1265        //const int * columnQuadratic = quadratic->getIndices();
     1266        //const int * columnQuadraticStart = quadratic->getVectorStarts();
     1267        const int * columnQuadraticLength = quadratic->getVectorLengths();
     1268        for (i=0;i<numberXColumns;i++) {
     1269          int length=columnQuadraticLength[i];
     1270#ifndef CORRECT_COLUMN_COUNTS
     1271          length=1;
     1272#endif
     1273          if (length)
     1274            numberQuadraticColumns++;
     1275        }
     1276        int numberXRows = numberRows-numberQuadraticColumns;
     1277        numberQuadraticColumns=0;
     1278        for (i=0;i<numberXColumns;i++) {
     1279          int length=columnQuadraticLength[i];
     1280#ifndef CORRECT_COLUMN_COUNTS
     1281          length=1;
     1282#endif
     1283          if (length) {
     1284            rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
     1285            numberQuadraticColumns++;
     1286          }
     1287        }   
     1288        int numberQuadraticRows=0;
     1289        for (i=0;i<numberXRows;i++) {
     1290          // See if any in row quadratic
     1291          int j;
     1292          int numberQ=0;
     1293          for (j=rowStart[i];j<rowStart[i+1];j++) {
     1294            int iColumn = column[j];
     1295            if (columnQuadraticLength[iColumn])
     1296              numberQ++;
     1297          }
     1298#ifndef CORRECT_ROW_COUNTS
     1299          numberQ=1;
     1300#endif
     1301          if (numberQ) {
     1302            columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
     1303            numberQuadraticRows++;
     1304          }
     1305        }
    12681306        // and make sure Sj okay
    1269         for (iColumn=numberRows;iColumn<numberColumns;iColumn++) {
     1307        for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
    12701308          CoinBigIndex j=columnStart[iColumn];
    12711309          assert(columnLength[iColumn]==1);
    … …  
    13301368  }
    13311369}
    1332 /* Unpacks a column into an CoinIndexedvector
    1333 ** in packed foramt
     1370/* Unpacks a column into a CoinIndexedVector
     1371** in packed format
    13341372Note that model is NOT const.  Bounds and objective could
    13351373be modified if doing column generation (just for this variable) */

  • branches/pre/ClpQuadraticObjective.cpp

    r199 r200  
    172172// Returns gradient
    173173double * 
    174 ClpQuadraticObjective::gradient(const double * solution, double & offset)
     174ClpQuadraticObjective::gradient(const double * solution, double & offset,bool refresh)
    175175{
    176176  offset=0.0;
    … …  
    178178    return objective_;
    179179  } else {
    180     if (!gradient_)
    181       gradient_ = new double[numberExtendedColumns_];
    182     const int * columnQuadratic = quadraticObjective_->getIndices();
    183     const int * columnQuadraticStart = quadraticObjective_->getVectorStarts();
    184     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
    185     const double * quadraticElement = quadraticObjective_->getElements();
    186     double offset=0.0;
    187     memcpy(gradient_,objective_,numberExtendedColumns_*sizeof(double));
    188     int iColumn;
    189     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    190       double valueI = solution[iColumn];
    191       int j;
    192       for (j=columnQuadraticStart[iColumn];
    193            j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    194         int jColumn = columnQuadratic[j];
    195         double valueJ = solution[jColumn];
    196         double elementValue = quadraticElement[j];
    197         elementValue *= 0.5;
    198         offset += valueI*valueJ*elementValue;
    199         double gradientI = valueJ*elementValue;
    200         double gradientJ = valueI*elementValue;
    201         offset -= gradientI*valueI;
    202         gradient_[iColumn] += gradientI;
    203         offset -= gradientJ*valueJ;
    204         gradient_[jColumn] += gradientJ;
     180    if (refresh||!gradient_) {
     181      if (!gradient_)
     182        gradient_ = new double[numberExtendedColumns_];
     183      const int * columnQuadratic = quadraticObjective_->getIndices();
     184      const int * columnQuadraticStart = quadraticObjective_->getVectorStarts();
     185      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
     186      const double * quadraticElement = quadraticObjective_->getElements();
     187      double offset=0.0;
     188      memcpy(gradient_,objective_,numberExtendedColumns_*sizeof(double));
     189      int iColumn;
     190      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     191        double valueI = solution[iColumn];
     192        int j;
     193        for (j=columnQuadraticStart[iColumn];
     194             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     195          int jColumn = columnQuadratic[j];
     196          double valueJ = solution[jColumn];
     197          double elementValue = quadraticElement[j];
     198          elementValue *= 0.5;
     199          offset += valueI*valueJ*elementValue;
     200          double gradientI = valueJ*elementValue;
     201          double gradientJ = valueI*elementValue;
     202          offset -= gradientI*valueI;
     203          gradient_[iColumn] += gradientI;
     204          offset -= gradientJ*valueJ;
     205          gradient_[jColumn] += gradientJ;
     206        }
    205207      }
    206208    }

  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r199 r200  
    490490  ClpQuadraticInfo info;
    491491  ClpSimplexPrimalQuadratic * model2 = makeQuadratic(info);
     492  if (!model2) {
     493    printf("** Not quadratic\n");
     494    primal(1);
     495    return 0;
     496  }
    492497#if 0
    493498  CoinMpsIO writer;
    … …  
    552557      const double * elementByColumn = quadratic->getElements();
    553558      double direction = optimizationDirection_;
     559      const int * which = info->quadraticSequence();
    554560      // direction is actually scale out not scale in
    555561      if (direction)
    … …  
    565571        for (j=0;j<number;j++) {
    566572          int iColumn2 = columnsInThisColumn[j];
     573          iColumn2 = which[iColumn2];
     574          assert (iColumn2>=0);
    567575          newElement[j] = elementsInThisColumn[j]*scale*rowScale_[iColumn2+numberXRows];
    568576        }
    … …  
    652660        int nextSequenceIn=-1;
    653661        for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    654           if (which[iColumn]>=0) {
    655             double value = solution_[iColumn];
    656             double lower = lower_[iColumn];
    657             double upper = upper_[iColumn];
    658             if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
    659                 ||getColumnStatus(iColumn)==superBasic) {
    660               if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    661                 if (phase!=2) {
    662                   phase=2;
    663                   nextSequenceIn=iColumn;
    664                 }
     662          double value = solution_[iColumn];
     663          double lower = lower_[iColumn];
     664          double upper = upper_[iColumn];
     665          if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
     666              ||getColumnStatus(iColumn)==superBasic) {
     667            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
     668              if (phase!=2) {
     669                phase=2;
     670                nextSequenceIn=iColumn;
    665671              }
    666672            }
     673          }
     674          if (which[iColumn]>=0) {
    667675            if (getColumnStatus(iColumn)==basic) {
    668676              int iS = which[iColumn]+numberRows_;
    … …  
    677685          }
    678686        }
    679         int offset=numberXColumns-numberColumns_;
     687        const int * whichRow = info->quadraticRowSequence();
    680688        for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;iColumn++) {
    681689          double value = solution_[iColumn];
    … …  
    691699            }
    692700          }
    693           if (getColumnStatus(iColumn)==basic) {
    694             int iS = iColumn+offset;
    695             assert (getColumnStatus(iS)!=basic);
    696             if(fabs(solution_[iS])>dualTolerance_) {
    697               if (phase==0) {
    698                 phase=1;
    699                 nextSequenceIn=iS;
     701          int iRow = iColumn-numberColumns_;
     702          if (whichRow[iRow]>=0) {
     703            if (getColumnStatus(iColumn)==basic) {
     704              int iS = whichRow[iRow]+numberXColumns;;
     705              assert (getColumnStatus(iS)!=basic);
     706              if(fabs(solution_[iS])>dualTolerance_) {
     707                if (phase==0) {
     708                  phase=1;
     709                  nextSequenceIn=iS;
     710                }
    700711              }
    701712            }
    … …  
    754765  double saveObjective = objectiveValue_;
    755766  int numberXColumns = info->numberXColumns();
    756   const int * which = info->quadraticSequence();
    757   const double * pi = solution_+numberXColumns;
     767  int numberXRows = info->numberXRows();
     768  int numberQuadraticRows = info->numberQuadraticRows();
     769  int numberQuadraticColumns = info->numberQuadraticColumns();
     770  const int * whichRow = info->quadraticRowSequence();
     771  const int * whichColumn = info->quadraticSequence();
     772  const int * backRow = info->backRowSequence();
     773  const int * backColumn = info->backSequence();
     774  //const double * pi = solution_+numberXColumns;
    758775  int nextSequenceIn=info->sequenceIn();
    759776  int oldSequenceIn=nextSequenceIn;
    … …  
    824841        createDjs(info,rowArray_[3],rowArray_[2]);
    825842        dualIn_=0.0; // so no updates
     843        rowArray_[1]->clear();
    826844        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
    827845                     columnArray_[0],columnArray_[1]);
    … …  
    861879        chosen=-1;
    862880        unpack(rowArray_[1]);
     881#if 1
    863882        // compute dj in case linear
    864883        {
    865           dualIn_=cost_[sequenceIn_];
     884          info->createGradient(this);
     885          double * gradient = info->gradient();
     886          dualIn_=gradient[sequenceIn_];
    866887          int j;
    867888          const double * element = rowArray_[1]->denseVector();
    … …  
    872893            for (j=0;j<numberNonZero; j++) {
    873894              int iRow = whichRow[j];
    874               dualIn_ -= element[j]*pi[iRow];
     895              dualIn_ -= element[j]*dual_[iRow];
    875896            }
    876897          } else {
    877898            for (j=0;j<numberNonZero; j++) {
    878899              int iRow = whichRow[j];
    879               dualIn_ -= element[iRow]*pi[iRow];
     900              dualIn_ -= element[iRow]*dual_[iRow];
    880901            }
    881902          }
    882903        }
     904#else
     905        dualIn_=dj_[sequenceIn_];
     906#endif
    883907        factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
    884908        if (cleanupIteration) {
    … …  
    935959          upperIn_=upper_[sequenceIn_];
    936960          if (sequenceIn_<numberXColumns) {
    937             int jSequence = which[sequenceIn_];
     961            int jSequence = whichColumn[sequenceIn_];
    938962            if (jSequence>=0) {
    939               dualIn_ = solution_[jSequence+numberRows_];
     963              dualIn_ = solution_[jSequence+numberXColumns+numberQuadraticRows];
    940964            } else {
    941965              // need coding
    … …  
    948972            }
    949973          } else {
    950             dualIn_ = solution_[sequenceIn_-numberColumns_+numberXColumns];
     974            int iRow = sequenceIn_-numberColumns_;
     975            assert (iRow>=0);
     976            if (whichRow[iRow]>=0)
     977              dualIn_ = solution_[numberXColumns+whichRow[iRow]];
    951978          }
    952979          valueIn_=solution_[sequenceIn_];
    … …  
    957984          if (sequenceIn_<numberColumns_) {
    958985            // Set dj as value of slack
    959             const int * which = info->quadraticSequence();
    960             crucialSj = which[sequenceIn_];
     986            crucialSj = whichColumn[sequenceIn_];
    961987            if (crucialSj>=0)
    962               crucialSj += numberRows_; // sj which should go to 0.0
     988              crucialSj += numberQuadraticRows+numberXColumns; // sj which should go to 0.0
    963989          } else {
    964990            // Set dj as value of pi
    965             crucialSj = sequenceIn_-numberColumns_+numberXColumns; // pi which should go to 0.0
    966             //dualIn_=solution_[crucialSj];
     991            int iRow = whichRow[sequenceIn_-numberColumns_];
     992            if (iRow>=0)
     993              crucialSj = iRow+numberXColumns; // pi which should go to 0.0
     994            else
     995              crucialSj=-1;
    967996          }
    968997          info->setCrucialSj(crucialSj);
    … …  
    12001229          int crucialSj=info->crucialSj();
    12011230          if (sequenceOut_<numberXColumns) {
    1202             chosen =sequenceOut_ + numberRows_; // sj variable
     1231            chosen =whichColumn[sequenceOut_] + numberQuadraticRows + numberXColumns; // sj variable
    12031232          } else if (sequenceOut_>=numberColumns_) {
    12041233            // Does this mean we can change pi
    12051234            int iRow = sequenceOut_-numberColumns_;
    1206             if (iRow<numberRows_-numberXColumns) {
    1207               int iPi = iRow+numberXColumns;
     1235            if (iRow<numberXRows) {
     1236              int jRow = whichRow[iRow];
     1237              assert (jRow>=0);
     1238              int iPi = jRow+numberXColumns;
    12081239              //printf("pi for row %d is %g\n",
    12091240              //     iRow,solution_[iPi]);
    … …  
    12131244              abort();
    12141245            }
    1215           } else if (sequenceOut_<numberRows_) {
     1246          } else if (sequenceOut_<numberQuadraticRows+numberXColumns) {
    12161247            // pi out
    1217             chosen = sequenceOut_-numberXColumns+numberColumns_;
     1248            chosen = backRow[sequenceOut_-numberXColumns]+numberColumns_;
    12181249          } else {
    12191250            // sj out
    1220             chosen = sequenceOut_-numberRows_;
     1251            chosen = backColumn[sequenceOut_-numberQuadraticRows-numberXColumns];
    12211252          }
    12221253          // But double check as original incoming might have gone out
    … …  
    12251256            break; // means original has gone out
    12261257          }
     1258          // In certain circumstances need to look further to see what to come in
     1259          // Think further
     1260          if (numberXColumns>numberQuadraticColumns) {
     1261            // this should be pre-computed
     1262            int iSjRow=-1;
     1263            int iRow;
     1264            for (iRow=0;iRow<numberRows_;iRow++) {
     1265              if (pivotVariable_[iRow]==crucialSj) {
     1266                iSjRow=iRow;
     1267                break;
     1268              }
     1269            }
     1270            assert (iSjRow>=0);
     1271            double direction;
     1272            if (solution_[crucialSj]>0)
     1273              direction=-1.0;
     1274            else
     1275              direction=1.0;
     1276            rowArray_[0]->createPacked(1,&iSjRow,&direction);
     1277            factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
     1278            // put row of tableau in rowArray[0] and columnArray[0]
     1279            matrix_->transposeTimes(this,-1.0,
     1280                                    rowArray_[0],rowArray_[3],columnArray_[0]);
     1281            int k;
     1282            // all packed
     1283            // Column
     1284            int number=columnArray_[0]->getNumElements();
     1285            int * which = columnArray_[0]->getIndices();
     1286            double * element = columnArray_[0]->denseVector();
     1287            int best=-1;
     1288            double bestAlpha=1.0e-6;
     1289            if (numberIterations_==29)
     1290              printf("iteration %d coming up\n",numberIterations_);
     1291            for (k=0;k<number;k++) {
     1292              int iSequence=which[k];
     1293              double value=element[k];
     1294              // choose chosen if possible
     1295              if (iSequence==chosen) {
     1296                if (value>0.0)
     1297                  printf("chosen %d going up alpha %g, sol %g\n",iSequence,value,solution_[crucialSj]);
     1298                else
     1299                  printf("chosen %d going down alpha %g, sol %g\n",iSequence,value,solution_[crucialSj]);
     1300                value = fabs(value)*1.0e5;
     1301              } else if (iSequence<numberXColumns) {
     1302                double djValue = dj_[iSequence];
     1303                if (value>0.0)
     1304                  printf("X%d going up alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
     1305                else
     1306                  printf("X%d going down alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
     1307                ClpSimplex::Status status = getColumnStatus(iSequence);
     1308               
     1309                switch(status) {
     1310                 
     1311                case ClpSimplex::basic:
     1312                  if (fabs(value)>1.0e-3)
     1313                    abort();
     1314                  break;
     1315                case ClpSimplex::isFixed:
     1316                  value=0.0;
     1317                  break;
     1318                case ClpSimplex::isFree:
     1319                case ClpSimplex::superBasic:
     1320                  value=fabs(value);
     1321                  break;
     1322                case ClpSimplex::atUpperBound:
     1323                  value = -value;
     1324                  break;
     1325                case ClpSimplex::atLowerBound:
     1326                  break;
     1327                }
     1328              } else {
     1329                value=0.0;
     1330              }
     1331              if (value>bestAlpha) {
     1332                bestAlpha=value;
     1333                best=iSequence;
     1334              }
     1335            }
     1336            // Row
     1337            number=rowArray_[0]->getNumElements();
     1338            which = rowArray_[0]->getIndices();
     1339            element = rowArray_[0]->denseVector();
     1340            for (k=0;k<number;k++) {
     1341              int iSequence=which[k];
     1342              double value=element[k];
     1343              // choose chosen if possible
     1344              if (iSequence==chosen) {
     1345                if (value>0.0)
     1346                  printf("chosen row %d going up alpha %g, sol %g\n",iSequence,value,solution_[crucialSj]);
     1347                else
     1348                  printf("chosen row %d going down alpha %g, sol %g\n",iSequence,value,solution_[crucialSj]);
     1349                value = fabs(value)*1.0e5;
     1350              } else if (iSequence<numberXRows) {
     1351                double djValue = dj_[iSequence+numberColumns_];
     1352                if (value>0.0)
     1353                  printf("R%d going up alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
     1354                else
     1355                  printf("R%d going down alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
     1356                ClpSimplex::Status status = getRowStatus(iSequence);
     1357               
     1358                switch(status) {
     1359                 
     1360                case ClpSimplex::basic:
     1361                  if (fabs(value)>1.0e-3)
     1362                    abort();
     1363                  break;
     1364                case ClpSimplex::isFixed:
     1365                  value=0.0;
     1366                  break;
     1367                case ClpSimplex::isFree:
     1368                case ClpSimplex::superBasic:
     1369                  value=fabs(value);
     1370                  break;
     1371                case ClpSimplex::atUpperBound:
     1372                  value = -value;
     1373                  break;
     1374                case ClpSimplex::atLowerBound:
     1375                  break;
     1376                }
     1377              } else {
     1378                value=0.0;
     1379              }
     1380              if (value>bestAlpha) {
     1381                bestAlpha=value;
     1382                best=iSequence+numberColumns_;
     1383              }
     1384            }
     1385            assert (best>=0);
     1386            if (best!=chosen)
     1387              printf("** chosen changed from %d to %d\n",chosen,best);
     1388            chosen=best;
     1389          }
     1390          rowArray_[0]->clear();
     1391          rowArray_[3]->checkClear();
     1392          columnArray_[0]->clear();
    12271393          if (returnCode==-2) {
    12281394            // factorization
    … …  
    13721538      crucialSj2 += numberRows_; // sj2 which should go to 0.0
    13731539  } else if (sequenceIn_>=numberColumns_) {
    1374     crucialSj2 = sequenceIn_-numberColumns_+numberXColumns; // pi which should go to 0.0
     1540    const int * which = info->quadraticRowSequence();
     1541    int iRow=sequenceIn_-numberColumns_;
     1542    crucialSj2 = which[iRow];
     1543    if (crucialSj2>=0)
     1544      crucialSj2 += numberXColumns; // sj2 which should go to 0.0
    13751545  }
    13761546  double tj = 0.0;
    … …  
    15261696  double upperTheta = maximumMovement;
    15271697  bool throwAway=false;
    1528   if (numberIterations_==1395) {
     1698  if (numberIterations_==1453) {
    15291699    printf("Bad iteration coming up after iteration %d\n",numberIterations_);
    15301700  }
    … …  
    23592529// Fills in reduced costs
    23602530void
    2361 ClpSimplexPrimalQuadratic::createDjs (const ClpQuadraticInfo * info,
     2531ClpSimplexPrimalQuadratic::createDjs (ClpQuadraticInfo * info,
    23622532                            CoinIndexedVector * array1, CoinIndexedVector * array2)
    23632533{
    23642534  const int * which = info->quadraticSequence();
     2535  const int * whichRow = info->quadraticRowSequence();
     2536  const int * backRow = info->backRowSequence();
     2537  int numberQuadraticRows = info->numberQuadraticRows();
     2538  //const int * back = info->backSequence();
     2539  //int numberQuadraticColumns = info->numberQuadraticColumns();
    23652540  int numberXColumns = info->numberXColumns();
    23662541  int numberXRows = info->numberXRows();
    … …  
    23682543  int iSequence;
    23692544  int start=numberXColumns+numberXRows;
    2370   int jSequence;
    23712545  const double * djWeight = info->djWeight();
    23722546  // Actually only need to do this after invert (use extra parameter to createDjs)
    … …  
    23932567    }
    23942568    // And make sure pi is zero if slack basic
    2395     for (iRow=0;iRow<numberXRows;iRow++) {
    2396       //if (getRowStatus(iRow)==basic)
    2397       //assert(!solution_[iRow+numberXColumns]);
    2398       if (getRowStatus(iRow)==basic)
     2569    for (iRow=0;iRow<numberQuadraticRows;iRow++) {
     2570      int jRow = backRow[iRow];
     2571      if (getRowStatus(jRow)==basic)
    23992572        solution_[iRow+numberXColumns]=0.0;
    24002573    }
    … …  
    24042577    // Now costs
    24052578    for (iSequence=0;iSequence<numberXColumns;iSequence++) {
    2406       iRow = iSequence+numberXRows;
    2407       double value=cost_[iSequence];
    2408       if (fabs(value)>1.0e5) {
    2409         assert (getColumnStatus(iSequence)==basic);
    2410       }
    2411       storedCost[iSequence]=value;
    2412       storedUpper[iSequence]=value;
    2413       storedValue[iSequence]=value;
    2414       value += modifiedCost[iRow];
    2415       if (value) {
    2416         modifiedCost[iRow]=value;
    2417         index[number++]=iRow;
    2418       } else {
    2419         modifiedCost[iRow]=0.0;
     2579      int jSequence = which[iSequence];
     2580      if (jSequence>=0) {
     2581        iRow = jSequence+numberXRows;
     2582        double value=cost_[iSequence];
     2583        if (fabs(value)>1.0e5) {
     2584          assert (getColumnStatus(iSequence)==basic);
     2585        }
     2586        storedCost[jSequence]=value;
     2587        storedUpper[jSequence]=value;
     2588        storedValue[jSequence]=value;
     2589        value += modifiedCost[iRow];
     2590        if (value) {
     2591          modifiedCost[iRow]=value;
     2592          index[number++]=iRow;
     2593        } else {
     2594          modifiedCost[iRow]=0.0;
     2595        }
    24202596      }
    24212597    }
    … …  
    24322608    // And slacks
    24332609    // sign? - or does any of this make sense
    2434     for (iRow=0;iRow<numberXRows;iRow++) {
    2435       int iSequence  = iRow + numberColumns_;
     2610    assert (numberQuadraticRows==numberXRows);
     2611    for (iRow=0;iRow<numberQuadraticRows;iRow++) {
     2612      int jRow = backRow[iRow];
     2613      int iSequence  = jRow + numberColumns_;
    24362614      double value = cost_[iSequence];
    24372615      // Its possible a fixed vector may have had this set
    … …  
    24592637    // Now modify pi values by slack costs to make djs
    24602638    // Later save matrix_->add results
    2461     for (iRow=0;iRow<numberXRows;iRow++) {
    2462       int iSequence  = iRow + numberColumns_;
     2639    for (iRow=0;iRow<numberQuadraticRows;iRow++) {
     2640      int jRow = backRow[iRow];
     2641      int iSequence  = jRow + numberColumns_;
    24632642      int iPi = iRow+numberXColumns;
    24642643      solution_[iPi]-=cost_[iSequence];
    … …  
    24672646    matrix_->times(1.0,solution_,modifiedCost,rowScale_,columnScale_);
    24682647    // Back to pi
    2469     for (iRow=0;iRow<numberXRows;iRow++) {
    2470       int iSequence  = iRow + numberColumns_;
     2648    for (iRow=0;iRow<numberQuadraticRows;iRow++) {
     2649      int jRow = backRow[iRow];
     2650      int iSequence  = jRow + numberColumns_;
    24712651      int iPi = iRow+numberXColumns;
    24722652      solution_[iPi]+=cost_[iSequence];
    … …  
    24812661        lower_[k+numberColumns_]=modifiedCost[k];
    24822662        upper_[k+numberColumns_]=modifiedCost[k];
     2663        //cost_[k+numberColumns_]=0.0;
    24832664      }
    24842665      modifiedCost[k]=0.0;
    24852666    }
    2486   }
    2487   // fill in linear ones
    2488   memcpy(dj_,cost_,numberXColumns*sizeof(double));
    2489   matrix_->transposeTimes(-1.0,pi,dj_,rowScale_,columnScale_);
    2490   memset(dj_+numberXColumns,0,(numberXRows+info->numberQuadraticColumns())*sizeof(double));
    2491   for (iSequence=0;iSequence<numberXColumns;iSequence++) {
    2492     int jSequence = which[iSequence];
    2493     double value;
    2494     if (jSequence>=0) {
    2495       jSequence += start;
    2496       value = solution_[jSequence];
     2667    memcpy(dj_,cost_,numberColumns_*sizeof(double));
     2668    matrix_->transposeTimes(-1.0,pi,dj_,rowScale_,columnScale_);
     2669    info->createGradient(this);
     2670    double * gradient = info->gradient();
     2671    // fill in as if linear
     2672    // will not be correct unless complementary - but that should not matter
     2673    number=0;
     2674    for (iRow=0;iRow<numberRows_;iRow++) {
     2675      int iPivot=pivotVariable_[iRow];
     2676      if (gradient[iPivot]) {
     2677        modifiedCost[iRow] = gradient[iPivot];
     2678        index[number++]=iRow;
     2679      }
     2680    }
     2681    array1->setNumElements(number);
     2682    factorization_->updateColumnTranspose(array2,array1);
     2683    memcpy(dual_,modifiedCost,numberRows_*sizeof(double));
     2684    matrix_->transposeTimes(-1.0,modifiedCost,gradient,rowScale_,columnScale_);
     2685   
     2686    array1->clear();
     2687    memset(dj_+numberXColumns,0,(numberXRows+info->numberQuadraticColumns())*sizeof(double));
     2688    for (iSequence=0;iSequence<numberXColumns;iSequence++) {
     2689      int jSequence = which[iSequence];
     2690      double value;
     2691      if (jSequence>=0) {
     2692        jSequence += start;
     2693        value = solution_[jSequence];
     2694      } else {
     2695        value=dj_[iSequence];
     2696        value=gradient[iSequence];
     2697      }
     2698      double value2 = value*djWeight[iSequence];
     2699      if (fabs(value2)>dualTolerance_)
     2700        value=value2;
     2701      else if (value<-dualTolerance_)
     2702        value = -1.001*dualTolerance_;
     2703      else if (value>dualTolerance_)
     2704        value = 1.001*dualTolerance_;
     2705      if (djWeight[iSequence]<1.0e-6)
     2706        value=value2;
     2707      dj_[iSequence]=value;
     2708    }
     2709#if 0
     2710    if (numberIterations_<20) {
     2711      for (iSequence=0;iSequence<min(20,numberXColumns);iSequence++)
     2712        printf("%d %g %g\n",iSequence,dj_[iSequence],gradient[iSequence]);
    24972713    } else {
    2498       value=dj_[iSequence];
    2499     }
    2500     double value2 = value*djWeight[iSequence];
    2501     if (fabs(value2)>dualTolerance_)
    2502       value=value2;
    2503     else if (value<-dualTolerance_)
    2504       value = -1.001*dualTolerance_;
    2505     else if (value>dualTolerance_)
    2506       value = 1.001*dualTolerance_;
    2507     if (djWeight[iSequence]<1.0e-6)
    2508       value=value2;
    2509     dj_[iSequence]=value;
    2510   }
    2511   // And slacks
    2512   // Value of sj is - value of pi
    2513   for (jSequence=0;jSequence<numberXRows;jSequence++) {
    2514     int iSequence  = jSequence + numberColumns_;
    2515     int iPi = jSequence+numberXColumns;
    2516     double value = solution_[iPi]+cost_[iSequence];
    2517     double value2 = value*djWeight[iSequence];
    2518     if (fabs(value2)>dualTolerance_)
    2519       value=value2;
    2520     else if (value<-dualTolerance_)
    2521       value = -1.001*dualTolerance_;
    2522     else if (value>dualTolerance_)
    2523       value = 1.001*dualTolerance_;
    2524     if (djWeight[iSequence]<1.0e-6)
    2525       value=value2;
    2526     dj_[iSequence]=value;
     2714      exit(22);
     2715    }
     2716#endif
     2717    // And slacks
     2718    // Value of sj is - value of pi
     2719    for (iSequence=numberColumns_;iSequence<numberColumns_+numberXRows;iSequence++) {
     2720      int jSequence = whichRow[iSequence-numberColumns_];
     2721      double value;
     2722      if (jSequence>=0) {
     2723        int iPi = jSequence+numberXColumns;
     2724        value = solution_[iPi]+cost_[iSequence];
     2725      } else {
     2726        value=dj_[iSequence];
     2727        value=gradient[iSequence];
     2728      }
     2729      double value2 = value*djWeight[iSequence];
     2730      if (fabs(value2)>dualTolerance_)
     2731        value=value2;
     2732      else if (value<-dualTolerance_)
     2733        value = -1.001*dualTolerance_;
     2734      else if (value>dualTolerance_)
     2735        value = 1.001*dualTolerance_;
     2736      if (djWeight[iSequence]<1.0e-6)
     2737        value=value2;
     2738      dj_[iSequence]=value;
     2739    }
     2740  }
     2741  if (info->crucialSj()<0) {
     2742    int i;
     2743    for (i=0;i<numberXRows;i++) {
     2744      double pi1 = dual_[i];
     2745      double pi2 = solution_[numberXColumns+i];
     2746      if (fabs(pi1-pi2)>1.0e-3)
     2747        printf("row %d, dual %g sol %g\n",i,pi1,pi2);
     2748    }
    25272749  }
    25282750}
    … …  
    27602982  // First workout how many rows extra
    27612983  info=ClpQuadraticInfo(this);
     2984  quadratic = info.quadraticObjective();
    27622985  int numberQuadratic = info.numberQuadraticColumns();
    27632986  int newNumberRows = numberRows+numberQuadratic;
    … …  
    28443067  newNumberColumns=numberColumns;
    28453068  // pi
     3069  int numberQuadraticRows = info.numberQuadraticRows();
    28463070  // Get a row copy in standard format
    28473071  CoinPackedMatrix copy;
    … …  
    28623086         j++) {
    28633087      double elementValue=elementByRow[j];
    2864       int jColumn = column[j];
    2865       elements2[numberElements]=elementValue;
    2866       row2[numberElements++]=jColumn+numberRows;
     3088      int jColumn = which[column[j]];
     3089      if (jColumn>=0) {
     3090        elements2[numberElements]=elementValue;
     3091        row2[numberElements++]=jColumn+numberRows;
     3092      }
    28673093      dj[jColumn]-= value*elementValue;
    28683094    }
     3095#ifndef CORRECT_ROW_COUNTS
    28693096    newNumberColumns++;
    28703097    start2[newNumberColumns]=numberElements;
     3098#else
     3099    if (numberElements>start2[newNumberColumns]) {
     3100      newNumberColumns++;
     3101      start2[newNumberColumns]=numberElements;
     3102    }
     3103#endif
    28713104  }
    28723105  // djs
    … …  
    28763109    solution2[newNumberColumns]=dj[iColumn];
    28773110    elements2[numberElements]=1.0;
    2878     row2[numberElements++]=back[iColumn]+numberRows;
     3111    row2[numberElements++]=iColumn+numberRows;
    28793112    newNumberColumns++;
    28803113    start2[newNumberColumns]=numberElements;
    … …  
    29363169  delete [] elements2;
    29373170  model2->allSlackBasis();
    2938   //model2->scaling(false);
     3171  model2->scaling(false);
     3172  printf("scaling off\n");
    29393173  model2->setLogLevel(this->logLevel());
    29403174  // Move solution across
    … …  
    29833217  }
    29843218#else
     3219  const int * whichQuadraticRow = info.quadraticRowSequence();
    29853220  for (iRow=0;iRow<numberRows;iRow++) {
    29863221    assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
    29873222    assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
    29883223    model2->setRowStatus(iRow,basic);
    2989     model2->setColumnStatus(iRow+numberColumns,superBasic);
    2990     solution2[iRow+numberColumns]=0.0;
    2991   }
    2992   for (iColumn=numberRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
     3224    int jRow = whichQuadraticRow[iRow];
     3225    if (jRow>=0) {
     3226      model2->setColumnStatus(jRow+numberColumns,superBasic);
     3227      solution2[jRow+numberColumns]=0.0;
     3228    }
     3229  }
     3230  for (iColumn=numberQuadraticRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
    29933231    model2->setColumnStatus(iColumn,basic);
    2994     model2->setRowStatus(iColumn-numberColumns,isFixed);
     3232    model2->setRowStatus(iColumn-numberQuadraticRows-numberColumns+numberRows,isFixed);
    29953233  }
    29963234#endif
    … …  
    29983236  model2->times(1.0,solution2,rowSolution2);
    29993237  for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
    3000     int iS = back[iColumn]+newNumberRows;
     3238    int iS = iColumn+numberColumns+numberQuadraticRows;
    30013239    int iRow = iColumn+numberRows;
    30023240    double value = rowSolution2[iRow];
    … …  
    30393277                   ClpQuadraticInfo & info)
    30403278{
    3041   memcpy(dualRowSolution(),quadraticModel->primalColumnSolution()+numberColumns_,numberRows_*sizeof(double));
     3279  memcpy(dualRowSolution(),quadraticModel->dualRowSolution(),numberRows_*sizeof(double));
    30423280  const double * solution2 = quadraticModel->primalColumnSolution();
    30433281  memcpy(primalColumnSolution(),solution2,numberColumns_*sizeof(double));
    … …  
    30463284  quadraticModel->times(1.0,quadraticModel->primalColumnSolution(),
    30473285               quadraticModel->primalRowSolution());
    3048   memcpy(dualColumnSolution(),
    3049          quadraticModel->primalColumnSolution()+numberRows_+numberColumns_,
    3050          numberColumns_*sizeof(double));
    30513286
    30523287  int iColumn;
    … …  
    30733308    const double * quadraticElement = quadratic->getElements();
    30743309    const int * which = info.quadraticSequence();
    3075     int start = numberRows_+numberColumns_;
     3310    int start = info.numberQuadraticRows()+numberColumns_;
    30763311    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    30773312      int jColumn = which[iColumn];
    … …  
    31343369    quadraticSequence_(NULL),
    31353370    backSequence_(NULL),
     3371    quadraticRowSequence_(NULL),
     3372    backRowSequence_(NULL),
    31363373    currentSequenceIn_(-1),
    31373374    crucialSj_(-1),
    … …  
    31433380    validSolution_(NULL),
    31443381    djWeight_(NULL),
     3382    gradient_(NULL),
    31453383    numberXRows_(-1),
    31463384    numberXColumns_(-1),
    31473385    numberQuadraticColumns_(0),
     3386    numberQuadraticRows_(0),
    31483387    infeasCost_(0.0)
    31493388{
    … …  
    31543393    quadraticSequence_(NULL),
    31553394    backSequence_(NULL),
     3395    quadraticRowSequence_(NULL),
     3396    backRowSequence_(NULL),
    31563397    currentSequenceIn_(-1),
    31573398    crucialSj_(-1),
    … …  
    31633404    validSolution_(NULL),
    31643405    djWeight_(NULL),
     3406    gradient_(NULL),
    31653407    numberXRows_(-1),
    31663408    numberXColumns_(-1),
    31673409    numberQuadraticColumns_(0),
     3410    numberQuadraticRows_(0),
    31683411    infeasCost_(0.0)
    31693412{
    … …  
    31713414    numberXRows_ = model->numberRows();
    31723415    numberXColumns_ = model->numberColumns();
     3416    //ClpQuadraticObjective *originalObjective = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
     3417    //assert (originalObjective);
    31733418    originalObjective_ = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
    31743419    assert (originalObjective_);
    31753420    quadraticSequence_ = new int[numberXColumns_];
    31763421    backSequence_ = new int[numberXColumns_];
     3422    quadraticRowSequence_ = new int[numberXRows_];
     3423    backRowSequence_ = new int[numberXRows_];
    31773424    int i;
    3178     numberQuadraticColumns_=numberXColumns_;
     3425    const CoinPackedMatrix * quadratic =
     3426      originalObjective_->quadraticObjective();
     3427    //const int * columnQuadratic = quadratic->getIndices();
     3428    //const int * columnQuadraticStart = quadratic->getVectorStarts();
     3429    const int * columnQuadraticLength = quadratic->getVectorLengths();
     3430    //const double * quadraticElement = quadratic->getElements();
     3431    memset(quadraticRowSequence_,0,numberXRows_*sizeof(int));
     3432    numberQuadraticColumns_=0;
     3433    assert (numberXColumns_==quadratic->getNumCols());
     3434    const int * row = model->matrix()->getIndices();
     3435    const int * columnStart = model->matrix()->getVectorStarts();
     3436    const int * columnLength = model->matrix()->getVectorLengths();
    31793437    for (i=0;i<numberXColumns_;i++) {
    3180       quadraticSequence_[i]=i;
    3181       backSequence_[i]=i;
    3182     }
    3183     int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     3438      int length = columnQuadraticLength[i];
     3439#ifndef CORRECT_COLUMN_COUNTS
     3440      length=1;
     3441#endif
     3442      if (length) {
     3443        int j;
     3444        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     3445          int iRow=row[j];
     3446          quadraticRowSequence_[iRow]++;
     3447        }
     3448        quadraticSequence_[i]=numberQuadraticColumns_;
     3449        backSequence_[numberQuadraticColumns_]=i;
     3450        numberQuadraticColumns_++;
     3451      } else {
     3452        quadraticSequence_[i]=-1;
     3453      }
     3454    }
     3455    // Do counts of quadratic rows
     3456    numberQuadraticRows_=0;
     3457    for (i=0;i<numberXRows_;i++) {
     3458#ifndef CORRECT_ROW_COUNTS
     3459      quadraticRowSequence_[i]=1;
     3460#endif
     3461      if (quadraticRowSequence_[i]>0) {
     3462        quadraticRowSequence_[i]=numberQuadraticRows_;
     3463        backRowSequence_[numberQuadraticRows_++]=i;
     3464      } else {
     3465        quadraticRowSequence_[i]=-1;
     3466      }
     3467    }
     3468#if 0
     3469    for (i=0;i<numberXColumns_;i++) {
     3470      int j;
     3471      int next = columnQuadraticStart[i]+columnQuadraticLength[i];
     3472      assert (next==columnQuadraticStart[i+1]);
     3473      for (j=columnQuadraticStart[i];j<next;j++) {
     3474        int iColumn = columnQuadratic[j];
     3475        iColumn = quadraticSequence_[iColumn];
     3476        assert (iColumn>=0);
     3477        newIndices[j]=iColumn;
     3478      }
     3479    }
     3480    originalObjective_ = new ClpQuadraticObjective(originalObjective->linearObjective(),
     3481                                                   originalObjective->numberColumns(),
     3482                                                   columnQuadraticStart,newIndices,quadraticElement,
     3483                                                   originalObjective->numberExtendedColumns());
     3484#endif
     3485    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
    31843486    int numberRows = numberXRows_+numberQuadraticColumns_;
    31853487    int size = numberRows+numberColumns;
    … …  
    31943496  delete [] quadraticSequence_;
    31953497  delete [] backSequence_;
     3498  delete [] quadraticRowSequence_;
     3499  delete [] backRowSequence_;
    31963500  delete [] currentSolution_;
    31973501  delete [] validSolution_;
    31983502  delete [] djWeight_;
     3503  delete [] gradient_;
    31993504}
    32003505// Copy
    … …  
    32033508    quadraticSequence_(NULL),
    32043509    backSequence_(NULL),
     3510    quadraticRowSequence_(NULL),
     3511    backRowSequence_(NULL),
    32053512    currentSequenceIn_(rhs.currentSequenceIn_),
    32063513    crucialSj_(rhs.crucialSj_),
    … …  
    32083515    validCrucialSj_(rhs.validCrucialSj_),
    32093516    currentPhase_(rhs.currentPhase_),
     3517    currentSolution_(NULL),
    32103518    validPhase_(rhs.validPhase_),
     3519    validSolution_(NULL),
     3520    djWeight_(NULL),
     3521    gradient_(NULL),
    32113522    numberXRows_(rhs.numberXRows_),
    32123523    numberXColumns_(rhs.numberXColumns_),
    32133524    numberQuadraticColumns_(rhs.numberQuadraticColumns_),
     3525    numberQuadraticRows_(rhs.numberQuadraticRows_),
    32143526    infeasCost_(rhs.infeasCost_)
    32153527{
    … …  
    32213533    memcpy(backSequence_,rhs.backSequence_,
    32223534           numberXColumns_*sizeof(int));
    3223     int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     3535    quadraticRowSequence_ = new int[numberXRows_];
     3536    memcpy(quadraticRowSequence_,rhs.quadraticRowSequence_,
     3537           numberXRows_*sizeof(int));
     3538    backRowSequence_ = new int[numberXRows_];
     3539    memcpy(backRowSequence_,rhs.backRowSequence_,
     3540           numberXRows_*sizeof(int));
     3541    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
    32243542    int numberRows = numberXRows_+numberQuadraticColumns_;
    32253543    int size = numberRows+numberColumns;
    … …  
    32413559    } else {
    32423560      djWeight_ = NULL;
     3561    }
     3562    if (rhs.gradient_) {
     3563      gradient_ = new double [size];
     3564      memcpy(gradient_,rhs.gradient_,size*sizeof(double));
     3565    } else {
     3566      gradient_ = NULL;
    32433567    }
    32443568  }
    … …  
    32523576    delete [] quadraticSequence_;
    32533577    delete [] backSequence_;
     3578    delete [] quadraticRowSequence_;
     3579    delete [] backRowSequence_;
    32543580    delete [] currentSolution_;
    32553581    delete [] validSolution_;
    32563582    delete [] djWeight_;
     3583    delete [] gradient_;
    32573584    currentSequenceIn_ = rhs.currentSequenceIn_;
    32583585    crucialSj_ = rhs.crucialSj_;
    … …  
    32653592    infeasCost_=rhs.infeasCost_;
    32663593    numberQuadraticColumns_=rhs.numberQuadraticColumns_;
     3594    numberQuadraticRows_=rhs.numberQuadraticRows_;
    32673595    if (numberXColumns_) {
    32683596      quadraticSequence_ = new int[numberXColumns_];
    … …  
    32723600      memcpy(backSequence_,rhs.backSequence_,
    32733601             numberXColumns_*sizeof(int));
     3602      quadraticRowSequence_ = new int[numberXRows_];
     3603      memcpy(quadraticRowSequence_,rhs.quadraticRowSequence_,
     3604             numberXRows_*sizeof(int));
     3605      backRowSequence_ = new int[numberXRows_];
     3606      memcpy(backRowSequence_,rhs.backRowSequence_,
     3607           numberXRows_*sizeof(int));
    32743608    } else {
    32753609      quadraticSequence_ = NULL;
    32763610      backSequence_ = NULL;
    3277     }
    3278     int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     3611      quadraticRowSequence_ = NULL;
     3612      backRowSequence_ = NULL;
     3613    }
     3614    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
    32793615    int numberRows = numberXRows_+numberQuadraticColumns_;
    32803616    int size = numberRows+numberColumns;
    … …  
    32963632    } else {
    32973633      djWeight_ = NULL;
     3634    }
     3635    if (rhs.gradient_) {
     3636      gradient_ = new double [size];
     3637      memcpy(gradient_,rhs.gradient_,size*sizeof(double));
     3638    } else {
     3639      gradient_ = NULL;
    32983640    }
    32993641  }
    … …  
    33073649  validCrucialSj_ = crucialSj_;
    33083650  validPhase_ = currentPhase_;
    3309   int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     3651  int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
    33103652  int numberRows = numberXRows_+numberQuadraticColumns_;
    33113653  int size = numberRows+numberColumns;
    … …  
    33313673{
    33323674  if (currentPhase_) {
    3333     int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     3675    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
    33343676    int numberRows = numberXRows_+numberQuadraticColumns_;
    33353677    int size = numberRows+numberColumns;
    … …  
    33543696  return originalObjective_->linearObjective();
    33553697}
     3698void
     3699ClpQuadraticInfo::createGradient(ClpSimplex * model)
     3700{
     3701  int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
     3702  int numberRows = numberXRows_+numberQuadraticColumns_;
     3703  int size = numberRows+numberColumns;
     3704  if (!gradient_)
     3705    gradient_= new double[size];
     3706  memcpy(gradient_,model->costRegion(),size*sizeof(double));
     3707  double * solution = model->solutionRegion();
     3708  const CoinPackedMatrix * quadratic = quadraticObjective();
     3709  const int * columnQuadratic = quadratic->getIndices();
     3710  const int * columnQuadraticStart = quadratic->getVectorStarts();
     3711  const int * columnQuadraticLength = quadratic->getVectorLengths();
     3712  const double * quadraticElement = quadratic->getElements();
     3713  // get current costs
     3714  int jSequence;
     3715  for (jSequence=0;jSequence<numberQuadraticColumns_;jSequence++) {
     3716    int iSequence = backSequence_[jSequence];
     3717    // get current gradient
     3718    double coeff1=gradient_[iSequence];
     3719    int j;
     3720    for (j=columnQuadraticStart[iSequence];
     3721         j<columnQuadraticStart[iSequence]+columnQuadraticLength[iSequence];j++) {
     3722      int jColumn = columnQuadratic[j];
     3723      double valueJ = solution[jColumn];
     3724      // maybe this is just if jColumn basic ??
     3725      double elementValue = quadraticElement[j];
     3726      double addValue = valueJ*elementValue;
     3727      coeff1 += addValue;
     3728    }
     3729    gradient_[iSequence]=coeff1;
     3730  }
     3731}

  • branches/pre/Makefile.Clp

    r199 r200  
    5050CXXFLAGS += -DCLP_IDIOT
    5151CXXFLAGS += -DUSE_PRESOLVE
     52CXXFLAGS += -DCORRECT_COLUMN_COUNTS
    5253ifeq ($(OptLevel),-g)
    5354#     CXXFLAGS += -DCLP_DEBUG

  • branches/pre/Samples/Makefile

    r185 r200  
    5959#       LDFLAGS += -lhistory -lreadline -ltermcap
    6060endif
    61 #LDFLAGS += -lefence
     61LDFLAGS += -lefence
    6262###############################################################################
    6363

  • branches/pre/Test/Makefile.test

    r182 r200  
    6262        LDFLAGS += -lhistory -lreadline -ltermcap
    6363endif
    64 #LDFLAGS += -lefence
     64LDFLAGS += -lefence
    6565
    6666###############################################################################

  • branches/pre/Test/unitTest.cpp

    r199 r200  
    11221122    int i;
    11231123    start[0]=0;
     1124    int n=0;
     1125    int kk=numberColumns-1;
     1126    int kk2=numberColumns-1;
    11241127    for (i=0;i<numberColumns;i++) {
    1125       start[i+1]=i+1;
    1126       column[i]=i;
    1127       element[i]=1.0e-1;
    1128       element[i]=0.0;
     1128      if (i>=kk) {
     1129        column[n]=i;
     1130        if (i>=kk2)
     1131          element[n]=1.0e-1;
     1132        else
     1133          element[n]=0.0;
     1134        n++;
     1135      }
     1136      start[i+1]=n;
    11291137    }
    11301138    // Load up objective
    … …  
    11401148    solution.quadraticPrimal(1);
    11411149    objValue = solution.getObjValue();
    1142     assert(eq(objValue,3.2368421));
     1150    //assert(eq(objValue,3.2368421));
    11431151    //exit(77);
    11441152  }
    … …  
    11591167    double * element = NULL;
    11601168    m.readQuadraticMps(NULL,start,column,element,2);
     1169    int column2[200];
     1170    double element2[200];
     1171    int start2[80];
     1172    int j;
     1173    start2[0]=0;
     1174    int nel=0;
     1175    bool good=false;
     1176    for (j=0;j<79;j++) {
     1177      if (start[j]==start[j+1]) {
     1178        column2[nel]=j;
     1179        element2[nel]=0.0;
     1180        nel++;
     1181      } else {
     1182        int i;
     1183        for (i=start[j];i<start[j+1];i++) {
     1184          column2[nel]=column[i];
     1185          element2[nel++]=element[i];
     1186        }
     1187      }
     1188      start2[j+1]=nel;
     1189    }
    11611190    // Load up objective
    1162     model.loadQuadraticObjective(model.numberColumns(),start,column,element);
     1191    if (good)
     1192      model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
     1193    else
     1194      model.loadQuadraticObjective(model.numberColumns(),start,column,element);
    11631195    delete [] start;
    11641196    delete [] column;
    11651197    delete [] element;
     1198    int numberColumns=model.numberColumns();
    11661199#if 0
    11671200    model.quadraticSLP(50,1.0e-4);
    … …  
    11691202    // Get feasible
    11701203    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
    1171     int numberColumns=model.numberColumns();
    11721204    ClpLinearObjective zeroObjective(NULL,numberColumns);
    11731205    model.setObjective(&zeroObjective);
    … …  
    11771209#endif
    11781210    model.setLogLevel(63);
     1211    //exit(77);
    11791212    model.quadraticPrimal(1);
    11801213    double objValue = model.getObjValue();
     1214    const double * solution = model.primalColumnSolution();
     1215    int i;
     1216    for (i=0;i<numberColumns;i++)
     1217      if (solution[i])
     1218        printf("%d %g\n",i,solution[i]);
    11811219    CoinRelFltEq eq(1.0e-4);
    11821220    assert(eq(objValue,-400.92));

  • branches/pre/include/ClpLinearObjective.hpp

    r180 r200  
    2121  /** Returns gradient.  If Linear then solution may be NULL,
    2222      also returns an offset (to be added to current one)
     23      If refresh is false then uses last solution
    2324  */
    24   virtual double * gradient(const double * solution, double & offset);
     25  virtual double * gradient(const double * solution, double & offset,bool refresh);
    2526  /// Resize objective
    2627  virtual void resize(int newNumberColumns) ;

  • branches/pre/include/ClpModel.hpp

    r199 r200  
    265265    if (objective_) {
    266266      double offset;
    267       return objective_->gradient(NULL,offset);
     267      return objective_->gradient(NULL,offset,false);
     268    } else {
     269      return NULL;
     270    }
     271  }
     272   inline double * objective(const double * solution, double & offset,bool refresh=true) const           
     273  {
     274    offset=0.0;
     275    if (objective_) {
     276      return objective_->gradient(solution,offset,refresh);
    268277    } else {
    269278      return NULL;
    … …  
    274283    if (objective_) {
    275284      double offset;
    276       return objective_->gradient(NULL,offset);
     285      return objective_->gradient(NULL,offset,false);
    277286    } else {
    278287      return NULL;

  • branches/pre/include/ClpObjective.hpp

    r180 r200  
    2222  /** Returns gradient.  If Linear then solution may be NULL,
    2323      also returns an offset (to be added to current one)
     24      If refresh is false then uses last solution
    2425  */
    25   virtual double * gradient(const double * solution, double & offset) = 0;
     26  virtual double * gradient(const double * solution, double & offset,bool refresh) = 0;
    2627  /// Resize objective
    2728  virtual void resize(int newNumberColumns) = 0;

  • branches/pre/include/ClpQuadraticObjective.hpp

    r199 r200  
    2222  /** Returns gradient.  If Quadratic then solution may be NULL,
    2323      also returns an offset (to be added to current one)
     24      If refresh is false then uses last solution
    2425  */
    25   virtual double * gradient(const double * solution, double & offset);
     26  virtual double * gradient(const double * solution, double & offset,bool refresh);
    2627  /// Resize objective
    2728  virtual void resize(int newNumberColumns) ;

  • branches/pre/include/ClpSimplexPrimalQuadratic.hpp

    r199 r200  
    6464                            CoinIndexedVector * array1, CoinIndexedVector * array2);
    6565  /// Fills in reduced costs
    66   void createDjs (const ClpQuadraticInfo * info,
     66  void createDjs (ClpQuadraticInfo * info,
    6767                  CoinIndexedVector * array1, CoinIndexedVector * array2);
    6868  /** Main part.
    … …  
    140140  inline int numberXRows() const
    141141  {return numberXRows_;};
     142  /// Number of Quadratic rows
     143  inline int numberQuadraticRows() const
     144  {return numberQuadraticRows_;};
    142145  /// Sequence number incoming
    143146  inline int sequenceIn() const
    … …  
    165168  inline const int * backSequence() const
    166169  {return backSequence_;};
     170  /// Row Quadratic sequence or -1 if linear
     171  inline const int * quadraticRowSequence() const
     172  {return quadraticRowSequence_;};
     173  /// Which rows are quadratic
     174  inline const int * backRowSequence() const
     175  {return backRowSequence_;};
    167176  /// Returns pointer to original objective
    168177  inline ClpQuadraticObjective * originalObjective() const
    … …  
    181190  inline double * djWeight() const
    182191  { return djWeight_;};
     192  /// create gradient
     193  void createGradient(ClpSimplex * model);
     194  /// Current gradient
     195  inline double * gradient() const
     196  { return gradient_;};
    183197  /// Infeas cost
    184198  inline double infeasCost() const
    … …  
    197211  /// Which ones are quadratic
    198212  int * backSequence_;
     213  /// Quadratic Row sequence
     214  int * quadraticRowSequence_;
     215  /// Which rows are quadratic
     216  int * backRowSequence_;
    199217  /// Current sequenceIn
    200218  int currentSequenceIn_;
    … …  
    215233  /// Dj weights to stop looping
    216234  double * djWeight_;
     235  /// Current gradient
     236  double * gradient_;
    217237  /// Number of original rows
    218238  int numberXRows_;
    … …  
    221241  /// Number of quadratic columns
    222242  int numberQuadraticColumns_;
     243  /// Number of quadratic rows
     244  int numberQuadraticRows_;
    223245  /// Infeasibility cost
    224246  double infeasCost_;
Note: See TracChangeset for help on using the changeset viewer.