Changeset 185 for branches


Ignore:
Timestamp:
Jul 23, 2003 4:28:06 PM (16 years ago)
Author:
forrest
Message:

Update to Quadratic

Location:
branches/pre
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpNonLinearCost.cpp

    r180 r185  
    4040   later may do dual analysis and even finding duplicate columns
    4141*/
    42 ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model,
    43                                      int numberOriginalColumns)
     42ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model)
    4443{
    4544  model_ = model;
     
    6867  double * cost = model_->costRegion();
    6968
    70   bool forQuadratic = (numberOriginalColumns>0);
    71 
    7269  // For quadratic we need -inf,0,0,+inf
    73   if (!forQuadratic) {
    74     for (iSequence=0;iSequence<numberTotal;iSequence++) {
    75       if (lower[iSequence]>-1.0e20)
    76         put++;
    77       if (upper[iSequence]<1.0e20)
    78         put++;
    79       put += 2;
    80     }
    81   } else {
    82     put = 4*numberTotal;
     70  for (iSequence=0;iSequence<numberTotal;iSequence++) {
     71    if (lower[iSequence]>-1.0e20)
     72      put++;
     73    if (upper[iSequence]<1.0e20)
     74      put++;
     75    put += 2;
    8376  }
    8477
     
    9487  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    9588
    96     if (!forQuadratic||iSequence<numberOriginalColumns) {
    97       if (lower[iSequence]>-COIN_DBL_MAX) {
    98         lower_[put] = -COIN_DBL_MAX;
    99         setInfeasible(put,true);
    100         cost_[put++] = cost[iSequence]-infeasibilityCost;
    101       }
    102       whichRange_[iSequence]=put;
    103       lower_[put] = lower[iSequence];
    104       cost_[put++] = cost[iSequence];
    105       lower_[put] = upper[iSequence];
    106       cost_[put++] = cost[iSequence]+infeasibilityCost;
    107       if (upper[iSequence]<1.0e20) {
    108         lower_[put] = COIN_DBL_MAX;
    109         setInfeasible(put-1,true);
    110         cost_[put++] = 1.0e50;
    111       }
    112       start_[iSequence+1]=put;
    113     } else {
    114       // quadratic  variable
     89    if (lower[iSequence]>-COIN_DBL_MAX) {
    11590      lower_[put] = -COIN_DBL_MAX;
    11691      setInfeasible(put,true);
    117       cost_[put++] = -infeasibilityCost;
    118       whichRange_[iSequence]=put;
    119       lower_[put] = max(-0.5*COIN_DBL_MAX,lower[iSequence]);
    120       cost_[put++] = 0.0;
    121       lower_[put] = min(0.5*COIN_DBL_MAX,upper[iSequence]);
    122       cost_[put++] = infeasibilityCost;
     92      cost_[put++] = cost[iSequence]-infeasibilityCost;
     93    }
     94    whichRange_[iSequence]=put;
     95    lower_[put] = lower[iSequence];
     96    cost_[put++] = cost[iSequence];
     97    lower_[put] = upper[iSequence];
     98    cost_[put++] = cost[iSequence]+infeasibilityCost;
     99    if (upper[iSequence]<1.0e20) {
    123100      lower_[put] = COIN_DBL_MAX;
    124101      setInfeasible(put-1,true);
    125       cost_[put++] = 0.0;
    126       start_[iSequence+1]=put;
    127     }
     102      cost_[put++] = 1.0e50;
     103    }
     104    start_[iSequence+1]=put;
    128105  }
    129106
     
    765742  return lower_[jRange];
    766743}
    767 // Sets inside bounds (i.e. non infinite - used in QP
    768 void
    769 ClpNonLinearCost::setBounds(int sequence, double lower, double upper)
    770 {
    771   int start = start_[sequence];
    772   assert(start_[sequence+1]==start+4);
    773   lower_[start+1]=lower;
    774   lower_[start+2]=upper;
    775 }
    776 
     744
  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r183 r185  
    477477ClpSimplexPrimalQuadratic::primalQuadratic(int phase)
    478478{
    479 #if 0
    480   // Dantzig's method looks easier - lets try that
    481 
    482   int numberColumns = this->numberColumns();
    483   double * columnLower = this->columnLower();
    484   double * columnUpper = this->columnUpper();
    485   double * objective = this->objective();
    486   double * solution = this->primalColumnSolution();
    487   double * dj = this->dualColumnSolution();
    488   double * pi = this->dualRowSolution();
    489 
    490   int numberRows = this->numberRows();
    491   double * rowLower = this->rowLower();
    492   double * rowUpper = this->rowUpper();
    493 
    494   // and elements
    495   CoinPackedMatrix * matrix = this->matrix();
    496   const int * row = matrix->getIndices();
    497   const int * columnStart = matrix->getVectorStarts();
    498   const double * element =  matrix->getElements();
    499   const int * columnLength = matrix->getVectorLengths();
    500 
    501   // Get list of non linear columns
    502   CoinPackedMatrix * quadratic = quadraticObjective();
    503   if (!quadratic||!quadratic->getNumElements()) {
    504     // no quadratic part
    505     return primal(1);
    506   }
    507 
    508   int iColumn;
    509   const int * columnQuadratic = quadratic->getIndices();
    510   const int * columnQuadraticStart = quadratic->getVectorStarts();
    511   const int * columnQuadraticLength = quadratic->getVectorLengths();
    512   const double * quadraticElement = quadratic->getElements();
    513 #if 0
    514   // deliberate bad solution
    515   // Change to use phase
    516   //double * saveO = new double[numberColumns];
    517   //memcpy(saveO,objective,numberColumns*sizeof(double));
    518   //memset(objective,0,numberColumns*sizeof(double));
    519   tempMessage messageHandler(this);;
    520   passInMessageHandler(&messageHandler);
    521   factorization()->maximumPivots(1);
    522   primal();
    523   CoinMessageHandler handler2;
    524   passInMessageHandler(&handler2);
    525   factorization()->maximumPivots(100);
    526   setMaximumIterations(1000);
    527 #endif
    528   //memcpy(objective,saveO,numberColumns*sizeof(double));
    529   //printf("For testing - deliberate bad solution\n");
    530   //columnUpper[0]=0.0;
    531   //columnLower[0]=0.0;
    532   //quadraticSLP(50,1.0e-4);
    533   //primal(1);
    534   //columnUpper[0]=COIN_DBL_MAX;
    535  
    536   // Create larger problem
    537   // First workout how many rows extra
    538   ClpQuadraticInfo info(this);
    539   int numberQuadratic = info.numberQuadraticColumns();
    540   int newNumberRows = numberRows+numberQuadratic;
    541   int newNumberColumns = numberColumns + numberRows + numberQuadratic;
    542   int numberElements = 2*matrix->getNumElements()
    543     +2*quadratic->getNumElements()
    544     + numberQuadratic;
    545   double * elements2 = new double[numberElements];
    546   int * start2 = new int[newNumberColumns+1];
    547   int * row2 = new int[numberElements];
    548   double * objective2 = new double[newNumberColumns];
    549   double * columnLower2 = new double[newNumberColumns];
    550   double * columnUpper2 = new double[newNumberColumns];
    551   double * rowLower2 = new double[newNumberRows];
    552   double * rowUpper2 = new double[newNumberRows];
    553   const int * which = info.quadraticSequence();
    554   const int * back = info.backSequence();
    555   memcpy(rowLower2,rowLower,numberRows*sizeof(double));
    556   memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
    557   int iRow;
    558   for (iRow=0;iRow<numberQuadratic;iRow++) {
    559     double cost = objective[back[iRow]];
    560     rowLower2[iRow+numberRows]=cost;
    561     rowUpper2[iRow+numberRows]=cost;
    562   }
    563   memset(objective2,0,newNumberColumns*sizeof(double));
    564   // Get a row copy of quadratic objective in standard format
    565   CoinPackedMatrix copyQ;
    566   copyQ.reverseOrderedCopyOf(*quadratic);
    567   const int * columnQ = copyQ.getIndices();
    568   const CoinBigIndex * rowStartQ = copyQ.getVectorStarts();
    569   const int * rowLengthQ = copyQ.getVectorLengths();
    570   const double * elementByRowQ = copyQ.getElements();
    571   // Move solution across
    572   double * solution2 = new double[newNumberColumns];
    573   memset(solution2,0,newNumberColumns*sizeof(double));
    574   newNumberColumns=0;
    575   numberElements=0;
    576   start2[0]=0;
    577   // x
    578   memcpy(dj,objective,numberColumns*sizeof(double));
    579   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    580     // Original matrix
    581     columnLower2[iColumn]=columnLower[iColumn];
    582     columnUpper2[iColumn]=columnUpper[iColumn];
    583     solution2[iColumn]=solution[iColumn];
    584     int j;
    585     for (j=columnStart[iColumn];
    586          j<columnStart[iColumn]+columnLength[iColumn];
    587          j++) {
    588       elements2[numberElements]=element[j];
    589       row2[numberElements++]=row[j];
    590     }
    591     if (which[iColumn]>=0) {
    592       // Quadratic and modify djs
    593       for (j=columnQuadraticStart[iColumn];
    594            j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
    595            j++) {
    596         int jColumn = columnQuadratic[j];
    597         double value = quadraticElement[j];
    598         if (iColumn!=jColumn)
    599           value *= 0.5;
    600         dj[jColumn] += solution[iColumn]*value;
    601         elements2[numberElements]=-value;
    602         row2[numberElements++]=which[jColumn]+numberRows;
    603       }
    604       for (j=rowStartQ[iColumn];
    605            j<rowStartQ[iColumn]+rowLengthQ[iColumn];
    606            j++) {
    607         int jColumn = columnQ[j];
    608         double value = elementByRowQ[j];
    609         if (iColumn!=jColumn) {
    610           value *= 0.5;
    611           dj[jColumn] += solution[iColumn]*value;
    612           elements2[numberElements]=-value;
    613           row2[numberElements++]=which[jColumn]+numberRows;
    614         }
    615       }
    616     }
    617     start2[iColumn+1]=numberElements;
    618   }
    619   newNumberColumns=numberColumns;
    620   // pi
    621   // Get a row copy in standard format
    622   CoinPackedMatrix copy;
    623   copy.reverseOrderedCopyOf(*(this->matrix()));
    624   // get matrix data pointers
    625   const int * column = copy.getIndices();
    626   const CoinBigIndex * rowStart = copy.getVectorStarts();
    627   const int * rowLength = copy.getVectorLengths();
    628   const double * elementByRow = copy.getElements();
    629   for (iRow=0;iRow<numberRows;iRow++) {
    630     solution2[newNumberColumns]=pi[iRow];
    631     double value = pi[iRow];
    632     columnLower2[newNumberColumns]=-COIN_DBL_MAX;
    633     columnUpper2[newNumberColumns]=COIN_DBL_MAX;
    634     int j;
    635     for (j=rowStart[iRow];
    636          j<rowStart[iRow]+rowLength[iRow];
    637          j++) {
    638       double elementValue=elementByRow[j];
    639       int jColumn = column[j];
    640       elements2[numberElements]=elementValue;
    641       row2[numberElements++]=jColumn+numberRows;
    642       dj[jColumn]-= value*elementValue;
    643     }
    644     newNumberColumns++;
    645     start2[newNumberColumns]=numberElements;
    646   }
    647   // djs
    648   for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
    649     columnLower2[newNumberColumns]=-COIN_DBL_MAX;
    650     columnUpper2[newNumberColumns]=COIN_DBL_MAX;
    651     solution2[newNumberColumns]=dj[iColumn];
    652     elements2[numberElements]=1.0;
    653     row2[numberElements++]=back[iColumn]+numberRows;
    654     newNumberColumns++;
    655     start2[newNumberColumns]=numberElements;
    656   }
    657   // Create model
    658   ClpSimplex model2(*this);
    659   model2.resize(0,0);
    660   model2.loadProblem(newNumberColumns,newNumberRows,
    661                      start2,row2, elements2,
    662                      columnLower2,columnUpper2,
    663                      objective2,
    664                      rowLower2,rowUpper2);
    665   delete [] objective2;
    666   delete [] rowLower2;
    667   delete [] rowUpper2;
    668   delete [] columnLower2;
    669   delete [] columnUpper2;
    670   // Now create expanded quadratic objective for use in primalRow
    671   // Later pack down in some way
    672   start2[0]=0;
    673   numberElements=0;
    674   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    675     // Quadratic
    676     int j;
    677     for (j=columnQuadraticStart[iColumn];
    678          j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
    679          j++) {
    680       int jColumn = columnQuadratic[j];
    681       double value = quadraticElement[j];
    682       if (iColumn!=jColumn)
    683         value *= 0.5;
    684       elements2[numberElements]=value;
    685       row2[numberElements++]=jColumn;
    686     }
    687     for (j=rowStartQ[iColumn];
    688          j<rowStartQ[iColumn]+rowLengthQ[iColumn];
    689          j++) {
    690       int jColumn = columnQ[j];
    691       double value = elementByRowQ[j];
    692       if (iColumn!=jColumn) {
    693         value *= 0.5;
    694         elements2[numberElements]=value;
    695         row2[numberElements++]=jColumn;
    696       }
    697     }
    698     start2[iColumn+1]=numberElements;
    699   }
    700   // and pad
    701   for (;iColumn<newNumberColumns;iColumn++)
    702     start2[iColumn+1]=numberElements;
    703   // Load up objective
    704   model2.loadQuadraticObjective(newNumberColumns,start2,row2,elements2);
    705   delete [] start2;
    706   delete [] row2;
    707   delete [] elements2;
    708   model2.allSlackBasis();
    709   model2.scaling(false);
    710   model2.setLogLevel(this->logLevel());
    711   // Move solution across
    712   memcpy(model2.primalColumnSolution(),solution2,
    713          newNumberColumns*sizeof(double));
    714   columnLower2 = model2.columnLower();
    715   columnUpper2 = model2.columnUpper();
    716   delete [] solution2;
    717   solution2 = model2.primalColumnSolution();
    718   // Compute row activities and check feasible
    719   double * rowSolution2 = model2.primalRowSolution();
    720   memset(rowSolution2,0,newNumberRows*sizeof(double));
    721   model2.times(1.0,solution2,rowSolution2);
    722   rowLower2 = model2.rowLower();
    723   rowUpper2 = model2.rowUpper();
    724 #if 0
    725   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    726     Status xStatus = getColumnStatus(iColumn);
    727     bool isSuperBasic;
    728     int iS = iColumn+newNumberRows;
    729     double value = solution2[iS];
    730     if (fabs(value)>dualTolerance_)
    731       isSuperBasic=true;
    732     else
    733       isSuperBasic=false;
    734     // For moment take all x out of basis
    735     // Does not seem right
    736     isSuperBasic=true;
    737     model2.setColumnStatus(iColumn,xStatus);
    738     if (xStatus==basic) {
    739       if (!isSuperBasic) {
    740         model2.setRowStatus(numberRows+iColumn,basic);
    741         model2.setColumnStatus(iS,superBasic);
    742       } else {
    743         model2.setRowStatus(numberRows+iColumn,isFixed);
    744         model2.setColumnStatus(iS,basic);
    745         model2.setColumnStatus(iColumn,superBasic);
    746       }
    747     } else {
    748       model2.setRowStatus(numberRows+iColumn,isFixed);
    749       model2.setColumnStatus(iS,basic);
    750     }
    751   }
    752   for (iRow=0;iRow<numberRows;iRow++) {
    753     Status rowStatus = getRowStatus(iRow);
    754     model2.setRowStatus(iRow,rowStatus);
    755     if (rowStatus!=basic) {
    756       model2.setColumnStatus(iRow+numberColumns,basic); // make dual basic
    757     }
    758     assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
    759     assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
    760   }
    761   // why ?? - take duals out and adjust
    762   for (iRow=0;iRow<numberRows;iRow++) {
    763     model2.setRowStatus(iRow,basic);
    764     model2.setColumnStatus(iRow+numberColumns,superBasic);
    765     solution2[iRow+numberColumns]=0.0;
    766   }
    767 #else
    768   for (iRow=0;iRow<numberRows;iRow++) {
    769     assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
    770     assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
    771     model2.setRowStatus(iRow,basic);
    772     model2.setColumnStatus(iRow+numberColumns,superBasic);
    773     solution2[iRow+numberColumns]=0.0;
    774   }
    775   for (iColumn=numberRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
    776     model2.setColumnStatus(iColumn,basic);
    777     model2.setRowStatus(iColumn-numberColumns,isFixed);
    778   }
    779 #endif
    780   memset(rowSolution2,0,newNumberRows*sizeof(double));
    781   model2.times(1.0,solution2,rowSolution2);
    782   for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
    783     int iS = back[iColumn]+newNumberRows;
    784     int iRow = iColumn+numberRows;
    785     double value = rowSolution2[iRow];
    786     if (value>rowUpper2[iRow]) {
    787       rowSolution2[iRow] = rowUpper2[iRow];
    788       solution2[iS]-=value-rowUpper2[iRow];
    789     } else {
    790       rowSolution2[iRow] = rowLower2[iRow];
    791       solution2[iS]-=value-rowLower2[iRow];
    792     }
    793   }
    794 
    795  
    796   // See if any s sub j have wrong sign and/or use djs from infeasibility objective
    797   double objectiveOffset;
    798   getDblParam(ClpObjOffset,objectiveOffset);
    799   double objValue = -objectiveOffset;
    800   for (iColumn=0;iColumn<numberColumns;iColumn++)
    801     objValue += objective[iColumn]*solution2[iColumn];
    802   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    803     double valueI = solution2[iColumn];
    804     int j;
    805     for (j=columnQuadraticStart[iColumn];
    806          j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    807       int jColumn = columnQuadratic[j];
    808       double valueJ = solution2[jColumn];
    809       double elementValue = quadraticElement[j];
    810       objValue += 0.5*valueI*valueJ*elementValue;
    811     }
    812   }
    813   printf("Objective value %g\n",objValue);
    814   for (iColumn=0;iColumn<newNumberColumns;iColumn++)
    815     printf("%d %g\n",iColumn,solution2[iColumn]);
    816 #if 1
    817   CoinMpsIO writer;
    818   writer.setMpsData(*model2.matrix(), COIN_DBL_MAX,
    819                     model2.getColLower(), model2.getColUpper(),
    820                     model2.getObjCoefficients(),
    821                     (const char*) 0 /*integrality*/,
    822                     model2.getRowLower(), model2.getRowUpper(),
    823                     NULL,NULL);
    824   writer.writeMps("xx.mps");
    825 #endif 
    826   // Now do quadratic
    827   // If we did not do Slp we should have primal feasible basic solution
    828   // Do safe cast as no data
    829   ClpSimplexPrimalQuadratic * modelPtr =
    830     (ClpSimplexPrimalQuadratic *) (&model2);
    831   ClpPrimalQuadraticDantzig dantzigP(modelPtr,numberRows);
    832   modelPtr->setPrimalColumnPivotAlgorithm(dantzigP);
    833   modelPtr->messageHandler()->setLogLevel(63);
    834   modelPtr->primalQuadratic2(this,phase);
    835   memcpy(dualRowSolution(),model2.primalColumnSolution()+numberColumns_,numberRows_*sizeof(double));
    836   memcpy(primalColumnSolution(),model2.primalColumnSolution(),numberColumns_*sizeof(double));
    837   memset(model2.primalRowSolution(),0,newNumberRows*sizeof(double));
    838   model2.times(1.0,model2.primalColumnSolution(),
    839                model2.primalRowSolution());
    840   memcpy(dualColumnSolution(),model2.primalColumnSolution()+numberRows_+numberColumns_,
    841          numberColumns_*sizeof(double));
    842   objValue = -objectiveOffset;
    843   for (iColumn=0;iColumn<numberColumns;iColumn++)
    844     objValue += objective[iColumn]*solution2[iColumn];
    845   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    846     double valueI = solution2[iColumn];
    847     if (fabs(valueI)>1.0e-5) {
    848       int djColumn = iColumn+numberRows+numberColumns;
    849       assert(solution2[djColumn]<1.0e-7);
    850     }
    851     int j;
    852     for (j=columnQuadraticStart[iColumn];
    853          j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    854       int jColumn = columnQuadratic[j];
    855       double valueJ = solution2[jColumn];
    856       double elementValue = quadraticElement[j];
    857       objValue += 0.5*valueI*valueJ*elementValue;
    858     }
    859   }
    860   printf("Objective value %g\n",objValue);
    861   objectiveValue_ = objValue + objectiveOffset;
    862   return 0;
    863 #else
    864479  // Get a feasible solution
    865480  if (numberPrimalInfeasibilities())
     
    887502  endQuadratic(model2,info);
    888503  return 0;
    889 #endif
    890504}
    891505int ClpSimplexPrimalQuadratic::primalQuadratic2 (ClpQuadraticInfo * info,
     
    906520  const int * which = info->quadraticSequence();
    907521  //const int * back = info->backSequence();
    908   // Save solution
    909   double * saveSolution = new double [numberRows_+numberColumns_];
    910522  assert (!scalingFlag_);
    911   memcpy(saveSolution,columnActivity_,numberColumns_*sizeof(double));
    912   memcpy(saveSolution+numberColumns_,rowActivity_,numberRows_*sizeof(double));
    913523  // initialize - values pass coding and algorithm_ is +1
    914524  if (!startup(1)) {
    915525
    916526    // Setup useful stuff
    917     //ClpQuadraticInfo info(originalModel);
    918     if (phase)
    919       memcpy(solution_,saveSolution,
    920              (numberRows_+numberColumns_)*sizeof(double));
     527    info->setCurrentPhase(phase);
    921528    int lastCleaned=0; // last time objective or bounds cleaned up
    922     int sequenceIn=-1;
     529    info->setSequenceIn(-1);
    923530    info->setCrucialSj(-1);
    924531   
    925     // special nonlinear cost
    926     delete nonLinearCost_;
    927     nonLinearCost_ = new ClpNonLinearCost(this,info->numberXColumns());
    928532    // Say no pivot has occurred (for steepest edge and updates)
    929533    pivotRow_=-2;
     
    959563      if (lastGoodIteration_==numberIterations_)
    960564        factorType=3;
    961       if (phase)
    962         memcpy(saveSolution,solution_,
    963                (numberRows_+numberColumns_)*sizeof(double));
    964       if (phase&&0) {
    965         // Clean solution
    966         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    967           if (getColumnStatus(iColumn)==isFree)
    968             solution_[iColumn]=0.0;
    969         }
    970       } 
     565      if (info->currentPhase())
     566        info->setCurrentSolution(solution_);
    971567      // may factorize, checks if problem finished
    972       statusOfProblemInPrimal(lastCleaned,factorType,progress_);
    973       if (phase)
    974         memcpy(solution_,saveSolution,
     568      statusOfProblemInPrimal(lastCleaned,factorType,progress_,info);
     569      if (info->currentPhase())
     570        memcpy(solution_,info->currentSolution(),
    975571               (numberRows_+numberColumns_)*sizeof(double));
    976572     
     
    1004600      // Check problem phase
    1005601      // We assume all X are feasible
    1006       phase=0;
    1007       if (saveSolution) {
    1008         sequenceIn=-1;
     602      if (info->currentSolution()&&info->sequenceIn()<0) {
     603        phase=0;
     604        int nextSequenceIn=-1;
    1009605        for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    1010606          if (which[iColumn]>=0) {
     
    1012608            double lower = lower_[iColumn];
    1013609            double upper = upper_[iColumn];
    1014             if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
     610            if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
     611                ||getColumnStatus(iColumn)==superBasic) {
    1015612              if (getColumnStatus(iColumn)!=basic) {
    1016613                if (phase!=2) {
    1017614                  phase=2;
    1018                   sequenceIn=iColumn;
     615                  nextSequenceIn=iColumn;
    1019616                }
    1020617              }
     
    1026623                if (phase==0) {
    1027624                  phase=1;
    1028                   sequenceIn=iS;
     625                  nextSequenceIn=iS;
    1029626                }
    1030627              }
     
    1037634          double lower = lower_[iColumn];
    1038635          double upper = upper_[iColumn];
    1039           if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
     636          if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
     637            ||getColumnStatus(iColumn)==superBasic) {
    1040638            if (getColumnStatus(iColumn)!=basic) {
    1041639              if (phase!=2) {
    1042640                phase=2;
    1043                 sequenceIn=iColumn;
     641                nextSequenceIn=iColumn;
    1044642              }
    1045643            }
     
    1051649              if (phase==0) {
    1052650                phase=1;
    1053                 sequenceIn=iS;
     651                nextSequenceIn=iS;
    1054652              }
    1055653            }
    1056654          }
    1057655        }
    1058       }
    1059       if (!phase) {
    1060         delete [] saveSolution;
    1061         saveSolution=NULL;
     656        info->setSequenceIn(nextSequenceIn);
     657        info->setCurrentPhase(phase);
    1062658      }
    1063659
    1064660      // exit if victory declared
    1065       if (!phase&&primalColumnPivot_->pivotColumn(rowArray_[1],
     661      if (!info->currentPhase()&&primalColumnPivot_->pivotColumn(rowArray_[1],
    1066662                                          rowArray_[2],rowArray_[3],
    1067663                                          columnArray_[0],
     
    1073669      // Iterate
    1074670      problemStatus_=-1;
    1075       whileIterating(sequenceIn,info,phase);
     671      whileIterating(info);
    1076672    }
    1077673  }
    1078674  // clean up
    1079   delete [] saveSolution;
    1080675  finish();
    1081676  restoreData(data);
     
    1093688*/
    1094689int
    1095 ClpSimplexPrimalQuadratic::whileIterating(int & sequenceIn,
    1096                       ClpQuadraticInfo * info,
    1097                       int phase)
     690ClpSimplexPrimalQuadratic::whileIterating(
     691                      ClpQuadraticInfo * info)
    1098692{
    1099   checkComplementary (info);
     693  checkComplementarity (info);
    1100694  int crucialSj = info->crucialSj();
    1101695  int returnCode=-1;
     696  int phase = info->currentPhase();
    1102697  double saveObjective = objectiveValue_;
    1103698  int numberXColumns = info->numberXColumns();
     
    1109704  const double * element =  matrix_->getElements();
    1110705  const int * columnLength = matrix_->getVectorLengths();
    1111   int oldSequenceIn=sequenceIn;
     706  int nextSequenceIn=info->sequenceIn();
     707  int oldSequenceIn=nextSequenceIn;
    1112708  // status stays at -1 while iterating, >=0 finished, -2 to invert
    1113709  // status -3 to go to top without an invert
     
    1136732    if (phase==2) {
    1137733      // values pass
    1138       if (sequenceIn<0) {
     734      if (nextSequenceIn<0) {
    1139735        // get next
    1140736        int iColumn;
     
    1146742          if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
    1147743            if (getColumnStatus(iColumn)!=basic) {
    1148               sequenceIn=iColumn;
     744              nextSequenceIn=iColumn;
    1149745              break;
    1150746            }
    1151747          }
    1152748        }
    1153         if (sequenceIn<0) {
     749        if (nextSequenceIn<0) {
    1154750          iStart=max(iStart,numberColumns_);
    1155751          int numberXRows = info->numberXRows();
     
    1161757            if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
    1162758              if (getColumnStatus(iColumn)!=basic) {
    1163                 sequenceIn=iColumn;
     759                nextSequenceIn=iColumn;
    1164760                break;
    1165761              }
     
    1168764        }
    1169765      }
    1170       oldSequenceIn=sequenceIn;
    1171       sequenceIn_ = sequenceIn;
     766      oldSequenceIn=nextSequenceIn;
     767      sequenceIn_ = nextSequenceIn;
    1172768      cleanupIteration=false;
    1173769      if (sequenceIn_<numberXColumns) {
     
    1185781        }
    1186782      } else {
    1187         dualIn_ = - solution_[sequenceIn_-numberColumns_+numberXColumns];
     783        dualIn_ = solution_[sequenceIn_-numberColumns_+numberXColumns];
    1188784      }
    1189785      valueIn_=solution_[sequenceIn_];
     
    1193789        directionIn_ = 1;
    1194790    } else {
    1195       if (sequenceIn<0) {
     791      if (nextSequenceIn<0) {
    1196792        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
    1197793                     columnArray_[0],columnArray_[1]);
    1198794        cleanupIteration=false;
    1199795      } else {
    1200         sequenceIn_ = sequenceIn;
     796        sequenceIn_ = nextSequenceIn;
    1201797        cleanupIteration=true;
    1202798      }
     
    1206802    rowArray_[1]->clear();
    1207803    if (sequenceIn_>=0) {
    1208       sequenceIn=-1;
     804      if (numberIterations_>=8796000)
     805        printf("loxx 4721 %g %g\n",lower_[4721],upper_[4721]);
     806      nextSequenceIn=-1;
    1209807      // we found a pivot column
    1210808      int chosen = sequenceIn_;
     
    1258856
    1259857          dualIn_=dj_[sequenceIn_];
    1260           if (sequenceIn_>numberXColumns&&
    1261               sequenceIn_<-numberColumns_) {
    1262             // We can let flip to 0.0
    1263             assert (cleanupIteration);
    1264             if (solution_[sequenceIn_]>0.0) {
    1265               nonLinearCost_->setBounds(sequenceIn_, 0.0,
    1266                                         0.5*COIN_DBL_MAX);
    1267               lower_[sequenceIn_]=0.0;
    1268               upper_[sequenceIn_]=0.5*COIN_DBL_MAX;
    1269             } else {
    1270               nonLinearCost_->setBounds(sequenceIn_, -0.5*COIN_DBL_MAX,
    1271                                         0.0);
    1272               lower_[sequenceIn_]=-0.5*COIN_DBL_MAX;
    1273               upper_[sequenceIn_]=0.0;
    1274             }
    1275           } else {
    1276             // Let go through bounds
    1277             nonLinearCost_->setBounds(sequenceIn_, -0.5*COIN_DBL_MAX,
    1278                                       0.5*COIN_DBL_MAX);
    1279             lower_[sequenceIn_]=-0.5*COIN_DBL_MAX;
    1280             upper_[sequenceIn_]=0.5*COIN_DBL_MAX;
    1281           }
    1282858          lowerIn_=lower_[sequenceIn_];
    1283859          upperIn_=upper_[sequenceIn_];
     
    1314890                                                           pivotRow_,
    1315891                                                           alpha_);
     892          if (result>=10) {
     893            updateStatus = max(updateStatus,result/10);
     894            result = result%10;
     895            if (updateStatus>1)
     896              alpha_=0.0; // force re-factorization
     897          }
    1316898          // if no pivots, bad update but reasonable alpha - take and invert
    1317899          if (updateStatus==2&&
     
    1334916            }
    1335917            // later we may need to unwind more e.g. fake bounds
    1336             if(lastGoodIteration_ != numberIterations_) {
     918            if(lastGoodIteration_ != numberIterations_||factorization_->pivots()) {
    1337919              rowArray_[1]->clear();
    1338920              pivotRow_=-1;
    1339921              returnCode=-4;
    1340922              // retry on return
    1341               sequenceIn = sequenceIn_;
     923              nextSequenceIn = sequenceIn_;
    1342924              break;
    1343925            } else {
     
    14241006          double lowerValue = lower_[sequenceOut_];
    14251007          double upperValue = upper_[sequenceOut_];
     1008          if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
     1009            // really free but going to zero
     1010            lowerValue=0.0;
     1011            upperValue=0.0;
     1012          }
    14261013          assert(valueOut_>=lowerValue-primalTolerance_&&
    14271014                 valueOut_<=upperValue+primalTolerance_);
     
    14481035        nonLinearCost_->setOne(sequenceIn_,valueIn_);
    14491036        int whatNext=housekeeping(objectiveChange);
    1450         checkComplementary (info);
     1037        if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
     1038          // really free but going to zero
     1039          setStatus(sequenceOut_,isFree);
     1040        }
     1041        checkComplementarity (info);
    14511042
    14521043        if (whatNext==1) {
     
    15091100    }
    15101101  }
     1102  info->setSequenceIn(nextSequenceIn);
    15111103  return returnCode;
    15121104}
     
    17071299  coeff2 *= 0.5;
    17081300  printf("coefficients %g %g - dualIn %g\n",coeff1,coeff2,dualIn_);
    1709   if (!cleanupIteration)
    1710     assert (fabs(way*coeff1-dualIn_)<1.0e-4);
     1301  int accuracyFlag=0;
     1302  if (!cleanupIteration) {
     1303    if (way*coeff1*dualIn_<0.0) {
     1304      // bad
     1305      accuracyFlag=2;
     1306    } else if (fabs(way*coeff1-dualIn_)>1.0e-2*(1.0e-3+fabs(dualIn_))) {
     1307      // bad
     1308      accuracyFlag=2;
     1309    } else if (fabs(way*coeff1-dualIn_)>1.0e-4*(1.0+fabs(dualIn_))) {
     1310      // not wondeful
     1311      accuracyFlag=1;
     1312    }
     1313    dualIn_ = way*coeff1;
     1314    if (crucialSj>=0) {
     1315      solution_[crucialSj]=dualIn_;
     1316    }
     1317  }
    17111318  // interesting places are where derivative zero or sj goes to zero
    17121319  double d1,d2=1.0e50;
    17131320  if (fabs(coeff2)>1.0e-9)
    17141321    d1 = - 0.5*coeff1/coeff2;
    1715   else if (coeff1<=1.0e-9)
     1322  else if (coeff1<=1.0e-6)
    17161323    d1 = maximumMovement;
    17171324  else
     
    20211628      upperOut_ = upperIn_;
    20221629      alpha_ = 0.0;
     1630      if (accuracyFlag<2)
     1631        accuracyFlag=0;
    20231632      if (way<0.0) {
    20241633        directionOut_=1;      // to lower bound
     
    20331642      result=0;
    20341643      assert (pivotRow_<0);
    2035       nonLinearCost_->setBounds(crucialSj, 0.0,0.0);
    2036       lower_[crucialSj]=0.0;
    2037       upper_[crucialSj]=0.0;
    2038       setColumnStatus(crucialSj,isFixed);
     1644      setColumnStatus(crucialSj,isFree);
    20391645      pivotRow_ = iSjRow;
    20401646      alpha_ = work[pivotRow_];
     
    20421648      sequenceOut_ = pivotVariable_[pivotRow_];
    20431649      valueOut_ = solution(sequenceOut_);
    2044       lowerOut_=lower_[sequenceOut_];
    2045       upperOut_=upper_[sequenceOut_];
    20461650      theta_ = d2;
    20471651      if (way<0.0)
    20481652        theta_ = - theta_;
    2049       double newValue = valueOut_ - theta_*alpha_;
    2050       if (alpha_*way<0.0) {
    2051         directionOut_=-1;      // to upper bound
    2052         if (fabs(theta_)>1.0e-6)
    2053           upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
    2054         else
    2055           upperOut_ = newValue;
    2056       } else {
    2057         directionOut_=1;      // to lower bound
    2058         if (fabs(theta_)>1.0e-6)
    2059           lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
    2060         else
    2061           lowerOut_ = newValue;
    2062       }
     1653      lowerOut_=0.0;
     1654      upperOut_=0.0;
    20631655      //????
    20641656      dualOut_ = reducedCost(sequenceOut_);
     
    21071699    printf("Objective value %g\n",objectiveValue());
    21081700  }
    2109   return result;
     1701  return result+10*accuracyFlag;
    21101702}
    2111 // Just for debug
     1703/* Checks if finished.  Updates status */
     1704void
     1705ClpSimplexPrimalQuadratic::statusOfProblemInPrimal(int & lastCleaned,int type,
     1706                                          ClpSimplexProgress * progress,
     1707                                                   ClpQuadraticInfo * info)
     1708{
     1709  if (type==2) {
     1710    // trouble - restore solution
     1711    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
     1712    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
     1713           numberRows_*sizeof(double));
     1714    memcpy(columnActivityWork_,savedSolution_ ,
     1715           numberColumns_*sizeof(double));
     1716    forceFactorization_=1; // a bit drastic but ..
     1717    pivotRow_=-1; // say no weights update
     1718    changeMade_++; // say change made
     1719    info->restoreStatus();
     1720  }
     1721  int saveThreshold = factorization_->sparseThreshold();
     1722  int tentativeStatus = problemStatus_;
     1723  if (problemStatus_>-3||problemStatus_==-4) {
     1724    // factorize
     1725    // later on we will need to recover from singularities
     1726    // also we could skip if first time
     1727    // do weights
     1728    // This may save pivotRow_ for use
     1729    primalColumnPivot_->saveWeights(this,1);
     1730
     1731    if (type) {
     1732      // is factorization okay?
     1733      if (internalFactorize(1)) {
     1734        if (solveType_==2) {
     1735          // say odd
     1736          problemStatus_=5;
     1737          return;
     1738        }
     1739
     1740        // no - restore previous basis
     1741        assert (type==1);
     1742        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
     1743        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
     1744               numberRows_*sizeof(double));
     1745        memcpy(columnActivityWork_,savedSolution_ ,
     1746               numberColumns_*sizeof(double));
     1747        forceFactorization_=1; // a bit drastic but ..
     1748        type = 2;
     1749        assert (internalFactorize(1)==0);
     1750        info->restoreStatus();
     1751        changeMade_++; // say change made
     1752      }
     1753    }
     1754    if (problemStatus_!=-4)
     1755      problemStatus_=-3;
     1756  }
     1757  // at this stage status is -3 or -5 if looks unbounded
     1758  // get primal and dual solutions
     1759  // put back original costs and then check
     1760  createRim(4);
     1761  gutsOfSolution(NULL,NULL);
     1762  checkComplementarity(info);
     1763  // Double check reduced costs if no action
     1764  if (progress->lastIterationNumber(0)==numberIterations_) {
     1765    if (primalColumnPivot_->looksOptimal()) {
     1766      numberDualInfeasibilities_ = 0;
     1767      sumDualInfeasibilities_ = 0.0;
     1768    }
     1769  }
     1770  // Check if looping
     1771  int loop;
     1772  if (type!=2)
     1773    loop = progress->looping();
     1774  else
     1775    loop=-1;
     1776  if (loop>=0) {
     1777    problemStatus_ = loop; //exit if in loop
     1778    return ;
     1779  } else if (loop<-1) {
     1780    // something may have changed
     1781    gutsOfSolution(NULL,NULL);
     1782    checkComplementarity(info);
     1783  }
     1784  // really for free variables in
     1785  //if((progressFlag_&2)!=0)
     1786  //problemStatus_=-1;;
     1787  progressFlag_ = 0; //reset progress flag
     1788
     1789  handler_->message(CLP_SIMPLEX_STATUS,messages_)
     1790    <<numberIterations_<<objectiveValue();
     1791  handler_->printing(sumPrimalInfeasibilities_>0.0)
     1792    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
     1793  handler_->printing(sumDualInfeasibilities_>0.0)
     1794    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
     1795  handler_->printing(numberDualInfeasibilitiesWithoutFree_
     1796                     <numberDualInfeasibilities_)
     1797                       <<numberDualInfeasibilitiesWithoutFree_;
     1798  handler_->message()<<CoinMessageEol;
     1799  assert (primalFeasible());
     1800  // we may wish to say it is optimal even if infeasible
     1801  bool alwaysOptimal = (specialOptions_&1)!=0;
     1802  // give code benefit of doubt
     1803  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
     1804      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
     1805    // say optimal (with these bounds etc)
     1806    numberDualInfeasibilities_ = 0;
     1807    sumDualInfeasibilities_ = 0.0;
     1808    numberPrimalInfeasibilities_ = 0;
     1809    sumPrimalInfeasibilities_ = 0.0;
     1810  }
     1811  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
     1812  if (dualFeasible()||problemStatus_==-4) {
     1813    if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
     1814      //may need infeasiblity cost changed
     1815      // we can see if we can construct a ray
     1816      // make up a new objective
     1817      double saveWeight = infeasibilityCost_;
     1818      // save nonlinear cost as we are going to switch off costs
     1819      ClpNonLinearCost * nonLinear = nonLinearCost_;
     1820      // do twice to make sure Primal solution has settled
     1821      // put non-basics to bounds in case tolerance moved
     1822      // put back original costs
     1823      createRim(4);
     1824      nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1825      gutsOfSolution(NULL,NULL);
     1826      checkComplementarity(info);
     1827
     1828      infeasibilityCost_=1.0e100;
     1829      // put back original costs
     1830      createRim(4);
     1831      nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1832      // may have fixed infeasibilities - double check
     1833      if (nonLinearCost_->numberInfeasibilities()==0) {
     1834        // carry on
     1835        problemStatus_ = -1;
     1836        infeasibilityCost_=saveWeight;
     1837      } else {
     1838        nonLinearCost_=NULL;
     1839        // scale
     1840        int i;
     1841        for (i=0;i<numberRows_+numberColumns_;i++)
     1842          cost_[i] *= 1.0e-100;
     1843        gutsOfSolution(NULL,NULL);
     1844        checkComplementarity(info);
     1845        nonLinearCost_=nonLinear;
     1846        infeasibilityCost_=saveWeight;
     1847        if ((infeasibilityCost_>=1.0e18||
     1848             numberDualInfeasibilities_==0)&&perturbation_==101) {
     1849          unPerturb(); // stop any further perturbation
     1850          numberDualInfeasibilities_=1; // carry on
     1851        }
     1852        if (infeasibilityCost_>=1.0e20||
     1853            numberDualInfeasibilities_==0) {
     1854          // we are infeasible - use as ray
     1855          delete [] ray_;
     1856          ray_ = new double [numberRows_];
     1857          memcpy(ray_,dual_,numberRows_*sizeof(double));
     1858          // and get feasible duals
     1859          infeasibilityCost_=0.0;
     1860          createRim(4);
     1861          nonLinearCost_->checkInfeasibilities(primalTolerance_);
     1862          gutsOfSolution(NULL,NULL);
     1863          checkComplementarity(info);
     1864          // so will exit
     1865          infeasibilityCost_=1.0e30;
     1866          // reset infeasibilities
     1867          sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
     1868          numberPrimalInfeasibilities_=
     1869            nonLinearCost_->numberInfeasibilities();
     1870        }
     1871        if (infeasibilityCost_<1.0e20) {
     1872          infeasibilityCost_ *= 5.0;
     1873          changeMade_++; // say change made
     1874          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
     1875            <<infeasibilityCost_
     1876            <<CoinMessageEol;
     1877          // put back original costs and then check
     1878          createRim(4);
     1879          nonLinearCost_->checkInfeasibilities(0.0);
     1880          gutsOfSolution(NULL,NULL);
     1881          checkComplementarity(info);
     1882          problemStatus_=-1; //continue
     1883        } else {
     1884          // say infeasible
     1885          problemStatus_ = 1;
     1886        }
     1887      }
     1888    } else {
     1889      // may be optimal
     1890      if (perturbation_==101) {
     1891        unPerturb(); // stop any further perturbation
     1892        lastCleaned=-1; // carry on
     1893      }
     1894      if ( lastCleaned!=numberIterations_||unflag()) {
     1895        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
     1896          <<primalTolerance_
     1897          <<CoinMessageEol;
     1898        if (numberTimesOptimal_<4) {
     1899          numberTimesOptimal_++;
     1900          changeMade_++; // say change made
     1901          if (numberTimesOptimal_==1) {
     1902            // better to have small tolerance even if slower
     1903            factorization_->zeroTolerance(1.0e-15);
     1904          }
     1905          lastCleaned=numberIterations_;
     1906          if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
     1907            handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
     1908              <<CoinMessageEol;
     1909          double oldTolerance = primalTolerance_;
     1910          primalTolerance_=dblParam_[ClpPrimalTolerance];
     1911          // put back original costs and then check
     1912          createRim(4);
     1913          nonLinearCost_->checkInfeasibilities(oldTolerance);
     1914          gutsOfSolution(NULL,NULL);
     1915          checkComplementarity(info);
     1916          problemStatus_ = -1;
     1917        } else {
     1918          problemStatus_=0; // optimal
     1919          if (lastCleaned<numberIterations_) {
     1920            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
     1921              <<CoinMessageEol;
     1922          }
     1923        }
     1924      } else {
     1925        problemStatus_=0; // optimal
     1926      }
     1927    }
     1928  } else {
     1929    // see if looks unbounded
     1930    if (problemStatus_==-5) {
     1931      if (nonLinearCost_->numberInfeasibilities()) {
     1932        if (infeasibilityCost_>1.0e18&&perturbation_==101) {
     1933          // back off weight
     1934          infeasibilityCost_ = 1.0e13;
     1935          unPerturb(); // stop any further perturbation
     1936        }
     1937        //we need infeasiblity cost changed
     1938        if (infeasibilityCost_<1.0e20) {
     1939          infeasibilityCost_ *= 5.0;
     1940          changeMade_++; // say change made
     1941          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
     1942            <<infeasibilityCost_
     1943            <<CoinMessageEol;
     1944          // put back original costs and then check
     1945          createRim(4);
     1946          gutsOfSolution(NULL,NULL);
     1947          checkComplementarity(info);
     1948          problemStatus_=-1; //continue
     1949        } else {
     1950          // say unbounded
     1951          problemStatus_ = 2;
     1952        }
     1953      } else {
     1954        // say unbounded
     1955        problemStatus_ = 2;
     1956      }
     1957    } else {
     1958      if(type==3&&problemStatus_!=-5)
     1959        unflag(); // odd
     1960      // carry on
     1961      problemStatus_ = -1;
     1962    }
     1963  }
     1964  if (type==0||type==1) {
     1965    if (type!=1||!saveStatus_) {
     1966      // create save arrays
     1967      delete [] saveStatus_;
     1968      delete [] savedSolution_;
     1969      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
     1970      savedSolution_ = new double [numberRows_+numberColumns_];
     1971    }
     1972    // save arrays
     1973    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
     1974    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
     1975           numberRows_*sizeof(double));
     1976    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
     1977    info->saveStatus();
     1978  }
     1979  // restore weights (if saved) - also recompute infeasibility list
     1980  if (tentativeStatus>-3)
     1981    primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
     1982  else
     1983    primalColumnPivot_->saveWeights(this,3);
     1984  if (problemStatus_<0&&!changeMade_) {
     1985    problemStatus_=4; // unknown
     1986  }
     1987  if (saveThreshold) {
     1988    // use default at present
     1989    factorization_->sparseThreshold(0);
     1990    factorization_->goSparse();
     1991  }
     1992  lastGoodIteration_ = numberIterations_;
     1993}
     1994// Fills in reduced costs
     1995void
     1996ClpSimplexPrimalQuadratic::createDjs (const ClpQuadraticInfo * info)
     1997{
     1998  const int * which = info->quadraticSequence();
     1999  // Matrix for linear stuff
     2000  const int * row = matrix_->getIndices();
     2001  const int * columnStart = matrix_->getVectorStarts();
     2002  const double * element =  matrix_->getElements();
     2003  const int * columnLength = matrix_->getVectorLengths();
     2004  int numberXColumns = info->numberXColumns();
     2005  int numberXRows = info->numberXRows();
     2006  const double * pi = solution_+numberXColumns;
     2007  int iSequence;
     2008  int start=numberXColumns+numberXRows;
     2009  for (iSequence=0;iSequence<numberXColumns;iSequence++) {
     2010    int jSequence = which[iSequence];
     2011    double value;
     2012    if (jSequence>=0) {
     2013      jSequence += start;
     2014      value = solution_[jSequence];
     2015    } else {
     2016      value=cost_[iSequence];
     2017      int j;
     2018      for (j=columnStart[iSequence];j<columnStart[iSequence]+columnLength[iSequence]; j++) {
     2019        int iRow = row[j];
     2020        value -= element[j]*pi[iRow];
     2021      }
     2022    }
     2023    dj_[iSequence]=value;
     2024  }
     2025  // And slacks
     2026  // Value of sj is - value of pi
     2027  int firstSlack = numberColumns_;
     2028  int jSequence;
     2029  for (jSequence=0;jSequence<numberXRows;jSequence++) {
     2030    int iSequence  = jSequence + firstSlack;
     2031    int iPi = jSequence+numberXColumns;
     2032    double value = solution_[iPi];
     2033    dj_[iSequence]=value;
     2034  }
     2035}
     2036
    21122037int
    2113 ClpSimplexPrimalQuadratic::checkComplementary (const
     2038ClpSimplexPrimalQuadratic::checkComplementarity (const
    21142039                                               ClpQuadraticInfo * info)
    21152040{
    2116   //const ClpSimplex * originalModel= info->originalModel();
     2041  createDjs(info);
    21172042  int numberXColumns = info->numberXColumns();
    2118   int i;
     2043  int numberXRows = info->numberXRows();
     2044  int iSequence;
     2045  int start=numberXColumns+numberXRows;
     2046  numberDualInfeasibilities_=0;
     2047  sumDualInfeasibilities_=0.0;
    21192048  const int * which = info->quadraticSequence();
    2120   for (i=0;i<numberXColumns;i++) {
    2121     if (which[i]>=0) {
    2122       int jSequence = i+ numberRows_;
    2123       if (getColumnStatus(i)==basic) {
     2049  for (iSequence=0;iSequence<numberXColumns;iSequence++) {
     2050    int jSequence = which[iSequence];
     2051    double value=dj_[iSequence];
     2052    ClpSimplex::Status status = getColumnStatus(iSequence);
     2053   
     2054    switch(status) {
     2055     
     2056    case ClpSimplex::basic:
     2057      if (jSequence>=0) {
     2058        jSequence += start;
    21242059        if (getColumnStatus(jSequence)==basic)
    21252060          printf("Struct %d (%g) and %d (%g) both basic\n",
    2126                  i,solution_[i],jSequence,solution_[jSequence]);
    2127       }
    2128     }
    2129   }
    2130   int originalNumberRows = info->numberXRows();
     2061                 iSequence,solution_[iSequence],jSequence,solution_[jSequence]);
     2062      }
     2063      break;
     2064    case ClpSimplex::isFixed:
     2065      break;
     2066    case ClpSimplex::isFree:
     2067    case ClpSimplex::superBasic:
     2068      if (fabs(value)>dualTolerance_) {
     2069        sumDualInfeasibilities_ += fabs(value)-dualTolerance_;
     2070        numberDualInfeasibilities_++;
     2071      }
     2072      break;
     2073    case ClpSimplex::atUpperBound:
     2074      if (value>dualTolerance_) {
     2075        sumDualInfeasibilities_ += value-dualTolerance_;
     2076        numberDualInfeasibilities_++;
     2077      }
     2078      break;
     2079    case ClpSimplex::atLowerBound:
     2080      if (value<-dualTolerance_) {
     2081        sumDualInfeasibilities_ -= value+dualTolerance_;
     2082        numberDualInfeasibilities_++;
     2083      }
     2084    }
     2085  }
    21312086  int offset = numberXColumns;
    2132   for (i=0;i<originalNumberRows;i++) {
    2133     int jSequence = i+offset;
    2134     if (getRowStatus(i)==basic) {
    2135       if (getColumnStatus(jSequence)==basic)
     2087  for (iSequence=0;iSequence<numberXRows;iSequence++) {
     2088    int jSequence = iSequence+offset;
     2089    double value=dj_[iSequence+numberColumns_];
     2090    ClpSimplex::Status status = getRowStatus(iSequence);
     2091   
     2092    switch(status) {
     2093     
     2094    case ClpSimplex::basic:
     2095      if (getRowStatus(jSequence)==basic)
    21362096        printf("Row %d (%g) and %d (%g) both basic\n",
    2137                i,solution_[i],jSequence,solution_[jSequence]);
    2138     }
    2139   }
    2140   return 0;
     2097               iSequence,solution_[iSequence],jSequence,solution_[jSequence]);
     2098      break;
     2099    case ClpSimplex::isFixed:
     2100      break;
     2101    case ClpSimplex::isFree:
     2102    case ClpSimplex::superBasic:
     2103      if (fabs(value)>dualTolerance_) {
     2104        sumDualInfeasibilities_ += fabs(value)-dualTolerance_;
     2105        numberDualInfeasibilities_++;
     2106      }
     2107      break;
     2108    case ClpSimplex::atUpperBound:
     2109      if (value>dualTolerance_) {
     2110        sumDualInfeasibilities_ += value-dualTolerance_;
     2111        numberDualInfeasibilities_++;
     2112      }
     2113      break;
     2114    case ClpSimplex::atLowerBound:
     2115      if (value<-dualTolerance_) {
     2116        sumDualInfeasibilities_ -= value+dualTolerance_;
     2117        numberDualInfeasibilities_++;
     2118      }
     2119    }
     2120  }
     2121  sumOfRelaxedDualInfeasibilities_ =sumDualInfeasibilities_;
     2122  return numberDualInfeasibilities_;
    21412123}
    21422124 
     
    23922374  rowUpper2 = model2->rowUpper();
    23932375#if 0
     2376  // redo as Dantzig says
    23942377  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    23952378    Status xStatus = getColumnStatus(iColumn);
    23962379    bool isSuperBasic;
    23972380    int iS = iColumn+newNumberRows;
    2398     double value = solution2[iS];
    2399     if (fabs(value)>dualTolerance_)
    2400       isSuperBasic=true;
    2401     else
    2402       isSuperBasic=false;
    2403     // For moment take all x out of basis
    2404     // Does not seem right
    2405     isSuperBasic=true;
    24062381    model2->setColumnStatus(iColumn,xStatus);
    24072382    if (xStatus==basic) {
    2408       if (!isSuperBasic) {
    2409         model2->setRowStatus(numberRows+iColumn,basic);
    2410         model2->setColumnStatus(iS,superBasic);
    2411       } else {
    2412         model2->setRowStatus(numberRows+iColumn,isFixed);
    2413         model2->setColumnStatus(iS,basic);
    2414         model2->setColumnStatus(iColumn,superBasic);
    2415       }
     2383      model2->setColumnStatus(iS,isFree);
     2384      solution2[iS]=0.0;
    24162385    } else {
    2417       model2->setRowStatus(numberRows+iColumn,isFixed);
    24182386      model2->setColumnStatus(iS,basic);
    24192387    }
     2388    // take slack out on equality rows
     2389    model2->setRowBasic(iColumn+numberRows,superBasic);
    24202390  }
    24212391  for (iRow=0;iRow<numberRows;iRow++) {
     
    25842554    quadraticSequence_(NULL),
    25852555    backSequence_(NULL),
     2556    currentSequenceIn_(-1),
    25862557    crucialSj_(-1),
     2558    validSequenceIn_(-1),
     2559    validCrucialSj_(-1),
     2560    currentPhase_(-1),
     2561    currentSolution_(NULL),
     2562    validPhase_(-1),
     2563    validSolution_(NULL),
    25872564    numberXRows_(-1),
    25882565    numberXColumns_(-1),
     
    25952572    quadraticSequence_(NULL),
    25962573    backSequence_(NULL),
     2574    currentSequenceIn_(-1),
    25972575    crucialSj_(-1),
     2576    validSequenceIn_(-1),
     2577    validCrucialSj_(-1),
     2578    currentPhase_(-1),
     2579    currentSolution_(NULL),
     2580    validPhase_(-1),
     2581    validSolution_(NULL),
    25982582    numberXRows_(-1),
    25992583    numberXColumns_(-1),
     
    26182602  delete [] quadraticSequence_;
    26192603  delete [] backSequence_;
     2604  delete [] currentSolution_;
     2605  delete [] validSolution_;
    26202606}
    26212607// Copy
     
    26242610    quadraticSequence_(NULL),
    26252611    backSequence_(NULL),
     2612    currentSequenceIn_(rhs.currentSequenceIn_),
    26262613    crucialSj_(rhs.crucialSj_),
     2614    validSequenceIn_(rhs.validSequenceIn_),
     2615    validCrucialSj_(rhs.validCrucialSj_),
     2616    currentPhase_(rhs.currentPhase_),
     2617    validPhase_(rhs.validPhase_),
    26272618    numberXRows_(rhs.numberXRows_),
    26282619    numberXColumns_(rhs.numberXColumns_),
     
    26362627    memcpy(backSequence_,rhs.backSequence_,
    26372628           numberXColumns_*sizeof(int));
     2629    int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     2630    int numberRows = numberXRows_+numberQuadraticColumns_;
     2631    int size = numberRows+numberColumns;
     2632    if (rhs.currentSolution_) {
     2633      currentSolution_ = new double [size];
     2634      memcpy(currentSolution_,rhs.currentSolution_,size*sizeof(double));
     2635    } else {
     2636      currentSolution_ = NULL;
     2637    }
     2638    if (rhs.validSolution_) {
     2639      validSolution_ = new double [size];
     2640      memcpy(validSolution_,rhs.validSolution_,size*sizeof(double));
     2641    } else {
     2642      validSolution_ = NULL;
     2643    }
    26382644  }
    26392645}
     
    26452651    originalModel_ = rhs.originalModel_;
    26462652    delete [] quadraticSequence_;
    2647     quadraticSequence_ = NULL;
    26482653    delete [] backSequence_;
    2649     backSequence_ = NULL;
     2654    delete [] currentSolution_;
     2655    delete [] validSolution_;
     2656    currentSequenceIn_ = rhs.currentSequenceIn_;
    26502657    crucialSj_ = rhs.crucialSj_;
     2658    validSequenceIn_ = rhs.validSequenceIn_;
     2659    validCrucialSj_ = rhs.validCrucialSj_;
     2660    currentPhase_ = rhs.currentPhase_;
     2661    validPhase_ = rhs.validPhase_;
    26512662    numberXRows_ = rhs.numberXRows_;
    26522663    numberXColumns_ = rhs.numberXColumns_;
     
    26592670      memcpy(backSequence_,rhs.backSequence_,
    26602671             numberXColumns_*sizeof(int));
     2672    } else {
     2673      quadraticSequence_ = NULL;
     2674      backSequence_ = NULL;
     2675    }
     2676    int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     2677    int numberRows = numberXRows_+numberQuadraticColumns_;
     2678    int size = numberRows+numberColumns;
     2679    if (rhs.currentSolution_) {
     2680      currentSolution_ = new double [size];
     2681      memcpy(currentSolution_,rhs.currentSolution_,size*sizeof(double));
     2682    } else {
     2683      currentSolution_ = NULL;
     2684    }
     2685    if (rhs.validSolution_) {
     2686      validSolution_ = new double [size];
     2687      memcpy(validSolution_,rhs.validSolution_,size*sizeof(double));
     2688    } else {
     2689      validSolution_ = NULL;
    26612690    }
    26622691  }
    26632692  return *this;
    26642693}
     2694// Save current In and Sj status
     2695 void
     2696ClpQuadraticInfo::saveStatus()
     2697{
     2698  validSequenceIn_ = currentSequenceIn_;
     2699  validCrucialSj_ = crucialSj_;
     2700  validPhase_ = currentPhase_;
     2701  int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     2702  int numberRows = numberXRows_+numberQuadraticColumns_;
     2703  int size = numberRows+numberColumns;
     2704  if (!validSolution_)
     2705    validSolution_ = new double [size];
     2706  memcpy(validSolution_,currentSolution_,size*sizeof(double));
     2707}
     2708// Restore previous
     2709void
     2710ClpQuadraticInfo::restoreStatus()
     2711{
     2712  currentSequenceIn_ = validSequenceIn_;
     2713  crucialSj_ = validCrucialSj_;
     2714  currentPhase_ = validPhase_;
     2715  delete [] currentSolution_;
     2716  currentSolution_ = validSolution_;
     2717  validSolution_=NULL;
     2718}
     2719void
     2720ClpQuadraticInfo::setCurrentSolution(const double * solution)
     2721{
     2722  if (currentPhase_) {
     2723    int numberColumns = numberXRows_+numberXColumns_+numberQuadraticColumns_;
     2724    int numberRows = numberXRows_+numberQuadraticColumns_;
     2725    int size = numberRows+numberColumns;
     2726    if (!currentSolution_)
     2727      currentSolution_ = new double [size];
     2728    memcpy(currentSolution_,solution,size*sizeof(double));
     2729  } else {
     2730    delete [] currentSolution_;
     2731    currentSolution_=NULL;
     2732  }
     2733}
  • branches/pre/Samples/Makefile

    r181 r185  
    5959#       LDFLAGS += -lhistory -lreadline -ltermcap
    6060endif
    61 LDFLAGS += -lefence
     61#LDFLAGS += -lefence
    6262###############################################################################
    6363
  • branches/pre/include/ClpNonLinearCost.hpp

    r180 r185  
    4141      This will just set up wasteful arrays for linear, but
    4242      later may do dual analysis and even finding duplicate columns .
    43       If for quadratic then free variables get extra 0,0 bits
    44       (flagged by numberOriginalColumns)
    45   */
    46   ClpNonLinearCost(ClpSimplex * model,int numberOriginalColumns=-1);
     43  */
     44  ClpNonLinearCost(ClpSimplex * model);
    4745  /** Constructor from simplex and list of non-linearities (columns only)
    4846      First lower of each column has to match real lower
     
    140138  inline double cost(int sequence) const
    141139  { return cost_[whichRange_[sequence]+offset_[sequence]];};
    142   /// Sets inside bounds (i.e. non infinite - used in QP
    143   void setBounds(int sequence, double lower, double upper);
    144140  //@}
    145141
  • branches/pre/include/ClpSimplexPrimalQuadratic.hpp

    r183 r185  
    5959  int endQuadratic(ClpSimplexPrimalQuadratic * quadraticModel,
    6060                   ClpQuadraticInfo & info);
    61   /// Just for debug
    62   int checkComplementary (const ClpQuadraticInfo * info);
     61  /// Checks complementarity and computes infeasibilities
     62  int checkComplementarity (const ClpQuadraticInfo * info);
     63  /// Fills in reduced costs
     64  void createDjs (const ClpQuadraticInfo * info);
    6365  /** Main part.
    6466      phase - 0 normal, 1 getting complementary solution,
    6567      2 getting basic solution. */
    66   int whileIterating (int & sequenceIn,
    67                       ClpQuadraticInfo * info,
    68                       int phase);
     68  int whileIterating (ClpQuadraticInfo * info);
    6969  /**
    7070      Row array has pivot column
     
    8484                ClpQuadraticInfo * info,
    8585                bool cleanupIteration);
     86  /**  Refactorizes if necessary
     87       Checks if finished.  Updates status.
     88       lastCleaned refers to iteration at which some objective/feasibility
     89       cleaning too place.
     90
     91       type - 0 initial so set up save arrays etc
     92            - 1 normal -if good update save
     93            - 2 restoring from saved
     94  */
     95  void statusOfProblemInPrimal(int & lastCleaned, int type,
     96                               ClpSimplexProgress * progress,
     97                               ClpQuadraticInfo * info);
    8698  //@}
    8799
     
    125137  inline int numberXRows() const
    126138  {return numberXRows_;};
     139  /// Sequence number incoming
     140  inline int sequenceIn() const
     141  {return currentSequenceIn_;};
     142  inline void setSequenceIn(int sequence)
     143  {currentSequenceIn_=sequence;};
    127144  /// Sequence number of binding Sj
    128145  inline int crucialSj() const
     
    130147  inline void setCrucialSj(int sequence)
    131148  {crucialSj_=sequence;};
     149  /// Current phase
     150  inline int currentPhase() const
     151  {return currentPhase_;};
     152  inline void setCurrentPhase(int phase)
     153  {currentPhase_=phase;};
     154  /// Current saved solution
     155  inline const double * currentSolution() const
     156  {return currentSolution_;};
     157  void setCurrentSolution(const double * solution);
    132158  /// Original objective
    133159  inline const double * originalObjective() const
     
    142168  inline const ClpSimplex * originalModel() const
    143169  { return originalModel_;};
     170  /// Save current In and Sj status
     171  void saveStatus();
     172  /// Restore previous
     173  void restoreStatus();
    144174  //@}
    145175   
     
    153183  /// Which ones are quadratic
    154184  int * backSequence_;
     185  /// Current sequenceIn
     186  int currentSequenceIn_;
    155187  /// Crucial Sj
    156188  int crucialSj_;
     189  /// Valid sequenceIn (for valid factorization)
     190  int validSequenceIn_;
     191  /// Valid crucial Sj
     192  int validCrucialSj_;
     193  /// Current phase
     194  int currentPhase_;
     195  /// Current saved solution
     196  double * currentSolution_;
     197  /// Valid phase
     198  int validPhase_;
     199  /// Valid saved solution
     200  double * validSolution_;
    157201  /// Number of original rows
    158202  int numberXRows_;
Note: See TracChangeset for help on using the changeset viewer.