Changeset 203 for branches


Ignore:
Timestamp:
Sep 2, 2003 4:10:27 PM (16 years ago)
Author:
forrest
Message:

Quadratic improvements

Location:
branches/pre
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpFactorization.cpp

    r196 r203  
    44#include "CoinPragma.hpp"
    55#include "ClpFactorization.hpp"
     6#include "ClpQuadraticObjective.hpp"
    67#include "CoinHelperFunctions.hpp"
    78#include "CoinIndexedVector.hpp"
     
    365366  }
    366367
     368  if (!status_) {
     369    // take out part if quadratic
     370    if (model->algorithm()==2) {
     371      ClpObjective * obj = model->objectiveAsObject();
     372      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
     373      assert (quadraticObj);
     374      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     375      int numberXColumns = quadratic->getNumCols();
     376      assert (numberXColumns<numberColumns);
     377      int base = numberColumns-numberXColumns;
     378      int * which = new int [numberXColumns];
     379      int * pivotVariable = model->pivotVariable();
     380      int * permute = pivotColumn();
     381      int i;
     382      int n=0;
     383      for (i=0;i<numberRows;i++) {
     384        int iSj = pivotVariable[i]-base;
     385        if (iSj>=0&&iSj<numberXColumns)
     386          which[n++]=permute[i];
     387      }
     388      if (n)
     389        emptyRows(n,which);
     390      delete [] which;
     391    }
     392  }
    367393  return status_;
    368394}
  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r202 r203  
    529529  int numberXColumns = info->numberXColumns();
    530530  int numberXRows = info->numberXRows();
    531   const int * which = info->quadraticSequence();
    532   //const int * back = info->backSequence();
    533531  // initialize - values pass coding and algorithm_ is +1
    534   // for moment just set objective for a bit
    535532  ClpObjective * saveObj = objectiveAsObject();
    536533  setObjectivePointer(info->originalObjective());
    537534  if (!startup(1)) {
    538535
    539     setObjectivePointer(saveObj);
     536    //setObjectivePointer(saveObj);
    540537    // Setup useful stuff
    541538    info->setCurrentPhase(phase);
     
    544541    info->setCrucialSj(-1);
    545542    bool deleteCosts=false;
    546     if (scalingFlag_) {
     543    if (scalingFlag_>0) {
    547544      // scale
    548545      CoinPackedMatrix * quadratic = info->quadraticObjective();
     
    557554      const double * elementByColumn = quadratic->getElements();
    558555      double direction = optimizationDirection_;
    559       const int * which = info->quadraticSequence();
    560556      // direction is actually scale out not scale in
    561557      if (direction)
     
    571567        for (j=0;j<number;j++) {
    572568          int iColumn2 = columnsInThisColumn[j];
    573           iColumn2 = which[iColumn2];
    574           assert (iColumn2>=0);
    575569          newElement[j] = elementsInThisColumn[j]*scale*rowScale_[iColumn2+numberXRows];
    576570        }
     
    660654        int nextSequenceIn=-1;
    661655        int numberQuadraticRows = info->numberQuadraticRows();
    662         const int * whichColumn = info->quadraticSequence();
    663656        for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    664657          double value = solution_[iColumn];
     
    668661              ||getColumnStatus(iColumn)==superBasic) {
    669662            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    670               int jSequence=whichColumn[iColumn];
    671663              if (fabs(dj_[iColumn])>10.0*dualTolerance_||
    672                        (jSequence>=0&&
    673                         getColumnStatus(jSequence+numberXColumns+numberQuadraticRows)==basic)) {
     664                  getColumnStatus(iColumn+numberXColumns+numberQuadraticRows)==basic) {
    674665                if (phase!=2) {
    675666                  phase=2;
     
    679670            }
    680671          }
    681 #if 0
    682           if (which[iColumn]>=0) {
    683             if (getColumnStatus(iColumn)==basic) {
    684               int iS = which[iColumn]+numberRows_;
    685               assert (getColumnStatus(iS)!=basic);
    686               if(fabs(solution_[iS])>dualTolerance_) {
    687                 if (phase==0) {
    688                   phase=1;
    689                   nextSequenceIn=iS;
    690                 }
    691               }
    692             }
    693           }
    694 #endif
    695         }
    696         const int * whichRow = info->quadraticRowSequence();
     672        }
    697673        for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;iColumn++) {
    698674          double value = solution_[iColumn];
     
    708684            }
    709685          }
    710 #if 0
    711           // may work later when whichRow working??
    712           int iRow = iColumn-numberColumns_;
    713           if (whichRow[iRow]>=0) {
    714             if (getColumnStatus(iColumn)==basic) {
    715               int iS = whichRow[iRow]+numberXColumns;;
    716               assert (getColumnStatus(iS)!=basic);
    717               if(fabs(solution_[iS])>dualTolerance_) {
    718                 if (phase==0) {
    719                   phase=1;
    720                   nextSequenceIn=iS;
    721                 }
    722               }
    723             }
    724           }
    725 #endif
    726686        }
    727687        info->setSequenceIn(nextSequenceIn);
     
    779739  int numberXRows = info->numberXRows();
    780740  int numberQuadraticRows = info->numberQuadraticRows();
    781   int numberQuadraticColumns = info->numberQuadraticColumns();
    782   const int * whichRow = info->quadraticRowSequence();
    783   const int * whichColumn = info->quadraticSequence();
    784   const int * backRow = info->backRowSequence();
    785   const int * backColumn = info->backSequence();
    786   //const double * pi = solution_+numberXColumns;
     741  // Make list of implied sj
     742  // And backward pointers to basic variables
     743  {
     744    int * impliedSj = info->impliedSj();
     745    int i;
     746    for (i=numberRows_;i<numberColumns_;i++) {
     747      if (getColumnStatus(i)==basic)
     748        impliedSj[i-numberRows_]=i;
     749      else
     750        impliedSj[i-numberRows_]=-1;
     751    }
     752    int * basicRow = info->basicRow();
     753    for (i=0;i<numberRows_+numberColumns_;i++)
     754      basicRow[i]=-1;
     755    for (i=0;i<numberRows_;i++)
     756      basicRow[pivotVariable_[i]]=i;
     757  }
    787758  int nextSequenceIn=info->sequenceIn();
    788759  int oldSequenceIn=nextSequenceIn;
     
    829800          if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
    830801            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    831               int jSequence=whichColumn[iColumn];
    832802              if (fabs(dj_[iColumn])>10.0*dualTolerance_||
    833                        (jSequence>=0&&
    834                         getColumnStatus(jSequence+numberXColumns+numberQuadraticRows)==basic)) {
     803                  getColumnStatus(iColumn+numberXColumns+numberQuadraticRows)==basic) {
    835804                nextSequenceIn=iColumn;
    836805                break;
     
    877846        cleanupIteration=1;
    878847      }
    879       if ((sequenceIn_>=numberXColumns&&sequenceIn_<numberColumns_)&&info->crucialSj()<0)
    880         cleanupIteration=2;
    881848    }
    882849    pivotRow_=-1;
     
    904871        chosen=-1;
    905872        unpack(rowArray_[1]);
    906 #if 1
    907873        // compute dj in case linear
    908874        {
     
    927893          }
    928894        }
    929 #else
    930         dualIn_=dj_[sequenceIn_];
    931 #endif
     895        if ((!cleanupIteration&&sequenceIn_<numberXColumns)||
     896            (cleanupIteration&&info->crucialSj()<numberColumns_&&info->crucialSj()>=numberRows_)) {
     897          int iSequence;
     898          if (!cleanupIteration)
     899            iSequence=sequenceIn_;
     900          else
     901            iSequence=info->crucialSj()-numberRows_;
     902          // may need to re-do row of basis
     903          int * impliedSj = info->impliedSj();
     904          if (impliedSj[iSequence]>=0) {
     905            // mark as valid
     906            impliedSj[iSequence]=-1;
     907            ClpPackedMatrix* rowCopy =
     908              dynamic_cast< ClpPackedMatrix*>(rowCopy_);
     909            assert(rowCopy);
     910            const int * column = rowCopy->getIndices();
     911            const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     912            const double * rowElement = rowCopy->getElements();
     913            int iRow=iSequence+numberXRows;
     914            int i;
     915            int * index = rowArray_[2]->getIndices();
     916            double * element = rowArray_[2]->denseVector();
     917            int n=0;
     918            int * basicRow = info->basicRow();
     919            int * permute = factorization_->pivotColumn();
     920            for (i=rowStart[iRow];i<rowStart[iRow+1];i++) {
     921              int iColumn = column[i];
     922              iColumn = basicRow[iColumn];
     923              if (iColumn>=0&&iColumn!=iRow) {
     924                index[n]=permute[iColumn];
     925                element[n++]=rowElement[i];
     926              }
     927            }
     928            factorization_->replaceRow(permute[iRow],n,index,element);
     929            memset(element,0,n*sizeof(double));
     930          }
     931        }
     932        // Take out elements in implied Sj rows
     933        if (1) {
     934          int n=rowArray_[1]->getNumElements();
     935          int * index = rowArray_[1]->getIndices();
     936          double * element = rowArray_[1]->denseVector();
     937          int i;
     938          int * impliedSj = info->impliedSj();
     939          int n2=0;
     940          for (i=0;i<n;i++) {
     941            int iRow=index[i]-numberXRows;
     942            if (iRow<0||impliedSj[iRow]<0) {
     943              index[n2++]=iRow+numberXRows;
     944            } else {
     945              element[iRow+numberXRows]=0.0;
     946            }
     947          }
     948          rowArray_[1]->setNumElements(n2);
     949        }
    932950        factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
    933951        if (cleanupIteration) {
     
    952970                if (fabs(solution_[crucialSj])<1.0e-8)
    953971                  printf("zero crucialSj - pivot out - get right way\n");
    954                 //if (sequenceIn_>=numberColumns_)
    955                 //d2=-d2; // Try both ways? - makes no difference!!!
    956                 //if (sequenceIn_>=numberXColumns&&sequenceIn_<numberColumns_)
    957                 //d2=-d2; // Try both ways? - makes no difference!!!
    958972                // see which way to go
    959973                if (d2>0)
     
    970984            //assert(tj);
    971985            }
    972           } else {
    973             assert (cleanupIteration==2);
    974             // get correct dj
    975             int iSequence = sequenceIn_-numberXColumns;
    976             if (iSequence<numberQuadraticRows) {
    977               int iRow=backRow[iSequence];
    978               dj_[sequenceIn_] = - dj_[numberColumns_+iRow];
    979             } else {
    980               iSequence = backColumn[iSequence-numberQuadraticRows];
    981               dj_[sequenceIn_]= - dj_[iSequence];
    982             }
    983986          }
    984987
     
    994997          lowerIn_=lower_[sequenceIn_];
    995998          upperIn_=upper_[sequenceIn_];
    996           if (sequenceIn_<numberXColumns) {
    997             int jSequence = whichColumn[sequenceIn_];
    998             if (jSequence>=0) {
    999               dualIn_ = solution_[jSequence+numberXColumns+numberQuadraticRows];
    1000               dualIn_=dj_[sequenceIn_];
    1001             } else {
    1002               // need coding
    1003               //dualIn_=cost_[sequenceIn_];
    1004               //int j;
    1005               //for (j=columnStart[sequenceIn_];j<columnStart[sequenceIn_]+columnLength[sequenceIn_]; j++) {
    1006               //int iRow = row[j];
    1007               //dualIn_ -= element[j]*pi[iRow];
    1008               //}
    1009             }
    1010           } else {
    1011             int iRow = sequenceIn_-numberColumns_;
    1012             assert (iRow>=0);
    1013             if (whichRow[iRow]>=0) {
    1014               dualIn_ = solution_[numberXColumns+whichRow[iRow]];
    1015               dualIn_=dj_[sequenceIn_];
    1016               // temporary?
    1017               const int * columnLength = matrix_->getVectorLengths();
    1018               if (!columnLength[numberXColumns+whichRow[iRow]])
    1019                 dualIn_ = dual_[iRow]+cost_[sequenceIn_];
    1020             }
    1021           }
     999          dualIn_=dj_[sequenceIn_];
    10221000          valueIn_=solution_[sequenceIn_];
    10231001          if (dualIn_>0.0)
     
    10271005          if (sequenceIn_<numberColumns_) {
    10281006            // Set dj as value of slack
    1029             crucialSj = whichColumn[sequenceIn_];
    1030             if (crucialSj>=0)
    1031               crucialSj += numberQuadraticRows+numberXColumns; // sj which should go to 0.0
     1007            crucialSj = sequenceIn_+ numberQuadraticRows+numberXColumns; // sj which should go to 0.0
    10321008          } else {
    10331009            // Set dj as value of pi
    1034             int iRow = whichRow[sequenceIn_-numberColumns_];
    1035             if (iRow>=0)
    1036               crucialSj = iRow+numberXColumns; // pi which should go to 0.0
    1037             else
    1038               crucialSj=-1;
    1039             // temporary?
    1040             const int * columnLength = matrix_->getVectorLengths();
    1041             if (!columnLength[numberXColumns+iRow])
    1042               crucialSj=-1;
     1010            int iRow = sequenceIn_-numberColumns_;
     1011            crucialSj = iRow+numberXColumns; // pi which should go to 0.0
    10431012          }
    10441013          if (crucialSj>=0&&getColumnStatus(crucialSj)!=basic)
     
    10471016        }
    10481017        double oldSj=1.0e30;
    1049         double oldObjectiveValue = objectiveValue_;
    10501018        if (info->crucialSj()>=0&&cleanupIteration)
    10511019          oldSj= solution_[info->crucialSj()];
     
    10711039        saveObjective = objectiveValue_;
    10721040        if (pivotRow_>=0) {
     1041          // If sj out AND not gone out since last invert then add back
    10731042          // if stable replace in basis
    10741043          int updateStatus = 0;
     
    12521221          setStatus(sequenceOut_,superBasic);
    12531222        }
     1223        // And update backward pointers to basic variables
     1224        if (sequenceIn_!=sequenceOut_) {
     1225          int * basicRow = info->basicRow();
     1226          basicRow[sequenceOut_]=-1;
     1227          basicRow[sequenceIn_]=pivotRow_;
     1228        }
    12541229        if (info->crucialSj()>=0) {
    12551230          assert(fabs(oldSj)>= fabs(solution_[info->crucialSj()]));
     
    12591234        printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
    12601235               objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost());
    1261         if (objectiveValue_>oldObjectiveValue) {
    1262           // make less likely this one will come in again
    1263           double multiplier = CoinDrand48();
    1264           int iSequence;
    1265           if (!cleanupIteration) {
    1266             iSequence = sequenceIn_;
    1267           } else {
    1268             iSequence = info->crucialSj();
    1269             if (iSequence>numberRows_)
    1270               iSequence -= numberRows_;
    1271             else
    1272               iSequence = numberColumns_ +(iSequence-numberXColumns);
    1273           }
    1274           info->djWeight()[iSequence] /= (5.0+multiplier);
    1275         }
    12761236        if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
    12771237          // really free but going to zero
     
    12831243          int crucialSj=info->crucialSj();
    12841244          if (crucialSj>=numberXColumns+numberQuadraticRows) {
    1285             if (sequenceOut_==backColumn[crucialSj-numberXColumns-numberQuadraticRows])
     1245            if (sequenceOut_==crucialSj-numberXColumns-numberQuadraticRows)
    12861246              info->setCrucialSj(-1);
    12871247          } else {
    1288             if (sequenceOut_-numberColumns_==whichRow[crucialSj-numberXColumns])
     1248            if (sequenceOut_-numberColumns_==crucialSj-numberXColumns)
    12891249              info->setCrucialSj(-1);
    12901250          }
     
    13081268          int crucialSj=info->crucialSj();
    13091269          if (sequenceOut_<numberXColumns) {
    1310             chosen =whichColumn[sequenceOut_] + numberQuadraticRows + numberXColumns; // sj variable
     1270            chosen =sequenceOut_ + numberQuadraticRows + numberXColumns; // sj variable
    13111271          } else if (sequenceOut_>=numberColumns_) {
    13121272            // Does this mean we can change pi
    13131273            int iRow = sequenceOut_-numberColumns_;
    13141274            if (iRow<numberXRows) {
    1315               int jRow = whichRow[iRow];
    1316               assert (jRow>=0);
    1317               int iPi = jRow+numberXColumns;
    1318               //printf("pi for row %d is %g\n",
    1319               //     iRow,solution_[iPi]);
    1320               chosen=iPi;
     1275              chosen=iRow+numberXColumns;
    13211276            } else {
    13221277              printf("Row %d is in column part\n",iRow);
     
    13251280          } else if (sequenceOut_<numberQuadraticRows+numberXColumns) {
    13261281            // pi out
    1327             chosen = backRow[sequenceOut_-numberXColumns]+numberColumns_;
     1282            chosen = sequenceOut_-numberXColumns+numberColumns_;
    13281283          } else {
    13291284            // sj out
    1330             chosen = backColumn[sequenceOut_-numberQuadraticRows-numberXColumns];
     1285            chosen = sequenceOut_-numberQuadraticRows-numberXColumns;
    13311286          }
    13321287          // But double check as original incoming might have gone out
     
    13351290            break; // means original has gone out
    13361291          }
    1337           // In certain circumstances need to look further to see what to come in
    1338           // Think further
    1339           if (numberXColumns>numberQuadraticColumns) {
    1340             // this should be pre-computed
    1341             int iSjRow=-1;
    1342             int iRow;
    1343             for (iRow=0;iRow<numberRows_;iRow++) {
    1344               if (pivotVariable_[iRow]==crucialSj) {
    1345                 iSjRow=iRow;
    1346                 break;
    1347               }
    1348             }
    1349             assert (iSjRow>=0);
    1350             double direction;
    1351             if (solution_[crucialSj]>0)
    1352               direction=-1.0;
    1353             else
    1354               direction=1.0;
    1355             rowArray_[0]->createPacked(1,&iSjRow,&direction);
    1356             factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
    1357             // put row of tableau in rowArray[0] and columnArray[0]
    1358             matrix_->transposeTimes(this,-1.0,
    1359                                     rowArray_[0],rowArray_[3],columnArray_[0]);
    1360             int k;
    1361             // all packed
    1362             // Column
    1363             int number=columnArray_[0]->getNumElements();
    1364             int * which = columnArray_[0]->getIndices();
    1365             double * element = columnArray_[0]->denseVector();
    1366             int best=-1;
    1367             double bestAlpha=1.0e-6;
    1368             if (numberIterations_==1180)
    1369               printf("iteration %d coming up\n",numberIterations_);
    1370             for (k=0;k<number;k++) {
    1371               int iSequence=which[k];
    1372               double value=element[k];
    1373               // choose chosen if possible
    1374               if (iSequence==chosen) {
    1375                 if (value>0.0)
    1376                   printf("chosen %d going up alpha %g, sol %g\n",iSequence,value,solution_[crucialSj]);
    1377                 else
    1378                   printf("chosen %d going down alpha %g, sol %g\n",iSequence,value,solution_[crucialSj]);
    1379                 value = fabs(value)*1.0e5;
    1380               } else if (iSequence<numberXColumns) {
    1381                 double djValue = dj_[iSequence];
    1382                 // take out if other one basic
    1383                 int sj = whichColumn[iSequence];
    1384                 if (sj>=0&&getColumnStatus(sj+numberXColumns+numberQuadraticRows)==basic)
    1385                   value=0.0;
    1386                 if (value>0.0)
    1387                   printf("X%d going up alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
    1388                 else if (value<0.0)
    1389                   printf("X%d going down alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
    1390                 ClpSimplex::Status status = getColumnStatus(iSequence);
    1391                
    1392                 switch(status) {
    1393                  
    1394                 case ClpSimplex::basic:
    1395                   if (fabs(value)>1.0e-3)
    1396                     abort();
    1397                   break;
    1398                 case ClpSimplex::isFixed:
    1399                   value=0.0;
    1400                   break;
    1401                 case ClpSimplex::isFree:
    1402                 case ClpSimplex::superBasic:
    1403                   value=fabs(value);
    1404                   break;
    1405                 case ClpSimplex::atUpperBound:
    1406                   if (djValue>=-1.0e-8)
    1407                     value = -value;
    1408                   else
    1409                     value = -1.0e-4*value;
    1410                   break;
    1411                 case ClpSimplex::atLowerBound:
    1412                   if (djValue>=1.0e-8)
    1413                     value = 1.0e-4*value;
    1414                   break;
    1415                 }
    1416               } else {
    1417                 value=0.0;
    1418               }
    1419               if (value>bestAlpha) {
    1420                 bestAlpha=value;
    1421                 best=iSequence;
    1422               }
    1423             }
    1424             // Row
    1425             number=rowArray_[0]->getNumElements();
    1426             which = rowArray_[0]->getIndices();
    1427             element = rowArray_[0]->denseVector();
    1428             for (k=0;k<number;k++) {
    1429               int iSequence=which[k];
    1430               double value=element[k];
    1431               // choose chosen if possible
    1432               if (iSequence==chosen) {
    1433                 if (value>0.0)
    1434                   printf("chosen row %d going up alpha %g, sol %g\n",iSequence,value,solution_[crucialSj]);
    1435                 else
    1436                   printf("chosen row %d going down alpha %g, sol %g\n",iSequence,value,solution_[crucialSj]);
    1437                 value = fabs(value)*1.0e5;
    1438               } else if (iSequence<numberXRows) {
    1439                 double djValue = dj_[iSequence+numberColumns_];
    1440                 // take out if other one basic
    1441                 int sj = whichRow[iSequence];
    1442                 if (sj>=0&&getColumnStatus(sj+numberXColumns)==basic)
    1443                   value=0.0;
    1444                 if (value>0.0)
    1445                   printf("R%d going up alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
    1446                 else if (value<0.0)
    1447                   printf("R%d going down alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
    1448                 ClpSimplex::Status status = getRowStatus(iSequence);
    1449                
    1450                 switch(status) {
    1451                  
    1452                 case ClpSimplex::basic:
    1453                   if (fabs(value)>1.0e-3)
    1454                     abort();
    1455                   break;
    1456                 case ClpSimplex::isFixed:
    1457                   value=0.0;
    1458                   break;
    1459                 case ClpSimplex::isFree:
    1460                 case ClpSimplex::superBasic:
    1461                   value=fabs(value);
    1462                   break;
    1463                 case ClpSimplex::atUpperBound:
    1464                   if (djValue>=-1.0e-8)
    1465                     value = -value;
    1466                   else
    1467                     value = -1.0e-4*value;
    1468                   break;
    1469                 case ClpSimplex::atLowerBound:
    1470                   if (djValue>=1.0e-8)
    1471                     value = 1.0e-4*value;
    1472                   break;
    1473                 }
    1474               } else {
    1475                 value=0.0;
    1476               }
    1477               if (value>bestAlpha) {
    1478                 bestAlpha=value;
    1479                 best=iSequence+numberColumns_;
    1480               }
    1481             }
    1482             assert (best>=0);
    1483             if (best!=chosen)
    1484               printf("** chosen changed from %d to %d\n",chosen,best);
    1485             chosen=best;
    1486           }
    14871292          rowArray_[0]->clear();
    1488           rowArray_[3]->checkClear();
    14891293          columnArray_[0]->clear();
    14901294          if (returnCode==-2) {
     
    14941298          }
    14951299        } else {
    1496           // may need to bring sj in
    1497           if (0&&(returnCode==-1||returnCode==-2||returnCode==-3)) {
    1498             assert(info->crucialSj()<0);
    1499             chosen=-1;
    1500             if (sequenceOut_<numberXColumns&&whichColumn[sequenceOut_]>=0) {
    1501               chosen =whichColumn[sequenceOut_] + numberQuadraticRows + numberXColumns; // sj variable
    1502               if (getColumnStatus(chosen)==basic)
    1503                 chosen=-1; // already in
    1504             } else if (sequenceOut_>=numberColumns_&&whichRow[sequenceOut_-numberColumns_]>=0) {
    1505               chosen = whichRow[sequenceOut_-numberColumns_]+numberXColumns;
    1506               if (getColumnStatus(chosen)==basic)
    1507                 chosen=-1; // already in
    1508               // temporary?
    1509               const int * columnLength = matrix_->getVectorLengths();
    1510               if (chosen>=0&&!columnLength[chosen])
    1511                 chosen=-1;
    1512             }
    1513             if (chosen>=0) {
    1514               printf("what now\n");
    1515               nextSequenceIn=chosen;
    1516             }
    1517           }
    15181300          break;
    15191301        }
     
    16521434  int iSjRow2=-1,crucialSj2=-1;
    16531435  if (sequenceIn_<numberXColumns) {
    1654     const int * which = info->quadraticSequence();
    1655     crucialSj2 = which[sequenceIn_];
    1656     if (crucialSj2>=0)
    1657       crucialSj2 += numberRows_; // sj2 which should go to 0.0
     1436    crucialSj2 = sequenceIn_;
     1437    crucialSj2 += numberRows_; // sj2 which should go to 0.0
    16581438  } else if (sequenceIn_>=numberColumns_) {
    1659     const int * which = info->quadraticRowSequence();
    16601439    int iRow=sequenceIn_-numberColumns_;
    1661     crucialSj2 = which[iRow];
    1662     if (crucialSj2>=0)
    1663       crucialSj2 += numberXColumns; // sj2 which should go to 0.0
     1440    crucialSj2 = iRow;
     1441    crucialSj2 += numberXColumns; // sj2 which should go to 0.0
    16641442  }
    16651443  double tj = 0.0;
     
    16721450    double alpha = -work[iRow]*way;
    16731451    int iPivot=pivotVariable_[iRow];
     1452    if (numberIterations_==17)
     1453      printf("row %d col %d alpha %g solution %g\n",iRow,iPivot,alpha,solution_[iPivot]);
    16741454    if (iPivot<numberXColumns) {
    16751455      index2[number2++]=iPivot;
     
    22792059        valueOut_ = solution(sequenceOut_);
    22802060        theta_ = fabs(valueOut_/alpha_);
    2281         assert (fabs(maximumMovement-theta_)<1.0e-3*(1.0+maximumMovement));
     2061        if (fabs(maximumMovement-theta_)>1.0e-3*(1.0+maximumMovement))
     2062          printf("maximumMovement %g mismatch with theta %g\n",maximumMovement,theta_);;
    22822063        if (way<0.0)
    22832064          theta_ = - theta_;
     
    27072488                            CoinIndexedVector * array1, CoinIndexedVector * array2)
    27082489{
    2709   const int * which = info->quadraticSequence();
    2710   const int * whichRow = info->quadraticRowSequence();
    2711   const int * backRow = info->backRowSequence();
    27122490  int numberQuadraticRows = info->numberQuadraticRows();
    2713   const int * back = info->backSequence();
    27142491  int numberQuadraticColumns = info->numberQuadraticColumns();
    27152492  int numberXColumns = info->numberXColumns();
    27162493  int numberXRows = info->numberXRows();
    2717   const double * pi = solution_+numberXColumns;
    2718   int iSequence;
    2719   int start=numberXColumns+numberXRows;
    2720   const double * djWeight = info->djWeight();
    27212494  // Actually only need to do this after invert (use extra parameter to createDjs)
    27222495  // or when infeasibilities change
    27232496  // recode to do this later and update rather than recompute
    27242497  // When it is "change" then we don't need b (original rhs)
    2725   if (nonLinearCost_->numberInfeasibilities()<-1) {
    2726   } else {
    2727     array1->checkClear();
    2728     array2->checkClear();
     2498  {
    27292499    // update costs
    27302500    double * storedCost = lower_+numberColumns_+numberXRows;
     
    27332503    int * index = array1->getIndices();
    27342504    double * modifiedCost = array1->denseVector();
    2735     //#define CHECK
    2736     //#ifdef CHECK
    2737 #if 0
    2738     // zero out basic values
    2739     int iRow;
    2740     for (iRow=0;iRow<numberRows_;iRow++) {
    2741       int iPivot=pivotVariable_[iRow];
    2742       // save
    2743       dj_[iPivot]=solution_[iPivot];
    2744       solution_[iPivot] = 0.0;
    2745     }
    2746     // And make sure pi is zero if slack basic
    2747     for (iRow=0;iRow<numberQuadraticRows;iRow++) {
    2748       int jRow = backRow[iRow];
    2749       if (getRowStatus(jRow)==basic)
    2750         solution_[iRow+numberXColumns]=0.0;
    2751     }
    2752     matrix_->times(-1.0,columnActivityWork_,modifiedCost,rowScale_,columnScale_);
    2753    
    2754     int number=0;
    2755     // Now costs
    2756     for (iSequence=0;iSequence<numberXColumns;iSequence++) {
    2757       int jSequence = which[iSequence];
    2758       if (jSequence>=0) {
    2759         iRow = jSequence+numberXRows;
    2760         double value=cost_[iSequence];
    2761         if (fabs(value)>1.0e5) {
    2762           assert (getColumnStatus(iSequence)==basic);
    2763         }
    2764         storedCost[jSequence]=value;
    2765         storedUpper[jSequence]=value;
    2766         storedValue[jSequence]=value;
    2767         value += modifiedCost[iRow];
    2768         if (value) {
    2769           modifiedCost[iRow]=value;
    2770           index[number++]=iRow;
    2771         } else {
    2772           modifiedCost[iRow]=0.0;
    2773         }
    2774       }
    2775     }
    2776     for (iRow=0;iRow<numberXRows;iRow++) {
    2777       double value = modifiedCost[iRow]+rowActivityWork_[iRow];;
    2778       if (value) {
    2779         modifiedCost[iRow]=value;
    2780         index[number++]=iRow;
    2781       } else {
    2782         modifiedCost[iRow]=0.0;
    2783       }
    2784     }
    2785     array1->setNumElements(number);
    2786     // And slacks
    2787     // sign? - or does any of this make sense
    2788     assert (numberQuadraticRows==numberXRows);
    2789     for (iRow=0;iRow<numberQuadraticRows;iRow++) {
    2790       int jRow = backRow[iRow];
    2791       int iSequence  = jRow + numberColumns_;
    2792       double value = cost_[iSequence];
    2793       // Its possible a fixed vector may have had this set
    2794       if (value&&getRowStatus(iRow)==basic) {
    2795         // add in pi column
    2796         matrix_->add(this,array1,iRow+numberXColumns,value);
    2797       }
    2798     }
    2799     // create pricing vector
    2800     int k;
    2801     factorization_->updateColumn(array2,array1);
    2802     int n=array1->getNumElements();
    2803     for (k=0;k<n;k++) {
    2804       int iRow = index[k];
    2805       double value=modifiedCost[iRow];
    2806       int iPivot = pivotVariable_[iRow];
    2807       if (iPivot>=numberXColumns&&iPivot<numberColumns_)
    2808         solution_[iPivot] = value;
    2809       else
    2810         solution_[iPivot] = dj_[iPivot];
    2811       //solution_[iPivot] = value;
    2812      
    2813     }
    2814     array1->clear();
    2815     // Now modify pi values by slack costs to make djs
    2816     // Later save matrix_->add results
    2817     for (iRow=0;iRow<numberQuadraticRows;iRow++) {
    2818       int jRow = backRow[iRow];
    2819       int iSequence  = jRow + numberColumns_;
    2820       int iPi = iRow+numberXColumns;
    2821       solution_[iPi]-=cost_[iSequence];
    2822     }
    2823     // check looks okay
    2824     matrix_->times(1.0,solution_,modifiedCost,rowScale_,columnScale_);
    2825     // Back to pi
    2826     for (iRow=0;iRow<numberQuadraticRows;iRow++) {
    2827       int jRow = backRow[iRow];
    2828       int iSequence  = jRow + numberColumns_;
    2829       int iPi = iRow+numberXColumns;
    2830       solution_[iPi]+=cost_[iSequence];
    2831     }
    2832     double * rowSolution = solution_+numberColumns_;
    2833     for (k=0;k<numberRows_;k++) {
    2834       if (k<numberXRows) {
    2835         //if (fabs(modifiedCost[k]-rowSolution[k])>1.0e-3*(1.0+fabs(rowSolution[k])))
    2836         //printf("bad row %d %g %g\n",k,modifiedCost[k],rowSolution[k]);
    2837       } else {
    2838         rowSolution[k]=modifiedCost[k];
    2839         lower_[k+numberColumns_]=modifiedCost[k];
    2840         upper_[k+numberColumns_]=modifiedCost[k];
    2841         //cost_[k+numberColumns_]=0.0;
    2842       }
    2843       modifiedCost[k]=0.0;
    2844     }
    2845     memcpy(dj_,cost_,numberColumns_*sizeof(double));
    2846     matrix_->transposeTimes(-1.0,pi,dj_,rowScale_,columnScale_);
     2505    // Compute duals
    28472506    info->createGradient(this);
    28482507    double * gradient = info->gradient();
    28492508    // fill in as if linear
    28502509    // will not be correct unless complementary - but that should not matter
    2851     number=0;
     2510    // Just save sj value in that case
     2511    //double saveValue=0.0;
     2512    //int crucialSj = info->crucialSj();
     2513    //if (crucialSj>=0)
     2514    int number=0;
     2515    int iRow;
    28522516    for (iRow=0;iRow<numberRows_;iRow++) {
    28532517      int iPivot=pivotVariable_[iRow];
     
    28602524    factorization_->updateColumnTranspose(array2,array1);
    28612525    memcpy(dual_,modifiedCost,numberRows_*sizeof(double));
    2862     matrix_->transposeTimes(-1.0,modifiedCost,gradient,rowScale_,columnScale_);
    2863    
    28642526    array1->clear();
    2865     memset(dj_+numberXColumns,0,(numberXRows+info->numberQuadraticColumns())*sizeof(double));
    2866     for (iSequence=0;iSequence<numberXColumns;iSequence++) {
    2867       int jSequence = which[iSequence];
    2868       double value;
    2869       if (jSequence>=0) {
    2870         jSequence += start;
    2871         //solution_[jSequence]=gradient[iSequence];
    2872         value = solution_[jSequence];
     2527    memcpy(dj_,gradient,(numberColumns_+numberRows_)*sizeof(double));
     2528    matrix_->transposeTimes(-1.0,dual_,dj_,rowScale_,columnScale_);
     2529    memset(dj_+numberXColumns,0,(numberXRows+numberQuadraticColumns)*sizeof(double));
     2530    // And djs
     2531    for (iRow=0;iRow<numberRows_;iRow++) {
     2532      dj_[numberColumns_+iRow]= cost_[numberColumns_+iRow]+dual_[iRow];
     2533    }
     2534    double saveSj=0.0;
     2535    if (info->crucialSj()>=0)
     2536      saveSj=solution_[info->crucialSj()];
     2537    // Set dual solution
     2538    for (iRow=0;iRow<numberQuadraticRows;iRow++) {
     2539      solution_[iRow+numberXColumns]=dual_[iRow];
     2540    }
     2541    // clear sj
     2542    int start = numberXColumns+numberQuadraticRows;
     2543    memset(solution_+start,0,numberQuadraticColumns*sizeof(double));
     2544    memset(dj_+numberXColumns,0,(numberQuadraticRows+numberQuadraticColumns)*sizeof(double));
     2545    matrix_->times(-1.0,columnActivityWork_,modifiedCost,rowScale_,columnScale_);
     2546    memset(modifiedCost,0,numberXRows*sizeof(double));
     2547    // Do costs as rhs and sj solution values
     2548    for (iRow=0;iRow<numberQuadraticColumns;iRow++) {
     2549      int jSequence = iRow;
     2550      double value=cost_[jSequence];
     2551      if (fabs(value)>1.0e5) {
     2552        assert (getColumnStatus(jSequence)==basic);
     2553      }
     2554      double value2=-modifiedCost[iRow+numberXRows];
     2555      modifiedCost[iRow+numberXRows]=0.0;
     2556      jSequence = iRow+start;
     2557      if (getColumnStatus(jSequence)!=basic) {
     2558        if (fabs(value2-value)>1.0e-3)
     2559          solution_[jSequence]=0.0;
     2560        value=value2;
    28732561      } else {
    2874         value=dj_[iSequence];
    2875         value=gradient[iSequence];
    2876       }
    2877       double value2 = value*djWeight[iSequence];
    2878       if (fabs(value2)>dualTolerance_)
    2879         value=value2;
    2880       else if (value<-dualTolerance_)
    2881         value = -1.001*dualTolerance_;
    2882       else if (value>dualTolerance_)
    2883         value = 1.001*dualTolerance_;
    2884       if (djWeight[iSequence]<1.0e-6)
    2885         value=value2;
    2886       dj_[iSequence]=value;
    2887     }
    2888 #if 0
    2889     if (numberIterations_<20) {
    2890       for (iSequence=0;iSequence<min(20,numberXColumns);iSequence++)
    2891         printf("%d %g %g\n",iSequence,dj_[iSequence],gradient[iSequence]);
    2892     } else {
    2893       exit(22);
    2894     }
    2895 #endif
    2896     // And slacks
    2897     // Value of sj is - value of pi
    2898     for (iSequence=numberColumns_;iSequence<numberColumns_+numberXRows;iSequence++) {
    2899       int iRow = iSequence-numberColumns_;
    2900       int jSequence = whichRow[iRow];
    2901       double value;
    2902       if (jSequence>=0) {
    2903         int iPi = jSequence+numberXColumns;
    2904         value = solution_[iPi]+cost_[iSequence];
    2905         const int * columnLength = matrix_->getVectorLengths();
    2906         if (!columnLength[iPi])
    2907           value = dual_[iRow]+cost_[iSequence];
    2908       } else {
    2909         value=dj_[iSequence];
    2910         value=dual_[iRow]+cost_[iSequence];
    2911       }
    2912       double value2 = value*djWeight[iSequence];
    2913       if (fabs(value2)>dualTolerance_)
    2914         value=value2;
    2915       else if (value<-dualTolerance_)
    2916         value = -1.001*dualTolerance_;
    2917       else if (value>dualTolerance_)
    2918         value = 1.001*dualTolerance_;
    2919       if (djWeight[iSequence]<1.0e-6)
    2920         value=value2;
    2921       dj_[iSequence]=value;
    2922     }
    2923 #if 0
    2924     if (info->crucialSj()<0) {
    2925       for (iSequence=numberXColumns+numberQuadraticRows;iSequence<numberColumns_;iSequence++)
    2926         if (fabs(solution_[iSequence])>1.0e-4)
    2927           printf("sj %d (%d) value %g\n",iSequence-numberXColumns-numberQuadraticRows,
    2928                  back[iSequence-numberXColumns-numberQuadraticRows],solution_[iSequence]);
    2929     }
    2930     if (info->crucialSj()<0) {
    2931       int i;
    2932       for (i=0;i<numberXRows;i++) {
    2933         double pi1 = dual_[i];
    2934         double pi2 = solution_[numberXColumns+i];
    2935         if (fabs(pi1-pi2)>1.0e-3)
    2936           printf("row %d, dual %g sol %g\n",i,pi1,pi2);
    2937       }
    2938     }
    2939 #endif
    2940 #endif
    2941 #if 1
    2942     {
    2943       //if (info->crucialSj()>=0)
    2944       //return;
    2945 #ifdef CHECK
    2946       double saveSol[2000];
    2947       memcpy(saveSol,solution_,(numberRows_+numberColumns_)*sizeof(double));
    2948       double saveDj[2000];
    2949       memcpy(saveDj,dj_,(numberRows_+numberColumns_)*sizeof(double));
    2950 #endif
    2951       // Compute duals
    2952       info->createGradient(this);
    2953       double * gradient = info->gradient();
    2954       // fill in as if linear
    2955       // will not be correct unless complementary - but that should not matter
    2956       // Just save sj value in that case
    2957       //double saveValue=0.0;
    2958       //int crucialSj = info->crucialSj();
    2959       //if (crucialSj>=0)
    2960       int number=0;
    2961       int iRow;
    2962       for (iRow=0;iRow<numberRows_;iRow++) {
    2963         int iPivot=pivotVariable_[iRow];
    2964         if (gradient[iPivot]) {
    2965           modifiedCost[iRow] = gradient[iPivot];
    2966           index[number++]=iRow;
    2967         }
    2968       }
    2969       array1->setNumElements(number);
    2970       factorization_->updateColumnTranspose(array2,array1);
    2971       memcpy(dual_,modifiedCost,numberRows_*sizeof(double));
    2972       array1->clear();
    2973       memcpy(dj_,gradient,(numberColumns_+numberRows_)*sizeof(double));
    2974       matrix_->transposeTimes(-1.0,dual_,dj_,rowScale_,columnScale_);
    2975       memset(dj_+numberXColumns,0,(numberXRows+numberQuadraticColumns)*sizeof(double));
    2976       // And djs
    2977       for (iRow=0;iRow<numberRows_;iRow++) {
    2978         dj_[numberColumns_+iRow]= cost_[numberColumns_+iRow]+dual_[iRow];
    2979       }
    2980       double saveSj=0.0;
    2981       if (info->crucialSj()>=0)
    2982         saveSj=solution_[info->crucialSj()];
    2983       if (info->crucialSj()>=0&&0) {
    2984 #ifdef CHECK
    2985         for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++)
    2986           if (fabs(solution_[iSequence]-saveSol[iSequence])>1.0e-4)
    2987             printf("Badxsol %d value %g true %g\n",iSequence,
    2988                    solution_[iSequence],saveSol[iSequence]);
    2989         for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++)
    2990           if (fabs(dj_[iSequence]-saveDj[iSequence])>1.0e-4)
    2991             printf("Badxdj %d value %g true %g\n",iSequence,
    2992                    dj_[iSequence],saveDj[iSequence]);
    2993         memcpy(solution_,saveSol,(numberRows_+numberColumns_)*sizeof(double));
    2994         memcpy(dj_,saveDj,(numberRows_+numberColumns_)*sizeof(double));
    2995 #endif
    2996         if (numberIterations_==-1289) {
    2997           int i;
    2998           printf ("returning on sj\n");
    2999           for (i=0;i<numberRows_+numberColumns_;i++) {
    3000             char x='N';
    3001             if (getColumnStatus(i)==basic)
    3002               x='B';
    3003             printf("CH %d %c sol %g dj %g\n",i,x,solution_[i],dj_[i]);
    3004           }
    3005         }
    3006         return;
    3007       }
    3008       // Set dual solution
    3009       for (iRow=0;iRow<numberQuadraticRows;iRow++) {
    3010         int jRow = backRow[iRow];
    3011         solution_[iRow+numberXColumns]=dual_[jRow];
    3012       }
    3013       // clear sj
    3014       int start = numberXColumns+numberQuadraticRows;
    3015       memset(solution_+start,0,numberQuadraticColumns*sizeof(double));
    3016       memset(dj_+numberXColumns,0,(numberQuadraticRows+numberQuadraticColumns)*sizeof(double));
    3017       matrix_->times(-1.0,columnActivityWork_,modifiedCost,rowScale_,columnScale_);
    3018       memset(modifiedCost,0,numberXRows*sizeof(double));
    3019       // Do costs as rhs and sj solution values
    3020       for (iRow=0;iRow<numberQuadraticColumns;iRow++) {
    3021         int jSequence = back[iRow];
    3022         double value=cost_[jSequence];
    3023         if (fabs(value)>1.0e5) {
    3024           assert (getColumnStatus(jSequence)==basic);
    3025         }
    3026         double value2=-modifiedCost[iRow+numberXRows];
    3027         modifiedCost[iRow+numberXRows]=0.0;
    3028         jSequence = iRow+start;
    3029         if (getColumnStatus(jSequence)!=basic) {
    3030           if (fabs(value2-value)>1.0e-3)
    3031             //printf("bad nonbasic match %d %g %g\n",iRow,value,value2);
    3032           // should adjust value?
    3033           solution_[jSequence]=0.0;
    3034           value=value2;
    3035         } else {
    3036           //if (fabs(value-value2-dj_[jSequence-start])>1.0e-3)
    3037           //printf("bad basic match %d %g %g\n",iRow,value,value2);
    3038           solution_[jSequence]=value-value2;
    3039         }
    3040         storedCost[iRow]=value;
    3041         storedUpper[iRow]=value;
    3042         storedValue[iRow]=value;
    3043       }
    3044       if (info->crucialSj()>=0)
    3045         solution_[info->crucialSj()]=saveSj;
    3046 #ifdef CHECK
    3047       for (iSequence=numberXColumns+numberQuadraticRows;iSequence<numberColumns_;iSequence++)
    3048         if (fabs(solution_[iSequence])>1.0e-4||fabs(saveSol[iSequence])>1.0e-4)
    3049           printf("sj %d (%d) value %g true %g\n",iSequence-numberXColumns-numberQuadraticRows,
    3050                  back[iSequence-numberXColumns-numberQuadraticRows],solution_[iSequence],saveSol[iSequence]);
    3051       for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++)
    3052         if (fabs(solution_[iSequence]-saveSol[iSequence])>1.0e-4)
    3053           printf("Badsol %d value %g true %g\n",iSequence,
    3054                  solution_[iSequence],saveSol[iSequence]);
    3055       for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++)
    3056         if (fabs(dj_[iSequence]-saveDj[iSequence])>1.0e-4)
    3057           printf("Baddj %d value %g true %g\n",iSequence,
    3058                  dj_[iSequence],saveDj[iSequence]);
    3059       memcpy(solution_,saveSol,(numberRows_+numberColumns_)*sizeof(double));
    3060       memcpy(dj_,saveDj,(numberRows_+numberColumns_)*sizeof(double));
    3061 #endif
    3062     }
    3063 #endif
     2562        solution_[jSequence]=value-value2;
     2563      }
     2564      storedCost[iRow]=value;
     2565      storedUpper[iRow]=value;
     2566      storedValue[iRow]=value;
     2567    }
     2568    if (info->crucialSj()>=0)
     2569      solution_[info->crucialSj()]=saveSj;
    30642570  }
    30652571  if (numberIterations_==-1289) {
     
    30852591  numberDualInfeasibilities_=0;
    30862592  sumDualInfeasibilities_=0.0;
    3087   const int * which = info->quadraticSequence();
    30882593  // Compute objective function from scratch
    30892594  const CoinPackedMatrix * quadratic =
     
    31292634      objectiveValue_ += valueI*valueJ*elementValue;
    31302635    }
    3131     int jSequence = which[iColumn];
    3132     if (jSequence>=0)
    3133       jSequence += start;
     2636    int jSequence = iColumn+start;
    31342637    double value=dj_[iColumn];
    31352638    ClpSimplex::Status status = getColumnStatus(iColumn);
     
    31442647     
    31452648    case ClpSimplex::basic:
    3146       if (jSequence>=0) {
    3147         if (getColumnStatus(jSequence)==basic) {
    3148           handler_->message(CLP_QUADRATIC_BOTH,messages_)
    3149             <<"Structural"<<iColumn
    3150             <<solution_[iColumn]<<jSequence<<solution_[jSequence]
    3151             <<CoinMessageEol;
    3152           number2++;
    3153           int crucialSj =info->crucialSj();
    3154           assert (crucialSj>=0);
    3155           assert (crucialSj==iColumn||crucialSj==jSequence);
    3156         } else {
    3157           number1++;
    3158         }
     2649      if (getColumnStatus(jSequence)==basic) {
     2650        handler_->message(CLP_QUADRATIC_BOTH,messages_)
     2651          <<"Structural"<<iColumn
     2652          <<solution_[iColumn]<<jSequence<<solution_[jSequence]
     2653          <<CoinMessageEol;
     2654        number2++;
     2655        assert (info->crucialSj()>=0);
     2656        assert (info->crucialSj()==iColumn||info->crucialSj()==jSequence);
     2657      } else {
     2658        number1++;
    31592659      }
    31602660      break;
     
    32132713               iRow,solution_[iColumn],jSequence,solution_[jSequence]);
    32142714          number2++;
    3215           int crucialSj =info->crucialSj();
    3216           assert (crucialSj>=0);
    3217           assert (crucialSj==iColumn||crucialSj==jSequence);
     2715          assert (info->crucialSj()>=0);
     2716          assert (info->crucialSj()==iColumn||info->crucialSj()==jSequence);
    32182717        } else {
    32192718          number1++;
     
    32422741    }
    32432742  }
    3244   printf("number 0 - %d, 1 - %d, 2 - %d\n",number0,number1,number2);
     2743  //printf("number 0 - %d, 1 - %d, 2 - %d\n",number0,number1,number2);
    32452744  objectiveValue_ += linearCost+infeasCost;
    32462745  assert (infeasCost>=0.0);
     
    33362835  double * rowLower2 = new double[newNumberRows];
    33372836  double * rowUpper2 = new double[newNumberRows];
    3338   const int * which = info.quadraticSequence();
    3339   const int * back = info.backSequence();
    33402837  memcpy(rowLower2,rowLower,numberRows*sizeof(double));
    33412838  memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
    33422839  int iRow;
    33432840  for (iRow=0;iRow<numberQuadratic;iRow++) {
    3344     double cost = objective[back[iRow]];
     2841    double cost = objective[iRow];
    33452842    rowLower2[iRow+numberRows]=cost;
    33462843    rowUpper2[iRow+numberRows]=cost;
     
    33762873      row2[numberElements++]=row[j];
    33772874    }
    3378     if (which[iColumn]>=0) {
    3379       // Quadratic and modify djs
    3380       for (j=columnQuadraticStart[iColumn];
    3381            j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
    3382            j++) {
    3383         int jColumn = columnQuadratic[j];
    3384         double value = quadraticElement[j];
    3385         if (iColumn!=jColumn)
    3386           value *= 0.5;
     2875    // Quadratic and modify djs
     2876    for (j=columnQuadraticStart[iColumn];
     2877         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
     2878         j++) {
     2879      int jColumn = columnQuadratic[j];
     2880      double value = quadraticElement[j];
     2881      if (iColumn!=jColumn)
     2882        value *= 0.5;
     2883      dj[jColumn] += solution[iColumn]*value;
     2884      elements2[numberElements]=-value;
     2885      row2[numberElements++]=jColumn+numberRows;
     2886    }
     2887    for (j=rowStartQ[iColumn];
     2888         j<rowStartQ[iColumn]+rowLengthQ[iColumn];
     2889         j++) {
     2890      int jColumn = columnQ[j];
     2891      double value = elementByRowQ[j];
     2892      if (iColumn!=jColumn) {
     2893        value *= 0.5;
    33872894        dj[jColumn] += solution[iColumn]*value;
    33882895        elements2[numberElements]=-value;
    3389         row2[numberElements++]=which[jColumn]+numberRows;
    3390       }
    3391       for (j=rowStartQ[iColumn];
    3392            j<rowStartQ[iColumn]+rowLengthQ[iColumn];
    3393            j++) {
    3394         int jColumn = columnQ[j];
    3395         double value = elementByRowQ[j];
    3396         if (iColumn!=jColumn) {
    3397           value *= 0.5;
    3398           dj[jColumn] += solution[iColumn]*value;
    3399           elements2[numberElements]=-value;
    3400           row2[numberElements++]=which[jColumn]+numberRows;
    3401         }
     2896        row2[numberElements++]=jColumn+numberRows;
    34022897      }
    34032898    }
     
    34252920         j++) {
    34262921      double elementValue=elementByRow[j];
    3427       int jColumn = which[column[j]];
     2922      int jColumn = column[j];
    34282923      if (jColumn>=0) {
    34292924        elements2[numberElements]=elementValue;
     
    35563051  }
    35573052#else
    3558   const int * whichQuadraticRow = info.quadraticRowSequence();
    35593053  for (iRow=0;iRow<numberRows;iRow++) {
    35603054    assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
    35613055    assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
    35623056    model2->setRowStatus(iRow,basic);
    3563     int jRow = whichQuadraticRow[iRow];
    3564     if (jRow>=0) {
    3565       model2->setColumnStatus(jRow+numberColumns,superBasic);
    3566       solution2[jRow+numberColumns]=0.0;
    3567     }
     3057    int jRow = iRow;
     3058    model2->setColumnStatus(jRow+numberColumns,superBasic);
     3059    solution2[jRow+numberColumns]=0.0;
    35683060  }
    35693061  for (iColumn=numberQuadraticRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
     
    36463138    const int * columnQuadraticLength = quadratic->getVectorLengths();
    36473139    const double * quadraticElement = quadratic->getElements();
    3648     const int * which = info.quadraticSequence();
    36493140    int start = info.numberQuadraticRows()+numberColumns_;
    36503141    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3651       int jColumn = which[iColumn];
     3142      int jColumn = iColumn;
    36523143      double valueI=solution2[iColumn];
    36533144      double value;
     
    37063197ClpQuadraticInfo::ClpQuadraticInfo()
    37073198  : originalObjective_(NULL),
    3708     quadraticSequence_(NULL),
    3709     backSequence_(NULL),
    3710     quadraticRowSequence_(NULL),
    3711     backRowSequence_(NULL),
     3199    basicRow_(NULL),
     3200    impliedSj_(NULL),
    37123201    currentSequenceIn_(-1),
    37133202    crucialSj_(-1),
     
    37303219ClpQuadraticInfo::ClpQuadraticInfo(const ClpSimplex * model)
    37313220  : originalObjective_(NULL),
    3732     quadraticSequence_(NULL),
    3733     backSequence_(NULL),
    3734     quadraticRowSequence_(NULL),
    3735     backRowSequence_(NULL),
     3221    basicRow_(NULL),
     3222    impliedSj_(NULL),
    37363223    currentSequenceIn_(-1),
    37373224    crucialSj_(-1),
     
    37573244    originalObjective_ = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
    37583245    assert (originalObjective_);
    3759     quadraticSequence_ = new int[numberXColumns_];
    3760     backSequence_ = new int[numberXColumns_];
    3761     quadraticRowSequence_ = new int[numberXRows_];
    3762     backRowSequence_ = new int[numberXRows_];
     3246    impliedSj_ = new int[numberXColumns_];
     3247    basicRow_ = new int[numberXColumns_];
    37633248    int i;
    3764     const CoinPackedMatrix * quadratic =
    3765       originalObjective_->quadraticObjective();
    3766     //const int * columnQuadratic = quadratic->getIndices();
    3767     //const int * columnQuadraticStart = quadratic->getVectorStarts();
    3768     const int * columnQuadraticLength = quadratic->getVectorLengths();
    3769     //const double * quadraticElement = quadratic->getElements();
    3770     memset(quadraticRowSequence_,0,numberXRows_*sizeof(int));
    3771     numberQuadraticColumns_=0;
    3772     assert (numberXColumns_==quadratic->getNumCols());
    3773     const int * row = model->matrix()->getIndices();
    3774     const int * columnStart = model->matrix()->getVectorStarts();
    3775     const int * columnLength = model->matrix()->getVectorLengths();
    3776     for (i=0;i<numberXColumns_;i++) {
    3777       int length = columnQuadraticLength[i];
    3778 #ifndef CORRECT_COLUMN_COUNTS
    3779       length=1;
    3780 #endif
    3781       if (length) {
    3782         int j;
    3783         for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
    3784           int iRow=row[j];
    3785           quadraticRowSequence_[iRow]++;
    3786         }
    3787         quadraticSequence_[i]=numberQuadraticColumns_;
    3788         backSequence_[numberQuadraticColumns_]=i;
    3789         numberQuadraticColumns_++;
    3790       } else {
    3791         quadraticSequence_[i]=-1;
    3792       }
    3793     }
    3794     // Do counts of quadratic rows
    3795     numberQuadraticRows_=0;
    3796     for (i=0;i<numberXRows_;i++) {
    3797 #ifndef CORRECT_ROW_COUNTS
    3798       quadraticRowSequence_[i]=1;
    3799 #endif
    3800       if (quadraticRowSequence_[i]>0) {
    3801         quadraticRowSequence_[i]=numberQuadraticRows_;
    3802         backRowSequence_[numberQuadraticRows_++]=i;
    3803       } else {
    3804         quadraticRowSequence_[i]=-1;
    3805       }
    3806     }
    3807 #if 0
    3808     for (i=0;i<numberXColumns_;i++) {
    3809       int j;
    3810       int next = columnQuadraticStart[i]+columnQuadraticLength[i];
    3811       assert (next==columnQuadraticStart[i+1]);
    3812       for (j=columnQuadraticStart[i];j<next;j++) {
    3813         int iColumn = columnQuadratic[j];
    3814         iColumn = quadraticSequence_[iColumn];
    3815         assert (iColumn>=0);
    3816         newIndices[j]=iColumn;
    3817       }
    3818     }
    3819     originalObjective_ = new ClpQuadraticObjective(originalObjective->linearObjective(),
    3820                                                    originalObjective->numberColumns(),
    3821                                                    columnQuadraticStart,newIndices,quadraticElement,
    3822                                                    originalObjective->numberExtendedColumns());
    3823 #endif
     3249    numberQuadraticColumns_=numberXColumns_;
     3250    numberQuadraticRows_=numberXRows_;
    38243251    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
    38253252    int numberRows = numberXRows_+numberQuadraticColumns_;
    38263253    int size = numberRows+numberColumns;
    38273254    djWeight_ = new double [size];
     3255    basicRow_ = new int[size];
    38283256    for (i=0;i<size;i++)
    38293257      djWeight_[i]=1.0;
     
    38333261ClpQuadraticInfo:: ~ClpQuadraticInfo()
    38343262{
    3835   delete [] quadraticSequence_;
    3836   delete [] backSequence_;
    3837   delete [] quadraticRowSequence_;
    3838   delete [] backRowSequence_;
     3263  delete [] impliedSj_;
     3264  delete [] basicRow_;
    38393265  delete [] currentSolution_;
    38403266  delete [] validSolution_;
     
    38453271ClpQuadraticInfo::ClpQuadraticInfo(const ClpQuadraticInfo& rhs)
    38463272  : originalObjective_(rhs.originalObjective_),
    3847     quadraticSequence_(NULL),
    3848     backSequence_(NULL),
    3849     quadraticRowSequence_(NULL),
    3850     backRowSequence_(NULL),
     3273    basicRow_(NULL),
     3274    impliedSj_(NULL),
    38513275    currentSequenceIn_(rhs.currentSequenceIn_),
    38523276    crucialSj_(rhs.crucialSj_),
     
    38663290{
    38673291  if (numberXColumns_) {
    3868     quadraticSequence_ = new int[numberXColumns_];
    3869     memcpy(quadraticSequence_,rhs.quadraticSequence_,
    3870            numberXColumns_*sizeof(int));
    3871     backSequence_ = new int[numberXColumns_];
    3872     memcpy(backSequence_,rhs.backSequence_,
    3873            numberXColumns_*sizeof(int));
    3874     quadraticRowSequence_ = new int[numberXRows_];
    3875     memcpy(quadraticRowSequence_,rhs.quadraticRowSequence_,
    3876            numberXRows_*sizeof(int));
    3877     backRowSequence_ = new int[numberXRows_];
    3878     memcpy(backRowSequence_,rhs.backRowSequence_,
    3879            numberXRows_*sizeof(int));
    38803292    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
    38813293    int numberRows = numberXRows_+numberQuadraticColumns_;
    38823294    int size = numberRows+numberColumns;
     3295    impliedSj_ = new int[numberXColumns_];
     3296    memcpy(impliedSj_,rhs.impliedSj_,numberXColumns_*sizeof(int));
     3297    basicRow_ = new int [size];
     3298    memcpy(basicRow_,rhs.basicRow_,size*sizeof(int));
    38833299    if (rhs.currentSolution_) {
    38843300      currentSolution_ = new double [size];
     
    39133329  if (this != &rhs) {
    39143330    originalObjective_ = rhs.originalObjective_;
    3915     delete [] quadraticSequence_;
    3916     delete [] backSequence_;
    3917     delete [] quadraticRowSequence_;
    3918     delete [] backRowSequence_;
     3331    delete [] impliedSj_;
     3332    delete [] basicRow_;
    39193333    delete [] currentSolution_;
    39203334    delete [] validSolution_;
     
    39323346    numberQuadraticColumns_=rhs.numberQuadraticColumns_;
    39333347    numberQuadraticRows_=rhs.numberQuadraticRows_;
    3934     if (numberXColumns_) {
    3935       quadraticSequence_ = new int[numberXColumns_];
    3936       memcpy(quadraticSequence_,rhs.quadraticSequence_,
    3937              numberXColumns_*sizeof(int));
    3938       backSequence_ = new int[numberXColumns_];
    3939       memcpy(backSequence_,rhs.backSequence_,
    3940              numberXColumns_*sizeof(int));
    3941       quadraticRowSequence_ = new int[numberXRows_];
    3942       memcpy(quadraticRowSequence_,rhs.quadraticRowSequence_,
    3943              numberXRows_*sizeof(int));
    3944       backRowSequence_ = new int[numberXRows_];
    3945       memcpy(backRowSequence_,rhs.backRowSequence_,
    3946            numberXRows_*sizeof(int));
    3947     } else {
    3948       quadraticSequence_ = NULL;
    3949       backSequence_ = NULL;
    3950       quadraticRowSequence_ = NULL;
    3951       backRowSequence_ = NULL;
    3952     }
    39533348    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
    39543349    int numberRows = numberXRows_+numberQuadraticColumns_;
    39553350    int size = numberRows+numberColumns;
     3351    impliedSj_ = new int[numberXColumns_];
     3352    memcpy(impliedSj_,rhs.impliedSj_,numberXColumns_*sizeof(int));
     3353    basicRow_ = new int [size];
     3354    memcpy(basicRow_,rhs.basicRow_,size*sizeof(int));
    39563355    if (rhs.currentSolution_) {
    39573356      currentSolution_ = new double [size];
     
    40533452  int jSequence;
    40543453  for (jSequence=0;jSequence<numberQuadraticColumns_;jSequence++) {
    4055     int iSequence = backSequence_[jSequence];
     3454    int iSequence = jSequence;
    40563455    // get current gradient
    40573456    double coeff1=gradient_[iSequence];
  • branches/pre/Makefile.Clp

    r200 r203  
    66# between then specify the exact level you want, e.g., -O1 or -O2
    77OptLevel := -g
    8 #OptLevel := -O3
     8OptLevel := -O3
    99
    1010
     
    5050CXXFLAGS += -DCLP_IDIOT
    5151CXXFLAGS += -DUSE_PRESOLVE
    52 CXXFLAGS += -DCORRECT_COLUMN_COUNTS
     52#CXXFLAGS += -DCORRECT_COLUMN_COUNTS
    5353ifeq ($(OptLevel),-g)
    5454#     CXXFLAGS += -DCLP_DEBUG
  • branches/pre/Test/Makefile.test

    r200 r203  
    6262        LDFLAGS += -lhistory -lreadline -ltermcap
    6363endif
    64 LDFLAGS += -lefence
     64#LDFLAGS += -lefence
    6565
    6666###############################################################################
  • branches/pre/Test/unitTest.cpp

    r201 r203  
    12101210    model.setLogLevel(63);
    12111211    //exit(77);
    1212     model.setFactorizationFrequency(1);
     1212    model.setFactorizationFrequency(10);
    12131213    model.quadraticPrimal(1);
    12141214    double objValue = model.getObjValue();
  • branches/pre/include/ClpSimplexPrimalQuadratic.hpp

    r201 r203  
    163163  {return currentSolution_;};
    164164  void setCurrentSolution(const double * solution);
    165   /// Quadratic sequence or -1 if linear
    166   inline const int * quadraticSequence() const
    167   {return quadraticSequence_;};
    168   /// Which ones are quadratic
    169   inline const int * backSequence() const
    170   {return backSequence_;};
    171   /// Row Quadratic sequence or -1 if linear
    172   inline const int * quadraticRowSequence() const
    173   {return quadraticRowSequence_;};
    174   /// Which rows are quadratic
    175   inline const int * backRowSequence() const
    176   {return backRowSequence_;};
    177165  /// Returns pointer to original objective
    178166  inline ClpQuadraticObjective * originalObjective() const
     
    201189  inline void setInfeasCost(double value)
    202190  { infeasCost_ = value;};
     191  /// Backward pointer to basis (inverse of pivotVariable_)
     192  inline int * basicRow() const
     193  { return basicRow_;};
     194  /// Set if Sj variable is implied
     195  inline int * impliedSj() const
     196  { return impliedSj_;};
    203197  //@}
    204198   
     
    208202  /// Objective
    209203  ClpQuadraticObjective * originalObjective_;
    210   /// Quadratic sequence
    211   int * quadraticSequence_;
    212   /// Which ones are quadratic
    213   int * backSequence_;
    214   /// Quadratic Row sequence
    215   int * quadraticRowSequence_;
    216   /// Which rows are quadratic
    217   int * backRowSequence_;
     204  /// Backward pointer to basis (inverse of pivotVariable_)
     205  int * basicRow_;
     206  /// Set if Sj variable is implied
     207  int * impliedSj_;
    218208  /// Current sequenceIn
    219209  int currentSequenceIn_;
Note: See TracChangeset for help on using the changeset viewer.