Changeset 166


Ignore:
Timestamp:
Apr 28, 2003 9:55:54 AM (17 years ago)
Author:
forrest
Message:

Some fixes - when starting primitive QP coding

Location:
trunk
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpNonLinearCost.cpp

    r162 r166  
    4040   later may do dual analysis and even finding duplicate columns
    4141*/
    42 ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model)
     42ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model,bool forQuadratic)
    4343{
    4444  model_ = model;
     
    6868
    6969  for (iSequence=0;iSequence<numberTotal;iSequence++) {
     70    if (lower[iSequence]>-1.0e20)
     71      put++;
    7072    if (upper[iSequence]<1.0e20)
    7173      put++;
    72     put += 3;
     74    put += 2;
     75  }
     76
     77  // For quadratic we need -inf,0,0,+inf
     78  if (forQuadratic) {
     79    for (iSequence=0;iSequence<numberTotal;iSequence++) {
     80      if (lower[iSequence]<-1.0e20&&upper[iSequence]>1.0e20)
     81        put+=2;
     82    }
    7383  }
    7484
     
    8494  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    8595
    86     lower_[put] = -COIN_DBL_MAX;
    87     setInfeasible(put,true);
    88 
    89     cost_[put++] = cost[iSequence]-infeasibilityCost;
    90     whichRange_[iSequence]=put;
    91     lower_[put] = lower[iSequence];
    92     cost_[put++] = cost[iSequence];
    93     lower_[put] = upper[iSequence];
    94     cost_[put++] = cost[iSequence]+infeasibilityCost;
    95     if (upper[iSequence]<1.0e20) {
     96    if (!forQuadratic||lower[iSequence]>-1.0e20||upper[iSequence]<1.0e20) {
     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 free variable
     115      lower_[put] = -COIN_DBL_MAX;
     116      setInfeasible(put,true);
     117      cost_[put++] = -infeasibilityCost;
     118      whichRange_[iSequence]=put;
     119      lower_[put] = -0.5*COIN_DBL_MAX;
     120      cost_[put++] = 0.0;
     121      lower_[put] = 0.5*COIN_DBL_MAX;
     122      cost_[put++] = infeasibilityCost;
    96123      lower_[put] = COIN_DBL_MAX;
    97124      setInfeasible(put-1,true);
    98       cost_[put++] = 1.0e50;
    99     }
    100     start_[iSequence+1]=put;
     125      cost_[put++] = 0.0;
     126      start_[iSequence+1]=put;
     127    }
    101128  }
    102129
     
    730757  return lower_[jRange];
    731758}
    732 
     759// Sets inside bounds (i.e. non infinite - used in QP
     760void
     761ClpNonLinearCost::setBounds(int sequence, double lower, double upper)
     762{
     763  int start = start_[sequence];
     764  assert(start_[sequence+1]==start+4);
     765  lower_[start+1]=lower;
     766  lower_[start+2]=upper;
     767}
     768
  • trunk/ClpSimplex.cpp

    r164 r166  
    24022402  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
    24032403}
     2404#ifdef QUADRATIC
    24042405#include "ClpSimplexPrimalQuadratic.hpp"
    24052406/* Solves quadratic problem using SLP - may be used as crash
     
    24112412  return ((ClpSimplexPrimalQuadratic *) this)->primalSLP(numberPasses,deltaTolerance);
    24122413}
    2413 // Solves quadratic using Beale's algorithm - primal
     2414// Solves quadratic using Wolfe's algorithm - primal
    24142415int
    2415 ClpSimplex::quadraticBeale()
    2416 {
    2417   return ((ClpSimplexPrimalQuadratic *) this)->primalBeale();
    2418 }
     2416ClpSimplex::quadraticWolfe()
     2417{
     2418  return ((ClpSimplexPrimalQuadratic *) this)->primalWolfe();
     2419}
     2420#endif
    24192421/* For strong branching.  On input lower and upper are new bounds
    24202422   while on output they are objective function values (>1.0e50 infeasible).
  • trunk/ClpSimplexDual.cpp

    r148 r166  
    19291929          sequenceIn_=-1;
    19301930          // now choose largest and sum all ones which will go through
    1931 #define MINIMUMTHETA 1.0e-12
    19321931          for (i=0;i<interesting[iFlip];i++) {
    19331932            int iSequence=index[i];
     
    20632062    }
    20642063   
     2064#define MINIMUMTHETA 1.0e-12
     2065    // Movement should be minimum for anti-degeneracy - unless
     2066    // fixed variable out
     2067    double minimumTheta;
     2068    if (upperOut_>lowerOut_)
     2069      minimumTheta=MINIMUMTHETA;
     2070    else
     2071      minimumTheta=0.0;
    20652072    if (sequenceIn_>=0) {
    20662073      // at this stage sequenceIn_ is just pointer into index array
     
    20742081      oldValue = dj_[sequenceIn_];
    20752082      theta_ = oldValue/alpha;
    2076       if (theta_<MINIMUMTHETA) {
     2083      if (theta_<minimumTheta) {
    20772084        // can't pivot to zero
    2078         if (oldValue-MINIMUMTHETA*alpha>=-dualTolerance_) {
    2079           theta_=MINIMUMTHETA;
    2080         } else if (oldValue-MINIMUMTHETA*alpha>=-newTolerance) {
    2081           theta_=MINIMUMTHETA;
     2085        if (oldValue-minimumTheta*alpha>=-dualTolerance_) {
     2086          theta_=minimumTheta;
     2087        } else if (oldValue-minimumTheta*alpha>=-newTolerance) {
     2088          theta_=minimumTheta;
    20822089          thisIncrease=true;
    20832090        } else {
  • trunk/ClpSimplexPrimal.cpp

    r162 r166  
    261261  }
    262262  // clean up
     263  unflag();
    263264  finish();
    264265  restoreData(data);
     
    637638      }
    638639    } else {
     640      if(type==3&&problemStatus_!=-5)
     641        unflag(); // odd
    639642      // carry on
    640643      problemStatus_ = -1;
     
    10021005  if (pivotRow_>=0) {
    10031006   
     1007    alpha_ = work[pivotRow_];
     1008    // translate to sequence
     1009    sequenceOut_ = pivotVariable_[pivotRow_];
     1010    valueOut_ = solution(sequenceOut_);
     1011    lowerOut_=lower_[sequenceOut_];
     1012    upperOut_=upper_[sequenceOut_];
    10041013#define MINIMUMTHETA 1.0e-12
     1014    // Movement should be minimum for anti-degeneracy - unless
     1015    // fixed variable out
     1016    double minimumTheta;
     1017    if (upperOut_>lowerOut_)
     1018      minimumTheta=MINIMUMTHETA;
     1019    else
     1020      minimumTheta=0.0;
    10051021    // will we need to increase tolerance
    10061022#ifdef CLP_DEBUG
     
    10081024#endif
    10091025    double largestInfeasibility = primalTolerance_;
    1010     if (theta_<MINIMUMTHETA) {
    1011       theta_=MINIMUMTHETA;
     1026    if (theta_<minimumTheta) {
     1027      theta_=minimumTheta;
    10121028      for (iIndex=0;iIndex<numberSwapped;iIndex++) {
    10131029        int iRow = indexSwapped[iIndex];
     
    10271043      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
    10281044    }
    1029     alpha_ = work[pivotRow_];
    1030     // translate to sequence
    1031     sequenceOut_ = pivotVariable_[pivotRow_];
    1032     valueOut_ = solution(sequenceOut_);
    1033     lowerOut_=lower_[sequenceOut_];
    1034     upperOut_=upper_[sequenceOut_];
    10351045
    10361046    if (way<0.0)
  • trunk/ClpSimplexPrimalQuadratic.cpp

    r163 r166  
    88#include "CoinHelperFunctions.hpp"
    99#include "ClpSimplexPrimalQuadratic.hpp"
     10#include "ClpPrimalQuadraticDantzig.hpp"
    1011#include "ClpFactorization.hpp"
    1112#include "ClpNonLinearCost.hpp"
     
    7778
    7879  // get feasible
    79   if (numberPrimalInfeasibilities())
     80  if (!this->status()||numberPrimalInfeasibilities())
    8081    primal(1);
    8182  // still infeasible
     
    9495    trueLower[jNon]=columnLower[iColumn];
    9596    trueUpper[jNon]=columnUpper[iColumn];
     97    if (solution[iColumn]<trueLower[jNon])
     98      solution[iColumn]=trueLower[jNon];
     99    else if (solution[iColumn]>trueUpper[jNon])
     100      solution[iColumn]=trueUpper[jNon];
    96101  }
    97102  int iPass;
     
    159164      printf("obj at saved %g, obj now %g zero deriv at %g - value %g\n",
    160165             coeff0+coeff1+coeff2,coeff0,lambda,lambdaValue);
     166      if (!iPass) lambda=0.0;
    161167      if (lambda>0.0&&lambda<=1.0) {
    162168        // update solution
     
    199205        double valueJ = solution[jColumn];
    200206        double elementValue = quadraticElement[j];
    201         objValue += 0.5*valueI*valueJ*elementValue;
    202         offset += 0.5*valueI*valueJ*elementValue;
     207        elementValue *= 0.5;
     208        objValue += valueI*valueJ*elementValue;
     209        offset += valueI*valueJ*elementValue;
    203210        double gradientI = valueJ*elementValue;
    204211        double gradientJ = valueI*elementValue;
     
    258265        }
    259266      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
    260         trust[jNon] *= 0.5;
     267        trust[jNon] *= 0.2;
    261268      }
    262269      maxGap = max(maxGap,trust[jNon]);
    263270    }
    264271    std::cout<<"largest gap is "<<maxGap<<std::endl;
     272    if (iPass>10000) {
     273      for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     274        trust[jNon] *=0.0001;
     275    }
    265276    if (goodMove>0) {
    266277      double drop = lastObjective-objValue;
     
    301312      memset(r,0,numberColumns*sizeof(double));
    302313    }
     314#if 0
    303315    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    304316      iColumn=listNonLinearColumn[jNon];
     
    324336      }
    325337    }
     338#endif
    326339    if (goodMove) {
    327340      memcpy(saveSolution,solution,numberColumns*sizeof(double));
     
    350363      this->scaling(false);
    351364      this->primal(1);
     365      if (this->status()) {
     366        CoinMpsIO writer;
     367        writer.setMpsData(*matrix(), COIN_DBL_MAX,
     368                          getColLower(), getColUpper(),
     369                          getObjCoefficients(),
     370                          (const char*) 0 /*integrality*/,
     371                          getRowLower(), getRowUpper(),
     372                          NULL,NULL);
     373        writer.writeMps("xx.mps");
     374      }
     375      assert (!this->status()); // must be feasible
    352376      goodMove=1;
    353377    } else {
     
    420444  return 0;
    421445}
    422 // Beale's method
     446// Wolfe's method
    423447int
    424 ClpSimplexPrimalQuadratic::primalBeale()
     448ClpSimplexPrimalQuadratic::primalWolfe()
    425449{
    426450  // Wolfe's method looks easier - lets try that
    427   // This is as a user would see
    428451
    429452  int numberColumns = this->numberColumns();
     
    464487  if (numberPrimalInfeasibilities())
    465488    return 0;
     489  //printf("For testing - deliberate bad solution\n");
     490  //columnUpper[0]=0.0;
     491  //columnLower[0]=0.0;
     492  //quadraticSLP(50,1.0e-4);
     493  //primal(1);
     494  //columnUpper[0]=COIN_DBL_MAX;
    466495 
    467496  // Create larger problem
     
    473502    +2*quadratic->getNumElements()
    474503    + 2*numberColumns;
    475   // type array
    476   // >=0 points to dj or x (whichever this isn't)
    477   // -2 row variable
    478   // -1 artificial (see previous)
    479   int * type = new int[newNumberColumns];
    480504  double * elements2 = new double[numberElements];
    481505  int * start2 = new int[newNumberColumns+1];
     
    500524  const int * rowLengthQ = copyQ.getVectorLengths();
    501525  const double * elementByRowQ = copyQ.getElements();
     526  // Move solution across
     527  double * solution2 = new double[newNumberColumns];
     528  memset(solution2,0,newNumberColumns*sizeof(double));
    502529  newNumberColumns=0;
    503530  numberElements=0;
    504531  start2[0]=0;
    505532  // x
     533  memcpy(dj,objective,numberColumns*sizeof(double));
    506534  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    507535    // Original matrix
    508536    columnLower2[iColumn]=columnLower[iColumn];
    509537    columnUpper2[iColumn]=columnUpper[iColumn];
     538    solution2[iColumn]=solution[iColumn];
    510539    int j;
    511540    for (j=columnStart[iColumn];
     
    523552      if (iColumn!=jColumn)
    524553        value *= 0.5;
    525       dj[iColumn] += solution[jColumn]*value;
     554      dj[jColumn] += solution[iColumn]*value;
    526555      elements2[numberElements]=-value;
    527556      row2[numberElements++]=jColumn+numberRows;
     
    534563      if (iColumn!=jColumn) {
    535564        value *= 0.5;
    536         dj[iColumn] += solution[jColumn]*value;
     565        dj[jColumn] += solution[iColumn]*value;
    537566        elements2[numberElements]=-value;
    538567        row2[numberElements++]=jColumn+numberRows;
     
    553582  const double * elementByRow = copy.getElements();
    554583  for (iRow=0;iRow<numberRows;iRow++) {
     584    solution2[newNumberColumns]=pi[iRow];
     585    double value = pi[iRow];
    555586    // should look at rows to get correct bounds
    556587    if (rowLower[iRow]==rowUpper[iRow]) {
     
    572603         j<rowStart[iRow]+rowLength[iRow];
    573604         j++) {
    574       elements2[numberElements]=elementByRow[j];
    575       row2[numberElements++]=column[j]+numberRows;
    576     }
    577     type[newNumberColumns]=-2;
     605      double elementValue=elementByRow[j];
     606      int jColumn = column[j];
     607      elements2[numberElements]=elementValue;
     608      row2[numberElements++]=jColumn+numberRows;
     609      dj[jColumn]-= value*elementValue;
     610    }
    578611    newNumberColumns++;
    579612    start2[newNumberColumns]=numberElements;
    580613  }
    581   // djs and artificials
     614  // djs
    582615  double tolerance = dualTolerance();
    583616  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    584     // dj
    585617    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
    586618    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
    587     elements2[numberElements]=1.0;
    588     row2[numberElements++]=iColumn+numberRows;
    589     type[newNumberColumns]=iColumn;
    590     type[iColumn]=newNumberColumns;
    591     newNumberColumns++;
    592     start2[newNumberColumns]=numberElements;
    593     // artificial (assuming no bounds)
     619    solution2[newNumberColumns]=dj[iColumn];
     620    // If basic then fix dj
    594621    if (getStatus(iColumn)==basic||
    595622        (solution[iColumn]>columnLower[iColumn]+1.0e-7&&
    596623         solution[iColumn]<columnUpper[iColumn]-1.0e-7)) {
    597       columnUpper2[newNumberColumns-1]=0.0; // fix for now
    598       columnLower2[newNumberColumns-1]=0.0; // fix for now
     624      columnUpper2[newNumberColumns]=0.0;
     625      columnLower2[newNumberColumns]=0.0;
     626      solution2[newNumberColumns]=0.0;
     627    } else {
     628      // let's try and encourage good djs
     629      if (solution[iColumn]<columnLower[iColumn]+1.0e-7) {
     630        columnUpper2[iColumn]=columnLower[iColumn]; // fix for now
     631        columnLower2[newNumberColumns]=0.0;
     632      } else {
     633        columnLower2[iColumn]=columnUpper[iColumn]; // fix for now
     634        columnUpper2[newNumberColumns]=0.0;
     635      }
     636    }
     637    elements2[numberElements]=1.0;
     638    row2[numberElements++]=iColumn+numberRows;
     639    newNumberColumns++;
     640    start2[newNumberColumns]=numberElements;
     641  }
     642  // Artificials and fix all non-basic x
     643  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     644    if (getStatus(iColumn)==basic||
     645        (solution[iColumn]>columnLower[iColumn]+1.0e-7&&
     646         solution[iColumn]<columnUpper[iColumn]-1.0e-7)) {
    599647      if (fabs(dj[iColumn])>tolerance) {
    600648        columnLower2[newNumberColumns]=0.0;
    601649        columnUpper2[newNumberColumns]=COIN_DBL_MAX;
    602650        objective2[newNumberColumns]=1.0;
     651        solution2[newNumberColumns]=fabs(dj[iColumn]);
    603652        if (dj[iColumn]>0.0)
    604653          elements2[numberElements]=1.0;
     
    611660    } else {
    612661      if (solution[iColumn]<columnLower[iColumn]+1.0e-7) {
    613         columnUpper2[iColumn]=columnLower[iColumn]; // fix for now
    614662        if (dj[iColumn]<-tolerance) {
    615663          columnLower2[newNumberColumns]=0.0;
    616664          columnUpper2[newNumberColumns]=COIN_DBL_MAX;
    617665          objective2[newNumberColumns]=1.0;
     666          solution2[newNumberColumns]=fabs(dj[iColumn]);
    618667          elements2[numberElements]=-1.0;
    619668          row2[numberElements++]=iColumn+numberRows;
     
    622671        }
    623672      } else {
    624         columnLower2[iColumn]=columnUpper[iColumn]; // fix for now
    625673        if (dj[iColumn]>tolerance) {
    626674          columnLower2[newNumberColumns]=0.0;
    627675          columnUpper2[newNumberColumns]=COIN_DBL_MAX;
    628676          objective2[newNumberColumns]=1.0;
     677          solution2[newNumberColumns]=fabs(dj[iColumn]);
    629678          elements2[numberElements]=1.0;
    630679          row2[numberElements++]=iColumn+numberRows;
     
    634683      }
    635684    }
    636     // temp
    637     //columnLower2[iColumn]=solution[iColumn];
    638     //columnUpper2[iColumn]=solution[iColumn];
    639685  }
    640686  // Create model
    641   ClpSimplex model2;
     687  ClpSimplex model2(*this);
     688  model2.resize(0,0);
    642689  model2.loadProblem(newNumberColumns,newNumberRows,
    643690                     start2,row2, elements2,
     
    645692                     objective2,
    646693                     rowLower2,rowUpper2);
     694  int numberColumns2 = newNumberColumns;
    647695  model2.allSlackBasis();
     696  //printf("Need to move stuff across to model2\n");
     697  //model2.setLogLevel(this->logLevel());
    648698  // Move solution across
    649   double * solution2 = model2.primalColumnSolution();
    650   memcpy(solution2,solution,numberColumns*sizeof(double));
    651   memcpy(solution2+numberColumns,pi,numberRows*sizeof(double));
     699  memcpy(model2.primalColumnSolution(),solution2,newNumberColumns*sizeof(double));
     700  delete [] solution2;
     701  solution2 = model2.primalColumnSolution();
    652702  for (iRow=0;iRow<numberRows;iRow++)
    653703    model2.setRowStatus(iRow,getRowStatus(iRow));
    654   // djs and artificials
     704  for (iColumn=0;iColumn<numberColumns;iColumn++)
     705    model2.setColumnStatus(iColumn,getColumnStatus(iColumn));
     706  // djs
    655707  newNumberColumns = numberRows+numberColumns;
    656708  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    657     model2.setStatus(iColumn,getStatus(iColumn));
    658     if (getStatus(iColumn)==basic||
    659         (solution[iColumn]>columnLower[iColumn]+1.0e-7&&
    660          solution[iColumn]<columnUpper[iColumn]-1.0e-7)) {
    661       solution2[newNumberColumns++]=0.0;
    662       if (fabs(dj[iColumn])>tolerance)
    663         solution2[newNumberColumns++]=fabs(dj[iColumn]);
    664     } else if (dj[iColumn]<-tolerance) {
    665       solution2[newNumberColumns++]=0.0;
    666       solution2[newNumberColumns++]=-dj[iColumn];
    667     } else {
    668       solution2[newNumberColumns++]=dj[iColumn];
    669     }
     709    if (solution2[newNumberColumns]) {
     710      model2.setRowStatus(numberRows+iColumn,isFixed);
     711      model2.setColumnStatus(newNumberColumns,basic);
     712    }
     713    newNumberColumns++;
     714  }
     715  for (iColumn=newNumberColumns;iColumn<numberColumns2;iColumn++) {
     716    model2.setStatus(iColumn,basic);
    670717  }
    671718  memset(model2.primalRowSolution(),0,newNumberRows*sizeof(double));
     
    673720               model2.primalRowSolution());
    674721  // solve
    675 #if 0
     722  model2.scaling(false);
     723  model2.primal(1);
     724  delete [] columnLower2;
     725  delete [] columnUpper2;
     726  columnLower2 = model2.columnLower();
     727  columnUpper2 = model2.columnUpper();
     728  // free up non basic djs
     729  int offset = numberRows+numberColumns;
     730  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     731    if (model2.getStatus(iColumn)!=basic) {
     732      if (solution2[iColumn]<columnLower[iColumn]+1.0e-7) {
     733        columnUpper2[iColumn]=columnLower[iColumn]; // fix for now
     734        columnLower2[offset+iColumn]=-COIN_DBL_MAX;
     735        columnUpper2[offset+iColumn]=COIN_DBL_MAX;
     736      }
     737    }
     738  }
     739  model2.primal(1);
     740  // Free up bounds
     741  // Drop artificials - pricing should check
     742  int jColumn;
     743  // use start2 for dropping columns
     744  for (jColumn=2*numberColumns+numberRows;jColumn<numberColumns2;jColumn++) {
     745    double elementValue=elements2[start2[jColumn]];
     746    int iColumn = row2[start2[jColumn]]-numberRows;
     747    int djColumn = iColumn+newNumberRows;
     748    double value = solution2[jColumn];
     749    double djSolution = solution2[djColumn]+
     750      elementValue*value;
     751    start2[jColumn-2*numberColumns-numberRows]=jColumn;
     752    // move solution across to dj
     753    solution2[djColumn]=djSolution;
     754    columnLower2[djColumn]=-COIN_DBL_MAX;
     755    columnUpper2[djColumn]=COIN_DBL_MAX;
     756    if (model2.getColumnStatus(jColumn)==basic)
     757      model2.setColumnStatus(djColumn,basic);
     758    else if (model2.getColumnStatus(djColumn)!=basic)
     759      model2.setColumnStatus(djColumn,isFree);
     760    if (model2.getColumnStatus(iColumn)==basic) {
     761      if (fabs(djSolution)>dualTolerance_)
     762        printf("bad dj/ basic combo %d\n",iColumn);
     763    }
     764  }
     765  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     766    if (model2.getColumnStatus(iColumn)!=basic) {
     767      if (columnUpper[iColumn]>columnLower[iColumn]) {
     768        double value = solution2[iColumn];
     769        columnLower2[iColumn]=columnLower[iColumn];
     770        columnUpper2[iColumn]=columnUpper[iColumn];
     771        if (value<columnLower[iColumn]+1.0e-8)
     772          model2.setColumnStatus(iColumn,atLowerBound);
     773        else if (value>columnUpper[iColumn]-1.0e-8)
     774          model2.setColumnStatus(iColumn,atUpperBound);
     775        else
     776          abort();
     777      }
     778    }
     779  }
     780  // delete artificials
     781  model2.deleteColumns(numberColumns2-2*numberColumns-numberRows,start2);
     782  delete [] start2;
     783  delete [] row2;
     784  delete [] elements2;
     785  delete [] objective2;
     786  delete [] rowLower2;
     787  delete [] rowUpper2;
     788  solution2 = model2.primalColumnSolution();
     789  columnLower2 = model2.columnLower();
     790  columnUpper2 = model2.columnUpper();
     791 
     792  // See if any s sub j have wrong sign and/or use djs from infeasibility objective
     793  double objectiveOffset;
     794  getDblParam(ClpObjOffset,objectiveOffset);
     795  double objValue = -objectiveOffset;
     796  for (iColumn=0;iColumn<numberColumns;iColumn++)
     797    objValue += objective[iColumn]*solution2[iColumn];
     798  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     799    double valueI = solution2[iColumn];
     800    if (fabs(valueI)>1.0e-5) {
     801      int djColumn = iColumn+numberRows+numberColumns;
     802      assert(solution2[djColumn]<1.0e-7);
     803    }
     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
    676817  CoinMpsIO writer;
    677818  writer.setMpsData(*model2.matrix(), COIN_DBL_MAX,
     
    683824  writer.writeMps("xx.mps");
    684825#endif 
    685   model2.scaling(false);
    686   model2.primal(1);
    687   // If objective 0.0 then we can drop artificials
    688   // See if any s sub j have wrong sign and/or use djs from infeasibility objective
    689   double objectiveOffset;
    690   getDblParam(ClpObjOffset,objectiveOffset);
    691   double objValue = -objectiveOffset;
     826  // Now do quadratic
     827  // Do safe cast as no data
     828  ClpSimplexPrimalQuadratic * modelPtr =
     829    (ClpSimplexPrimalQuadratic *) (&model2);
     830  ClpPrimalQuadraticDantzig dantzigP(modelPtr,numberRows);;
     831  modelPtr->setPrimalColumnPivotAlgorithm(dantzigP);
     832  modelPtr->primalWolfe2(this);
     833  memcpy(dualRowSolution(),model2.primalColumnSolution()+numberColumns_,numberRows_*sizeof(double));
     834  memcpy(primalColumnSolution(),model2.primalColumnSolution(),numberColumns_*sizeof(double));
     835  memset(model2.primalRowSolution(),0,newNumberRows*sizeof(double));
     836  model2.times(1.0,model2.primalColumnSolution(),
     837               model2.primalRowSolution());
     838  memcpy(dualColumnSolution(),model2.primalColumnSolution()+numberRows_+numberColumns_,
     839         numberColumns_*sizeof(double));
     840  objValue = -objectiveOffset;
    692841  for (iColumn=0;iColumn<numberColumns;iColumn++)
    693842    objValue += objective[iColumn]*solution2[iColumn];
    694843  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    695844    double valueI = solution2[iColumn];
    696     if (fabs(valueI)>1.0e-5)
    697       assert(solution2[type[iColumn]]<1.0e-7);
     845    if (fabs(valueI)>1.0e-5) {
     846      int djColumn = iColumn+numberRows+numberColumns;
     847      assert(solution2[djColumn]<1.0e-7);
     848    }
    698849    int j;
    699850    for (j=columnQuadraticStart[iColumn];
     
    706857  }
    707858  printf("Objective value %g\n",objValue);
    708   for (iColumn=0;iColumn<newNumberColumns;iColumn++)
    709     printf("%d %g\n",iColumn,solution2[iColumn]);
    710   delete [] type;
     859  objectiveValue_ = objValue + objectiveOffset;
    711860  return 0;
    712861}
     862int ClpSimplexPrimalQuadratic::primalWolfe2 (const ClpSimplexPrimalQuadratic * originalModel)
     863{
     864
     865  algorithm_ = +2;
     866
     867  // save data
     868  ClpDataSave data = saveData();
    713869 
    714 
     870 
     871  // initialize - maybe values pass and algorithm_ is +1
     872  if (!startup(0)) {
     873   
     874    int lastCleaned=0; // last time objective or bounds cleaned up
     875    int sequenceIn=-1;
     876    int crucialSj=-1;
     877   
     878    // special nonlinear cost
     879    delete nonLinearCost_;
     880    nonLinearCost_ = new ClpNonLinearCost(this,true);
     881    // Progress indicator
     882    ClpSimplexProgress progress(this);
     883   
     884    // Say no pivot has occurred (for steepest edge and updates)
     885    pivotRow_=-2;
     886   
     887    // This says whether to restore things etc
     888    int factorType=0;
     889    /*
     890      Status of problem:
     891      0 - optimal
     892      1 - infeasible
     893      2 - unbounded
     894      -1 - iterating
     895      -2 - factorization wanted
     896      -3 - redo checking without factorization
     897      -4 - looks infeasible
     898      -5 - looks unbounded
     899    */
     900    while (problemStatus_<0) {
     901      int iRow,iColumn;
     902      // clear
     903      for (iRow=0;iRow<4;iRow++) {
     904        rowArray_[iRow]->clear();
     905      }   
     906     
     907      for (iColumn=0;iColumn<2;iColumn++) {
     908        columnArray_[iColumn]->clear();
     909      }   
     910     
     911      // give matrix (and model costs and bounds a chance to be
     912      // refreshed (normally null)
     913      matrix_->refresh(this);
     914      // If we have done no iterations - special
     915      if (lastGoodIteration_==numberIterations_)
     916        factorType=3;
     917      // may factorize, checks if problem finished
     918      statusOfProblemInPrimal(lastCleaned,factorType,progress);
     919     
     920      // Compute objective function from scratch
     921      const CoinPackedMatrix * quadratic = originalModel->quadraticObjective();
     922      const int * columnQuadratic = quadratic->getIndices();
     923      const int * columnQuadraticStart = quadratic->getVectorStarts();
     924      const int * columnQuadraticLength = quadratic->getVectorLengths();
     925      const double * quadraticElement = quadratic->getElements();
     926      const double * originalCost = originalModel->objective();
     927      objectiveValue_=0.0;
     928      int originalNumberColumns = originalModel->numberColumns();
     929      for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
     930        double valueI = solution_[iColumn];
     931        objectiveValue_ += valueI*originalCost[iColumn];
     932        int j;
     933        for (j=columnQuadraticStart[iColumn];
     934             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     935          int jColumn = columnQuadratic[j];
     936          double valueJ = solution_[jColumn];
     937          double elementValue = 0.5*quadraticElement[j];
     938          objectiveValue_ += valueI*valueJ*elementValue;
     939        }
     940      }
     941
     942      // Say good factorization
     943      factorType=1;
     944     
     945      // Say no pivot has occurred (for steepest edge and updates)
     946      pivotRow_=-2;
     947
     948      // exit if victory declared
     949      if (primalColumnPivot_->pivotColumn(rowArray_[1],
     950                                          rowArray_[2],rowArray_[3],
     951                                          columnArray_[0],
     952                                          columnArray_[1]) < 0) {
     953        problemStatus_=0;
     954        break;
     955      }
     956     
     957      // Iterate
     958      problemStatus_=-1;
     959      whileIterating(originalModel,sequenceIn,crucialSj);
     960    }
     961  }
     962  // clean up
     963  finish();
     964  restoreData(data);
     965  return problemStatus_;
     966}
     967/*
     968  Reasons to come out:
     969  -1 iterations etc
     970  -2 inaccuracy
     971  -3 slight inaccuracy (and done iterations)
     972  -4 end of values pass and done iterations
     973  +0 looks optimal (might be infeasible - but we will investigate)
     974  +2 looks unbounded
     975  +3 max iterations
     976*/
     977int
     978ClpSimplexPrimalQuadratic::whileIterating(
     979                      const ClpSimplexPrimalQuadratic * originalModel,
     980                      int & sequenceIn,int & crucialSj)
     981{
     982
     983  int returnCode=-1;
     984  double saveObjective = objectiveValue_;
     985  int originalNumberColumns = originalModel->numberColumns();
     986  // status stays at -1 while iterating, >=0 finished, -2 to invert
     987  // status -3 to go to top without an invert
     988  while (problemStatus_==-1) {
     989#ifdef CLP_DEBUG
     990    {
     991      int i;
     992      // not [1] as has information
     993      for (i=0;i<4;i++) {
     994        if (i!=1)
     995          rowArray_[i]->checkClear();
     996      }   
     997      for (i=0;i<2;i++) {
     998        columnArray_[i]->checkClear();
     999      }   
     1000    }     
     1001#endif
     1002    // choose column to come in
     1003    // can use pivotRow_ to update weights
     1004    // pass in list of cost changes so can do row updates (rowArray_[1])
     1005    // NOTE rowArray_[0] is used by computeDuals which is a
     1006    // slow way of getting duals but might be used
     1007    // Initially Dantzig and look at s variables
     1008    // Only do if one not already chosen
     1009    if (sequenceIn<0)
     1010      primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
     1011                   columnArray_[0],columnArray_[1]);
     1012    else
     1013      sequenceIn_ = sequenceIn;
     1014    pivotRow_=-1;
     1015    sequenceOut_=-1;
     1016    rowArray_[1]->clear();
     1017    if (sequenceIn_>=0) {
     1018      sequenceIn=-1;
     1019      // we found a pivot column
     1020      int chosen = sequenceIn_;
     1021      // do second half of iteration
     1022      while (chosen>=0) {
     1023        objectiveValue_=saveObjective;
     1024        returnCode=-1;
     1025        pivotRow_=-1;
     1026        sequenceOut_=-1;
     1027        rowArray_[1]->clear();
     1028        // we found a pivot column
     1029        // update the incoming column
     1030        sequenceIn_=chosen;
     1031        chosen=-1;
     1032        if (sequenceIn_>=originalNumberColumns) {
     1033          // move back to a version of primalColumn?
     1034          valueIn_=solution_[sequenceIn_];
     1035          dualIn_=dj_[sequenceIn_];
     1036          nonLinearCost_->setBounds(sequenceIn_, -0.5*COIN_DBL_MAX,
     1037                                    0.5*COIN_DBL_MAX);
     1038          lower_[sequenceIn_]=-0.5*COIN_DBL_MAX;
     1039          upper_[sequenceIn_]=0.5*COIN_DBL_MAX;
     1040          lowerIn_=lower_[sequenceIn_];
     1041          upperIn_=upper_[sequenceIn_];
     1042          if (dualIn_>0.0)
     1043            directionIn_ = -1;
     1044          else
     1045            directionIn_ = 1;
     1046        } else {
     1047          // Set dj as value of slack
     1048          dualIn_=solution_[sequenceIn_+numberRows_];
     1049          crucialSj = sequenceIn_+numberRows_; // sj which should go to 0.0
     1050        }
     1051        unpack(rowArray_[1]);
     1052        // save reduced cost
     1053        //double saveDj = dualIn_;
     1054        factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
     1055        // do ratio test and re-compute dj
     1056        // Note second parameter long enough for columns
     1057        int result=primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
     1058                             originalModel,crucialSj);
     1059        saveObjective = objectiveValue_;
     1060        if (pivotRow_>=0) {
     1061          // if stable replace in basis
     1062          int updateStatus = factorization_->replaceColumn(rowArray_[2],
     1063                                                           pivotRow_,
     1064                                                           alpha_);
     1065          // if no pivots, bad update but reasonable alpha - take and invert
     1066          if (updateStatus==2&&
     1067              lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
     1068            updateStatus=4;
     1069          if (updateStatus==1||updateStatus==4) {
     1070            // slight error
     1071            if (factorization_->pivots()>5||updateStatus==4) {
     1072              returnCode=-3;
     1073            }
     1074          } else if (updateStatus==2) {
     1075            // major error
     1076            // better to have small tolerance even if slower
     1077            factorization_->zeroTolerance(1.0e-15);
     1078            int maxFactor = factorization_->maximumPivots();
     1079            if (maxFactor>10) {
     1080              if (forceFactorization_<0)
     1081                forceFactorization_= maxFactor;
     1082              forceFactorization_ = max (1,(forceFactorization_>>1));
     1083            }
     1084            // later we may need to unwind more e.g. fake bounds
     1085            if(lastGoodIteration_ != numberIterations_) {
     1086              rowArray_[1]->clear();
     1087              pivotRow_=-1;
     1088              returnCode=-4;
     1089              // retry on return
     1090              sequenceIn = sequenceIn_;
     1091              break;
     1092            } else {
     1093              // need to reject something
     1094              char x = isColumn(sequenceIn_) ? 'C' :'R';
     1095              handler_->message(CLP_SIMPLEX_FLAG,messages_)
     1096                <<x<<sequenceWithin(sequenceIn_)
     1097                <<CoinMessageEol;
     1098              setFlagged(sequenceIn_);
     1099              lastBadIteration_ = numberIterations_; // say be more cautious
     1100              rowArray_[1]->clear();
     1101              pivotRow_=-1;
     1102              returnCode = -5;
     1103              break;
     1104             
     1105            }
     1106          } else if (updateStatus==3) {
     1107            // out of memory
     1108            // increase space if not many iterations
     1109            if (factorization_->pivots()<
     1110                0.5*factorization_->maximumPivots()&&
     1111                factorization_->pivots()<200)
     1112              factorization_->areaFactor(
     1113                                         factorization_->areaFactor() * 1.1);
     1114            returnCode =-2; // factorize now
     1115          }
     1116          // here do part of steepest - ready for next iteration
     1117          primalColumnPivot_->updateWeights(rowArray_[1]);
     1118        } else {
     1119          if (pivotRow_==-1) {
     1120            // no outgoing row is valid
     1121            rowArray_[0]->clear();
     1122            if (!factorization_->pivots()) {
     1123              returnCode = 2; //say looks unbounded
     1124              // do ray
     1125              primalRay(rowArray_[1]);
     1126            } else {
     1127              returnCode = 4; //say looks unbounded but has iterated
     1128            }
     1129            break;
     1130          } else {
     1131            // flipping from bound to bound
     1132          }
     1133        }
     1134       
     1135       
     1136        // update primal solution
     1137       
     1138        double objectiveChange=0.0;
     1139        // Cost on pivot row may change - may need to change dualIn
     1140        double oldCost=0.0;
     1141        if (pivotRow_>=0)
     1142          oldCost = cost(pivotVariable_[pivotRow_]);
     1143        // rowArray_[1] is not empty - used to update djs
     1144        updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
     1145        if (pivotRow_>=0)
     1146          dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
     1147        double oldValue = valueIn_;
     1148        if (directionIn_==-1) {
     1149          // as if from upper bound
     1150          if (sequenceIn_!=sequenceOut_) {
     1151            // variable becoming basic
     1152            valueIn_ -= fabs(theta_);
     1153          } else {
     1154            valueIn_=lowerIn_;
     1155          }
     1156        } else {
     1157          // as if from lower bound
     1158          if (sequenceIn_!=sequenceOut_) {
     1159            // variable becoming basic
     1160            valueIn_ += fabs(theta_);
     1161          } else {
     1162            valueIn_=upperIn_;
     1163          }
     1164        }
     1165        objectiveChange += dualIn_*(valueIn_-oldValue);
     1166        // outgoing
     1167        if (sequenceIn_!=sequenceOut_) {
     1168          if (directionOut_>0) {
     1169            valueOut_ = lowerOut_;
     1170          } else {
     1171            valueOut_ = upperOut_;
     1172          }
     1173          double lowerValue = lower_[sequenceOut_];
     1174          double upperValue = upper_[sequenceOut_];
     1175          assert(valueOut_>=lowerValue-primalTolerance_&&
     1176                 valueOut_<=upperValue+primalTolerance_);
     1177          // may not be exactly at bound and bounds may have changed
     1178          if (valueOut_<=lowerValue+primalTolerance_) {
     1179            directionOut_=1;
     1180          } else if (valueOut_>=upperValue-primalTolerance_) {
     1181            directionOut_=-1;
     1182          } else {
     1183            printf("*** variable wandered off bound %g %g %g!\n",
     1184                   lowerValue,valueOut_,upperValue);
     1185            if (valueOut_-lowerValue<=upperValue-valueOut_) {
     1186              valueOut_=lowerValue;
     1187              directionOut_=1;
     1188            } else {
     1189              valueOut_=upperValue;
     1190              directionOut_=-1;
     1191            }
     1192          }
     1193          solution_[sequenceOut_]=valueOut_;
     1194          nonLinearCost_->setOne(sequenceOut_,valueOut_);
     1195        }
     1196        // change cost and bounds on incoming if primal
     1197        nonLinearCost_->setOne(sequenceIn_,valueIn_);
     1198        int whatNext=housekeeping(objectiveChange);
     1199       
     1200        if (whatNext==1) {
     1201          returnCode =-2; // refactorize
     1202        } else if (whatNext==2) {
     1203          // maximum iterations or equivalent
     1204          returnCode=3;
     1205        } else if(numberIterations_ == lastGoodIteration_
     1206                  + 2 * factorization_->maximumPivots()) {
     1207          // done a lot of flips - be safe
     1208          returnCode =-2; // refactorize
     1209        }
     1210        // may need to go round again
     1211        if (result) {
     1212          assert(sequenceOut_<originalNumberColumns||
     1213                 sequenceOut_>=numberColumns_);
     1214          if (sequenceOut_<originalNumberColumns) {
     1215            chosen =sequenceOut_ + numberRows_; // sj variable
     1216          } else {
     1217            // Does this mean we can change pi
     1218            int iRow = sequenceOut_-numberColumns_;
     1219            if (iRow<numberRows_-originalNumberColumns) {
     1220              int iPi = iRow+originalNumberColumns;
     1221              printf("pi for row %d is %g\n",
     1222                     iRow,solution_[iPi]);
     1223            } else {
     1224              printf("Row %d is in column part\n",iRow);
     1225            }
     1226            abort();
     1227          }
     1228        } else {
     1229          break;
     1230        }
     1231      }
     1232      if (returnCode<-1&&returnCode>-5) {
     1233        problemStatus_=-2; //
     1234      } else if (returnCode==-5) {
     1235        // something flagged - continue;
     1236      } else if (returnCode==2) {
     1237        problemStatus_=-5; // looks unbounded
     1238      } else if (returnCode==4) {
     1239        problemStatus_=-2; // looks unbounded but has iterated
     1240      } else if (returnCode!=-1) {
     1241        assert(returnCode==3);
     1242        problemStatus_=3;
     1243      }
     1244    } else {
     1245      // no pivot column
     1246#ifdef CLP_DEBUG
     1247      if (handler_->logLevel()&32)
     1248        printf("** no column pivot\n");
     1249#endif
     1250      if (nonLinearCost_->numberInfeasibilities())
     1251        problemStatus_=-4; // might be infeasible
     1252      returnCode=0;
     1253      break;
     1254    }
     1255  }
     1256  return returnCode;
     1257}
     1258/*
     1259   Row array has pivot column
     1260   This chooses pivot row.
     1261   For speed, we may need to go to a bucket approach when many
     1262   variables go through bounds
     1263   On exit rhsArray will have changes in costs of basic variables
     1264*/
     1265int
     1266ClpSimplexPrimalQuadratic::primalRow(CoinIndexedVector * rowArray,
     1267                                     CoinIndexedVector * rhsArray,
     1268                                     CoinIndexedVector * spareArray,
     1269                                     CoinIndexedVector * spareArray2,
     1270                                     const ClpSimplexPrimalQuadratic * originalModel,
     1271                                     int crucialSj)
     1272{
     1273  int result=1;
     1274  // sequence stays as row number until end
     1275  pivotRow_=-1;
     1276  int numberSwapped=0;
     1277  int numberRemaining=0;
     1278
     1279  int numberThru =0; // number gone thru a barrier
     1280  int lastThru =0; // number gone thru a barrier on last time
     1281 
     1282  double totalThru=0.0; // for when variables flip
     1283  double acceptablePivot=1.0e-7;
     1284  if (factorization_->pivots())
     1285    acceptablePivot=1.0e-5; // if we have iterated be more strict
     1286  double bestEverPivot=acceptablePivot;
     1287  int lastPivotRow = -1;
     1288  double lastPivot=0.0;
     1289  double lastTheta=1.0e50;
     1290  int lastNumberSwapped=0;
     1291
     1292  // use spareArrays to put ones looked at in
     1293  // First one is list of candidates
     1294  // We could compress if we really know we won't need any more
     1295  // Second array has current set of pivot candidates
     1296  // with a backup list saved in double * part of indexed vector
     1297
     1298  // for zeroing out arrays after
     1299  int maximumSwapped=0;
     1300  // pivot elements
     1301  double * spare;
     1302  // indices
     1303  int * index, * indexSwapped;
     1304  int * saveSwapped;
     1305  spareArray->clear();
     1306  spareArray2->clear();
     1307  spare = spareArray->denseVector();
     1308  index = spareArray->getIndices();
     1309  saveSwapped = (int *) spareArray2->denseVector();
     1310  indexSwapped = spareArray2->getIndices();
     1311
     1312  // we also need somewhere for effective rhs
     1313  double * rhs=rhsArray->denseVector();
     1314
     1315  /*
     1316    First we get a list of possible pivots.  We can also see if the
     1317    problem looks unbounded.
     1318
     1319    At first we increase theta and see what happens.  We start
     1320    theta at a reasonable guess.  If in right area then we do bit by bit.
     1321    We save possible pivot candidates
     1322
     1323   */
     1324
     1325  // do first pass to get possibles
     1326  // We can also see if unbounded
     1327
     1328  double * work=rowArray->denseVector();
     1329  int number=rowArray->getNumElements();
     1330  int * which=rowArray->getIndices();
     1331
     1332  // we need to swap sign if coming in from ub
     1333  double way = directionIn_;
     1334  double maximumMovement;
     1335  if (way>0.0)
     1336    maximumMovement = min(1.0e30,upperIn_-valueIn_);
     1337  else
     1338    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
     1339
     1340  int iIndex;
     1341
     1342  // Work out coefficients for quadratic term
     1343  const CoinPackedMatrix * quadratic = originalModel->quadraticObjective();
     1344  const int * columnQuadratic = quadratic->getIndices();
     1345  const int * columnQuadraticStart = quadratic->getVectorStarts();
     1346  const int * columnQuadraticLength = quadratic->getVectorLengths();
     1347  const double * quadraticElement = quadratic->getElements();
     1348  const double * originalCost = originalModel->objective();
     1349  // Use rhsArray
     1350  rhsArray->clear();
     1351  int * index2 = rhsArray->getIndices();
     1352  int originalNumberColumns = originalModel->numberColumns();
     1353  int number2=0;
     1354  //int numberOriginalRows = originalModel->numberRows();
     1355  // sj
     1356  int iSjRow=-1;
     1357  double tj = -1.0;
     1358  for (iIndex=0;iIndex<number;iIndex++) {
     1359    int iRow = which[iIndex];
     1360    double alpha = -work[iRow]*way;
     1361    int iPivot=pivotVariable_[iRow];
     1362    if (iPivot<originalNumberColumns) {
     1363      index2[number2++]=iPivot;
     1364      rhs[iPivot]=alpha;
     1365      //printf("col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
     1366    } else {
     1367      //printf("new col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
     1368      if (iPivot==crucialSj) {
     1369        tj = alpha;
     1370        iSjRow = iRow;
     1371      }
     1372    }
     1373  }
     1374  // Incoming
     1375  if (sequenceIn_<originalNumberColumns) {
     1376    index2[number2++]=sequenceIn_;
     1377    rhs[sequenceIn_]=way;
     1378    printf("incoming col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
     1379  } else {
     1380    printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
     1381  }
     1382  rhsArray->setNumElements(number2);
     1383  // Change in objective will be theta*coeff1 + theta*theta*coeff2
     1384  double coeff1 = 0.0;
     1385  double coeff2 = 0.0;
     1386 
     1387  for (iIndex=0;iIndex<number2;iIndex++) {
     1388    int iColumn=index2[iIndex];
     1389    double valueI = solution_[iColumn];
     1390    double alphaI = rhs[iColumn];
     1391    coeff1 += alphaI*originalCost[iColumn];
     1392    int j;
     1393    for (j=columnQuadraticStart[iColumn];
     1394         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1395      int jColumn = columnQuadratic[j];
     1396      double valueJ = solution_[jColumn];
     1397      double alphaJ = rhs[jColumn];
     1398      double elementValue = 0.5*quadraticElement[j];
     1399      coeff1 += (valueI*alphaJ+alphaI*valueJ)*elementValue;
     1400      coeff2 += (alphaI*alphaJ)*elementValue;
     1401    }
     1402  }
     1403  printf("coefficients %g %g\n",coeff1,coeff2);
     1404  // interesting places are where derivative zero or sj goes to zero
     1405  double d1,d2=1.0e50;
     1406  if (fabs(coeff2)>1.0e-9)
     1407    d1 = - 0.5*coeff1/coeff2;
     1408  else if (coeff1<=1.0e-9)
     1409    d1 = maximumMovement;
     1410  else
     1411    d1 = 0.0;
     1412  if (tj<=0.0) {
     1413    if (sequenceIn_<originalNumberColumns)
     1414      printf("column %d is basically linear\n",sequenceIn_);
     1415    //assert(!columnQuadraticLength[sequenceIn_]);
     1416  } else {
     1417    d2 = fabs(solution_[crucialSj]/tj);
     1418  }
     1419  printf("derivative zero at %g, sj zero at %g\n",d1,d2);
     1420  if (d1>1.0e10&&d2>1.0e10) {
     1421    // linear variable entering
     1422    // maybe we should have done dual iteration to force sj to 0.0
     1423  }
     1424  maximumMovement = min(maximumMovement,d1);
     1425  maximumMovement = min(maximumMovement,d2);
     1426   
     1427  rhsArray->clear();
     1428  double tentativeTheta = maximumMovement;
     1429  double upperTheta = maximumMovement;
     1430
     1431
     1432  for (iIndex=0;iIndex<number;iIndex++) {
     1433
     1434    int iRow = which[iIndex];
     1435    double alpha = work[iRow];
     1436    int iPivot=pivotVariable_[iRow];
     1437    alpha *= way;
     1438    double oldValue = solution(iPivot);
     1439    // get where in bound sequence
     1440    if (alpha>0.0) {
     1441      // basic variable going towards lower bound
     1442      double bound = lower(iPivot);
     1443      oldValue -= bound;
     1444    } else if (alpha<0.0) {
     1445      // basic variable going towards upper bound
     1446      double bound = upper(iPivot);
     1447      oldValue = bound-oldValue;
     1448    }
     1449    double value = oldValue-tentativeTheta*fabs(alpha);
     1450    assert (oldValue>=-primalTolerance_*1.0001);
     1451    if (value<-primalTolerance_) {
     1452      // add to list
     1453      spare[numberRemaining]=alpha;
     1454      rhs[iRow]=oldValue;
     1455      index[numberRemaining++]=iRow;
     1456      double value=oldValue-upperTheta*fabs(alpha);
     1457      if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot)
     1458        upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
     1459    }
     1460  }
     1461
     1462  // we need to keep where rhs non-zeros are
     1463  int numberInRhs=numberRemaining;
     1464  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
     1465  rhsArray->setNumElements(numberInRhs);
     1466
     1467  theta_=maximumMovement;
     1468
     1469  bool goBackOne = false;
     1470
     1471  if (numberRemaining) {
     1472
     1473    // looks like pivoting
     1474    // now try until reasonable theta
     1475    tentativeTheta = max(10.0*upperTheta,1.0e-7);
     1476    tentativeTheta = min(tentativeTheta,maximumMovement);
     1477   
     1478    // loops increasing tentative theta until can't go through
     1479   
     1480    while (tentativeTheta <= maximumMovement) {
     1481      double thruThis = 0.0;
     1482     
     1483      double bestPivot=acceptablePivot;
     1484      pivotRow_ = -1;
     1485     
     1486      numberSwapped = 0;
     1487     
     1488      upperTheta = maximumMovement;
     1489     
     1490      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
     1491
     1492        int iRow = index[iIndex];
     1493        double alpha = spare[iIndex];
     1494        double oldValue = rhs[iRow];
     1495        double value = oldValue-tentativeTheta*fabs(alpha);
     1496
     1497        if (value<-primalTolerance_) {
     1498          // how much would it cost to go thru
     1499          thruThis += alpha*
     1500            nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
     1501          // goes on swapped list (also means candidates if too many)
     1502          indexSwapped[numberSwapped++]=iRow;
     1503          if (fabs(alpha)>bestPivot) {
     1504            bestPivot=fabs(alpha);
     1505            pivotRow_ = iRow;
     1506            theta_ = oldValue/bestPivot;
     1507          }
     1508        } else {
     1509          value = oldValue-upperTheta*fabs(alpha);
     1510          if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
     1511            upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
     1512        }
     1513      }
     1514     
     1515      maximumSwapped = max(maximumSwapped,numberSwapped);
     1516
     1517      double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 99999999;
     1518      // but make a bit more pessimistic
     1519      dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
     1520      if (totalThru+thruThis>=dualCheck) {
     1521        // We should be pivoting in this batch
     1522        // so compress down to this lot
     1523
     1524        int saveNumber = numberRemaining;
     1525        numberRemaining=0;
     1526        for (iIndex=0;iIndex<numberSwapped;iIndex++) {
     1527          int iRow = indexSwapped[iIndex];
     1528          spare[numberRemaining]=way*work[iRow];
     1529          index[numberRemaining++]=iRow;
     1530        }
     1531        memset(spare+numberRemaining,0,
     1532               (saveNumber-numberRemaining)*sizeof(double));
     1533        int iTry;
     1534#define MAXTRY 100
     1535        // first get ratio with tolerance
     1536        for (iTry=0;iTry<MAXTRY;iTry++) {
     1537         
     1538          upperTheta=maximumMovement;
     1539          numberSwapped = 0;
     1540         
     1541          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
     1542           
     1543            int iRow = index[iIndex];
     1544            double alpha = fabs(spare[iIndex]);
     1545            double oldValue = rhs[iRow];
     1546            double value = oldValue-upperTheta*alpha;
     1547           
     1548            if (value<-primalTolerance_ && alpha>=acceptablePivot)
     1549              upperTheta = (oldValue+primalTolerance_)/alpha;
     1550           
     1551          }
     1552         
     1553          // now look at best in this lot
     1554          bestPivot=acceptablePivot;
     1555          pivotRow_=-1;
     1556          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
     1557           
     1558            int iRow = index[iIndex];
     1559            double alpha = spare[iIndex];
     1560            double oldValue = rhs[iRow];
     1561            double value = oldValue-upperTheta*fabs(alpha);
     1562           
     1563            if (value<=0) {
     1564              // how much would it cost to go thru
     1565              totalThru += alpha*
     1566                nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
     1567              // goes on swapped list (also means candidates if too many)
     1568              indexSwapped[numberSwapped++]=iRow;
     1569              if (fabs(alpha)>bestPivot) {
     1570                bestPivot=fabs(alpha);
     1571                theta_ = oldValue/bestPivot;
     1572                pivotRow_=iRow;
     1573              }
     1574            } else {
     1575              value = oldValue-upperTheta*fabs(alpha);
     1576              if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
     1577                upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
     1578            }
     1579          }
     1580         
     1581          maximumSwapped = max(maximumSwapped,numberSwapped);
     1582          if (bestPivot<0.1*bestEverPivot&&
     1583              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
     1584            // back to previous one
     1585            goBackOne = true;
     1586            break;
     1587          } else if (pivotRow_==-1&&upperTheta>largeValue_) {
     1588            if (lastPivot>acceptablePivot) {
     1589              // back to previous one
     1590              goBackOne = true;
     1591              //break;
     1592            } else {
     1593              // can only get here if all pivots so far too small
     1594            }
     1595            break;
     1596          } else {
     1597            dualCheck = - 2.0*coeff2*theta_ - coeff1-9999999;
     1598            if (totalThru>=dualCheck) {
     1599              break; // no point trying another loop
     1600            } else {
     1601              // skip this lot
     1602              nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
     1603              lastPivotRow=pivotRow_;
     1604              lastTheta = theta_;
     1605              lastThru = numberThru;
     1606              numberThru += numberSwapped;
     1607              lastNumberSwapped = numberSwapped;
     1608              memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
     1609              if (bestPivot>bestEverPivot)
     1610                bestEverPivot=bestPivot;
     1611            }
     1612          }
     1613        }
     1614        break;
     1615      } else {
     1616        // skip this lot
     1617        nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
     1618        lastPivotRow=pivotRow_;
     1619        lastTheta = theta_;
     1620        lastThru = numberThru;
     1621        numberThru += numberSwapped;
     1622        lastNumberSwapped = numberSwapped;
     1623        memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
     1624        if (bestPivot>bestEverPivot)
     1625          bestEverPivot=bestPivot;
     1626        totalThru += thruThis;
     1627        tentativeTheta = 2.0*upperTheta;
     1628      }
     1629    }
     1630    // can get here without pivotRow_ set but with lastPivotRow
     1631    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
     1632      // back to previous one
     1633      pivotRow_=lastPivotRow;
     1634      theta_ = lastTheta;
     1635            // undo this lot
     1636      nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
     1637      memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
     1638      numberSwapped = lastNumberSwapped;
     1639    }
     1640  }
     1641
     1642  if (pivotRow_>=0) {
     1643   
     1644#define MINIMUMTHETA 1.0e-12
     1645    // will we need to increase tolerance
     1646#ifdef CLP_DEBUG
     1647    bool found=false;
     1648#endif
     1649    double largestInfeasibility = primalTolerance_;
     1650    if (theta_<MINIMUMTHETA) {
     1651      theta_=MINIMUMTHETA;
     1652      for (iIndex=0;iIndex<numberSwapped;iIndex++) {
     1653        int iRow = indexSwapped[iIndex];
     1654#ifdef CLP_DEBUG
     1655        if (iRow==pivotRow_)
     1656          found=true;
     1657#endif
     1658        largestInfeasibility = max (largestInfeasibility,
     1659                                    -(rhs[iRow]-fabs(work[iRow])*theta_));
     1660      }
     1661#ifdef CLP_DEBUG
     1662      assert(found);
     1663      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32))
     1664        printf("Primal tolerance increased from %g to %g\n",
     1665               primalTolerance_,largestInfeasibility);
     1666#endif
     1667      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
     1668    }
     1669    alpha_ = work[pivotRow_];
     1670    // translate to sequence
     1671    sequenceOut_ = pivotVariable_[pivotRow_];
     1672    valueOut_ = solution(sequenceOut_);
     1673    lowerOut_=lower_[sequenceOut_];
     1674    upperOut_=upper_[sequenceOut_];
     1675
     1676    if (way<0.0)
     1677      theta_ = - theta_;
     1678    double newValue = valueOut_ - theta_*alpha_;
     1679    if (alpha_*way<0.0) {
     1680      directionOut_=-1;      // to upper bound
     1681      if (fabs(theta_)>1.0e-6)
     1682        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
     1683      else
     1684        upperOut_ = newValue;
     1685    } else {
     1686      directionOut_=1;      // to lower bound
     1687      if (fabs(theta_)>1.0e-6)
     1688        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
     1689      else
     1690        lowerOut_ = newValue;
     1691    }
     1692    dualOut_ = reducedCost(sequenceOut_);
     1693  } else {
     1694    double trueMaximumMovement;
     1695    if (way>0.0)
     1696      trueMaximumMovement = min(1.0e30,upperIn_-valueIn_);
     1697    else
     1698      trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_);
     1699    if (maximumMovement<1.0e20&&maximumMovement==trueMaximumMovement) {
     1700      // flip
     1701      theta_ = maximumMovement;
     1702      pivotRow_ = -2; // so we can tell its a flip
     1703      result=0;
     1704      sequenceOut_ = sequenceIn_;
     1705      valueOut_ = valueIn_;
     1706      dualOut_ = dualIn_;
     1707      lowerOut_ = lowerIn_;
     1708      upperOut_ = upperIn_;
     1709      alpha_ = 0.0;
     1710      if (way<0.0) {
     1711        directionOut_=1;      // to lower bound
     1712        theta_ = lowerOut_ - valueOut_;
     1713      } else {
     1714        directionOut_=-1;      // to upper bound
     1715        theta_ = upperOut_ - valueOut_;
     1716      }
     1717    } else if (fabs(maximumMovement-d2)<dualTolerance_) {
     1718      // sj going to zero
     1719      result=0;
     1720      assert (pivotRow_<0);
     1721      nonLinearCost_->setBounds(crucialSj, 0.0,0.0);
     1722      lower_[crucialSj]=0.0;
     1723      upper_[crucialSj]=0.0;
     1724      setColumnStatus(crucialSj,isFixed);
     1725      pivotRow_ = iSjRow;
     1726      alpha_ = work[pivotRow_];
     1727      // translate to sequence
     1728      sequenceOut_ = pivotVariable_[pivotRow_];
     1729      valueOut_ = solution(sequenceOut_);
     1730      lowerOut_=lower_[sequenceOut_];
     1731      upperOut_=upper_[sequenceOut_];
     1732      theta_ = d2;
     1733      if (way<0.0)
     1734        theta_ = - theta_;
     1735      double newValue = valueOut_ - theta_*alpha_;
     1736      if (alpha_*way<0.0) {
     1737        directionOut_=-1;      // to upper bound
     1738        if (fabs(theta_)>1.0e-6)
     1739          upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
     1740        else
     1741          upperOut_ = newValue;
     1742      } else {
     1743        directionOut_=1;      // to lower bound
     1744        if (fabs(theta_)>1.0e-6)
     1745          lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
     1746        else
     1747          lowerOut_ = newValue;
     1748      }
     1749      //????
     1750      dualOut_ = reducedCost(sequenceOut_);
     1751    } else {
     1752      // need to do something
     1753      abort();
     1754    }
     1755  }
     1756
     1757  // clear arrays
     1758
     1759  memset(spare,0,numberRemaining*sizeof(double));
     1760  memset(saveSwapped,0,maximumSwapped*sizeof(int));
     1761
     1762  // put back original bounds etc
     1763  nonLinearCost_->goBackAll(rhsArray);
     1764
     1765  rhsArray->clear();
     1766  // Change in objective will be theta*coeff1 + theta*theta*coeff2
     1767  objectiveValue_ += theta_*coeff1+theta_*theta_*coeff2;
     1768  printf("New objective value %g\n",objectiveValue_);
     1769  return result;
     1770}
     1771 
     1772
  • trunk/Idiot.cpp

    r76 r166  
    153153  printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
    154154  solve2();
     155#ifdef CLP_IDIOT
    155156  crossOver(3); // make basis ?
     157#endif
    156158}
    157159void
  • trunk/Test/unitTest.cpp

    r152 r166  
    793793    }
    794794  }
    795 #if 0 
    796795  // test network
    797796  {   
     
    804803    if (!fp) {
    805804      fprintf(stderr,"Unable to open file input.130 in mpsDir directory\n");
    806       exit(1);
    807     }
    808     int problem;
    809     char temp[100];
    810     // read and skip
    811     fscanf(fp,"%s",temp);
    812     assert (!strcmp(temp,"BEGIN"));
    813     fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows,
    814            &numberColumns);
    815     // scan down to SUPPLY
    816     while (fgets(temp,100,fp)) {
    817       if (!strncmp(temp,"SUPPLY",6))
    818         break;
    819     }
    820     if (strncmp(temp,"SUPPLY",6)) {
    821       fprintf(stderr,"Unable to find SUPPLY\n");
    822       exit(2);
    823     }
    824     // get space for rhs
    825     double * lower = new double[numberRows];
    826     double * upper = new double[numberRows];
     805    } else {
     806      int problem;
     807      char temp[100];
     808      // read and skip
     809      fscanf(fp,"%s",temp);
     810      assert (!strcmp(temp,"BEGIN"));
     811      fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows,
     812             &numberColumns);
     813      // scan down to SUPPLY
     814      while (fgets(temp,100,fp)) {
     815        if (!strncmp(temp,"SUPPLY",6))
     816          break;
     817      }
     818      if (strncmp(temp,"SUPPLY",6)) {
     819        fprintf(stderr,"Unable to find SUPPLY\n");
     820        exit(2);
     821      }
     822      // get space for rhs
     823      double * lower = new double[numberRows];
     824      double * upper = new double[numberRows];
     825      int i;
     826      for (i=0;i<numberRows;i++) {
     827        lower[i]=0.0;
     828        upper[i]=0.0;
     829      }
     830      // ***** Remember to convert to C notation
     831      while (fgets(temp,100,fp)) {
     832        int row;
     833        int value;
     834        if (!strncmp(temp,"ARCS",4))
     835          break;
     836        sscanf(temp,"%d %d",&row,&value);
     837        upper[row-1]=-value;
     838        lower[row-1]=-value;
     839      }
     840      if (strncmp(temp,"ARCS",4)) {
     841        fprintf(stderr,"Unable to find ARCS\n");
     842        exit(2);
     843      }
     844      // number of columns may be underestimate
     845      int * head = new int[2*numberColumns];
     846      int * tail = new int[2*numberColumns];
     847      int * ub = new int[2*numberColumns];
     848      int * cost = new int[2*numberColumns];
     849      // ***** Remember to convert to C notation
     850      numberColumns=0;
     851      while (fgets(temp,100,fp)) {
     852        int iHead;
     853        int iTail;
     854        int iUb;
     855        int iCost;
     856        if (!strncmp(temp,"DEMAND",6))
     857          break;
     858        sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
     859        iHead--;
     860        iTail--;
     861        head[numberColumns]=iHead;
     862        tail[numberColumns]=iTail;
     863        ub[numberColumns]=iUb;
     864        cost[numberColumns]=iCost;
     865        numberColumns++;
     866      }
     867      if (strncmp(temp,"DEMAND",6)) {
     868        fprintf(stderr,"Unable to find DEMAND\n");
     869        exit(2);
     870      }
     871      // ***** Remember to convert to C notation
     872      while (fgets(temp,100,fp)) {
     873        int row;
     874        int value;
     875        if (!strncmp(temp,"END",3))
     876          break;
     877        sscanf(temp,"%d %d",&row,&value);
     878        upper[row-1]=value;
     879        lower[row-1]=value;
     880      }
     881      if (strncmp(temp,"END",3)) {
     882        fprintf(stderr,"Unable to find END\n");
     883        exit(2);
     884      }
     885      printf("Problem %d has %d rows and %d columns\n",problem,
     886             numberRows,numberColumns);
     887      fclose(fp);
     888      ClpSimplex  model;
     889      // now build model
     890     
     891      double * objective =new double[numberColumns];
     892      double * lowerColumn = new double[numberColumns];
     893      double * upperColumn = new double[numberColumns];
     894     
     895      double * element = new double [2*numberColumns];
     896      int * start = new int[numberColumns+1];
     897      int * row = new int[2*numberColumns];
     898      start[numberColumns]=2*numberColumns;
     899      for (i=0;i<numberColumns;i++) {
     900        start[i]=2*i;
     901        element[2*i]=-1.0;
     902        element[2*i+1]=1.0;
     903        row[2*i]=head[i];
     904        row[2*i+1]=tail[i];
     905        lowerColumn[i]=0.0;
     906        upperColumn[i]=ub[i];
     907        objective[i]=cost[i];
     908      }
     909      // Create Packed Matrix
     910      CoinPackedMatrix matrix;
     911      int * lengths = NULL;
     912      matrix.assignMatrix(true,numberRows,numberColumns,
     913                          2*numberColumns,element,row,start,lengths);
     914      // load model
     915     
     916      model.loadProblem(matrix,
     917                        lowerColumn,upperColumn,objective,
     918                        lower,upper);
     919     
     920      model.factorization()->maximumPivots(200+model.numberRows()/100);
     921      model.createStatus();
     922      double time1 = cpuTime();
     923      model.dual();
     924      std::cout<<"Network problem, ClpPackedMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
     925      ClpPlusMinusOneMatrix plusMinus(matrix);
     926      assert (plusMinus.getIndices()); // would be zero if not +- one
     927      model.loadProblem(plusMinus,
     928                        lowerColumn,upperColumn,objective,
     929                        lower,upper);
     930     
     931      model.factorization()->maximumPivots(200+model.numberRows()/100);
     932      model.createStatus();
     933      time1 = cpuTime();
     934      model.dual();
     935      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
     936      ClpNetworkMatrix network(numberColumns,head,tail);
     937      model.loadProblem(network,
     938                        lowerColumn,upperColumn,objective,
     939                        lower,upper);
     940     
     941      model.factorization()->maximumPivots(200+model.numberRows()/100);
     942      model.createStatus();
     943      time1 = cpuTime();
     944      model.dual();
     945      std::cout<<"Network problem, ClpNetworkMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
     946      delete [] lower;
     947      delete [] upper;
     948      delete [] head;
     949      delete [] tail;
     950      delete [] ub;
     951      delete [] cost;
     952      delete [] objective;
     953      delete [] lowerColumn;
     954      delete [] upperColumn;
     955    }
     956  }
     957#ifdef QUADRATIC
     958  // Test quadratic to solve linear
     959  if (1) {   
     960    CoinMpsIO m;
     961    std::string fn = "tiny";
     962    m.readMps(fn.c_str(),"mps");
     963    ClpSimplex solution;
     964    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
     965                         m.getObjCoefficients(),
     966                         m.getRowLower(),m.getRowUpper());
     967    solution.dual();
     968    // get quadratic part
     969    int numberColumns=solution.numberColumns();
     970    int * start=new int [numberColumns+1];
     971    int * column = new int[numberColumns];
     972    double * element = new double[numberColumns];
    827973    int i;
    828     for (i=0;i<numberRows;i++) {
    829       lower[i]=0.0;
    830       upper[i]=0.0;
    831     }
    832     // ***** Remember to convert to C notation
    833     while (fgets(temp,100,fp)) {
    834       int row;
    835       int value;
    836       if (!strncmp(temp,"ARCS",4))
    837         break;
    838       sscanf(temp,"%d %d",&row,&value);
    839       upper[row-1]=-value;
    840       lower[row-1]=-value;
    841     }
    842     if (strncmp(temp,"ARCS",4)) {
    843       fprintf(stderr,"Unable to find ARCS\n");
    844       exit(2);
    845     }
    846     // number of columns may be underestimate
    847     int * head = new int[2*numberColumns];
    848     int * tail = new int[2*numberColumns];
    849     int * ub = new int[2*numberColumns];
    850     int * cost = new int[2*numberColumns];
    851     // ***** Remember to convert to C notation
    852     numberColumns=0;
    853     while (fgets(temp,100,fp)) {
    854       int iHead;
    855       int iTail;
    856       int iUb;
    857       int iCost;
    858       if (!strncmp(temp,"DEMAND",6))
    859         break;
    860       sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
    861       iHead--;
    862       iTail--;
    863       head[numberColumns]=iHead;
    864       tail[numberColumns]=iTail;
    865       ub[numberColumns]=iUb;
    866       cost[numberColumns]=iCost;
    867       numberColumns++;
    868     }
    869     if (strncmp(temp,"DEMAND",6)) {
    870       fprintf(stderr,"Unable to find DEMAND\n");
    871       exit(2);
    872     }
    873     // ***** Remember to convert to C notation
    874     while (fgets(temp,100,fp)) {
    875       int row;
    876       int value;
    877       if (!strncmp(temp,"END",3))
    878         break;
    879       sscanf(temp,"%d %d",&row,&value);
    880       upper[row-1]=value;
    881       lower[row-1]=value;
    882     }
    883     if (strncmp(temp,"END",3)) {
    884       fprintf(stderr,"Unable to find END\n");
    885       exit(2);
    886     }
    887     printf("Problem %d has %d rows and %d columns\n",problem,
    888            numberRows,numberColumns);
    889     fclose(fp);
    890     ClpSimplex  model;
    891     // now build model
    892    
    893     double * objective =new double[numberColumns];
    894     double * lowerColumn = new double[numberColumns];
    895     double * upperColumn = new double[numberColumns];
    896    
    897     double * element = new double [2*numberColumns];
    898     int * start = new int[numberColumns+1];
    899     int * row = new int[2*numberColumns];
    900     start[numberColumns]=2*numberColumns;
     974    start[0]=0;
    901975    for (i=0;i<numberColumns;i++) {
    902       start[i]=2*i;
    903       element[2*i]=-1.0;
    904       element[2*i+1]=1.0;
    905       row[2*i]=head[i];
    906       row[2*i+1]=tail[i];
    907       lowerColumn[i]=0.0;
    908       upperColumn[i]=ub[i];
    909       objective[i]=cost[i];
    910     }
    911     // Create Packed Matrix
    912     CoinPackedMatrix matrix;
    913     int * lengths = NULL;
    914     matrix.assignMatrix(true,numberRows,numberColumns,
    915                         2*numberColumns,element,row,start,lengths);
    916     // load model
    917    
    918     model.loadProblem(matrix,
    919                       lowerColumn,upperColumn,objective,
    920                       lower,upper);
    921    
    922     model.factorization()->maximumPivots(200+model.numberRows()/100);
    923     model.createStatus();
    924     double time1 = cpuTime();
    925     model.dual();
    926     std::cout<<"Network problem, ClpPackedMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
    927     ClpPlusMinusOneMatrix plusMinus(matrix);
    928     assert (plusMinus.getIndices()); // would be zero if not +- one
    929     model.loadProblem(plusMinus,
    930                       lowerColumn,upperColumn,objective,
    931                       lower,upper);
    932    
    933     model.factorization()->maximumPivots(200+model.numberRows()/100);
    934     model.createStatus();
    935     time1 = cpuTime();
    936     model.dual();
    937     std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
    938     ClpNetworkMatrix network(numberColumns,head,tail);
    939     model.loadProblem(network,
    940                       lowerColumn,upperColumn,objective,
    941                       lower,upper);
    942    
    943     model.factorization()->maximumPivots(200+model.numberRows()/100);
    944     model.createStatus();
    945     time1 = cpuTime();
    946     model.dual();
    947     std::cout<<"Network problem, ClpNetworkMatrix took "<<cpuTime()-time1<<" seconds"<<std::endl;
    948     delete [] lower;
    949     delete [] upper;
    950     delete [] head;
    951     delete [] tail;
    952     delete [] ub;
    953     delete [] cost;
    954     delete [] objective;
    955     delete [] lowerColumn;
    956     delete [] upperColumn;
    957   }
    958 #endif
     976      start[i+1]=i+1;
     977      column[i]=i;
     978      element[i]=1.0e-6;
     979      element[i]=0.0;
     980    }
     981    // Load up objective
     982    solution.loadQuadraticObjective(numberColumns,start,column,element);
     983    delete [] start;
     984    delete [] column;
     985    delete [] element;
     986    solution.quadraticSLP(50,1.0e-4);
     987    double objValue = solution.getObjValue();
     988    CoinRelFltEq eq(1.0e-4);
     989    assert(eq(objValue,-8.42857));
     990    solution.setLogLevel(63);
     991    solution.quadraticWolfe();
     992    objValue = solution.getObjValue();
     993    assert(eq(objValue,-8.42857));
     994    //exit(77);
     995  }
    959996  // Test quadratic
    960   {   
     997  if (1) {   
    961998    CoinMpsIO m;
    962999    std::string fn = mpsDir+"share2qp";
     1000    fn = "share2qpb";
    9631001    m.readMps(fn.c_str(),"mps");
    9641002    ClpSimplex solution;
     
    9771015    delete [] column;
    9781016    delete [] element;
     1017    //solution.quadraticSLP(50,1.0e-4);
     1018    solution.setLogLevel(63);
     1019    solution.quadraticWolfe();
     1020    double objValue = solution.getObjValue();
     1021    CoinRelFltEq eq(1.0e-4);
     1022    assert(eq(objValue,-400.92));
     1023  }
     1024  if (1) {   
     1025    CoinMpsIO m;
     1026    std::string fn = "./beale";
     1027    //fn = "./jensen";
     1028    m.readMps(fn.c_str(),"mps");
     1029    ClpSimplex solution;
     1030    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
     1031                         m.getObjCoefficients(),
     1032                         m.getRowLower(),m.getRowUpper());
     1033    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
     1034    solution.dual();
     1035    // get quadratic part
     1036    int * start=NULL;
     1037    int * column = NULL;
     1038    double * element = NULL;
     1039    m.readQuadraticMps(NULL,start,column,element,2);
     1040    // Load up objective
     1041    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
     1042    delete [] start;
     1043    delete [] column;
     1044    delete [] element;
     1045    solution.quadraticWolfe();
    9791046    solution.quadraticSLP(50,1.0e-4);
    9801047    double objValue = solution.getObjValue();
    981     CoinRelFltEq eq(1.0e-2);
    982     assert(eq(objValue,-400.92));
    983   }
    984  
     1048    CoinRelFltEq eq(1.0e-4);
     1049    assert(eq(objValue,0.5));
     1050    solution.quadraticWolfe();
     1051    objValue = solution.getObjValue();
     1052    assert(eq(objValue,0.5));
     1053  }
     1054#endif 
    9851055}
  • trunk/include/ClpNonLinearCost.hpp

    r162 r166  
    4040  /** Constructor from simplex.
    4141      This will just set up wasteful arrays for linear, but
    42       later may do dual analysis and even finding duplicate columns
    43   */
    44   ClpNonLinearCost(ClpSimplex * model);
     42      later may do dual analysis and even finding duplicate columns .
     43      If for quadratic then free varaibles get extra 0,0 bits
     44  */
     45  ClpNonLinearCost(ClpSimplex * model,bool forQuadratic=false);
    4546  /** Constructor from simplex and list of non-linearities (columns only)
    4647      First lower of each column has to match real lower
     
    137138  inline double cost(int sequence) const
    138139  { return cost_[whichRange_[sequence]+offset_[sequence]];};
     140  /// Sets inside bounds (i.e. non infinite - used in QP
     141  void setBounds(int sequence, double lower, double upper);
    139142  //@}
    140143
  • trunk/include/ClpSimplex.hpp

    r162 r166  
    144144  */
    145145  int quadraticSLP(int numberPasses,double deltaTolerance);
    146   /// Solves quadratic using Beale's algorithm - primal
    147   int quadraticBeale();
     146  /// Solves quadratic using Wolfe's algorithm - primal
     147  int quadraticWolfe();
    148148  /// Passes in factorization
    149149  void setFactorization( ClpFactorization & factorization);
  • trunk/include/ClpSimplexPrimalQuadratic.hpp

    r152 r166  
    1515/** This solves LPs using the primal simplex method
    1616
    17     It inherits from ClpSimplex.  It has no data of its own and
    18     is never created - only cast from a ClpSimplex object at algorithm time.
     17    It inherits from ClpSimplexPrimal.  It has no data of its own and
     18    is never created - only cast from a ClpSimplexPrimal object at algorithm time.
     19    If needed create new class and pass around
    1920
    2021*/
     
    3738  /// A sequential LP method
    3839  int primalSLP(int numberPasses, double deltaTolerance);
    39   /// Beale's method
    40   int primalBeale();
     40  /// Wolfe's method (actually a mixture with Jensen and King)
     41  int primalWolfe();
     42  /// This is done after first pass
     43  int primalWolfe2 (const ClpSimplexPrimalQuadratic * originalModel);
     44  /// Main part
     45  int whileIterating (const ClpSimplexPrimalQuadratic * originalModel,
     46                      int & sequenceIn,
     47                      int & crucialSj);
     48  /**
     49      Row array has pivot column
     50      This chooses pivot row.
     51      Rhs array is used for distance to next bound (for speed)
     52      For speed, we may need to go to a bucket approach when many
     53      variables go through bounds
     54      On exit rhsArray will have changes in costs of basic variables
     55      Initially no go thru
     56      Returns 0 - can do normal iteration
     57      1 - losing complementarity
     58  */
     59  int primalRow(CoinIndexedVector * rowArray,
     60                CoinIndexedVector * rhsArray,
     61                CoinIndexedVector * spareArray,
     62                CoinIndexedVector * spareArray2,
     63                const ClpSimplexPrimalQuadratic * originalModel,
     64                int crucialSj);
    4165  //@}
    4266
Note: See TracChangeset for help on using the changeset viewer.