Changeset 196 for branches


Ignore:
Timestamp:
Aug 13, 2003 6:43:02 AM (16 years ago)
Author:
forrest
Message:

64 bit, Sbb etc

Location:
branches/pre
Files:
16 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpDualRowSteepest.cpp

    r180 r196  
    263263    double * array = alternateWeights_->denseVector();
    264264    int * which = alternateWeights_->getIndices();
     265    int i;
    265266    for (i=0;i<numberRows;i++) {
    266267      double value=0.0;
     
    280281      if (fabs(weights_[i]-value)>1.0e-4)
    281282        printf("%d old %g, true %g\n",i,weights_[i],value);
     283      //else
     284      //printf("%d matches %g\n",i,value);
    282285    }
    283286    delete temp;
     
    312315    }
    313316    spare->setNumElements ( numberNonZero );
     317    // Only one array active as already permuted
    314318    model_->factorization()->updateColumn(spare,spare,true);
    315319    // permute back
     
    317321    // alternateWeights_ should still be empty
    318322    int pivotRow = model_->pivotRow();
     323#ifdef CLP_DEBUG
     324    if ( model_->logLevel (  ) >4  &&
     325         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm))
     326      printf("on row %d, true weight %g, old %g\n",
     327             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
     328#endif
    319329    // could re-initialize here (could be expensive)
    320330    norm /= model_->alpha() * model_->alpha();
  • branches/pre/ClpFactorization.cpp

    r180 r196  
    450450{
    451451#ifdef CLP_DEBUG
    452   regionSparse->checkClear();
     452  if (!noPermute)
     453    regionSparse->checkClear();
    453454#endif
    454455  if (!networkBasis_) {
  • branches/pre/ClpModel.cpp

    r183 r196  
    2020#include "ClpMessage.hpp"
    2121#include "ClpLinearObjective.hpp"
     22#include "ClpQuadraticObjective.hpp"
    2223
    2324//#############################################################################
     
    4243  matrix_(NULL),
    4344  rowCopy_(NULL),
    44   quadraticObjective_(NULL),
    4545  ray_(NULL),
    4646  status_(NULL),
     
    106106  delete rowCopy_;
    107107  rowCopy_=NULL;
    108   delete quadraticObjective_;
    109   quadraticObjective_ = NULL;
    110108  delete [] ray_;
    111109  ray_ = NULL;
     
    373371      rowCopy_=NULL;
    374372    }
    375     if (rhs.quadraticObjective_) {
    376       quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
    377     } else {
    378       quadraticObjective_=NULL;
    379     }
    380373    matrix_=NULL;
    381374    if (rhs.matrix_) {
     
    395388    matrix_ = rhs.matrix_;
    396389    rowCopy_ = NULL;
    397     quadraticObjective_ = rhs.quadraticObjective_;
    398390    ray_ = rhs.ray_;
    399391    lengthNames_ = 0;
     
    441433  matrix_ = NULL;
    442434  rowCopy_ = NULL;
    443   quadraticObjective_=NULL,
    444435  delete [] otherModel.ray_;
    445436  otherModel.ray_ = ray_;
     
    677668      which[i-newNumberColumns]=i;
    678669    matrix_->deleteCols(numberColumns_-newNumberColumns,which);
    679     if (quadraticObjective_) {
    680       quadraticObjective_->deleteCols(numberColumns_-newNumberColumns,which);
    681       quadraticObjective_->deleteRows(numberColumns_-newNumberColumns,which);
    682     }
    683670    delete [] which;
    684671  }
     
    752739  matrix_->deleteCols(number,which);
    753740  //matrix_->removeGaps();
    754   if (quadraticObjective_) {
    755     quadraticObjective_->deleteCols(number,which);
    756     quadraticObjective_->deleteRows(number,which);
    757   }
    758741  // status
    759742  if (status_) {
     
    12091192{
    12101193  assert (numberColumns==numberColumns_);
    1211   quadraticObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
    1212                                              start[numberColumns],element,column,start,NULL);
     1194  assert ((dynamic_cast< ClpLinearObjective*>(objective_)));
     1195  double offset;
     1196  ClpObjective * obj = new ClpQuadraticObjective(objective_->gradient(NULL,offset),numberColumns,
     1197                                             start,column,element);
     1198  delete objective_;
     1199  objective_ = obj;
     1200
    12131201}
    12141202void
     
    12161204{
    12171205  assert (matrix.getNumCols()==numberColumns_);
    1218   quadraticObjective_ = new CoinPackedMatrix(matrix);
     1206  assert ((dynamic_cast< ClpLinearObjective*>(objective_)));
     1207  double offset;
     1208  ClpQuadraticObjective * obj =
     1209    new ClpQuadraticObjective(objective_->gradient(NULL,offset),numberColumns_,
     1210                                                 NULL,NULL,NULL);
     1211  delete objective_;
     1212  objective_ = obj;
     1213  obj->loadQuadraticObjective(matrix);
    12191214}
    12201215// Get rid of quadratic objective
     
    12221217ClpModel::deleteQuadraticObjective()
    12231218{
    1224   delete quadraticObjective_;
    1225   quadraticObjective_ = NULL;
     1219  ClpQuadraticObjective * obj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     1220  if (obj)
     1221    obj->deleteQuadraticObjective();
    12261222}
    12271223void
     
    14171413    rowCopy_=NULL;
    14181414  }
    1419   if (rhs->quadraticObjective_) {
    1420     quadraticObjective_ = new CoinPackedMatrix(*rhs->quadraticObjective_,
    1421                                                numberColumns,whichColumn,
    1422                                                numberColumns,whichColumn);
    1423   } else {
    1424     quadraticObjective_=NULL;
    1425   }
    14261415  matrix_=NULL;
    14271416  if (rhs->matrix_) {
  • branches/pre/ClpPrimalColumnDantzig.cpp

    r180 r196  
    8585    anyUpdates=false;
    8686  }
    87 
     87  if (!updates->getNumElements())
     88    anyUpdates=false; // in case dj 0.0 (for quadratic)
    8889  if (anyUpdates) {
    8990    model_->factorization()->updateColumnTranspose(spareRow2,updates);
  • branches/pre/ClpPrimalColumnSteepest.cpp

    r180 r196  
    414414  }
    415415
    416 #ifdef RECOMPUTE_ALL_DJS
    417   for (iSection=0;iSection<2;iSection++) {
     416  // If in quadratic re-compute all
     417  if (model_->algorithm()==2) {
     418    for (iSection=0;iSection<2;iSection++) {
    418419     
    419     reducedCost=model_->djRegion(iSection);
    420     int addSequence;
    421     int iSequence;
     420      reducedCost=model_->djRegion(iSection);
     421      int addSequence;
     422      int iSequence;
    422423     
    423     if (!iSection) {
    424       number = model_->numberRows();
    425       addSequence = model_->numberColumns();
    426     } else {
    427       number = model_->numberColumns();
    428       addSequence = 0;
    429     }
    430 
    431     if (!model_->nonLinearCost()->lookBothWays()) {
    432       for (iSequence=0;iSequence<number;iSequence++) {
    433         double value = reducedCost[iSequence];
    434         ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    435 
    436         switch(status) {
     424      if (!iSection) {
     425        number = model_->numberRows();
     426        addSequence = model_->numberColumns();
     427      } else {
     428        number = model_->numberColumns();
     429        addSequence = 0;
     430      }
     431     
     432      if (!model_->nonLinearCost()->lookBothWays()) {
     433        for (iSequence=0;iSequence<number;iSequence++) {
     434          double value = reducedCost[iSequence];
     435          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    437436         
    438         case ClpSimplex::basic:
    439           infeasible_->zero(iSequence+addSequence);
    440         case ClpSimplex::isFixed:
    441           break;
    442         case ClpSimplex::isFree:
    443           if (fabs(value)>tolerance) {
    444             // we are going to bias towards free (but only if reasonable)
    445             value *= FREE_BIAS;
    446             // store square in list
    447             if (infeas[iSequence+addSequence])
    448               infeas[iSequence+addSequence] = value*value; // already there
    449             else
    450               infeasible_->quickAdd(iSequence+addSequence,value*value);
    451           } else {
     437          switch(status) {
     438           
     439          case ClpSimplex::basic:
    452440            infeasible_->zero(iSequence+addSequence);
    453           }
    454           break;
    455         case ClpSimplex::atUpperBound:
    456           if (value>tolerance) {
    457             // store square in list
    458             if (infeas[iSequence+addSequence])
    459               infeas[iSequence+addSequence] = value*value; // already there
    460             else
    461               infeasible_->quickAdd(iSequence+addSequence,value*value);
    462           } else {
    463             infeasible_->zero(iSequence+addSequence);
    464           }
    465           break;
    466         case ClpSimplex::atLowerBound:
    467           if (value<-tolerance) {
    468             // store square in list
    469             if (infeas[iSequence+addSequence])
    470               infeas[iSequence+addSequence] = value*value; // already there
    471             else
    472               infeasible_->quickAdd(iSequence+addSequence,value*value);
    473           } else {
    474             infeasible_->zero(iSequence+addSequence);
    475           }
    476         }
    477       }
    478     } else {
    479       // we can go both ways
    480       ClpNonLinearCost * nonLinear = model_->nonLinearCost();
    481       for (iSequence=0;iSequence<number;iSequence++) {
    482         double value = reducedCost[iSequence];
    483         ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    484 
    485         switch(status) {
    486          
    487         case ClpSimplex::basic:
    488           infeasible_->zero(iSequence+addSequence);
    489         case ClpSimplex::isFixed:
    490           break;
    491         case ClpSimplex::isFree:
    492           if (fabs(value)>tolerance) {
    493             // we are going to bias towards free (but only if reasonable)
    494             value *= FREE_BIAS;
    495             // store square in list
    496             if (infeas[iSequence+addSequence])
    497               infeas[iSequence+addSequence] = value*value; // already there
    498             else
    499               infeasible_->quickAdd(iSequence+addSequence,value*value);
    500           } else {
    501             infeasible_->zero(iSequence+addSequence);
    502           }
    503           break;
    504         case ClpSimplex::atUpperBound:
    505           if (value>tolerance) {
    506             // store square in list
    507             if (infeas[iSequence+addSequence])
    508               infeas[iSequence+addSequence] = value*value; // already there
    509             else
    510               infeasible_->quickAdd(iSequence+addSequence,value*value);
    511           } else {
    512             // look other way - change up should be negative
    513             value -= nonLinear->changeUpInCost(iSequence+addSequence);
    514             if (value>-tolerance) {
    515               infeasible_->zero(iSequence+addSequence);
    516             } else {
     441          case ClpSimplex::isFixed:
     442            break;
     443          case ClpSimplex::isFree:
     444          case ClpSimplex::superBasic:
     445            if (fabs(value)>tolerance) {
     446              // we are going to bias towards free (but only if reasonable)
     447              value *= FREE_BIAS;
    517448              // store square in list
    518449              if (infeas[iSequence+addSequence])
     
    520451              else
    521452                infeasible_->quickAdd(iSequence+addSequence,value*value);
    522             }
    523           }
    524           break;
    525         case ClpSimplex::atLowerBound:
    526           if (value<-tolerance) {
    527             // store square in list
    528             if (infeas[iSequence+addSequence])
    529               infeas[iSequence+addSequence] = value*value; // already there
    530             else
    531               infeasible_->quickAdd(iSequence+addSequence,value*value);
    532           } else {
    533             // look other way - change down should be positive
    534             value -= nonLinear->changeDownInCost(iSequence+addSequence);
    535             if (value<tolerance) {
     453            } else {
    536454              infeasible_->zero(iSequence+addSequence);
    537             } else {
     455            }
     456            break;
     457          case ClpSimplex::atUpperBound:
     458            if (value>tolerance) {
    538459              // store square in list
    539460              if (infeas[iSequence+addSequence])
     
    541462              else
    542463                infeasible_->quickAdd(iSequence+addSequence,value*value);
     464            } else {
     465              infeasible_->zero(iSequence+addSequence);
     466            }
     467            break;
     468          case ClpSimplex::atLowerBound:
     469            if (value<-tolerance) {
     470              // store square in list
     471              if (infeas[iSequence+addSequence])
     472                infeas[iSequence+addSequence] = value*value; // already there
     473              else
     474                infeasible_->quickAdd(iSequence+addSequence,value*value);
     475            } else {
     476              infeasible_->zero(iSequence+addSequence);
    543477            }
    544478          }
    545479        }
    546       }
    547     }
    548   }
    549 #endif
     480      } else {
     481        // we can go both ways
     482        ClpNonLinearCost * nonLinear = model_->nonLinearCost();
     483        for (iSequence=0;iSequence<number;iSequence++) {
     484          double value = reducedCost[iSequence];
     485          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     486         
     487          switch(status) {
     488           
     489          case ClpSimplex::basic:
     490            infeasible_->zero(iSequence+addSequence);
     491          case ClpSimplex::isFixed:
     492            break;
     493          case ClpSimplex::isFree:
     494          case ClpSimplex::superBasic:
     495            if (fabs(value)>tolerance) {
     496              // we are going to bias towards free (but only if reasonable)
     497              value *= FREE_BIAS;
     498              // store square in list
     499              if (infeas[iSequence+addSequence])
     500                infeas[iSequence+addSequence] = value*value; // already there
     501              else
     502                infeasible_->quickAdd(iSequence+addSequence,value*value);
     503            } else {
     504              infeasible_->zero(iSequence+addSequence);
     505            }
     506            break;
     507          case ClpSimplex::atUpperBound:
     508            if (value>tolerance) {
     509              // store square in list
     510              if (infeas[iSequence+addSequence])
     511                infeas[iSequence+addSequence] = value*value; // already there
     512              else
     513                infeasible_->quickAdd(iSequence+addSequence,value*value);
     514            } else {
     515              // look other way - change up should be negative
     516              value -= nonLinear->changeUpInCost(iSequence+addSequence);
     517              if (value>-tolerance) {
     518                infeasible_->zero(iSequence+addSequence);
     519              } else {
     520                // store square in list
     521                if (infeas[iSequence+addSequence])
     522                  infeas[iSequence+addSequence] = value*value; // already there
     523                else
     524                  infeasible_->quickAdd(iSequence+addSequence,value*value);
     525              }
     526            }
     527            break;
     528          case ClpSimplex::atLowerBound:
     529            if (value<-tolerance) {
     530              // store square in list
     531              if (infeas[iSequence+addSequence])
     532                infeas[iSequence+addSequence] = value*value; // already there
     533              else
     534                infeasible_->quickAdd(iSequence+addSequence,value*value);
     535            } else {
     536              // look other way - change down should be positive
     537              value -= nonLinear->changeDownInCost(iSequence+addSequence);
     538              if (value<tolerance) {
     539                infeasible_->zero(iSequence+addSequence);
     540              } else {
     541                // store square in list
     542                if (infeas[iSequence+addSequence])
     543                  infeas[iSequence+addSequence] = value*value; // already there
     544                else
     545                  infeasible_->quickAdd(iSequence+addSequence,value*value);
     546              }
     547            }
     548          }
     549        }
     550      }
     551    }
     552  }
    550553  // for weights update we use pivotSequence
    551554  pivotRow = pivotSequence_;
     
    704707    //weight=1.0;
    705708    if (value>bestDj*weight&&value>tolerance) {
    706       // check flagged variable
     709      // check flagged variable and correct dj
    707710      if (!model_->flagged(iSequence)) {
    708711        /*if (model_->numberIterations()%100==0)
     
    11061109      int nWords = (number+31)>>5;
    11071110      reference_ = new unsigned int[nWords];
    1108       assert (sizeof(unsigned int)==4);
    11091111      // tiny overhead to zero out (stops valgrind complaining)
    11101112      memset(reference_,0,nWords*sizeof(int));
  • branches/pre/ClpPrimalQuadraticDantzig.cpp

    r190 r196  
    9090  //double tolerance=model_->currentPrimalTolerance();
    9191
    92   double bestDj = model_->dualTolerance();
    93   double bestInfeasibleDj = 10.0*model_->dualTolerance();
     92  double dualTolerance = model_->dualTolerance()*1000.0;
     93  double bestDj = dualTolerance;
     94  double bestInfeasibleDj = 10.0*dualTolerance;
    9495  int bestSequence=-1;
    9596  double dualIn=0.0;
     
    105106  const int * columnLength = matrix->getVectorLengths();
    106107  const int * which = quadraticInfo_->quadraticSequence();
    107   const double * objective = quadraticInfo_->originalObjective();
     108  const double * objective = quadraticInfo_->linearObjective();
     109  const double * djWeight = quadraticInfo_->djWeight();
    108110  for (iSequence=0;iSequence<originalNumberColumns;iSequence++) {
    109111    if (model_->flagged(iSequence))
     
    122124      }
    123125    }
     126    double value2 = value*djWeight[iSequence];
     127    if (fabs(value2)>dualTolerance)
     128      value=value2;
     129    else if (value<-dualTolerance)
     130      value = -1.001*dualTolerance;
     131    else if (value>dualTolerance)
     132      value = 1.001*dualTolerance;
     133    if (djWeight[iSequence]<1.0e-6)
     134      value=value2;
    124135    ClpSimplex::Status status = model_->getStatus(iSequence);
    125136   
     
    133144        bestSequence = jSequence;
    134145        dualIn = value;
     146        abort();
    135147      }
    136148      break;
     
    171183    // for L rows pi negative okay so choose if positive
    172184    double value = solution[iPi];
     185    double value2 = value*djWeight[iSequence];
     186    if (fabs(value2)>dualTolerance)
     187      value=value2;
     188    else if (value<-dualTolerance)
     189      value = -1.001*dualTolerance;
     190    else if (value>dualTolerance)
     191      value = 1.001*dualTolerance;
     192    if (djWeight[iSequence]<1.0e-6)
     193      value=value2;
    173194    ClpSimplex::Status status = model_->getStatus(iSequence);
    174195    switch(status) {
     
    208229  if (bestSequence>=0) {
    209230    model_->djRegion()[bestSequence] = dualIn;
    210 #if 0
    211     // free up depending which bound
    212     if (solution[iSequence]<trueLower[iSequence]+tolerance) {
    213       model_->setColumnStatus(iSequence,ClpSimplex::atLowerBound);
    214       upper[iSequence]=trueUpper[iSequence];
    215     } else if (solution[iSequence]>trueUpper[iSequence]-tolerance) {
    216       model_->setColumnStatus(iSequence,ClpSimplex::atUpperBound);
    217       lower[iSequence]=trueLower[iSequence];
    218     } else {
    219       abort();
    220     }
    221 #endif
    222231  }
    223232  return bestSequence;
  • branches/pre/ClpSimplex.cpp

    r190 r196  
    26812681int ClpSimplex::strongBranching(int numberVariables,const int * variables,
    26822682                                double * newLower, double * newUpper,
     2683                                double ** outputSolution,
     2684                                int * outputStatus, int * outputIterations,
    26832685                                bool stopOnFirstInfeasible,
    26842686                                bool alwaysFinish)
    26852687{
    26862688  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
    2687                                                     newLower,  newUpper,
     2689                                                    newLower,  newUpper,outputSolution,
     2690                                                    outputStatus, outputIterations,
    26882691                                                    stopOnFirstInfeasible,
    26892692                                                    alwaysFinish);
     
    44844487static bool equalDouble(double value1, double value2)
    44854488{
    4486   assert(sizeof(unsigned int)*2==sizeof(double));
    44874489  unsigned int *i1 = (unsigned int *) &value1;
    44884490  unsigned int *i2 = (unsigned int *) &value2;
    4489   return (i1[0]==i2[0]&&i1[1]==i2[1]);
     4491  if (sizeof(unsigned int)*2==sizeof(double))
     4492    return (i1[0]==i2[0]&&i1[1]==i2[1]);
     4493  else
     4494    return (i1[0]==i2[0]);
    44904495}
    44914496int
     
    45804585      return -2;
    45814586    } else {
    4582       model_->messageHandler()->message(CLP_LOOP,model_->messages())
    4583         <<CoinMessageEol;
     4587      // look at solution and maybe declare victory
     4588      if (infeasibility<1.0e-4) {
     4589        return 0;
     4590      } else {
     4591        model_->messageHandler()->message(CLP_LOOP,model_->messages())
     4592          <<CoinMessageEol;
    45844593#ifndef NDEBUG
    4585       abort();
     4594        abort();
    45864595#endif
    4587       // look at solution and maybe declare victory
    4588       if (infeasibility<1.0e-4)
    4589         return 0;
    4590       else
    45914596        return 3;
     4597      }
    45924598    }
    45934599  }
  • branches/pre/ClpSimplexDual.cpp

    r195 r196  
    29432943int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
    29442944                                    double * newLower, double * newUpper,
     2945                                    double ** outputSolution,
     2946                                    int * outputStatus, int * outputIterations,
    29452947                                    bool stopOnFirstInfeasible,
    29462948                                    bool alwaysFinish)
     
    29692971  factorization_->slackValue(-1.0);
    29702972  factorization_->zeroTolerance(1.0e-13);
    2971   // save if sparse factorization wanted
    2972   int saveSparse = factorization_->sparseThreshold();
    29732973
    29742974  int factorizationStatus = internalFactorize(0);
     
    29792979      <<factorizationStatus
    29802980      <<CoinMessageEol;
    2981   if (saveSparse) {
    2982     // use default at present
    2983     factorization_->sparseThreshold(0);
    2984     factorization_->goSparse();
    2985   }
    29862981 
    29872982  // save stuff
     
    30083003  // need to save/restore weights.
    30093004
     3005  int iSolution = 0;
    30103006  for (i=0;i<numberVariables;i++) {
    30113007    int iColumn = variables[i];
     
    30243020    // Start of fast iterations
    30253021    int status = fastDual(alwaysFinish);
    3026 
     3022    if (status) {
     3023      // not finished - might be optimal
     3024      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
     3025      double limit = 0.0;
     3026      getDblParam(ClpDualObjectiveLimit, limit);
     3027      if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
     3028          optimizationDirection_*limit) {
     3029        problemStatus_=0;
     3030      }
     3031      status=problemStatus_;
     3032    }
     3033
     3034    if (scalingFlag_<=0) {
     3035      memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
     3036    } else {
     3037      int j;
     3038      double * sol = outputSolution[iSolution];
     3039      for (j=0;j<numberColumns_;j++)
     3040        sol[j] = solution_[j]*columnScale_[j];
     3041    }
     3042    outputStatus[iSolution]=status;
     3043    outputIterations[iSolution]=numberIterations_;
     3044    iSolution++;
    30273045    // restore
    30283046    memcpy(solution_,saveSolution,
     
    30613079    // Start of fast iterations
    30623080    status = fastDual(alwaysFinish);
     3081    if (status) {
     3082      // not finished - might be optimal
     3083      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
     3084      double limit = 0.0;
     3085      getDblParam(ClpDualObjectiveLimit, limit);
     3086      if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
     3087          optimizationDirection_*limit) {
     3088        problemStatus_=0;
     3089      }
     3090      status=problemStatus_;
     3091    }
     3092    if (scalingFlag_<=0) {
     3093      memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
     3094    } else {
     3095      int j;
     3096      double * sol = outputSolution[iSolution];
     3097      for (j=0;j<numberColumns_;j++)
     3098        sol[j] = solution_[j]*columnScale_[j];
     3099    }
     3100    outputStatus[iSolution]=status;
     3101    outputIterations[iSolution]=numberIterations_;
     3102    iSolution++;
    30633103
    30643104    // restore
     
    31213161  delete [] savePivot;
    31223162
    3123   factorization_->sparseThreshold(saveSparse);
    31243163  // Get rid of some arrays and empty factorization
    31253164  deleteRim();
     
    31383177  }
    31393178  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
     3179  int iSolution =0;
    31403180  for (i=0;i<numberVariables;i++) {
    31413181    int iColumn = variables[i];
     
    31463186    double saveUpper = columnUpper_[iColumn];
    31473187    columnUpper_[iColumn] = newUpper[i];
    3148     dual();
     3188    int status=dual(0);
     3189    memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
     3190    outputStatus[iSolution]=status;
     3191    outputIterations[iSolution]=numberIterations_;
     3192    iSolution++;
    31493193
    31503194    // restore
     
    31683212    double saveLower = columnLower_[iColumn];
    31693213    columnLower_[iColumn] = newLower[i];
    3170     dual();
     3214    status=dual(0);
     3215    memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
     3216    outputStatus[iSolution]=status;
     3217    outputIterations[iSolution]=numberIterations_;
     3218    iSolution++;
    31713219
    31723220    // restore
     
    32243272{
    32253273  algorithm_ = -1;
     3274  // save data
     3275  ClpDataSave data = saveData();
    32263276  dualTolerance_=dblParam_[ClpDualTolerance];
    32273277  primalTolerance_=dblParam_[ClpPrimalTolerance];
     
    32393289  problemStatus_ = -1;
    32403290  numberIterations_=0;
     3291  factorization_->sparseThreshold(0);
     3292  factorization_->goSparse();
    32413293
    32423294  int lastCleaned=0; // last time objective or bounds cleaned up
     
    32703322  int returnCode = 0;
    32713323
     3324  int iRow,iColumn;
    32723325  while (problemStatus_<0) {
    3273     int iRow,iColumn;
    32743326    // clear
    32753327    for (iRow=0;iRow<4;iRow++) {
     
    32953347      double * givenPi=NULL;
    32963348      returnCode = whileIterating(givenPi);
    3297       if (!alwaysFinish&&returnCode<1) {
     3349      if (!alwaysFinish&&(returnCode<1||returnCode==3)) {
    32983350        double limit = 0.0;
    32993351        getDblParam(ClpDualObjectiveLimit, limit);
     
    33013353           optimizationDirection_*limit||
    33023354           numberAtFakeBound()) {
     3355          returnCode=1;
    33033356          // can't say anything interesting - might as well return
    33043357#ifdef CLP_DEBUG
     
    33063359                 numberIterations_,returnCode);
    33073360#endif
    3308           returnCode=1;
    33093361          break;
    33103362        }
     
    33143366  }
    33153367
     3368  // clear
     3369  for (iRow=0;iRow<4;iRow++) {
     3370    rowArray_[iRow]->clear();
     3371  }   
     3372 
     3373  for (iColumn=0;iColumn<2;iColumn++) {
     3374    columnArray_[iColumn]->clear();
     3375  }   
    33163376  assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
     3377  // Restore any saved stuff
     3378  restoreData(data);
    33173379  dualBound_ = saveDualBound;
    33183380  return returnCode;
  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r191 r196  
    99#include "ClpSimplexPrimalQuadratic.hpp"
    1010#include "ClpPrimalQuadraticDantzig.hpp"
     11#include "ClpPrimalColumnDantzig.hpp"
     12#include "ClpQuadraticObjective.hpp"
    1113#include "ClpFactorization.hpp"
    1214#include "ClpNonLinearCost.hpp"
     
    7375
    7476  // Get list of non linear columns
    75   CoinPackedMatrix * quadratic = quadraticObjective();
     77  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     78  CoinPackedMatrix * quadratic = NULL;
     79  if (quadraticObj)
     80    quadratic = quadraticObj->quadraticObjective();
    7681  if (!quadratic) {
    7782    // no quadratic part
     
    496501#endif 
    497502  // Now do quadratic
    498   ClpPrimalQuadraticDantzig dantzigP(model2,&info);
    499   model2->setPrimalColumnPivotAlgorithm(dantzigP);
     503  //ClpPrimalQuadraticDantzig dantzigP(model2,&info);
     504  //ClpPrimalColumnDantzig dantzigP;
     505  //model2->setPrimalColumnPivotAlgorithm(dantzigP);
    500506  model2->messageHandler()->setLogLevel(63);
    501507  model2->primalQuadratic2(&info,phase);
     
    520526  const int * which = info->quadraticSequence();
    521527  //const int * back = info->backSequence();
    522   assert (!scalingFlag_);
    523528  // initialize - values pass coding and algorithm_ is +1
    524529  if (!startup(1)) {
     
    529534    info->setSequenceIn(-1);
    530535    info->setCrucialSj(-1);
     536    bool deleteCosts=false;
     537    if (scalingFlag_) {
     538      // get own copy of original objective and scale
     539      ClpQuadraticObjective * quadraticObj = new ClpQuadraticObjective(*info->originalObjective());
     540      info->setOriginalObjective(quadraticObj);
     541      CoinPackedMatrix * quadratic = info->quadraticObjective();
     542      double * objective = info->linearObjective();
     543      // replace column by column
     544      double * newElement = new double[numberXColumns];
     545      // scale column copy
     546      // get matrix data pointers
     547      const int * row = quadratic->getIndices();
     548      const CoinBigIndex * columnStart = quadratic->getVectorStarts();
     549      const int * columnLength = quadratic->getVectorLengths();
     550      const double * elementByColumn = quadratic->getElements();
     551      double direction = optimizationDirection_;
     552      // direction is actually scale out not scale in
     553      if (direction)
     554        direction = 1.0/direction;
     555      int iColumn;
     556      for (iColumn=0;iColumn<numberXColumns;iColumn++) {
     557        int j;
     558        double scale = columnScale_[iColumn]*direction;
     559        const double * elementsInThisColumn = elementByColumn + columnStart[iColumn];
     560        const int * columnsInThisColumn = row + columnStart[iColumn];
     561        int number = columnLength[iColumn];
     562        assert (number<=numberXColumns);
     563        for (j=0;j<number;j++) {
     564          int iColumn2 = columnsInThisColumn[j];
     565          newElement[j] = elementsInThisColumn[j]*scale*columnScale_[iColumn2];
     566        }
     567        quadratic->replaceVector(iColumn,number,newElement);
     568        objective[iColumn] *= direction*columnScale_[iColumn];
     569      }
     570      delete [] newElement;
     571      deleteCosts=true;
     572    } else if (optimizationDirection_!=1.0) {
     573      // get own copy of original objective and scale
     574      ClpQuadraticObjective * quadraticObj = new ClpQuadraticObjective(*info->originalObjective());
     575      info->setOriginalObjective(quadraticObj);
     576      CoinPackedMatrix * quadratic = info->quadraticObjective();
     577      double * objective = info->linearObjective();
     578      // replace column by column
     579      double * newElement = new double[numberXColumns];
     580      // get matrix data pointers
     581      const CoinBigIndex * columnStart = quadratic->getVectorStarts();
     582      const int * columnLength = quadratic->getVectorLengths();
     583      const double * elementByColumn = quadratic->getElements();
     584      double direction = optimizationDirection_;
     585      // direction is actually scale out not scale in
     586      if (direction)
     587        direction = 1.0/direction;
     588      int iColumn;
     589      for (iColumn=0;iColumn<numberXColumns;iColumn++) {
     590        int j;
     591        const double * elementsInThisColumn = elementByColumn + columnStart[iColumn];
     592        int number = columnLength[iColumn];
     593        assert (number<=numberXColumns);
     594        for (j=0;j<number;j++) {
     595          newElement[j] = elementsInThisColumn[j]*direction;
     596        }
     597        quadratic->replaceVector(iColumn,number,newElement);
     598        objective[iColumn] *= direction;
     599      }
     600      delete [] newElement;
     601      deleteCosts=true;
     602    }
    531603   
    532604    // Say no pivot has occurred (for steepest edge and updates)
     
    568640        break; // declare victory
    569641     
    570       // Compute objective function from scratch
    571       const CoinPackedMatrix * quadratic =
    572         info->originalModel()->quadraticObjective();
    573       const int * columnQuadratic = quadratic->getIndices();
    574       const int * columnQuadraticStart = quadratic->getVectorStarts();
    575       const int * columnQuadraticLength = quadratic->getVectorLengths();
    576       const double * quadraticElement = quadratic->getElements();
    577       const double * originalCost = info->originalObjective();
    578       objectiveValue_=0.0;
    579       for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    580         double valueI = solution_[iColumn];
    581         objectiveValue_ += valueI*originalCost[iColumn];
    582         int j;
    583         for (j=columnQuadraticStart[iColumn];
    584              j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    585           int jColumn = columnQuadratic[j];
    586           double valueJ = solution_[jColumn];
    587           double elementValue = 0.5*quadraticElement[j];
    588           objectiveValue_ += valueI*valueJ*elementValue;
    589         }
    590       }
     642      checkComplementarity (info,rowArray_[0],rowArray_[1]);
    591643
    592644      // Say good factorization
     
    656708
    657709      // exit if victory declared
     710      dualIn_=0.0; // so no updates
    658711      if (!info->currentPhase()&&info->sequenceIn()<0&&primalColumnPivot_->pivotColumn(rowArray_[1],
    659712                                          rowArray_[2],rowArray_[3],
     
    667720      problemStatus_=-1;
    668721      whileIterating(info);
     722    }
     723    if (deleteCosts) {
     724      // delete scaled copy of objective
     725      delete info->originalObjective();
    669726    }
    670727  }
     
    697754  double saveObjective = objectiveValue_;
    698755  int numberXColumns = info->numberXColumns();
    699   int numberXRows = info->numberXRows();
    700756  const int * which = info->quadraticSequence();
    701757  const double * pi = solution_+numberXColumns;
    702   // Matrix for linear stuff
    703   const int * row = matrix_->getIndices();
    704   const int * columnStart = matrix_->getVectorStarts();
    705   const double * element =  matrix_->getElements();
    706   const int * columnLength = matrix_->getVectorLengths();
    707758  int nextSequenceIn=info->sequenceIn();
    708759  int oldSequenceIn=nextSequenceIn;
     
    771822      } else {
    772823        saveSequenceIn=sequenceIn_;
     824        createDjs(info,rowArray_[1],rowArray_[2]);
     825        dualIn_=0.0; // so no updates
    773826        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
    774827                     columnArray_[0],columnArray_[1]);
     
    777830      saveSequenceIn=sequenceIn_;
    778831      sequenceIn_ = nextSequenceIn;
    779       if (sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)
    780         cleanupIteration=false;
    781       else
     832      if (phase==2) {
     833        if (sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)
     834          cleanupIteration=false;
     835        else
     836          cleanupIteration=true;
     837      } else {
    782838        cleanupIteration=true;
     839      }
    783840    }
    784841    pivotRow_=-1;
     
    786843    rowArray_[1]->clear();
    787844    if (sequenceIn_>=0) {
    788       if (numberIterations_>=8796000)
    789         printf("loxx 4721 %g %g\n",lower_[4721],upper_[4721]);
    790845      nextSequenceIn=-1;
    791846      // we found a pivot column
     
    793848      // do second half of iteration
    794849      while (chosen>=0) {
    795         {
    796           // Compute objective function from scratch
    797           const CoinPackedMatrix * quadratic =
    798             info->originalModel()->quadraticObjective();
    799           const int * columnQuadratic = quadratic->getIndices();
    800           const int * columnQuadraticStart = quadratic->getVectorStarts();
    801           const int * columnQuadraticLength = quadratic->getVectorLengths();
    802           const double * quadraticElement = quadratic->getElements();
    803           const double * originalCost = info->originalObjective();
    804           int iColumn;
    805           objectiveValue_=0.0;
    806           double infeasCost=0.0;
    807           double linearCost=0.0;
    808           for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    809             double valueI = solution_[iColumn];
    810             linearCost += valueI*originalCost[iColumn];
    811             double diff =cost_[iColumn]-originalCost[iColumn];
    812             // WE are always feasible so look at low not up!
    813             if (diff>0.0) {
    814               double above = valueI-lower_[iColumn];
    815               assert(above>=0.0);
    816               infeasCost += diff*above;
    817             } else if (diff<0.0) {
    818               double below = upper_[iColumn]-valueI;
    819               assert(below>=0.0);
    820               infeasCost -= diff*below;
    821             }
    822             int j;
    823             for (j=columnQuadraticStart[iColumn];
    824                  j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    825               int jColumn = columnQuadratic[j];
    826               double valueJ = solution_[jColumn];
    827               double elementValue = 0.5*quadraticElement[j];
    828               objectiveValue_ += valueI*valueJ*elementValue;
    829             }
    830           }
    831           int iRow;
    832           for (iRow=0;iRow<numberXRows;iRow++) {
    833             int iColumn = iRow + numberColumns_;
    834             double diff =cost_[iColumn];
    835             double valueI = solution_[iColumn];
    836             if (diff>0.0) {
    837               double above = valueI-lower_[iColumn];
    838               assert(above>=0.0);
    839               infeasCost += diff*above;
    840             } else if (diff<0.0) {
    841               double below = upper_[iColumn]-valueI;
    842               assert(below>=0.0);
    843               infeasCost -= diff*below;
    844             }
    845           }
    846           objectiveValue_ += linearCost;
    847           assert (infeasCost>=0.0);
    848           printf("True objective is %g, infeas cost %g, sum %g\n",
    849                  objectiveValue_,infeasCost,objectiveValue_+infeasCost);
    850         }
     850        checkComplementarity (info,rowArray_[3],rowArray_[1]);
     851        printf("True objective is %g, infeas cost %g, sum %g\n",
     852               objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost());
    851853        objectiveValue_=saveObjective;
    852854        returnCode=-1;
     
    859861        chosen=-1;
    860862        unpack(rowArray_[1]);
     863        // compute dj in case linear
     864        {
     865          dualIn_=cost_[sequenceIn_];
     866          int j;
     867          const double * element = rowArray_[1]->denseVector();
     868          int numberNonZero=rowArray_[1]->getNumElements();
     869          const int * whichRow = rowArray_[1]->getIndices();
     870          bool packed = rowArray_[1]->packedMode();
     871          if (packed) {
     872            for (j=0;j<numberNonZero; j++) {
     873              int iRow = whichRow[j];
     874              dualIn_ -= element[j]*pi[iRow];
     875            }
     876          } else {
     877            for (j=0;j<numberNonZero; j++) {
     878              int iRow = whichRow[j];
     879              dualIn_ -= element[iRow]*pi[iRow];
     880            }
     881          }
     882        }
    861883        factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
    862884        if (cleanupIteration) {
     
    918940            } else {
    919941              // need coding
    920               dualIn_=cost_[sequenceIn_];
    921               int j;
    922               for (j=columnStart[sequenceIn_];j<columnStart[sequenceIn_]+columnLength[sequenceIn_]; j++) {
    923                 int iRow = row[j];
    924                 dualIn_ -= element[j]*pi[iRow];
    925               }
     942              //dualIn_=cost_[sequenceIn_];
     943              //int j;
     944              //for (j=columnStart[sequenceIn_];j<columnStart[sequenceIn_]+columnLength[sequenceIn_]; j++) {
     945              //int iRow = row[j];
     946              //dualIn_ -= element[j]*pi[iRow];
     947              //}
    926948            }
    927949          } else {
     
    957979                             info,
    958980                             cleanupIteration);
    959         if (pivotRow_==-1&&phase==2&&fabs(dualIn_)<1.0e-3) {
     981        if (pivotRow_==-1&&(phase==2||cleanupIteration)&&fabs(dualIn_)<1.0e-3) {
    960982          // try other way
    961983          dualIn_=-dualIn_;
     
    963985          if (info->crucialSj()>=0)
    964986            setColumnStatus(info->crucialSj(),basic);
    965           rowArray_[3]->checkClear();
    966           rowArray_[2]->checkClear();
    967           rowArray_[0]->checkClear();
    968987          result=primalRow(rowArray_[1],rowArray_[3],
    969988                           rowArray_[2],rowArray_[0],
     
    11391158          printf("%d Sj value went from %g to %g\n",crucialSj,oldSj,solution_[info->crucialSj()]);
    11401159        }
    1141         {
    1142           double oldValue = objectiveValue_;
    1143           // Compute objective function from scratch
    1144           const CoinPackedMatrix * quadratic =
    1145             info->originalModel()->quadraticObjective();
    1146           const int * columnQuadratic = quadratic->getIndices();
    1147           const int * columnQuadraticStart = quadratic->getVectorStarts();
    1148           const int * columnQuadraticLength = quadratic->getVectorLengths();
    1149           const double * quadraticElement = quadratic->getElements();
    1150           const double * originalCost = info->originalObjective();
    1151           int iColumn;
    1152           objectiveValue_=0.0;
    1153           double infeasCost=0.0;
    1154           double linearCost=0.0;
    1155           for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    1156             double valueI = solution_[iColumn];
    1157             linearCost += valueI*originalCost[iColumn];
    1158             double diff =cost_[iColumn]-originalCost[iColumn];
    1159             // WE are always feasible so look at low not up!
    1160             if (diff>0.0) {
    1161               double above = valueI-lower_[iColumn];
    1162               assert(above>=0.0);
    1163               infeasCost += diff*above;
    1164             } else if (diff<0.0) {
    1165               double below = upper_[iColumn]-valueI;
    1166               assert(below>=0.0);
    1167               infeasCost -= diff*below;
    1168             }
    1169             int j;
    1170             for (j=columnQuadraticStart[iColumn];
    1171                  j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    1172               int jColumn = columnQuadratic[j];
    1173               double valueJ = solution_[jColumn];
    1174               double elementValue = 0.5*quadraticElement[j];
    1175               objectiveValue_ += valueI*valueJ*elementValue;
    1176             }
     1160        double oldObjectiveValue = objectiveValue_;
     1161        checkComplementarity (info,rowArray_[2],rowArray_[3]);
     1162        printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
     1163               objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost());
     1164        if (objectiveValue_>oldObjectiveValue) {
     1165          // make less likely this one will come in again
     1166          double multiplier = CoinDrand48();
     1167          int iSequence;
     1168          if (!cleanupIteration) {
     1169            iSequence = sequenceIn_;
     1170          } else {
     1171            iSequence = info->crucialSj();
     1172            if (iSequence>numberRows_)
     1173              iSequence -= numberRows_;
     1174            else
     1175              iSequence = numberColumns_ +(iSequence-numberXColumns);
    11771176          }
    1178           int iRow;
    1179           for (iRow=0;iRow<numberXRows;iRow++) {
    1180             int iColumn = iRow + numberColumns_;
    1181             double diff =cost_[iColumn];
    1182             double valueI = solution_[iColumn];
    1183             if (diff>0.0) {
    1184               double above = valueI-lower_[iColumn];
    1185               assert(above>=0.0);
    1186               infeasCost += diff*above;
    1187             } else if (diff<0.0) {
    1188               double below = upper_[iColumn]-valueI;
    1189               assert(below>=0.0);
    1190               infeasCost -= diff*below;
    1191             }
    1192           }
    1193           objectiveValue_ += linearCost;
    1194           assert (infeasCost>=0.0);
    1195           printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
    1196                  objectiveValue_,infeasCost,objectiveValue_+infeasCost);
    1197           if (objectiveValue_>oldValue) {
    1198             // make less likely this one will come in again
    1199             double multiplier = CoinDrand48();
    1200             int iSequence;
    1201             if (!cleanupIteration) {
    1202               iSequence = sequenceIn_;
    1203             } else {
    1204               int iSequence = info->crucialSj();
    1205               if (iSequence>numberRows_)
    1206                 iSequence -= numberRows_;
    1207               else
    1208                 iSequence = numberColumns_ +(iSequence-numberXColumns);
    1209             }
    1210             info->djWeight()[iSequence] /= (5.0+multiplier);
    1211           }
     1177          info->djWeight()[iSequence] /= (5.0+multiplier);
    12121178        }
    12131179        if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
     
    12171183            info->setCrucialSj(-1);
    12181184        }
    1219         checkComplementarity (info,rowArray_[2],rowArray_[3]);
    12201185
    12211186        if (whatNext==1) {
     
    12331198        // may not be correct on second time
    12341199        if (result&&(returnCode==-1||returnCode==-2||returnCode==-3)) {
    1235           assert(sequenceOut_<numberXColumns||
    1236                  sequenceOut_>=numberColumns_);
    12371200          int crucialSj=info->crucialSj();
    12381201          if (sequenceOut_<numberXColumns) {
    12391202            chosen =sequenceOut_ + numberRows_; // sj variable
    1240           } else {
     1203          } else if (sequenceOut_>=numberColumns_) {
    12411204            // Does this mean we can change pi
    12421205            int iRow = sequenceOut_-numberColumns_;
     
    12501213              abort();
    12511214            }
     1215          } else if (sequenceOut_<numberRows_) {
     1216            // pi out
     1217            chosen = sequenceOut_-numberXColumns+numberColumns_;
     1218          } else {
     1219            // sj out
     1220            chosen = sequenceOut_-numberRows_;
    12521221          }
    12531222          // But double check as original incoming might have gone out
     
    13811350  // Work out coefficients for quadratic term
    13821351  // This is expanded one
    1383   const CoinPackedMatrix * quadratic = quadraticObjective();
     1352  const CoinPackedMatrix * quadratic = info->quadraticObjective();
    13841353  const int * columnQuadratic = quadratic->getIndices();
    13851354  const int * columnQuadraticStart = quadratic->getVectorStarts();
     
    14201389      //printf("col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
    14211390    } else {
    1422       coeffSlack += alpha*cost_[iPivot];
    1423       //printf("new col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
    1424       if (iPivot==crucialSj) {
    1425         tj = alpha;
    1426         iSjRow = iRow;
     1391      if (iPivot>=numberColumns_) {
     1392        // ? do we go through column of pi
     1393        assert (iPivot!=crucialSj);
     1394        coeffSlack += alpha*cost_[iPivot];
     1395      } else if (iPivot<numberRows_) {
     1396        // ? do we go through column of pi
     1397        if (iPivot==crucialSj) {
     1398          tj = alpha;
     1399          iSjRow = iRow;
     1400        }
     1401      } else {
     1402        if (iPivot==crucialSj) {
     1403          tj = alpha;
     1404          iSjRow = iRow;
     1405        }
    14271406      }
    14281407    }
     
    14411420  }
    14421421  rhsArray->setNumElements(number2);
    1443   if (numberIterations_!=-701) {  //if (numberIterations_>=0||cleanupIteration) {
    1444     for (iIndex=0;iIndex<number2;iIndex++) {
    1445       int iColumn=index2[iIndex];
    1446       //double valueI = solution_[iColumn];
    1447       double alphaI = rhs[iColumn];
    1448       int j;
    1449       for (j=columnQuadraticStart[iColumn];
    1450            j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    1451         int jColumn = columnQuadratic[j];
    1452         double valueJ = solution_[jColumn];
    1453         double alphaJ = rhs[jColumn];
    1454         double elementValue = quadraticElement[j];
    1455         coeff1 += (valueJ*alphaI)*elementValue;
    1456         coeff2 += (alphaI*alphaJ)*elementValue;
    1457       }
    1458     }
    1459   } else {
    1460     const int * row = matrix_->getIndices();
    1461     const int * columnStart = matrix_->getVectorStarts();
    1462     const int * columnLength = matrix_->getVectorLengths();
    1463     const double * element = matrix_->getElements();
     1422  for (iIndex=0;iIndex<number2;iIndex++) {
     1423    int iColumn=index2[iIndex];
     1424    //double valueI = solution_[iColumn];
     1425    double alphaI = rhs[iColumn];
    14641426    int j;
    1465     int jRow = sequenceIn_+info->numberXRows();
    1466     printf("sequence in %d, cost %g\n",sequenceIn_,
    1467            upper_[jRow+numberColumns_]);
    1468     double xx=0.0;
    1469     for (j=0;j<numberColumns_;j++) {
    1470       int k;
    1471       for (k=columnStart[j];k<columnStart[j]+columnLength[j];k++) {
    1472         if (row[k]==jRow) {
    1473           printf ("col %d, el %g, sol %g, contr %g\n",
    1474                   j,element[k],solution_[j],element[k]*solution_[j]);
    1475           xx+= element[k]*solution_[j];
    1476         }
    1477       }
    1478     }
    1479     printf("sum %g\n",xx);
    1480     for (iIndex=0;iIndex<number2;iIndex++) {
    1481       int iColumn=index2[iIndex];
    1482       double valueI = solution_[iColumn];
    1483       double alphaI = rhs[iColumn];
    1484       //rhs[iColumn]=0.0;
    1485       printf("Column %d, alpha %g sol %g cost %g, contr %g\n",
    1486              iColumn,alphaI,valueI,cost_[iColumn],
    1487              alphaI*cost_[iColumn]);
    1488       int j;
    1489       for (j=columnQuadraticStart[iColumn];
    1490            j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
    1491         int jColumn = columnQuadratic[j];
    1492         double valueJ = solution_[jColumn];
    1493         double alphaJ = rhs[jColumn];
    1494         double elementValue = quadraticElement[j];
    1495         printf("j %d alphaJ %g solJ %g el %g, contr %g\n",
    1496                jColumn,alphaJ,valueJ,elementValue,
    1497                (valueJ*alphaI)*elementValue);
    1498         coeff1 += (valueJ*alphaI)*elementValue;
    1499         coeff2 += (alphaI*alphaJ)*elementValue;
    1500       }
     1427    for (j=columnQuadraticStart[iColumn];
     1428         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1429      int jColumn = columnQuadratic[j];
     1430      double valueJ = solution_[jColumn];
     1431      double alphaJ = rhs[jColumn];
     1432      double elementValue = quadraticElement[j];
     1433      coeff1 += (valueJ*alphaI)*elementValue;
     1434      coeff2 += (alphaI*alphaJ)*elementValue;
    15011435    }
    15021436  }
     
    15081442  int accuracyFlag=0;
    15091443  if (!cleanupIteration) {
    1510     //assert (fabs(way*coeff1-dualIn_)<1.0e-4*(1.0+fabs(dualIn_)));
     1444    if (fabs(dualIn_)<dualTolerance_&&fabs(coeff1)<dualTolerance_) {
     1445      dualIn_=0.0;
     1446      coeff1=0.0;
     1447    }
     1448    //assert (fabs(way*coeff1-dualIn_)<1.0e-1*(1.0+fabs(dualIn_)));
    15111449    //assert (way*coeff1*dualIn_>=0.0);
    15121450    if (way*coeff1*dualIn_<0.0) {
    15131451      // bad
    15141452      accuracyFlag=2;
    1515     } else if (fabs(way*coeff1-dualIn_)>1.0e-4*(1.0+fabs(dualIn_))) {
     1453    } else if (fabs(way*coeff1-dualIn_)>1.0e-3*(1.0+fabs(dualIn_))) {
    15161454      // not wondeful
    15171455      accuracyFlag=1;
     
    15541492    <<coeff1<<coeff2<<coeffSlack<<dualIn_<<d1<<d2
    15551493    <<CoinMessageEol;
    1556   if (!cleanupIteration)
    1557     dualIn_ = way*coeff1;
     1494  //if (!cleanupIteration)
     1495  //dualIn_ = way*coeff1;
    15581496  double mind1d2=1.0e50;
    15591497  if (cleanupIteration) {
     
    15701508  }
    15711509  maximumMovement = min(maximumMovement,mind1d2);
     1510  double trueMaximumMovement;
     1511  if (way>0.0)
     1512    trueMaximumMovement = min(1.0e30,upperIn_-valueIn_);
     1513  else
     1514    trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_);
    15721515  rhsArray->clear();
    15731516  double tentativeTheta = maximumMovement;
    15741517  double upperTheta = maximumMovement;
    15751518  bool throwAway=false;
    1576   if (numberIterations_>=2007) {
     1519  if (numberIterations_==126) {
    15771520    printf("Bad iteration coming up after iteration %d\n",numberIterations_);
    15781521  }
    15791522
     1523  double minimumAny=1.0e50;
     1524  int whichMinimum=-1;
     1525  double minimumAlpha=0.0;
    15801526  for (iIndex=0;iIndex<number;iIndex++) {
    15811527
     
    15951541      oldValue = bound-oldValue;
    15961542    }
     1543    if (iPivot>=numberXColumns&&iPivot<numberColumns_) {
     1544      double bound=0.0;
     1545      double oldValue2 = solution(iPivot);
     1546      if (alpha>0.0) {
     1547        // basic variable going towards lower bound
     1548        oldValue2 -= bound;
     1549      } else if (alpha<0.0) {
     1550        // basic variable going towards upper bound
     1551        oldValue2 = bound-oldValue;
     1552      }
     1553      double value = oldValue2-minimumAny*fabs(alpha);
     1554      if (value*oldValue2<0.0) {
     1555        double ratio = fabs(oldValue2/alpha);
     1556        if (ratio<minimumAny&&fabs(alpha)>1.0e2*acceptablePivot) {
     1557          minimumAny = ratio;
     1558          whichMinimum = iRow;
     1559          minimumAlpha= fabs(alpha);
     1560        }
     1561      }
     1562    }
    15971563    double value = oldValue-tentativeTheta*fabs(alpha);
    15981564    assert (oldValue>=-primalTolerance_*1.0001);
     
    16161582
    16171583  bool goBackOne = false;
     1584  double fake=1.0e-2;
    16181585
    16191586  if (numberRemaining) {
     
    16691636      //dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
    16701637      //dualCheck=0.0;
    1671       if (totalThru+thruThis>=dualCheck) {
     1638      if ((totalThru+thruThis>=dualCheck&&bestPivot>acceptablePivot)||fake*bestPivot>1.0e-3) {
    16721639        // We should be pivoting in this batch
    16731640        // so compress down to this lot
     
    17481715          } else {
    17491716            dualCheck = - 2.0*coeff2*theta_ - coeff1;
    1750             //dualCheck =0.0;
    1751             if (totalThru>=dualCheck&&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
     1717            if ((totalThru>=dualCheck||fake*bestPivot>1.0e-3)
     1718                &&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
    17521719              if (!cleanupIteration) {
    17531720                // We can pivot out sj
    17541721                if (upperTheta==maximumMovement) {
    17551722                  printf("figures\n");
    1756                 } else {
    1757                   if (lastThru>=dualCheck)
     1723                } else if (iSjRow2>=0) {
     1724                  if (lastThru>=dualCheck&&theta_<trueMaximumMovement) {
    17581725                    printf("totalThru %g, lastThru %g - taking sj to zero?\n",totalThru,lastThru);
     1726                    // make sure will take correct path
     1727                    mind1d2=1.0;
     1728                    maximumMovement=1.0;
     1729                    d2=0.0;
     1730                    pivotRow_=-1;
     1731                  }
    17591732                }
    17601733              } else {
     
    17631736              }
    17641737              break; // no point trying another loop
    1765             } else if (totalThru>=dualCheck||upperTheta==maximumMovement) {
     1738            } else if (totalThru>=dualCheck||fake*bestPivot>1.0e-3||upperTheta==maximumMovement) {
    17661739              break; // no point trying another loop
    17671740            } else {
     
    18071780    }
    18081781  }
     1782  if (0&&minimumAny<theta_&&cleanupIteration&&bestEverPivot<1.0e-2
     1783      &&minimumAlpha>bestEverPivot*1.0e2&&minimumAny>1.0e-1&&info->currentPhase()==0) {
     1784    printf("Possible other pivot %d %g %g - alpha %g\n",
     1785           whichMinimum,minimumAny,theta_,minimumAlpha);
     1786    pivotRow_ = whichMinimum;
     1787    alpha_ = work[pivotRow_];
     1788    // translate to sequence
     1789    sequenceOut_ = pivotVariable_[pivotRow_];
     1790    valueOut_ = solution(sequenceOut_);
     1791    lowerOut_=0.0;
     1792    upperOut_=0.0;
     1793    theta_= fabs(valueOut_/alpha_);
     1794
     1795    if (way<0.0)
     1796      theta_ = - theta_;
     1797    if (alpha_*way<0.0) {
     1798      directionOut_=-1;      // to upper bound
     1799    } else {
     1800      directionOut_=1;      // to lower bound
     1801    }
     1802    dualOut_ = reducedCost(sequenceOut_);
     1803  } else {
    18091804
    18101805  if (pivotRow_>=0) {
    1811    
     1806    if (0) {
     1807      double move = coeff1*theta_+coeff2*theta_*theta_;
     1808      printf("Predicted change in obj is %g %g %g\n",
     1809             move,objectiveValue_-move,objectiveValue_+move);
     1810    }
    18121811#define MINIMUMTHETA 1.0e-12
    18131812    // will we need to increase tolerance
     
    18651864      accuracyFlag=2;
    18661865    }
    1867     double trueMaximumMovement;
    1868     if (way>0.0)
    1869       trueMaximumMovement = min(1.0e30,upperIn_-valueIn_);
    1870     else
    1871       trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_);
    18721866    if (maximumMovement<1.0e20&&maximumMovement==trueMaximumMovement) {
    18731867      // flip
     
    18981892      } else {
    18991893        // incoming dj going to zero
    1900         assert(!cleanupIteration);
    1901         if (!cleanupIteration)
     1894        if (!cleanupIteration) {
    19021895          result=0;
    1903         iSjRow=iSjRow2;
    1904         crucialSj = pivotVariable_[iSjRow];
    1905       }
    1906       assert (pivotRow_<0);
    1907       setColumnStatus(crucialSj,isFree);
    1908       pivotRow_ = iSjRow;
    1909       alpha_ = work[pivotRow_];
    1910       // translate to sequence
    1911       sequenceOut_ = pivotVariable_[pivotRow_];
    1912       assert (sequenceOut_==crucialSj);
    1913       valueOut_ = solution(sequenceOut_);
    1914       theta_ = maximumMovement;
    1915       if (way<0.0)
    1916         theta_ = - theta_;
    1917       lowerOut_=0.0;
    1918       upperOut_=0.0;
    1919       //????
    1920       dualOut_ = reducedCost(sequenceOut_);
     1896          iSjRow=iSjRow2;
     1897          crucialSj = pivotVariable_[iSjRow];
     1898          assert (pivotRow_<0);
     1899        } else {
     1900          assert (fabs(dualIn_)<1.0e-3);
     1901          result=1;
     1902          if (minimumAlpha>0.0) {
     1903            // could do this in other places if pivot looks good
     1904            printf("Should take minimum\n");
     1905            pivotRow_ = whichMinimum;
     1906            alpha_ = work[pivotRow_];
     1907            // translate to sequence
     1908            sequenceOut_ = pivotVariable_[pivotRow_];
     1909            valueOut_ = solution(sequenceOut_);
     1910            theta_ = minimumAny;
     1911            if (way<0.0)
     1912              theta_ = - theta_;
     1913            lowerOut_=0.0;
     1914            upperOut_=0.0;
     1915            //????
     1916            dualOut_ = reducedCost(sequenceOut_);
     1917          }
     1918        }
     1919      }
     1920      if (!result) {
     1921        setColumnStatus(crucialSj,isFree);
     1922        pivotRow_ = iSjRow;
     1923        alpha_ = work[pivotRow_];
     1924        // translate to sequence
     1925        sequenceOut_ = pivotVariable_[pivotRow_];
     1926        assert (sequenceOut_==crucialSj);
     1927        valueOut_ = solution(sequenceOut_);
     1928        theta_ = fabs(valueOut_/alpha_);
     1929        assert (fabs(maximumMovement-theta_)<1.0e-3*1.0+maximumMovement);
     1930        if (way<0.0)
     1931          theta_ = - theta_;
     1932        lowerOut_=0.0;
     1933        upperOut_=0.0;
     1934        //????
     1935        dualOut_ = reducedCost(sequenceOut_);
     1936      }
    19211937    } else {
    19221938      // need to do something
     
    19241940    }
    19251941  }
     1942  }
     1943
    19261944
    19271945  // clear arrays
     
    19401958    objectiveValue_ =0.0;
    19411959    CoinPackedMatrix * quadratic =
    1942       info->originalModel()->quadraticObjective();
     1960      info->originalObjective()->quadraticObjective();
    19431961    const int * columnQuadratic = quadratic->getIndices();
    19441962    const int * columnQuadraticStart = quadratic->getVectorStarts();
     
    19461964    const double * quadraticElement = quadratic->getElements();
    19471965    int numberColumns = info->numberXColumns();
    1948     const double * objective = info->originalObjective();
     1966    const double * objective = info->linearObjective();
    19491967    for (iColumn=0;iColumn<numberColumns;iColumn++)
    19501968      objectiveValue_ += objective[iColumn]*solution_[iColumn];
     
    20732091  }
    20742092  // Check if looping
    2075   // Take out for now
    20762093  int loop;
    2077   if (type!=2&&0)
    2078     loop = progress->looping();
     2094  if (type!=2)
     2095    loop = progress->looping(); // saves iteration number
    20792096  else
    20802097    loop=-1;
     2098  //loop=-1;
    20812099  if (loop>=0) {
    20822100    problemStatus_ = loop; //exit if in loop
     
    21072125  // give code benefit of doubt
    21082126  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
    2109       sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
     2127      sumOfRelaxedPrimalInfeasibilities_ == 0.0&&info->sequenceIn()<0) {
    21102128    // say optimal (with these bounds etc)
    21112129    numberDualInfeasibilities_ = 0;
     
    23322350{
    23332351  const int * which = info->quadraticSequence();
    2334   // Matrix for linear stuff
    2335   const int * row = matrix_->getIndices();
    2336   const int * columnStart = matrix_->getVectorStarts();
    2337   const double * element =  matrix_->getElements();
    2338   const int * columnLength = matrix_->getVectorLengths();
    23392352  int numberXColumns = info->numberXColumns();
    23402353  int numberXRows = info->numberXRows();
     
    23432356  int start=numberXColumns+numberXRows;
    23442357  int jSequence;
     2358  const double * djWeight = info->djWeight();
    23452359  // Actually only need to do this after invert (use extra parameter to createDjs)
    23462360  // or when infeasibilities change
     
    24582472    }
    24592473  }
     2474  // fill in linear ones
     2475  memcpy(dj_,cost_,numberXColumns*sizeof(double));
     2476  if (!rowScale_) {
     2477    matrix_->transposeTimes(-1.0,pi,dj_);
     2478  } else {
     2479    matrix_->transposeTimes(-1.0,pi,dj_,rowScale_,columnScale_);
     2480  }
     2481  memset(dj_+numberXColumns,0,(numberXRows+info->numberQuadraticColumns())*sizeof(double));
    24602482  for (iSequence=0;iSequence<numberXColumns;iSequence++) {
    24612483    int jSequence = which[iSequence];
     
    24652487      value = solution_[jSequence];
    24662488    } else {
    2467       value=0.0;
    2468       //value=cost_[iSequence];
    2469       int j;
    2470       for (j=columnStart[iSequence];j<columnStart[iSequence]+columnLength[iSequence]; j++) {
    2471         int iRow = row[j];
    2472         value -= element[j]*pi[iRow];
    2473       }
    2474     }
     2489      value=dj_[iSequence];
     2490    }
     2491    double value2 = value*djWeight[iSequence];
     2492    if (fabs(value2)>dualTolerance_)
     2493      value=value2;
     2494    else if (value<-dualTolerance_)
     2495      value = -1.001*dualTolerance_;
     2496    else if (value>dualTolerance_)
     2497      value = 1.001*dualTolerance_;
     2498    if (djWeight[iSequence]<1.0e-6)
     2499      value=value2;
    24752500    dj_[iSequence]=value;
    24762501  }
     
    24802505    int iSequence  = jSequence + numberColumns_;
    24812506    int iPi = jSequence+numberXColumns;
    2482     double value = solution_[iPi];
    2483     dj_[iSequence]=value+cost_[iSequence];
     2507    double value = solution_[iPi]+cost_[iSequence];
     2508    double value2 = value*djWeight[iSequence];
     2509    if (fabs(value2)>dualTolerance_)
     2510      value=value2;
     2511    else if (value<-dualTolerance_)
     2512      value = -1.001*dualTolerance_;
     2513    else if (value>dualTolerance_)
     2514      value = 1.001*dualTolerance_;
     2515    if (djWeight[iSequence]<1.0e-6)
     2516      value=value2;
     2517    dj_[iSequence]=value;
    24842518  }
    24852519}
    24862520
    24872521int
    2488 ClpSimplexPrimalQuadratic::checkComplementarity (const  ClpQuadraticInfo * info,
     2522ClpSimplexPrimalQuadratic::checkComplementarity (ClpQuadraticInfo * info,
    24892523                            CoinIndexedVector * array1, CoinIndexedVector * array2)
    24902524{
     
    24982532  // Compute objective function from scratch
    24992533  const CoinPackedMatrix * quadratic =
    2500     info->originalModel()->quadraticObjective();
     2534    info->quadraticObjective();
    25012535  const int * columnQuadratic = quadratic->getIndices();
    25022536  const int * columnQuadraticStart = quadratic->getVectorStarts();
    25032537  const int * columnQuadraticLength = quadratic->getVectorLengths();
    25042538  const double * quadraticElement = quadratic->getElements();
    2505   const double * originalCost = info->originalObjective();
     2539  const double * originalCost = info->linearObjective();
    25062540  int iColumn;
    25072541  objectiveValue_=0.0;
    25082542  double infeasCost=0.0;
    25092543  double linearCost=0.0;
     2544  int number0=0,number1=0,number2=0;
    25102545  for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    25112546    double valueI = solution_[iColumn];
     
    25312566    }
    25322567    int jSequence = which[iColumn];
     2568    if (jSequence>=0)
     2569      jSequence += start;
    25332570    double value=dj_[iColumn];
    25342571    ClpSimplex::Status status = getColumnStatus(iColumn);
     2572    if (status!=basic&&jSequence>=0) {
     2573      if (getColumnStatus(jSequence)==basic)
     2574        number1++;
     2575      else
     2576        number0++;
     2577    }
    25352578   
    25362579    switch(status) {
     
    25382581    case ClpSimplex::basic:
    25392582      if (jSequence>=0) {
    2540         jSequence += start;
    2541         if (getColumnStatus(jSequence)==basic)
     2583        if (getColumnStatus(jSequence)==basic) {
    25422584          handler_->message(CLP_QUADRATIC_BOTH,messages_)
    25432585            <<"Structural"<<iColumn
    25442586            <<solution_[iColumn]<<jSequence<<solution_[jSequence]
    25452587            <<CoinMessageEol;
     2588          number2++;
     2589        } else {
     2590          number1++;
     2591        }
    25462592      }
    25472593      break;
     
    25862632    double value=dj_[iRow+numberColumns_];
    25872633    ClpSimplex::Status status = getRowStatus(iRow);
     2634    if (status!=basic) {
     2635      if (getColumnStatus(jSequence)==basic)
     2636        number1++;
     2637      else
     2638        number0++;
     2639    }
    25882640   
    25892641    switch(status) {
    25902642     
    25912643    case ClpSimplex::basic:
    2592       if (getColumnStatus(jSequence)==basic)
     2644      if (getColumnStatus(jSequence)==basic) {
    25932645        printf("Row %d (%g) and %d (%g) both basic\n",
    25942646               iRow,solution_[iColumn],jSequence,solution_[jSequence]);
     2647          number2++;
     2648        } else {
     2649          number1++;
     2650        }
    25952651      break;
    25962652    case ClpSimplex::isFixed:
     
    26162672    }
    26172673  }
     2674  //printf("number 0 - %d, 1 - %d, 2 - %d\n",number0,number1,number2);
    26182675  objectiveValue_ += linearCost+infeasCost;
    26192676  assert (infeasCost>=0.0);
    26202677  sumOfRelaxedDualInfeasibilities_ =sumDualInfeasibilities_;
     2678  // But not zero if cleanup iteration
     2679  if (info->sequenceIn()>=0&&!numberDualInfeasibilities_)
     2680    numberDualInfeasibilities_=1;
    26212681  numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_;
     2682  info->setInfeasCost(infeasCost);
    26222683  return numberDualInfeasibilities_;
    26232684}
     
    26312692
    26322693  // Get list of non linear columns
    2633   CoinPackedMatrix * quadratic = quadraticObjective();
     2694  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     2695  assert (quadraticObj);
     2696  info.setOriginalObjective(quadraticObj);
     2697  CoinPackedMatrix * quadratic = info.quadraticObjective();
    26342698  if (!quadratic||!quadratic->getNumElements()) {
    26352699    // no quadratic part
     
    26402704  double * columnLower = this->columnLower();
    26412705  double * columnUpper = this->columnUpper();
    2642   double * objective = this->objective();
     2706  double * objective = info.linearObjective();
    26432707  double * solution = this->primalColumnSolution();
    26442708  double * dj = this->dualColumnSolution();
     
    29823046  for (iColumn=0;iColumn<numberColumns_;iColumn++)
    29833047    objValue += objective[iColumn]*solution2[iColumn];
    2984   CoinPackedMatrix * quadratic = quadraticObjective();
     3048  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
     3049  assert (quadraticObj);
     3050  CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    29853051  double * dj = dualColumnSolution();
    29863052  // Matrix for linear stuff
     
    30533119/// Default constructor.
    30543120ClpQuadraticInfo::ClpQuadraticInfo()
    3055   : originalModel_(NULL),
     3121  : originalObjective_(NULL),
    30563122    quadraticSequence_(NULL),
    30573123    backSequence_(NULL),
     
    30673133    numberXRows_(-1),
    30683134    numberXColumns_(-1),
    3069     numberQuadraticColumns_(0)
     3135    numberQuadraticColumns_(0),
     3136    infeasCost_(0.0)
    30703137{
    30713138}
    30723139// Constructor from original model
    30733140ClpQuadraticInfo::ClpQuadraticInfo(const ClpSimplex * model)
    3074   : originalModel_(model),
     3141  : originalObjective_(NULL),
    30753142    quadraticSequence_(NULL),
    30763143    backSequence_(NULL),
     
    30863153    numberXRows_(-1),
    30873154    numberXColumns_(-1),
    3088     numberQuadraticColumns_(0)
     3155    numberQuadraticColumns_(0),
     3156    infeasCost_(0.0)
    30893157{
    3090   if (originalModel_) {
    3091     numberXRows_ = originalModel_->numberRows();
    3092     numberXColumns_ = originalModel_->numberColumns();
     3158  if (model) {
     3159    numberXRows_ = model->numberRows();
     3160    numberXColumns_ = model->numberColumns();
     3161    originalObjective_ = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
     3162    assert (originalObjective_);
    30933163    quadraticSequence_ = new int[numberXColumns_];
    30943164    backSequence_ = new int[numberXColumns_];
     
    31183188// Copy
    31193189ClpQuadraticInfo::ClpQuadraticInfo(const ClpQuadraticInfo& rhs)
    3120   : originalModel_(rhs.originalModel_),
     3190  : originalObjective_(rhs.originalObjective_),
    31213191    quadraticSequence_(NULL),
    31223192    backSequence_(NULL),
     
    31293199    numberXRows_(rhs.numberXRows_),
    31303200    numberXColumns_(rhs.numberXColumns_),
    3131     numberQuadraticColumns_(rhs.numberQuadraticColumns_)
     3201    numberQuadraticColumns_(rhs.numberQuadraticColumns_),
     3202    infeasCost_(rhs.infeasCost_)
    31323203{
    31333204  if (numberXColumns_) {
     
    31663237{
    31673238  if (this != &rhs) {
    3168     originalModel_ = rhs.originalModel_;
     3239    originalObjective_ = rhs.originalObjective_;
    31693240    delete [] quadraticSequence_;
    31703241    delete [] backSequence_;
     
    31803251    numberXRows_ = rhs.numberXRows_;
    31813252    numberXColumns_ = rhs.numberXColumns_;
     3253    infeasCost_=rhs.infeasCost_;
    31823254    numberQuadraticColumns_=rhs.numberQuadraticColumns_;
    31833255    if (numberXColumns_) {
     
    32583330  }
    32593331}
     3332// Quadratic objective
     3333CoinPackedMatrix *
     3334ClpQuadraticInfo::quadraticObjective() const     
     3335{
     3336  return originalObjective_->quadraticObjective();
     3337}
     3338// Linear objective
     3339double *
     3340ClpQuadraticInfo::linearObjective() const     
     3341{
     3342  return originalObjective_->linearObjective();
     3343}
  • branches/pre/Makefile.Clp

    r183 r196  
    66# between then specify the exact level you want, e.g., -O1 or -O2
    77OptLevel := -g
    8 #OptLevel := -O3
     8OptLevel := -O3
    99
    1010
     
    2828LIBSRC += ClpPrimalColumnPivot.cpp
    2929LIBSRC += ClpPrimalColumnSteepest.cpp
     30LIBSRC += ClpQuadraticObjective.cpp
    3031LIBSRC += ClpSimplex.cpp
    3132LIBSRC += ClpSimplexDual.cpp
    3233LIBSRC += ClpSimplexPrimal.cpp
    3334LIBSRC += ClpSimplexPrimalQuadratic.cpp
    34 LIBSRC += ClpPrimalQuadraticDantzig.cpp
     35#LIBSRC += ClpPrimalQuadraticDantzig.cpp
    3536# and Presolve stuff
    3637LIBSRC += ClpPresolve.cpp
  • branches/pre/Test/ClpMain.cpp

    r180 r196  
    17031703                double * primalRowSolution =
    17041704                  models[iModel].primalRowSolution();
     1705                double * rowLower = models[iModel].rowLower();
     1706                double * rowUpper = models[iModel].rowUpper();
     1707                double primalTolerance = models[iModel].primalTolerance();
    17051708                char format[6];
    17061709                sprintf(format,"%%-%ds",max(lengthName,8));
    17071710                for (iRow=0;iRow<numberRows;iRow++) {
     1711                  if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance||
     1712                      primalRowSolution[iRow]<rowLower[iRow]-primalTolerance)
     1713                    fprintf(fp,"** ");
    17081714                  fprintf(fp,"%7d ",iRow);
    17091715                  if (lengthName)
     
    17181724                double * primalColumnSolution =
    17191725                  models[iModel].primalColumnSolution();
     1726                double * columnLower = models[iModel].columnLower();
     1727                double * columnUpper = models[iModel].columnUpper();
    17201728                for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1729                  if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance||
     1730                      primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance)
     1731                    fprintf(fp,"** ");
    17211732                  fprintf(fp,"%7d ",iColumn);
    17221733                  if (lengthName)
  • branches/pre/Test/unitTest.cpp

    r192 r196  
    936936  }
    937937  // test network
    938 #define QUADRATIC
     938  //#define QUADRATIC
    939939#ifndef QUADRATIC
    940940  if (1) {   
  • branches/pre/include/ClpModel.hpp

    r183 r196  
    308308   /// Clp Matrix
    309309   inline ClpMatrixBase * clpMatrix() const     { return matrix_; }
    310    /// Quadratic objective
    311    inline CoinPackedMatrix * quadraticObjective() const     { return quadraticObjective_; }
    312310   /// Objective value
    313311   inline double objectiveValue() const {
     
    488486  /// Row copy if wanted
    489487  ClpMatrixBase * rowCopy_;
    490   /// Quadratic objective if any
    491   CoinPackedMatrix * quadraticObjective_;
    492488  /// Infeasible/unbounded ray
    493489  double * ray_;
  • branches/pre/include/ClpSimplex.hpp

    r180 r196  
    195195      Return code is 0 if nothing interesting, -1 if infeasible both
    196196      ways and +1 if infeasible one way (check values to see which one(s))
     197      Solutions are filled in as well - even down, odd up - also
     198      status and number of iterations
    197199  */
    198200  int strongBranching(int numberVariables,const int * variables,
    199201                      double * newLower, double * newUpper,
     202                      double ** outputSolution,
     203                      int * outputStatus, int * outputIterations,
    200204                      bool stopOnFirstInfeasible=true,
    201205                      bool alwaysFinish=false);
  • branches/pre/include/ClpSimplexDual.hpp

    r180 r196  
    120120      Return code is 0 if nothing interesting, -1 if infeasible both
    121121      ways and +1 if infeasible one way (check values to see which one(s))
     122      Solutions are filled in as well - even down, odd up - also
     123      status and number of iterations
    122124  */
    123125  int strongBranching(int numberVariables,const int * variables,
    124126                      double * newLower, double * newUpper,
     127                      double ** outputSolution,
     128                      int * outputStatus, int * outputIterations,
    125129                      bool stopOnFirstInfeasible=true,
    126130                      bool alwaysFinish=false);
  • branches/pre/include/ClpSimplexPrimalQuadratic.hpp

    r191 r196  
    1212
    1313class ClpQuadraticInfo;
     14class ClpQuadraticObjective;
    1415
    1516#include "ClpSimplexPrimal.hpp"
     
    6061                   ClpQuadraticInfo & info);
    6162  /// Checks complementarity and computes infeasibilities
    62   int checkComplementarity (const ClpQuadraticInfo * info,
     63  int checkComplementarity (ClpQuadraticInfo * info,
    6364                            CoinIndexedVector * array1, CoinIndexedVector * array2);
    6465  /// Fills in reduced costs
     
    158159  {return currentSolution_;};
    159160  void setCurrentSolution(const double * solution);
    160   /// Original objective
    161   inline const double * originalObjective() const
    162   {return originalModel_->objective();};
    163161  /// Quadratic sequence or -1 if linear
    164162  inline const int * quadraticSequence() const
     
    167165  inline const int * backSequence() const
    168166  {return backSequence_;};
    169   /// Returns pointer to original model
    170   inline const ClpSimplex * originalModel() const
    171   { return originalModel_;};
     167  /// Returns pointer to original objective
     168  inline const ClpQuadraticObjective * originalObjective() const
     169  { return originalObjective_;};
     170  inline void setOriginalObjective( ClpQuadraticObjective * obj)
     171  { originalObjective_ = obj;};
     172  /// Quadratic objective
     173  CoinPackedMatrix * quadraticObjective() const;
     174  /// Linear objective
     175  double * linearObjective() const;
    172176  /// Save current In and Sj status
    173177  void saveStatus();
     
    177181  inline double * djWeight() const
    178182  { return djWeight_;};
     183  /// Infeas cost
     184  inline double infeasCost() const
     185  { return infeasCost_;};
     186  inline void setInfeasCost(double value)
     187  { infeasCost_ = value;};
    179188  //@}
    180189   
     
    182191  /**@name Data members */
    183192  //@{
    184   /// Model
    185   const ClpSimplex * originalModel_;
     193  /// Objective
     194  const ClpQuadraticObjective * originalObjective_;
    186195  /// Quadratic sequence
    187196  int * quadraticSequence_;
     
    212221  /// Number of quadratic columns
    213222  int numberQuadraticColumns_;
     223  /// Infeasibility cost
     224  double infeasCost_;
    214225  //@}
    215226};
Note: See TracChangeset for help on using the changeset viewer.