ClpSimplexPrimalQuadratic Class Reference

#include <ClpSimplexPrimalQuadratic.hpp>

Inheritance diagram for ClpSimplexPrimalQuadratic:

List of all members.

Public Member Functions
Description of algorithm
int	primalSLP (int numberPasses, double deltaTolerance)
	A sequential LP method.
int	primalQuadratic (int phase=2)
int	primalQuadratic2 (ClpQuadraticInfo *info, int phase=2)
ClpSimplexPrimalQuadratic *	makeQuadratic (ClpQuadraticInfo &info)
int	endQuadratic (ClpSimplexPrimalQuadratic *quadraticModel, ClpQuadraticInfo &info)
	This moves solution back.
int	checkComplementarity (ClpQuadraticInfo *info, CoinIndexedVector *array1, CoinIndexedVector *array2)
	Checks complementarity and computes infeasibilities.
void	createDjs (ClpQuadraticInfo *info, CoinIndexedVector *array1, CoinIndexedVector *array2)
	Fills in reduced costs.
int	whileIterating (ClpQuadraticInfo *info)
int	primalRow (CoinIndexedVector *rowArray, CoinIndexedVector *rhsArray, CoinIndexedVector *spareArray, CoinIndexedVector *spareArray2, ClpQuadraticInfo *info, int cleanupIteration)
void	statusOfProblemInPrimal (int &lastCleaned, int type, ClpSimplexProgress *progress, ClpQuadraticInfo *info)

---

Detailed Description

This solves Quadratic LPs using the primal simplex method

It inherits from ClpSimplexPrimal. It has no data of its own and is never created - only cast from a ClpSimplexPrimal object at algorithm time. If needed create new class and pass around

Definition at line 26 of file ClpSimplexPrimalQuadratic.hpp.

---

Member Function Documentation

ClpSimplexPrimalQuadratic * ClpSimplexPrimalQuadratic::makeQuadratic (  ClpQuadraticInfo &  info  )

This creates the large version of QP and fills in quadratic information. Returns NULL if no quadratic information Definition at line 2770 of file ClpSimplexPrimalQuadratic.cpp. References ClpSimplex::allSlackBasis(), ClpSimplex::ClpSimplex(), ClpModel::columnLower(), ClpModel::columnUpper(), ClpModel::dualColumnSolution(), ClpModel::dualRowSolution(), ClpSimplex::factorization(), ClpSimplex::loadProblem(), ClpModel::matrix(), ClpModel::messageHandler(), ClpQuadraticInfo::numberQuadraticColumns(), ClpQuadraticInfo::numberQuadraticRows(), ClpModel::numberRows(), ClpModel::objective(), ClpModel::passInMessageHandler(), ClpSimplexPrimal::primal(), ClpModel::primalColumnSolution(), ClpModel::primalRowSolution(), ClpQuadraticInfo::quadraticObjective(), ClpQuadraticObjective::quadraticObjective(), ClpModel::resize(), ClpModel::rowLower(), ClpModel::rowUpper(), ClpSimplex::setColumnStatus(), ClpModel::setLogLevel(), ClpQuadraticInfo::setOriginalObjective(), ClpSimplex::setRowStatus(), and ClpSimplex::times(). Referenced by primalQuadratic(). { // Get list of non linear columns ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_)); assert (quadraticObj); CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective(); if (!quadratic||!quadratic->getNumElements()) { // no quadratic part return NULL; } int numberColumns = this->numberColumns(); double * columnLower = this->columnLower(); double * columnUpper = this->columnUpper(); double * objective = this->objective(); double * solution = this->primalColumnSolution(); double * dj = this->dualColumnSolution(); double * pi = this->dualRowSolution(); int numberRows = this->numberRows(); double * rowLower = this->rowLower(); double * rowUpper = this->rowUpper(); // and elements CoinPackedMatrix * matrix = this->matrix(); const int * row = matrix->getIndices(); const int * columnStart = matrix->getVectorStarts(); const double * element = matrix->getElements(); const int * columnLength = matrix->getVectorLengths(); int iColumn; const int * columnQuadratic = quadratic->getIndices(); const int * columnQuadraticStart = quadratic->getVectorStarts(); const int * columnQuadraticLength = quadratic->getVectorLengths(); const double * quadraticElement = quadratic->getElements(); #if 0 // deliberate bad solution // Change to use phase //double * saveO = new double[numberColumns]; //memcpy(saveO,objective,numberColumns*sizeof(double)); //memset(objective,0,numberColumns*sizeof(double)); tempMessage messageHandler(this);; passInMessageHandler(&messageHandler); factorization()->maximumPivots(1); primal(); CoinMessageHandler handler2; passInMessageHandler(&handler2); factorization()->maximumPivots(100); setMaximumIterations(1000); #endif //memcpy(objective,saveO,numberColumns*sizeof(double)); // Get a feasible solution //printf("For testing - deliberate bad solution\n"); //columnUpper[0]=0.0; //columnLower[0]=0.0; //quadraticSLP(50,1.0e-4); //primal(1); //columnUpper[0]=COIN_DBL_MAX; // Create larger problem // First workout how many rows extra info=ClpQuadraticInfo(this); quadratic = info.quadraticObjective(); int numberQuadratic = info.numberQuadraticColumns(); int newNumberRows = numberRows+numberQuadratic; int newNumberColumns = numberColumns + numberRows + numberQuadratic; int numberElements = 2*matrix->getNumElements() +2*quadratic->getNumElements() + numberQuadratic; double * elements2 = new double[numberElements]; int * start2 = new int[newNumberColumns+1]; int * row2 = new int[numberElements]; double * objective2 = new double[newNumberColumns]; double * columnLower2 = new double[newNumberColumns]; double * columnUpper2 = new double[newNumberColumns]; double * rowLower2 = new double[newNumberRows]; double * rowUpper2 = new double[newNumberRows]; memcpy(rowLower2,rowLower,numberRows*sizeof(double)); memcpy(rowUpper2,rowUpper,numberRows*sizeof(double)); int iRow; for (iRow=0;iRow<numberQuadratic;iRow++) { double cost = objective[iRow]; rowLower2[iRow+numberRows]=cost; rowUpper2[iRow+numberRows]=cost; } memset(objective2,0,newNumberColumns*sizeof(double)); // Get a row copy of quadratic objective in standard format CoinPackedMatrix copyQ; copyQ.reverseOrderedCopyOf(*quadratic); const int * columnQ = copyQ.getIndices(); const CoinBigIndex * rowStartQ = copyQ.getVectorStarts(); const int * rowLengthQ = copyQ.getVectorLengths(); const double * elementByRowQ = copyQ.getElements(); // Move solution across double * solution2 = new double[newNumberColumns]; memset(solution2,0,newNumberColumns*sizeof(double)); newNumberColumns=0; numberElements=0; start2[0]=0; // x memcpy(dj,objective,numberColumns*sizeof(double)); for (iColumn=0;iColumn<numberColumns;iColumn++) { // Original matrix columnLower2[iColumn]=columnLower[iColumn]; columnUpper2[iColumn]=columnUpper[iColumn]; solution2[iColumn]=solution[iColumn]; // Put in cost so we know it objective2[iColumn]=objective[iColumn]; int j; for (j=columnStart[iColumn]; j<columnStart[iColumn]+columnLength[iColumn]; j++) { elements2[numberElements]=element[j]; row2[numberElements++]=row[j]; } // Quadratic and modify djs for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn]; j++) { int jColumn = columnQuadratic[j]; double value = quadraticElement[j]; if (iColumn!=jColumn) value *= 0.5; dj[jColumn] += solution[iColumn]*value; elements2[numberElements]=-value; row2[numberElements++]=jColumn+numberRows; } for (j=rowStartQ[iColumn]; j<rowStartQ[iColumn]+rowLengthQ[iColumn]; j++) { int jColumn = columnQ[j]; double value = elementByRowQ[j]; if (iColumn!=jColumn) { value *= 0.5; dj[jColumn] += solution[iColumn]*value; elements2[numberElements]=-value; row2[numberElements++]=jColumn+numberRows; } } start2[iColumn+1]=numberElements; } newNumberColumns=numberColumns; // pi int numberQuadraticRows = info.numberQuadraticRows(); // Get a row copy in standard format CoinPackedMatrix copy; copy.reverseOrderedCopyOf(*(this->matrix())); // get matrix data pointers const int * column = copy.getIndices(); const CoinBigIndex * rowStart = copy.getVectorStarts(); const int * rowLength = copy.getVectorLengths(); const double * elementByRow = copy.getElements(); for (iRow=0;iRow<numberRows;iRow++) { solution2[newNumberColumns]=pi[iRow]; double value = pi[iRow]; columnLower2[newNumberColumns]=-COIN_DBL_MAX; columnUpper2[newNumberColumns]=COIN_DBL_MAX; int j; for (j=rowStart[iRow]; j<rowStart[iRow]+rowLength[iRow]; j++) { double elementValue=elementByRow[j]; int jColumn = column[j]; if (jColumn>=0) { elements2[numberElements]=elementValue; row2[numberElements++]=jColumn+numberRows; } dj[jColumn]-= value*elementValue; } #ifndef CORRECT_ROW_COUNTS newNumberColumns++; start2[newNumberColumns]=numberElements; #else if (numberElements>start2[newNumberColumns]) { newNumberColumns++; start2[newNumberColumns]=numberElements; } #endif } // djs for (iColumn=0;iColumn<numberQuadratic;iColumn++) { columnLower2[newNumberColumns]=-COIN_DBL_MAX; columnUpper2[newNumberColumns]=COIN_DBL_MAX; solution2[newNumberColumns]=dj[iColumn]; elements2[numberElements]=1.0; row2[numberElements++]=iColumn+numberRows; newNumberColumns++; start2[newNumberColumns]=numberElements; } // Create model ClpSimplex * model2 = new ClpSimplex(*this); model2->resize(0,0); model2->loadProblem(newNumberColumns,newNumberRows, start2,row2, elements2, columnLower2,columnUpper2, objective2, rowLower2,rowUpper2); delete [] rowLower2; delete [] rowUpper2; delete [] columnLower2; delete [] columnUpper2; // Now create expanded quadratic objective for use in primalRow // Later pack down in some way start2[0]=0; numberElements=0; for (iColumn=0;iColumn<numberColumns;iColumn++) { // Quadratic int j; for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn]; j++) { int jColumn = columnQuadratic[j]; double value = quadraticElement[j]; if (iColumn!=jColumn) value *= 0.5; elements2[numberElements]=value; row2[numberElements++]=jColumn; } for (j=rowStartQ[iColumn]; j<rowStartQ[iColumn]+rowLengthQ[iColumn]; j++) { int jColumn = columnQ[j]; double value = elementByRowQ[j]; if (iColumn!=jColumn) { value *= 0.5; elements2[numberElements]=value; row2[numberElements++]=jColumn; } } start2[iColumn+1]=numberElements; } // and pad //for (;iColumn<newNumberColumns;iColumn++) //start2[iColumn+1]=numberElements; // Load up objective with expanded linear ClpQuadraticObjective * obj = new ClpQuadraticObjective(objective2,numberColumns, start2,row2,elements2,newNumberColumns); delete [] objective2; info.setOriginalObjective(obj); //model2->loadQuadraticObjective(newNumberColumns,start2,row2,elements2); delete [] start2; delete [] row2; delete [] elements2; model2->allSlackBasis(); //model2->scaling(false); //printf("scaling off\n"); model2->setLogLevel(this->logLevel()); // Move solution across memcpy(model2->primalColumnSolution(),solution2, newNumberColumns*sizeof(double)); columnLower2 = model2->columnLower(); columnUpper2 = model2->columnUpper(); delete [] solution2; solution2 = model2->primalColumnSolution(); // Compute row activities and check feasible double * rowSolution2 = model2->primalRowSolution(); memset(rowSolution2,0,newNumberRows*sizeof(double)); model2->times(1.0,solution2,rowSolution2); rowLower2 = model2->rowLower(); rowUpper2 = model2->rowUpper(); #if 0 // redo as Dantzig says for (iColumn=0;iColumn<numberColumns;iColumn++) { Status xStatus = getColumnStatus(iColumn); bool isSuperBasic; int iS = iColumn+newNumberRows; model2->setColumnStatus(iColumn,xStatus); if (xStatus==basic) { model2->setColumnStatus(iS,isFree); solution2[iS]=0.0; } else { model2->setColumnStatus(iS,basic); } // take slack out on equality rows model2->setRowBasic(iColumn+numberRows,superBasic); } for (iRow=0;iRow<numberRows;iRow++) { Status rowStatus = getRowStatus(iRow); model2->setRowStatus(iRow,rowStatus); if (rowStatus!=basic) { model2->setColumnStatus(iRow+numberColumns,basic); // make dual basic } assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_); assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_); } // why ?? - take duals out and adjust for (iRow=0;iRow<numberRows;iRow++) { model2->setRowStatus(iRow,basic); model2->setColumnStatus(iRow+numberColumns,superBasic); solution2[iRow+numberColumns]=0.0; } #else for (iRow=0;iRow<numberRows;iRow++) { assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_); assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_); model2->setRowStatus(iRow,basic); int jRow = iRow; model2->setColumnStatus(jRow+numberColumns,superBasic); solution2[jRow+numberColumns]=0.0; } for (iColumn=numberQuadraticRows+numberColumns;iColumn<newNumberColumns;iColumn++) { model2->setColumnStatus(iColumn,basic); model2->setRowStatus(iColumn-numberQuadraticRows-numberColumns+numberRows,isFixed); } #endif memset(rowSolution2,0,newNumberRows*sizeof(double)); model2->times(1.0,solution2,rowSolution2); for (iColumn=0;iColumn<numberQuadratic;iColumn++) { int iS = iColumn+numberColumns+numberQuadraticRows; int iRow = iColumn+numberRows; double value = rowSolution2[iRow]; if (value>rowUpper2[iRow]) { rowSolution2[iRow] = rowUpper2[iRow]; solution2[iS]-=value-rowUpper2[iRow]; } else { rowSolution2[iRow] = rowLower2[iRow]; solution2[iS]-=value-rowLower2[iRow]; } } // See if any s sub j have wrong sign and/or use djs from infeasibility objective double objectiveOffset; getDblParam(ClpObjOffset,objectiveOffset); double objValue = -objectiveOffset; for (iColumn=0;iColumn<numberColumns;iColumn++) objValue += objective[iColumn]*solution2[iColumn]; for (iColumn=0;iColumn<numberColumns;iColumn++) { double valueI = solution2[iColumn]; int j; for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { int jColumn = columnQuadratic[j]; double valueJ = solution2[jColumn]; double elementValue = quadraticElement[j]; objValue += 0.5*valueI*valueJ*elementValue; } } printf("Objective value %g\n",objValue); //for (iColumn=0;iColumn<newNumberColumns;iColumn++) //printf("%d %g\n",iColumn,solution2[iColumn]); return (ClpSimplexPrimalQuadratic *) model2; }

int ClpSimplexPrimalQuadratic::primalQuadratic (  int  phase = 2  )

Dantzig's method (actually a mixture with Jensen and King) phase - 0 normal, 1 getting complementary solution, 2 getting basic solution. Returns 0 if okay, 1 if LP infeasible. Definition at line 487 of file ClpSimplexPrimalQuadratic.cpp. References endQuadratic(), ClpModel::getColLower(), ClpModel::getColUpper(), ClpModel::getObjCoefficients(), ClpModel::getRowLower(), ClpModel::getRowUpper(), makeQuadratic(), ClpModel::matrix(), ClpModel::messageHandler(), ClpSimplex::numberPrimalInfeasibilities(), ClpSimplexPrimal::primal(), primalQuadratic2(), CoinMessageHandler::setLogLevel(), CoinMpsIO::setMpsData(), ClpSimplex::setPrimalColumnPivotAlgorithm(), and CoinMpsIO::writeMps(). { // Get a feasible solution if (numberPrimalInfeasibilities()) primal(1); // still infeasible if (numberPrimalInfeasibilities()) return 1; ClpQuadraticInfo info; ClpSimplexPrimalQuadratic * model2 = makeQuadratic(info); if (!model2) { printf("** Not quadratic\n"); primal(1); return 0; } #if 0 CoinMpsIO writer; writer.setMpsData(*model2->matrix(), COIN_DBL_MAX, model2->getColLower(), model2->getColUpper(), model2->getObjCoefficients(), (const char*) 0 /*integrality*/, model2->getRowLower(), model2->getRowUpper(), NULL,NULL); writer.writeMps("xx.mps"); #endif // Now do quadratic //ClpPrimalQuadraticDantzig dantzigP(model2,&info); ClpPrimalColumnDantzig dantzigP; model2->setPrimalColumnPivotAlgorithm(dantzigP); model2->messageHandler()->setLogLevel(63); model2->primalQuadratic2(&info,phase); endQuadratic(model2,info); return 0; }

int ClpSimplexPrimalQuadratic::primalQuadratic2 (  ClpQuadraticInfo *  info, int  phase = 2 )

This is done after first pass phase - 0 normal, 1 getting complementary solution, 2 getting basic solution. Definition at line 521 of file ClpSimplexPrimalQuadratic.cpp. References checkComplementarity(), CoinIndexedVector::clear(), ClpQuadraticInfo::currentPhase(), ClpQuadraticInfo::currentSolution(), ClpQuadraticInfo::linearObjective(), ClpQuadraticInfo::numberQuadraticRows(), ClpQuadraticInfo::numberXColumns(), ClpQuadraticInfo::numberXRows(), ClpModel::objectiveAsObject(), ClpQuadraticInfo::originalObjective(), ClpPrimalColumnPivot::pivotColumn(), ClpQuadraticInfo::quadraticObjective(), ClpMatrixBase::refresh(), ClpSimplex::restoreData(), ClpSimplex::saveData(), ClpQuadraticInfo::sequenceIn(), ClpQuadraticInfo::setCrucialSj(), ClpQuadraticInfo::setCurrentPhase(), ClpQuadraticInfo::setSequenceIn(), ClpSimplex::startup(), statusOfProblemInPrimal(), and whileIterating(). Referenced by primalQuadratic(). { algorithm_ = +2; // save data ClpDataSave data = saveData(); // Only ClpPackedMatrix at present assert ((dynamic_cast< ClpPackedMatrix*>(matrix_))); // Assume problem is feasible // Stuff below will be moved into a class int numberXColumns = info->numberXColumns(); int numberXRows = info->numberXRows(); // initialize - values pass coding and algorithm_ is +1 ClpObjective * saveObj = objectiveAsObject(); setObjectivePointer(info->originalObjective()); factorization_->setBiasLU(0); if (!startup(1)) { //setObjectivePointer(saveObj); // Setup useful stuff info->setCurrentPhase(phase); int lastCleaned=0; // last time objective or bounds cleaned up info->setSequenceIn(-1); info->setCrucialSj(-1); bool deleteCosts=false; if (scalingFlag_>0) { // scale CoinPackedMatrix * quadratic = info->quadraticObjective(); double * objective = info->linearObjective(); // replace column by column double * newElement = new double[numberXColumns]; // scale column copy // get matrix data pointers const int * row = quadratic->getIndices(); const CoinBigIndex * columnStart = quadratic->getVectorStarts(); const int * columnLength = quadratic->getVectorLengths(); const double * elementByColumn = quadratic->getElements(); double direction = optimizationDirection_; // direction is actually scale out not scale in if (direction) direction = 1.0/direction; int iColumn; for (iColumn=0;iColumn<numberXColumns;iColumn++) { int j; double scale = columnScale_[iColumn]*direction; const double * elementsInThisColumn = elementByColumn + columnStart[iColumn]; const int * columnsInThisColumn = row + columnStart[iColumn]; int number = columnLength[iColumn]; assert (number<=numberXColumns); for (j=0;j<number;j++) { int iColumn2 = columnsInThisColumn[j]; newElement[j] = elementsInThisColumn[j]*scale*rowScale_[iColumn2+numberXRows]; } quadratic->replaceVector(iColumn,number,newElement); objective[iColumn] *= direction*columnScale_[iColumn]; } delete [] newElement; deleteCosts=true; } else if (optimizationDirection_!=1.0) { CoinPackedMatrix * quadratic = info->quadraticObjective(); double * objective = info->linearObjective(); // replace column by column double * newElement = new double[numberXColumns]; // get matrix data pointers const CoinBigIndex * columnStart = quadratic->getVectorStarts(); const int * columnLength = quadratic->getVectorLengths(); const double * elementByColumn = quadratic->getElements(); double direction = optimizationDirection_; // direction is actually scale out not scale in if (direction) direction = 1.0/direction; int iColumn; for (iColumn=0;iColumn<numberXColumns;iColumn++) { int j; const double * elementsInThisColumn = elementByColumn + columnStart[iColumn]; int number = columnLength[iColumn]; assert (number<=numberXColumns); for (j=0;j<number;j++) { newElement[j] = elementsInThisColumn[j]*direction; } quadratic->replaceVector(iColumn,number,newElement); objective[iColumn] *= direction; } delete [] newElement; deleteCosts=true; } // Say no pivot has occurred (for steepest edge and updates) pivotRow_=-2; // This says whether to restore things etc int factorType=0; /* Status of problem: 0 - optimal 1 - infeasible 2 - unbounded -1 - iterating -2 - factorization wanted -3 - redo checking without factorization -4 - looks infeasible -5 - looks unbounded */ while (problemStatus_<0) { int iRow,iColumn; // clear for (iRow=0;iRow<4;iRow++) { rowArray_[iRow]->clear(); } for (iColumn=0;iColumn<2;iColumn++) { columnArray_[iColumn]->clear(); } // give matrix (and model costs and bounds a chance to be // refreshed (normally null) matrix_->refresh(this); // If we have done no iterations - special if (lastGoodIteration_==numberIterations_) factorType=3; // may factorize, checks if problem finished statusOfProblemInPrimal(lastCleaned,factorType,progress_,info); if (problemStatus_>=0) break; // declare victory checkComplementarity (info,rowArray_[0],rowArray_[1]); // Say good factorization factorType=1; // Say no pivot has occurred (for steepest edge and updates) pivotRow_=-2; // Check problem phase // We assume all X are feasible if (info->currentSolution()&&info->sequenceIn()<0) { phase=0; int nextSequenceIn=-1; int numberQuadraticRows = info->numberQuadraticRows(); for (iColumn=0;iColumn<numberXColumns;iColumn++) { double value = solution_[iColumn]; double lower = lower_[iColumn]; double upper = upper_[iColumn]; if ((value>lower+primalTolerance_&&value<upper-primalTolerance_) ||getColumnStatus(iColumn)==superBasic) { if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) { if (fabs(dj_[iColumn])>10.0*dualTolerance_|| getColumnStatus(iColumn+numberXColumns+numberQuadraticRows)==basic) { if (phase!=2) { phase=2; nextSequenceIn=iColumn; } } } } } for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;iColumn++) { double value = solution_[iColumn]; double lower = lower_[iColumn]; double upper = upper_[iColumn]; if ((value>lower+primalTolerance_&&value<upper-primalTolerance_) ||getColumnStatus(iColumn)==superBasic) { if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)&&fabs(dj_[iColumn])>10.0*dualTolerance_) { if (phase!=2) { phase=2; nextSequenceIn=iColumn; } } } } info->setSequenceIn(nextSequenceIn); info->setCurrentPhase(phase); } // exit if victory declared dualIn_=0.0; // so no updates if (!info->currentPhase()&&info->sequenceIn()<0&&primalColumnPivot_->pivotColumn(rowArray_[1], rowArray_[2],rowArray_[3], columnArray_[0], columnArray_[1]) < 0) { problemStatus_=0; break; } // Iterate problemStatus_=-1; whileIterating(info); } if (deleteCosts) { // delete scaled copy of objective delete info->originalObjective(); } } setObjectivePointer(saveObj); // clean up finish(); restoreData(data); return problemStatus_; }

int ClpSimplexPrimalQuadratic::primalRow (  CoinIndexedVector *  rowArray, CoinIndexedVector *  rhsArray, CoinIndexedVector *  spareArray, CoinIndexedVector *  spareArray2, ClpQuadraticInfo *  info, int  cleanupIteration )

Row array has pivot column This chooses pivot row. Rhs array is used for distance to next bound (for speed) For speed, we may need to go to a bucket approach when many variables go through bounds On exit rhsArray will have changes in costs of basic variables Initially no go thru Returns 0 - can do normal iteration 1 - losing complementarity cleanupIteration - 0 no, 1 yes, 2 restoring one of x/s in basis Definition at line 1349 of file ClpSimplexPrimalQuadratic.cpp. References ClpNonLinearCost::changeInCost(), CoinIndexedVector::clear(), ClpQuadraticInfo::crucialSj(), ClpQuadraticInfo::currentPhase(), CoinIndexedVector::denseVector(), CoinIndexedVector::getIndices(), CoinIndexedVector::getNumElements(), ClpNonLinearCost::goBack(), ClpNonLinearCost::goBackAll(), ClpNonLinearCost::goThru(), ClpQuadraticInfo::linearObjective(), CoinMessageHandler::logLevel(), CoinMessageHandler::message(), ClpNonLinearCost::nearest(), ClpQuadraticInfo::numberXColumns(), ClpQuadraticInfo::originalObjective(), ClpQuadraticObjective::quadraticObjective(), ClpQuadraticInfo::quadraticObjective(), CoinIndexedVector::setNumElements(), and ClpSimplex::solution(). Referenced by whileIterating(). { int result=1; // sequence stays as row number until end pivotRow_=-1; int numberSwapped=0; int numberRemaining=0; int crucialSj = info->crucialSj(); if (crucialSj<0) result=0; // linear int numberThru =0; // number gone thru a barrier int lastThru =0; // number gone thru a barrier on last time double totalThru=0.0; // for when variables flip double acceptablePivot=1.0e-7; if (factorization_->pivots()) acceptablePivot=1.0e-5; // if we have iterated be more strict double bestEverPivot=acceptablePivot; int lastPivotRow = -1; double lastPivot=0.0; double lastTheta=1.0e50; int lastNumberSwapped=0; // use spareArrays to put ones looked at in // First one is list of candidates // We could compress if we really know we won't need any more // Second array has current set of pivot candidates // with a backup list saved in double * part of indexed vector // for zeroing out arrays after int maximumSwapped=0; // pivot elements double * spare; // indices int * index, * indexSwapped; int * saveSwapped; spareArray->clear(); spareArray2->clear(); spare = spareArray->denseVector(); index = spareArray->getIndices(); saveSwapped = (int *) spareArray2->denseVector(); indexSwapped = spareArray2->getIndices(); // we also need somewhere for effective rhs double * rhs=rhsArray->denseVector(); /* First we get a list of possible pivots. We can also see if the problem looks unbounded. At first we increase theta and see what happens. We start theta at a reasonable guess. If in right area then we do bit by bit. We save possible pivot candidates */ // do first pass to get possibles // We can also see if unbounded double * work=rowArray->denseVector(); int number=rowArray->getNumElements(); int * which=rowArray->getIndices(); // we need to swap sign if coming in from ub double way = directionIn_; double maximumMovement; if (way>0.0) maximumMovement = min(1.0e30,upperIn_-valueIn_); else maximumMovement = min(1.0e30,valueIn_-lowerIn_); int iIndex; // Work out coefficients for quadratic term // This is expanded one const CoinPackedMatrix * quadratic = info->quadraticObjective(); const int * columnQuadratic = quadratic->getIndices(); const int * columnQuadraticStart = quadratic->getVectorStarts(); const int * columnQuadraticLength = quadratic->getVectorLengths(); const double * quadraticElement = quadratic->getElements(); //const double * originalCost = info->originalObjective(); // Use rhsArray rhsArray->clear(); int * index2 = rhsArray->getIndices(); int numberXColumns = info->numberXColumns(); int number2=0; //int numberOriginalRows = info->numberXRows(); // sj int iSjRow=-1; // sj for incoming int iSjRow2=-1,crucialSj2=-1; if (sequenceIn_<numberXColumns) { crucialSj2 = sequenceIn_; crucialSj2 += numberRows_; // sj2 which should go to 0.0 } else if (sequenceIn_>=numberColumns_) { int iRow=sequenceIn_-numberColumns_; crucialSj2 = iRow; crucialSj2 += numberXColumns; // sj2 which should go to 0.0 } double tj = 0.0; // Change in objective will be theta*coeff1 + theta*theta*coeff2 double coeff1 = 0.0; double coeff2 = 0.0; double coeffSlack=0.0; for (iIndex=0;iIndex<number;iIndex++) { int iRow = which[iIndex]; double alpha = -work[iRow]*way; int iPivot=pivotVariable_[iRow]; if (numberIterations_==17) printf("row %d col %d alpha %g solution %g\n",iRow,iPivot,alpha,solution_[iPivot]); if (iPivot<numberXColumns) { index2[number2++]=iPivot; rhs[iPivot]=alpha; coeff1 += alpha*cost_[iPivot]; //printf("col %d alpha %g solution %g cost %g scale %g\n",iPivot,alpha,solution_[iPivot], // cost_[iPivot],columnScale_[iPivot]); } else { if (iPivot>=numberColumns_) { // ? do we go through column of pi assert (iPivot!=crucialSj); coeffSlack += alpha*cost_[iPivot]; } else if (iPivot<numberRows_) { // ? do we go through column of pi if (iPivot==crucialSj) { tj = alpha; iSjRow = iRow; } } else { if (iPivot==crucialSj) { tj = alpha; iSjRow = iRow; } } } if (iPivot==crucialSj2) iSjRow2=iRow; } // Incoming if (sequenceIn_<numberXColumns) { index2[number2++]=sequenceIn_; rhs[sequenceIn_]=way; //printf("incoming col %d alpha %g solution %g cost %g scale %g\n",sequenceIn_,way,valueIn_, // cost_[sequenceIn_],columnScale_[sequenceIn_]); coeff1 += way*cost_[sequenceIn_]; } else { //printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_); coeffSlack += way*cost_[sequenceIn_]; } printf("coeff1 now %g\n",coeff1); rhsArray->setNumElements(number2); double largestCoeff1=1.0e-20; for (iIndex=0;iIndex<number2;iIndex++) { int iColumn=index2[iIndex]; //double valueI = solution_[iColumn]; double alphaI = rhs[iColumn]; int j; for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { int jColumn = columnQuadratic[j]; double valueJ = solution_[jColumn]; double alphaJ = rhs[jColumn]; double elementValue = quadraticElement[j]; double addValue = (valueJ*alphaI)*elementValue; largestCoeff1 = max(largestCoeff1,fabs(addValue)); coeff1 += addValue; coeff2 += (alphaI*alphaJ)*elementValue; } } coeff2 *= 0.5; if (coeffSlack) printf("slack has non-zero cost - trouble?\n"); coeff1 += coeffSlack; //coeff1 = way*dualIn_; int accuracyFlag=0; if (!cleanupIteration) { if (fabs(dualIn_)<dualTolerance_&&fabs(coeff1)<dualTolerance_) { dualIn_=0.0; coeff1=0.0; } if (fabs(way*coeff1-dualIn_)>1.0e-2*(1.0+fabs(dualIn_))) printf("primal error %g, largest %g, coeff1 %g, dualin %g\n", largestPrimalError_,largestCoeff1,way*coeff1,dualIn_); assert (fabs(way*coeff1-dualIn_)<1.0e-1*(1.0+fabs(dualIn_))); assert (way*coeff1*dualIn_>=0.0); if (way*coeff1*dualIn_<0.0) { // bad accuracyFlag=2; } else if (fabs(way*coeff1-dualIn_)>1.0e-3*(1.0+fabs(dualIn_))) { // not wonderful accuracyFlag=1; } if (crucialSj>=0) { solution_[crucialSj]=way*coeff1; } } else if (dualIn_<0.0&&way*coeff1>1.0e-2) { printf("bad pair\n"); accuracyFlag=1; } else if (dualIn_>0.0&&way*coeff1<-1.0e-2) { accuracyFlag=1; printf("bad pair\n"); } // interesting places are where derivative zero or sj goes to zero double d1,d2=1.0e50; if (fabs(coeff2)>1.0e-9) d1 = - 0.5*coeff1/coeff2; else if (coeff1<=1.0e-6) d1 = maximumMovement; else d1 = 0.0; if (fabs(tj)<1.0e-7||!cleanupIteration) { //if (sequenceIn_<numberXColumns) //printf("column %d is basically linear\n",sequenceIn_); //assert(!columnQuadraticLength[sequenceIn_]); } else { d2 = -solution_[crucialSj]/tj; if (d2<0.0) { printf("d2 would be negative at %g\n",d2); d2=1.0e50; } } if (d1>1.0e10&&d2>1.0e10) { // linear variable entering // maybe we should have done dual iteration to force sj to 0.0 //printf("linear variable\n"); } handler_->message(CLP_QUADRATIC_PRIMAL_DETAILS,messages_) <<coeff1<<coeff2<<coeffSlack<<dualIn_<<d1<<d2 <<CoinMessageEol; // reality check { if (crucialSj2>=0) { // see if works out if (getColumnStatus(crucialSj2)==basic) { if (iSjRow2>=0) { //corresponding alpha nonzero double alpha=work[iSjRow2]; printf("Sj alpha %g sol %g ratio %g\n",alpha,solution_[crucialSj2],solution_[crucialSj2]/alpha); if (fabs(fabs(d1)-fabs(solution_[crucialSj2]/alpha))>1.0e-3*(1.0+fabs(d1))) { printf("bad test\n"); if (factorization_->pivots()) accuracyFlag=2; } d1=fabs(solution_[crucialSj2]/alpha); } else { printf("Sj basic but no alpha\n"); } } else { printf("Sj not basic\n"); } } } //if (!cleanupIteration) //dualIn_ = way*coeff1; double mind1d2=1.0e50; if (cleanupIteration) { // we are only interested in d1 if smaller than d2 mind1d2 = min(maximumMovement,d2); //if (d1>1.0e-8&&d1<0.999*d2) //mind1d2=d1; } else { // There is no d2 - d1 refers to crucialSj if (d1>1.0e-5) { mind1d2 = min(maximumMovement,d1); //assert (d1>=0.0); } } maximumMovement = min(maximumMovement,mind1d2); double trueMaximumMovement; if (way>0.0) trueMaximumMovement = min(1.0e30,upperIn_-valueIn_); else trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_); rhsArray->clear(); double tentativeTheta = maximumMovement; double upperTheta = maximumMovement; bool throwAway=false; if (numberIterations_==1750) { printf("Bad iteration coming up after iteration %d\n",numberIterations_); } double minimumAny=1.0e50; int whichMinimum=-1; double minimumAlpha=0.0; for (iIndex=0;iIndex<number;iIndex++) { int iRow = which[iIndex]; double alpha = work[iRow]; int iPivot=pivotVariable_[iRow]; alpha *= way; double oldValue = solution(iPivot); // get where in bound sequence if (alpha>0.0) { // basic variable going towards lower bound double bound = lower(iPivot); oldValue -= bound; } else if (alpha<0.0) { // basic variable going towards upper bound double bound = upper(iPivot); oldValue = bound-oldValue; } if (iPivot>=numberXColumns&&iPivot<numberColumns_) { double bound=0.0; double oldValue2 = solution(iPivot); if (alpha>0.0) { // basic variable going towards lower bound oldValue2 -= bound; } else if (alpha<0.0) { // basic variable going towards upper bound oldValue2 = bound-oldValue; } double value = oldValue2-minimumAny*fabs(alpha); if (value*oldValue2<0.0) { double ratio = fabs(oldValue2/alpha); if (ratio<minimumAny&&fabs(alpha)>1.0e2*acceptablePivot) { minimumAny = ratio; whichMinimum = iRow; minimumAlpha= fabs(alpha); } } } double value = oldValue-tentativeTheta*fabs(alpha); assert (oldValue>=-primalTolerance_*1.0001); if (value<-primalTolerance_) { // add to list spare[numberRemaining]=alpha; rhs[iRow]=oldValue; index[numberRemaining++]=iRow; double value=oldValue-upperTheta*fabs(alpha); if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot) upperTheta = (oldValue+primalTolerance_)/fabs(alpha); } } // we need to keep where rhs non-zeros are int numberInRhs=numberRemaining; memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int)); rhsArray->setNumElements(numberInRhs); theta_=maximumMovement; bool goBackOne = false; double fake=1.0e-2; if (numberRemaining) { // looks like pivoting // now try until reasonable theta tentativeTheta = max(10.0*upperTheta,1.0e-7); tentativeTheta = min(tentativeTheta,maximumMovement); // loops increasing tentative theta until can't go through while (tentativeTheta <= maximumMovement) { double bestPivotBeforeInteresting=0.0; double thruThis = 0.0; double bestPivot=acceptablePivot; pivotRow_ = -1; numberSwapped = 0; upperTheta = maximumMovement; for (iIndex=0;iIndex<numberRemaining;iIndex++) { int iRow = index[iIndex]; double alpha = spare[iIndex]; double oldValue = rhs[iRow]; double value = oldValue-tentativeTheta*fabs(alpha); if (value<-primalTolerance_) { // how much would it cost to go thru thruThis += alpha* nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha); // goes on swapped list (also means candidates if too many) indexSwapped[numberSwapped++]=iRow; if (fabs(alpha)>bestPivot) { bestPivot=fabs(alpha); pivotRow_ = iRow; theta_ = oldValue/bestPivot; } } else { value = oldValue-upperTheta*fabs(alpha); if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) upperTheta = (oldValue+primalTolerance_)/fabs(alpha); } } maximumSwapped = max(maximumSwapped,numberSwapped); bestPivotBeforeInteresting=bestPivot; //double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 0.1*infeasibilityCost_; double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1; // but make a bit more pessimistic if pivot //if (bestPivot>acceptablePivot) //dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck); //dualCheck=0.0; if ((totalThru+thruThis>=dualCheck&&bestPivot>acceptablePivot)||fake*bestPivot>1.0e-3) { // We should be pivoting in this batch // so compress down to this lot int saveNumber = numberRemaining; numberRemaining=0; for (iIndex=0;iIndex<numberSwapped;iIndex++) { int iRow = indexSwapped[iIndex]; spare[numberRemaining]=way*work[iRow]; index[numberRemaining++]=iRow; } memset(spare+numberRemaining,0, (saveNumber-numberRemaining)*sizeof(double)); double lastThru = totalThru; int iTry; #define MAXTRY 100 // first get ratio with tolerance for (iTry=0;iTry<MAXTRY;iTry++) { upperTheta=maximumMovement; numberSwapped = 0; for (iIndex=0;iIndex<numberRemaining;iIndex++) { int iRow = index[iIndex]; double alpha = fabs(spare[iIndex]); double oldValue = rhs[iRow]; double value = oldValue-upperTheta*alpha; if (value<-primalTolerance_ && alpha>=acceptablePivot) upperTheta = (oldValue+primalTolerance_)/alpha; } // now look at best in this lot bestPivot=acceptablePivot; pivotRow_=-1; for (iIndex=0;iIndex<numberRemaining;iIndex++) { int iRow = index[iIndex]; double alpha = spare[iIndex]; double oldValue = rhs[iRow]; double value = oldValue-upperTheta*fabs(alpha); if (value<=0) { // how much would it cost to go thru totalThru += alpha* nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha); // goes on swapped list (also means candidates if too many) indexSwapped[numberSwapped++]=iRow; if (fabs(alpha)>bestPivot) { bestPivot=fabs(alpha); theta_ = oldValue/bestPivot; pivotRow_=iRow; } } else { value = oldValue-upperTheta*fabs(alpha); if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) upperTheta = (oldValue+primalTolerance_)/fabs(alpha); } } maximumSwapped = max(maximumSwapped,numberSwapped); if (bestPivot<0.1*bestEverPivot&& bestEverPivot>1.0e-6&&bestPivot<1.0e-3) { // back to previous one goBackOne = true; break; } else if (pivotRow_==-1&&upperTheta>largeValue_) { if (lastPivot>acceptablePivot) { // back to previous one goBackOne = true; //break; } else { // can only get here if all pivots so far too small } break; } else { dualCheck = - 2.0*coeff2*theta_ - coeff1; if (bestPivotBeforeInteresting>1.0e-4&&bestPivot<1.0e-6) dualCheck=1.0e7; if ((totalThru>=dualCheck||fake*bestPivot>1.0e-3) &&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) { if (!cleanupIteration) { // We can pivot out sj if (upperTheta==maximumMovement) { printf("figures\n"); } else if (iSjRow2>=0) { if (lastThru>=dualCheck&&theta_<trueMaximumMovement) { printf("totalThru %g, lastThru %g - taking sj to zero?\n",totalThru,lastThru); // make sure will take correct path mind1d2=1.0; maximumMovement=1.0; d2=0.0; pivotRow_=-1; } } } else { if (pivotRow_<0&&lastPivotRow<0) throwAway=true; } break; // no point trying another loop } else if (totalThru>=dualCheck||fake*bestPivot>1.0e-3||upperTheta==maximumMovement) { break; // no point trying another loop } else { // skip this lot nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs); lastPivotRow=pivotRow_; lastTheta = theta_; lastThru = numberThru; numberThru += numberSwapped; lastNumberSwapped = numberSwapped; memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int)); if (bestPivot>bestEverPivot) bestEverPivot=bestPivot; lastThru = totalThru; } } } break; } else { // skip this lot nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs); lastPivotRow=pivotRow_; lastTheta = theta_; lastThru = numberThru; numberThru += numberSwapped; lastNumberSwapped = numberSwapped; memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int)); if (bestPivot>bestEverPivot) bestEverPivot=bestPivot; totalThru += thruThis; tentativeTheta = 2.0*upperTheta; } } // can get here without pivotRow_ set but with lastPivotRow if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) { // back to previous one pivotRow_=lastPivotRow; theta_ = lastTheta; // undo this lot nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs); memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int)); numberSwapped = lastNumberSwapped; } } if (0&&minimumAny<theta_&&cleanupIteration&&bestEverPivot<1.0e-2 &&minimumAlpha>bestEverPivot*1.0e2&&minimumAny>1.0e-1&&info->currentPhase()==0) { printf("Possible other pivot %d %g %g - alpha %g\n", whichMinimum,minimumAny,theta_,minimumAlpha); pivotRow_ = whichMinimum; alpha_ = work[pivotRow_]; // translate to sequence sequenceOut_ = pivotVariable_[pivotRow_]; valueOut_ = solution(sequenceOut_); lowerOut_=0.0; upperOut_=0.0; theta_= fabs(valueOut_/alpha_); if (way<0.0) theta_ = - theta_; if (alpha_*way<0.0) { directionOut_=-1; // to upper bound } else { directionOut_=1; // to lower bound } dualOut_ = reducedCost(sequenceOut_); } else { if (pivotRow_>=0) { if (0) { double move = coeff1*theta_+coeff2*theta_*theta_; printf("Predicted change in obj is %g %g %g\n", move,objectiveValue_-move,objectiveValue_+move); } #define MINIMUMTHETA 1.0e-12 // will we need to increase tolerance #ifdef CLP_DEBUG bool found=false; #endif double largestInfeasibility = primalTolerance_; if (theta_<MINIMUMTHETA) { theta_=MINIMUMTHETA; for (iIndex=0;iIndex<numberSwapped;iIndex++) { int iRow = indexSwapped[iIndex]; #ifdef CLP_DEBUG if (iRow==pivotRow_) found=true; #endif largestInfeasibility = max (largestInfeasibility, -(rhs[iRow]-fabs(work[iRow])*theta_)); } #ifdef CLP_DEBUG assert(found); if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32)) printf("Primal tolerance increased from %g to %g\n", primalTolerance_,largestInfeasibility); #endif primalTolerance_ = max(primalTolerance_,largestInfeasibility); } alpha_ = work[pivotRow_]; // translate to sequence sequenceOut_ = pivotVariable_[pivotRow_]; valueOut_ = solution(sequenceOut_); lowerOut_=lower_[sequenceOut_]; upperOut_=upper_[sequenceOut_]; if (way<0.0) theta_ = - theta_; double newValue = valueOut_ - theta_*alpha_; if (alpha_*way<0.0) { directionOut_=-1; // to upper bound if (fabs(theta_)>1.0e-6) upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue); else upperOut_ = newValue; } else { directionOut_=1; // to lower bound if (fabs(theta_)>1.0e-6) lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue); else lowerOut_ = newValue; } dualOut_ = reducedCost(sequenceOut_); } else { // If no pivot but bad move then throw away if (throwAway&&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) { assert (pivotRow_<0); accuracyFlag=2; } if (maximumMovement<1.0e20&&maximumMovement==trueMaximumMovement) { // flip theta_ = maximumMovement; pivotRow_ = -2; // so we can tell its a flip result=0; sequenceOut_ = sequenceIn_; valueOut_ = valueIn_; dualOut_ = dualIn_; lowerOut_ = lowerIn_; upperOut_ = upperIn_; alpha_ = 0.0; if (accuracyFlag<2) accuracyFlag=0; if (way<0.0) { directionOut_=1; // to lower bound theta_ = lowerOut_ - valueOut_; } else { directionOut_=-1; // to upper bound theta_ = upperOut_ - valueOut_; } // we may still have sj to get rid of } else if (fabs(maximumMovement-mind1d2)<dualTolerance_) { // crucialSj local copy i.e. dj going to zero if (maximumMovement==d2) { // sj going to zero result=0; } else { // incoming dj going to zero if (!cleanupIteration) { result=0; iSjRow=iSjRow2; if (iSjRow>=0) crucialSj = pivotVariable_[iSjRow]; else crucialSj=-1; assert (pivotRow_<0); // If crucialSj <0 then make superbasic? if (crucialSj<0) { printf("**ouch\n"); theta_ = maximumMovement; pivotRow_ = -2; // so we can tell its a flip result=0; sequenceOut_ = sequenceIn_; valueOut_ = valueIn_; dualOut_ = dualIn_; lowerOut_ = lowerIn_; upperOut_ = upperIn_; alpha_ = 0.0; if (accuracyFlag<2) accuracyFlag=0; if (way<0.0) { directionOut_=1; // to lower bound lowerIn_ = valueOut_-theta_; theta_ = lowerIn_ - valueOut_; } else { directionOut_=-1; // to upper bound upperIn_ = valueOut_+theta_; theta_ = upperIn_ - valueOut_; } } } else { assert (fabs(dualIn_)<1.0e-3); result=1; if (minimumAlpha>0.0) { // could do this in other places if pivot looks good printf("Should take minimum\n"); pivotRow_ = whichMinimum; alpha_ = work[pivotRow_]; // translate to sequence sequenceOut_ = pivotVariable_[pivotRow_]; valueOut_ = solution(sequenceOut_); theta_ = minimumAny; if (way<0.0) theta_ = - theta_; lowerOut_=0.0; upperOut_=0.0; //???? dualOut_ = reducedCost(sequenceOut_); } } } if (!result&&crucialSj>=0) { setColumnStatus(crucialSj,isFree); pivotRow_ = iSjRow; alpha_ = work[pivotRow_]; // translate to sequence sequenceOut_ = pivotVariable_[pivotRow_]; assert (sequenceOut_==crucialSj); valueOut_ = solution(sequenceOut_); theta_ = fabs(valueOut_/alpha_); if (fabs(maximumMovement-theta_)>1.0e-3*(1.0+maximumMovement)) printf("maximumMovement %g mismatch with theta %g\n",maximumMovement,theta_);; if (way<0.0) theta_ = - theta_; lowerOut_=0.0; upperOut_=0.0; //???? dualOut_ = reducedCost(sequenceOut_); } } else { // need to do something accuracyFlag=2; } } } // clear arrays memset(spare,0,numberRemaining*sizeof(double)); memset(saveSwapped,0,maximumSwapped*sizeof(int)); // put back original bounds etc nonLinearCost_->goBackAll(rhsArray); rhsArray->clear(); // Change in objective will be theta*coeff1 + theta*theta*coeff2 objectiveValue_ += theta_*coeff1+theta_*theta_*coeff2; { int iColumn; objectiveValue_ =0.0; CoinPackedMatrix * quadratic = info->originalObjective()->quadraticObjective(); const int * columnQuadratic = quadratic->getIndices(); const int * columnQuadraticStart = quadratic->getVectorStarts(); const int * columnQuadraticLength = quadratic->getVectorLengths(); const double * quadraticElement = quadratic->getElements(); int numberColumns = info->numberXColumns(); const double * objective = info->linearObjective(); for (iColumn=0;iColumn<numberColumns;iColumn++) objectiveValue_ += objective[iColumn]*solution_[iColumn]; for (iColumn=0;iColumn<numberColumns;iColumn++) { double valueI = solution_[iColumn]; int j; for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { int jColumn = columnQuadratic[j]; double valueJ = solution_[jColumn]; double elementValue = quadraticElement[j]; objectiveValue_ += 0.5*valueI*valueJ*elementValue; } } } if (accuracyFlag<2) { return result+10*accuracyFlag; } else { pivotRow_=1; // make sure correct path return 20; } }

int ClpSimplexPrimalQuadratic::primalSLP (  int  numberPasses, double  deltaTolerance )

A sequential LP method. Primal algorithms for quadratic At present we have two algorithms: a) Dantzig's algorithm b) Using a semi-trust region approach as for pooling problem This is in because I have it lying around Definition at line 58 of file ClpSimplexPrimalQuadratic.cpp. References ClpModel::columnLower(), ClpModel::columnUpper(), ClpModel::dualColumnSolution(), ClpModel::matrix(), ClpSimplex::numberPrimalInfeasibilities(), ClpModel::numberRows(), ClpModel::objective(), ClpModel::optimizationDirection(), ClpSimplexPrimal::primal(), ClpModel::primalColumnSolution(), ClpQuadraticObjective::quadraticObjective(), ClpSimplex::scaling(), ClpModel::setDblParam(), ClpModel::setLogLevel(), CoinMpsIO::setMpsData(), ClpModel::status(), and CoinMpsIO::writeMps(). { // Are we minimizing or maximizing double whichWay=optimizationDirection(); if (whichWay<0.0) whichWay=-1.0; else if (whichWay>0.0) whichWay=1.0; // This is as a user would see int numberColumns = this->numberColumns(); int numberRows = this->numberRows(); double * columnLower = this->columnLower(); double * columnUpper = this->columnUpper(); double * objective = this->objective(); double * solution = this->primalColumnSolution(); // Save objective double * saveObjective = new double [numberColumns]; memcpy(saveObjective,objective,numberColumns*sizeof(double)); // Get list of non linear columns ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_)); CoinPackedMatrix * quadratic = NULL; if (quadraticObj) quadratic = quadraticObj->quadraticObjective(); if (!quadratic) { // no quadratic part return primal(0); } int numberNonLinearColumns = 0; int iColumn; int * listNonLinearColumn = new int[numberColumns]; memset(listNonLinearColumn,0,numberColumns*sizeof(int)); const int * columnQuadratic = quadratic->getIndices(); const int * columnQuadraticStart = quadratic->getVectorStarts(); const int * columnQuadraticLength = quadratic->getVectorLengths(); const double * quadraticElement = quadratic->getElements(); for (iColumn=0;iColumn<numberColumns;iColumn++) { int j; for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { int jColumn = columnQuadratic[j]; listNonLinearColumn[jColumn]=1; listNonLinearColumn[iColumn]=1; } } for (iColumn=0;iColumn<numberColumns;iColumn++) { if(listNonLinearColumn[iColumn]) listNonLinearColumn[numberNonLinearColumns++]=iColumn; } if (!numberNonLinearColumns) { delete [] listNonLinearColumn; // no quadratic part return primal(0); } // get feasible if (!this->status()||numberPrimalInfeasibilities()) primal(1); // still infeasible if (numberPrimalInfeasibilities()) return 0; int jNon; int * last[3]; double * trust = new double[numberNonLinearColumns]; double * trueLower = new double[numberNonLinearColumns]; double * trueUpper = new double[numberNonLinearColumns]; for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; trust[jNon]=0.5; trueLower[jNon]=columnLower[iColumn]; trueUpper[jNon]=columnUpper[iColumn]; if (solution[iColumn]<trueLower[jNon]) solution[iColumn]=trueLower[jNon]; else if (solution[iColumn]>trueUpper[jNon]) solution[iColumn]=trueUpper[jNon]; } int iPass; double lastObjective=1.0e31; double * saveSolution = new double [numberColumns]; double * savePi = new double [numberRows]; unsigned char * saveStatus = new unsigned char[numberRows+numberColumns]; double targetDrop=1.0e31; double objectiveOffset; getDblParam(ClpObjOffset,objectiveOffset); // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change for (iPass=0;iPass<3;iPass++) { last[iPass]=new int[numberNonLinearColumns]; for (jNon=0;jNon<numberNonLinearColumns;jNon++) last[iPass][jNon]=0; } // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing int goodMove=-2; char * statusCheck = new char[numberColumns]; for (iPass=0;iPass<numberPasses;iPass++) { // redo objective double offset=0.0; double objValue=-objectiveOffset; double lambda=-1.0; if (goodMove>=0) { // get best interpolation double coeff0=-objectiveOffset,coeff1=0.0,coeff2=0.0; for (iColumn=0;iColumn<numberColumns;iColumn++) { coeff0 += saveObjective[iColumn]*solution[iColumn]; coeff1 += saveObjective[iColumn]*(saveSolution[iColumn]-solution[iColumn]); } for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; double valueI = solution[iColumn]; double valueISave = saveSolution[iColumn]; int j; for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { int jColumn = columnQuadratic[j]; double valueJ = solution[jColumn]; double valueJSave = saveSolution[jColumn]; double elementValue = 0.5*quadraticElement[j]; coeff0 += valueI*valueJ*elementValue; coeff1 += (valueI*valueJSave+valueISave*valueJ-2.0*valueI*valueJ)*elementValue; coeff2 += (valueISave*valueJSave+valueI*valueJ-valueISave*valueJ-valueI*valueJSave)*elementValue; } } double lambdaValue; if (fabs(coeff2)<1.0e-9) { if (coeff1+coeff2>=0.0) lambda = 0.0; else lambda = 1.0; } else { lambda = -(0.5*coeff1)/coeff2; if (lambda>1.0||lambda<0.0) { if (coeff1+coeff2>=0.0) lambda = 0.0; else lambda = 1.0; } } lambdaValue = lambda*lambda*coeff2+lambda*coeff1+coeff0; printf("coeffs %g %g %g - lastobj %g\n",coeff0,coeff1,coeff2,lastObjective); printf("obj at saved %g, obj now %g zero deriv at %g - value %g\n", coeff0+coeff1+coeff2,coeff0,lambda,lambdaValue); if (!iPass) lambda=0.0; if (lambda>0.0&&lambda<=1.0) { // update solution for (iColumn=0;iColumn<numberColumns;iColumn++) solution[iColumn] = lambda * saveSolution[iColumn] + (1.0-lambda) * solution[iColumn]; if (lambda>0.999) { memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double)); memcpy(status_,saveStatus,numberRows+numberColumns); } if (lambda>0.99999&&fabs(coeff1+coeff2)>1.0e-2) { // tighten all goodMove=-1; } } } memcpy(objective,saveObjective,numberColumns*sizeof(double)); for (iColumn=0;iColumn<numberColumns;iColumn++) objValue += objective[iColumn]*solution[iColumn]; for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; if (getColumnStatus(iColumn)==basic) { if (solution[iColumn]<columnLower[iColumn]+1.0e-8) statusCheck[iColumn]='l'; else if (solution[iColumn]>columnUpper[iColumn]-1.0e-8) statusCheck[iColumn]='u'; else statusCheck[iColumn]='B'; } else { if (solution[iColumn]<columnLower[iColumn]+1.0e-8) statusCheck[iColumn]='L'; else statusCheck[iColumn]='U'; } double valueI = solution[iColumn]; int j; for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { int jColumn = columnQuadratic[j]; double valueJ = solution[jColumn]; double elementValue = quadraticElement[j]; elementValue *= 0.5; objValue += valueI*valueJ*elementValue; offset += valueI*valueJ*elementValue; double gradientI = valueJ*elementValue; double gradientJ = valueI*elementValue; offset -= gradientI*valueI; objective[iColumn] += gradientI; offset -= gradientJ*valueJ; objective[jColumn] += gradientJ; } } printf("objective %g, objective offset %g\n",objValue,offset); setDblParam(ClpObjOffset,objectiveOffset-offset); objValue *= whichWay; int * temp=last[2]; last[2]=last[1]; last[1]=last[0]; last[0]=temp; for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; double change = solution[iColumn]-saveSolution[iColumn]; if (change<-1.0e-5) { if (fabs(change+trust[jNon])<1.0e-5) temp[jNon]=-1; else temp[jNon]=-2; } else if(change>1.0e-5) { if (fabs(change-trust[jNon])<1.0e-5) temp[jNon]=1; else temp[jNon]=2; } else { temp[jNon]=0; } } // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing double maxDelta=0.0; if (goodMove>=0) { if (objValue<=lastObjective) goodMove=1; else goodMove=0; } else { maxDelta=1.0e10; } double maxGap=0.0; for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; maxDelta = max(maxDelta, fabs(solution[iColumn]-saveSolution[iColumn])); if (goodMove>0) { if (last[0][jNon]*last[1][jNon]<0) { // halve trust[jNon] *= 0.5; } else { if (last[0][jNon]==last[1][jNon]&& last[0][jNon]==last[2][jNon]) trust[jNon] *= 1.5; } } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) { trust[jNon] *= 0.2; } maxGap = max(maxGap,trust[jNon]); } std::cout<<"largest gap is "<<maxGap<<std::endl; if (iPass>10000) { for (jNon=0;jNon<numberNonLinearColumns;jNon++) trust[jNon] *=0.0001; } if (goodMove>0) { double drop = lastObjective-objValue; std::cout<<"True drop was "<<drop<<std::endl; std::cout<<"largest delta is "<<maxDelta<<std::endl; if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&lambda<0.99999) { std::cout<<"Exiting"<<std::endl; break; } } if (!iPass) goodMove=1; targetDrop=0.0; double * r = this->dualColumnSolution(); for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; columnLower[iColumn]=max(solution[iColumn] -trust[jNon], trueLower[jNon]); columnUpper[iColumn]=min(solution[iColumn] +trust[jNon], trueUpper[jNon]); } if (iPass) { // get reduced costs this->matrix()->transposeTimes(savePi, this->dualColumnSolution()); for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; double dj = objective[iColumn]-r[iColumn]; r[iColumn]=dj; if (dj<0.0) targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]); else targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]); } } else { memset(r,0,numberColumns*sizeof(double)); } #if 0 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) { columnLower[iColumn]=max(solution[iColumn], trueLower[jNon]); columnUpper[iColumn]=min(solution[iColumn] +trust[jNon], trueUpper[jNon]); } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) { columnLower[iColumn]=max(solution[iColumn] -trust[jNon], trueLower[jNon]); columnUpper[iColumn]=min(solution[iColumn], trueUpper[jNon]); } else { columnLower[iColumn]=max(solution[iColumn] -trust[jNon], trueLower[jNon]); columnUpper[iColumn]=min(solution[iColumn] +trust[jNon], trueUpper[jNon]); } } #endif if (goodMove) { memcpy(saveSolution,solution,numberColumns*sizeof(double)); memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double)); memcpy(saveStatus,status_,numberRows+numberColumns); std::cout<<"Pass - "<<iPass <<", target drop is "<<targetDrop <<std::endl; lastObjective = objValue; if (targetDrop<1.0e-5&&goodMove&&iPass) { printf("Exiting on target drop %g\n",targetDrop); break; } { double * r = this->dualColumnSolution(); for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n", jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn], r[iColumn],statusCheck[iColumn],columnLower[iColumn], columnUpper[iColumn]); } } setLogLevel(63); this->scaling(false); this->primal(1); if (this->status()) { CoinMpsIO writer; writer.setMpsData(*matrix(), COIN_DBL_MAX, getColLower(), getColUpper(), getObjCoefficients(), (const char*) 0 /*integrality*/, getRowLower(), getRowUpper(), NULL,NULL); writer.writeMps("xx.mps"); } assert (!this->status()); // must be feasible goodMove=1; } else { // bad pass - restore solution printf("Backtracking\n"); memcpy(solution,saveSolution,numberColumns*sizeof(double)); memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double)); memcpy(status_,saveStatus,numberRows+numberColumns); iPass--; goodMove=-1; } } // restore solution memcpy(solution,saveSolution,numberColumns*sizeof(double)); for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; columnLower[iColumn]=max(solution[iColumn], trueLower[jNon]); columnUpper[iColumn]=min(solution[iColumn], trueUpper[jNon]); } delete [] statusCheck; delete [] savePi; delete [] saveStatus; // redo objective double offset=0.0; double objValue=-objectiveOffset; memcpy(objective,saveObjective,numberColumns*sizeof(double)); for (iColumn=0;iColumn<numberColumns;iColumn++) objValue += objective[iColumn]*solution[iColumn]; for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; double valueI = solution[iColumn]; int j; for (j=columnQuadraticStart[iColumn]; j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { int jColumn = columnQuadratic[j]; double valueJ = solution[jColumn]; double elementValue = quadraticElement[j]; objValue += 0.5*valueI*valueJ*elementValue; offset += 0.5*valueI*valueJ*elementValue; double gradientI = valueJ*elementValue; double gradientJ = valueI*elementValue; offset -= gradientI*valueI; objective[iColumn] += gradientI; offset -= gradientJ*valueJ; objective[jColumn] += gradientJ; } } printf("objective %g, objective offset %g\n",objValue,offset); setDblParam(ClpObjOffset,objectiveOffset-offset); this->primal(1); // redo values setDblParam(ClpObjOffset,objectiveOffset); objectiveValue_ += whichWay*offset; for (jNon=0;jNon<numberNonLinearColumns;jNon++) { iColumn=listNonLinearColumn[jNon]; columnLower[iColumn]= trueLower[jNon]; columnUpper[iColumn]= trueUpper[jNon]; } memcpy(objective,saveObjective,numberColumns*sizeof(double)); delete [] saveSolution; for (iPass=0;iPass<3;iPass++) delete [] last[iPass]; delete [] trust; delete [] trueUpper; delete [] trueLower; delete [] saveObjective; delete [] listNonLinearColumn; return 0; }

void ClpSimplexPrimalQuadratic::statusOfProblemInPrimal (  int &  lastCleaned, int  type, ClpSimplexProgress *  progress, ClpQuadraticInfo *  info )

Refactorizes if necessary Checks if finished. Updates status. lastCleaned refers to iteration at which some objective/feasibility cleaning too place. type - 0 initial so set up save arrays etc 1 normal -if good update save 2 restoring from saved Definition at line 2134 of file ClpSimplexPrimalQuadratic.cpp. References checkComplementarity(), ClpNonLinearCost::checkInfeasibilities(), ClpSimplex::createRim(), ClpQuadraticInfo::crucialSj(), ClpQuadraticInfo::currentPhase(), ClpQuadraticInfo::currentSolution(), ClpSimplex::dualFeasible(), ClpSimplex::gutsOfSolution(), ClpSimplexProgress::lastIterationNumber(), ClpPrimalColumnPivot::looksOptimal(), ClpSimplexProgress::looping(), CoinMessageHandler::message(), ClpNonLinearCost::numberInfeasibilities(), ClpQuadraticInfo::numberXColumns(), ClpModel::objectiveValue(), ClpSimplex::primalFeasible(), CoinMessageHandler::printing(), ClpQuadraticInfo::restoreStatus(), ClpQuadraticInfo::saveStatus(), ClpPrimalColumnPivot::saveWeights(), ClpQuadraticInfo::sequenceIn(), ClpQuadraticInfo::setCurrentSolution(), ClpSimplex::setFlagged(), ClpQuadraticInfo::setSequenceIn(), ClpNonLinearCost::sumInfeasibilities(), ClpSimplexPrimal::unflag(), and ClpSimplexPrimal::unPerturb(). Referenced by primalQuadratic2(). { if (info->currentPhase()) info->setCurrentSolution(solution_); //info->saveStatus(); if (type==2) { // trouble - restore solution memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char)); memcpy(rowActivityWork_,savedSolution_+numberColumns_ , numberRows_*sizeof(double)); memcpy(columnActivityWork_,savedSolution_ , numberColumns_*sizeof(double)); forceFactorization_=1; // a bit drastic but .. pivotRow_=-1; // say no weights update changeMade_++; // say change made info->restoreStatus(); } int saveThreshold = factorization_->sparseThreshold(); int tentativeStatus = problemStatus_; if (problemStatus_>-3||problemStatus_==-4) { // factorize // later on we will need to recover from singularities // also we could skip if first time // do weights // This may save pivotRow_ for use primalColumnPivot_->saveWeights(this,1); if (type) { // is factorization okay? int numberPivots = factorization_->pivots(); if (internalFactorize(1)) { if (solveType_==2) { // say odd problemStatus_=5; return; } numberIterations_ = progress_->lastIterationNumber(0); // no - restore previous basis assert (info->crucialSj()<0); // need to work out how to unwind assert (type==1); memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char)); memcpy(rowActivityWork_,savedSolution_+numberColumns_ , numberRows_*sizeof(double)); memcpy(columnActivityWork_,savedSolution_ , numberColumns_*sizeof(double)); forceFactorization_=1; // a bit drastic but .. type = 2; assert (internalFactorize(1)==0); info->restoreStatus(); // flag incoming if (numberPivots==1) { int crucialSj = info->crucialSj(); if (crucialSj<0) { setFlagged(sequenceIn_); } else { int iSequence; int numberXColumns = info->numberXColumns(); if (crucialSj<numberRows_) { // row iSequence = crucialSj-numberXColumns+numberColumns_; } else { // column iSequence = crucialSj-numberRows_; } if (getColumnStatus(iSequence)!=basic) { setFlagged(iSequence); info->setSequenceIn(-1); } else { printf("**** %d is basic ?\n",iSequence); } } } changeMade_++; // say change made } } if (problemStatus_!=-4) problemStatus_=-3; } if (info->crucialSj()<0) { // at this stage status is -3 or -5 if looks unbounded // get primal and dual solutions // put back original costs and then check createRim(4); if (info->currentPhase()) { memcpy(solution_,info->currentSolution(), (numberRows_+numberColumns_)*sizeof(double)); } gutsOfSolution(NULL,NULL); if (info->currentPhase()) { memcpy(solution_,info->currentSolution(), (numberRows_+numberColumns_)*sizeof(double)); nonLinearCost_->checkInfeasibilities(primalTolerance_); } checkComplementarity(info,rowArray_[2],rowArray_[3]); // Double check reduced costs if no action if (progress->lastIterationNumber(0)==numberIterations_) { if (primalColumnPivot_->looksOptimal()) { numberDualInfeasibilities_ = 0; sumDualInfeasibilities_ = 0.0; } } // Check if looping int loop; if (type!=2) loop = progress->looping(); // saves iteration number else loop=-1; //loop=-1; if (loop>=0) { problemStatus_ = loop; //exit if in loop return ; } else if (loop<-1) { // something may have changed gutsOfSolution(NULL,NULL); checkComplementarity(info,rowArray_[2],rowArray_[3]); } // really for free variables in //if((progressFlag_&2)!=0) //problemStatus_=-1;; progressFlag_ = 0; //reset progress flag handler_->message(CLP_SIMPLEX_STATUS,messages_) <<numberIterations_<<objectiveValue(); handler_->printing(sumPrimalInfeasibilities_>0.0) <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_; handler_->printing(sumDualInfeasibilities_>0.0) <<sumDualInfeasibilities_<<numberDualInfeasibilities_; handler_->printing(numberDualInfeasibilitiesWithoutFree_ <numberDualInfeasibilities_) <<numberDualInfeasibilitiesWithoutFree_; handler_->message()<<CoinMessageEol; assert (primalFeasible()); // we may wish to say it is optimal even if infeasible bool alwaysOptimal = (specialOptions_&1)!=0; // give code benefit of doubt if (sumOfRelaxedDualInfeasibilities_ == 0.0&& sumOfRelaxedPrimalInfeasibilities_ == 0.0&&info->sequenceIn()<0) { // say optimal (with these bounds etc) numberDualInfeasibilities_ = 0; sumDualInfeasibilities_ = 0.0; numberPrimalInfeasibilities_ = 0; sumPrimalInfeasibilities_ = 0.0; } // had ||(type==3&&problemStatus_!=-5) -- ??? why ???? if (dualFeasible()||problemStatus_==-4) { if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) { //may need infeasiblity cost changed // we can see if we can construct a ray // make up a new objective double saveWeight = infeasibilityCost_; // save nonlinear cost as we are going to switch off costs ClpNonLinearCost * nonLinear = nonLinearCost_; // do twice to make sure Primal solution has settled // put non-basics to bounds in case tolerance moved // put back original costs createRim(4); // put all non-basic variables to bounds int iSequence; for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) { Status status = getStatus(iSequence); if (status==atLowerBound||status==isFixed) solution_[iSequence]=lower_[iSequence]; else if (status==atUpperBound) solution_[iSequence]=upper_[iSequence]; } nonLinearCost_->checkInfeasibilities(primalTolerance_); gutsOfSolution(NULL,NULL); nonLinearCost_->checkInfeasibilities(primalTolerance_); checkComplementarity(info,rowArray_[2],rowArray_[3]); infeasibilityCost_=1.0e100; // put back original costs createRim(4); nonLinearCost_->checkInfeasibilities(primalTolerance_); // may have fixed infeasibilities - double check if (nonLinearCost_->numberInfeasibilities()==0) { // carry on problemStatus_ = -1; infeasibilityCost_=saveWeight; nonLinearCost_->checkInfeasibilities(primalTolerance_); checkComplementarity(info,rowArray_[2],rowArray_[3]); } else { nonLinearCost_=NULL; // scale int i; for (i=0;i<numberRows_+numberColumns_;i++) cost_[i] *= 1.0e-100; gutsOfSolution(NULL,NULL); nonLinearCost_=nonLinear; infeasibilityCost_=saveWeight; if ((infeasibilityCost_>=1.0e18|| numberDualInfeasibilities_==0)&&perturbation_==101) { unPerturb(); // stop any further perturbation numberDualInfeasibilities_=1; // carry on } if (infeasibilityCost_>=1.0e20|| numberDualInfeasibilities_==0) { // we are infeasible - use as ray delete [] ray_; ray_ = new double [numberRows_]; memcpy(ray_,dual_,numberRows_*sizeof(double)); // and get feasible duals infeasibilityCost_=0.0; createRim(4); // put all non-basic variables to bounds int iSequence; for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) { Status status = getStatus(iSequence); if (status==atLowerBound||status==isFixed) solution_[iSequence]=lower_[iSequence]; else if (status==atUpperBound) solution_[iSequence]=upper_[iSequence]; } nonLinearCost_->checkInfeasibilities(primalTolerance_); gutsOfSolution(NULL,NULL); checkComplementarity(info,rowArray_[2],rowArray_[3]); // so will exit infeasibilityCost_=1.0e30; // reset infeasibilities sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();; numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities(); } if (infeasibilityCost_<1.0e20) { infeasibilityCost_ *= 5.0; changeMade_++; // say change made handler_->message(CLP_PRIMAL_WEIGHT,messages_) <<infeasibilityCost_ <<CoinMessageEol; // put back original costs and then check createRim(4); nonLinearCost_->checkInfeasibilities(0.0); gutsOfSolution(NULL,NULL); checkComplementarity(info,rowArray_[2],rowArray_[3]); problemStatus_=-1; //continue } else { // say infeasible problemStatus_ = 1; } } } else { // may be optimal if (perturbation_==101) { unPerturb(); // stop any further perturbation lastCleaned=-1; // carry on } if ( lastCleaned!=numberIterations_||unflag()) { handler_->message(CLP_PRIMAL_OPTIMAL,messages_) <<primalTolerance_ <<CoinMessageEol; if (numberTimesOptimal_<4) { numberTimesOptimal_++; changeMade_++; // say change made if (numberTimesOptimal_==1) { // better to have small tolerance even if slower factorization_->zeroTolerance(1.0e-15); } lastCleaned=numberIterations_; if (primalTolerance_!=dblParam_[ClpPrimalTolerance]) handler_->message(CLP_PRIMAL_ORIGINAL,messages_) <<CoinMessageEol; double oldTolerance = primalTolerance_; primalTolerance_=dblParam_[ClpPrimalTolerance]; // put back original costs and then check createRim(4); // put all non-basic variables to bounds int iSequence; for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) { Status status = getStatus(iSequence); if (status==atLowerBound||status==isFixed) solution_[iSequence]=lower_[iSequence]; else if (status==atUpperBound) solution_[iSequence]=upper_[iSequence]; } nonLinearCost_->checkInfeasibilities(oldTolerance); gutsOfSolution(NULL,NULL); checkComplementarity(info,rowArray_[2],rowArray_[3]); problemStatus_ = -1; } else { problemStatus_=0; // optimal if (lastCleaned<numberIterations_) { handler_->message(CLP_SIMPLEX_GIVINGUP,messages_) <<CoinMessageEol; } } } else { problemStatus_=0; // optimal } } } else { // see if looks unbounded if (problemStatus_==-5) { if (nonLinearCost_->numberInfeasibilities()) { if (infeasibilityCost_>1.0e18&&perturbation_==101) { // back off weight infeasibilityCost_ = 1.0e13; unPerturb(); // stop any further perturbation } //we need infeasiblity cost changed if (infeasibilityCost_<1.0e20) { infeasibilityCost_ *= 5.0; changeMade_++; // say change made handler_->message(CLP_PRIMAL_WEIGHT,messages_) <<infeasibilityCost_ <<CoinMessageEol; // put back original costs and then check createRim(4); gutsOfSolution(NULL,NULL); checkComplementarity(info,rowArray_[2],rowArray_[3]); problemStatus_=-1; //continue } else { // say unbounded problemStatus_ = 2; } } else { // say unbounded problemStatus_ = 2; } } else { if(type==3&&problemStatus_!=-5) unflag(); // odd // carry on problemStatus_ = -1; } } } else { // don't update solution problemStatus_=-1; } if (type==0||type==1) { if (type!=1||!saveStatus_) { // create save arrays delete [] saveStatus_; delete [] savedSolution_; saveStatus_ = new unsigned char [numberRows_+numberColumns_]; savedSolution_ = new double [numberRows_+numberColumns_]; } // save arrays memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char)); memcpy(savedSolution_+numberColumns_ ,rowActivityWork_, numberRows_*sizeof(double)); memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double)); info->saveStatus(); } // restore weights (if saved) - also recompute infeasibility list if (tentativeStatus>-3) primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4); else primalColumnPivot_->saveWeights(this,3); if (problemStatus_<0&&!changeMade_) { problemStatus_=4; // unknown } if (saveThreshold) { // use default at present factorization_->sparseThreshold(0); factorization_->goSparse(); } lastGoodIteration_ = numberIterations_; }

int ClpSimplexPrimalQuadratic::whileIterating (  ClpQuadraticInfo *  info  )

Main part. phase - 0 normal, 1 getting complementary solution, 2 getting basic solution. Definition at line 733 of file ClpSimplexPrimalQuadratic.cpp. References ClpQuadraticInfo::basicRow(), CoinIndexedVector::checkClear(), checkComplementarity(), CoinIndexedVector::clear(), createDjs(), ClpQuadraticInfo::createGradient(), ClpQuadraticInfo::crucialSj(), ClpQuadraticInfo::currentPhase(), CoinIndexedVector::denseVector(), ClpPackedMatrix::getElements(), CoinIndexedVector::getIndices(), ClpPackedMatrix::getIndices(), CoinIndexedVector::getNumElements(), ClpPackedMatrix::getVectorStarts(), ClpQuadraticInfo::gradient(), ClpSimplex::housekeeping(), ClpQuadraticInfo::impliedSj(), ClpQuadraticInfo::infeasCost(), ClpSimplex::isColumn(), CoinMessageHandler::logLevel(), CoinMessageHandler::message(), ClpNonLinearCost::numberInfeasibilities(), ClpQuadraticInfo::numberQuadraticRows(), ClpQuadraticInfo::numberXColumns(), ClpQuadraticInfo::numberXRows(), ClpSimplexPrimal::primalColumn(), ClpSimplexPrimal::primalRay(), primalRow(), ClpQuadraticInfo::sequenceIn(), ClpSimplex::sequenceWithin(), ClpQuadraticInfo::setCrucialSj(), ClpSimplex::setFlagged(), ClpNonLinearCost::setOne(), ClpNonLinearCost::setOneOutgoing(), ClpQuadraticInfo::setSequenceIn(), ClpSimplex::unpack(), ClpSimplexPrimal::updatePrimalsInPrimal(), and ClpPrimalColumnPivot::updateWeights(). Referenced by primalQuadratic2(). { checkComplementarity (info,rowArray_[0],rowArray_[1]); int crucialSj = info->crucialSj(); if (info->crucialSj()>=0) { printf("after inv %d Sj value %g\n",crucialSj,solution_[info->crucialSj()]); } int returnCode=-1; int phase = info->currentPhase(); double saveObjective = objectiveValue_; int numberXColumns = info->numberXColumns(); int numberXRows = info->numberXRows(); int numberQuadraticRows = info->numberQuadraticRows(); // Make list of implied sj // And backward pointers to basic variables { int * impliedSj = info->impliedSj(); int i; for (i=numberRows_;i<numberColumns_;i++) { if (getColumnStatus(i)==basic) impliedSj[i-numberRows_]=i; else impliedSj[i-numberRows_]=-1; } int * basicRow = info->basicRow(); for (i=0;i<numberRows_+numberColumns_;i++) basicRow[i]=-1; for (i=0;i<numberRows_;i++) basicRow[pivotVariable_[i]]=i; } int nextSequenceIn=info->sequenceIn(); int oldSequenceIn=nextSequenceIn; int saveSequenceIn = sequenceIn_; // status stays at -1 while iterating, >=0 finished, -2 to invert // status -3 to go to top without an invert while (problemStatus_==-1) { if (info->crucialSj()<0&&factorization_->pivots()>=10) { returnCode =-2; // refactorize break; } #ifdef CLP_DEBUG { int i; // not [1] as has information for (i=0;i<4;i++) { if (i!=1) rowArray_[i]->checkClear(); } for (i=0;i<2;i++) { columnArray_[i]->checkClear(); } } #endif // choose column to come in // can use pivotRow_ to update weights // pass in list of cost changes so can do row updates (rowArray_[1]) // NOTE rowArray_[0] is used by computeDuals which is a // slow way of getting duals but might be used // Initially Dantzig and look at s variables // Only do if one not already chosen int cleanupIteration; if (nextSequenceIn<0) { cleanupIteration=0; if (phase==2) { // values pass // get next int iColumn; int iStart = oldSequenceIn+1; for (iColumn=iStart;iColumn<numberXColumns;iColumn++) { double value = solution_[iColumn]; double lower = lower_[iColumn]; double upper = upper_[iColumn]; if (value>lower+primalTolerance_&&value<upper-primalTolerance_) { if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) { if (fabs(dj_[iColumn])>10.0*dualTolerance_|| getColumnStatus(iColumn+numberXColumns+numberQuadraticRows)==basic) { nextSequenceIn=iColumn; break; } } } } if (nextSequenceIn<0) { iStart=max(iStart,numberColumns_); int numberXRows = info->numberXRows(); for (iColumn=iStart;iColumn<numberColumns_+numberXRows; iColumn++) { double value = solution_[iColumn]; double lower = lower_[iColumn]; double upper = upper_[iColumn]; if (value>lower+primalTolerance_&&value<upper-primalTolerance_) { if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)&&fabs(dj_[iColumn])>10.0*dualTolerance_) { nextSequenceIn=iColumn; break; } } } } oldSequenceIn=nextSequenceIn; saveSequenceIn=sequenceIn_; sequenceIn_ = nextSequenceIn; } else { saveSequenceIn=sequenceIn_; createDjs(info,rowArray_[3],rowArray_[2]); dualIn_=0.0; // so no updates rowArray_[1]->clear(); primalColumn(rowArray_[1],rowArray_[2],rowArray_[3], columnArray_[0],columnArray_[1]); } } else { saveSequenceIn=sequenceIn_; sequenceIn_ = nextSequenceIn; if (phase==2) { if ((sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)&&info->crucialSj()<0) cleanupIteration=0; else cleanupIteration=1; } else { cleanupIteration=1; } } pivotRow_=-1; sequenceOut_=-1; rowArray_[1]->clear(); if (sequenceIn_>=0) { nextSequenceIn=-1; // we found a pivot column int chosen = sequenceIn_; // do second half of iteration while (chosen>=0) { int saveSequenceInInfo=chosen; int saveCrucialSjInfo=info->crucialSj(); rowArray_[1]->clear(); checkComplementarity (info,rowArray_[3],rowArray_[1]); printf("True objective is %g, infeas cost %g, sum %g\n", objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost()); objectiveValue_=saveObjective; returnCode=-1; pivotRow_=-1; sequenceOut_=-1; // we found a pivot column // update the incoming column sequenceIn_=chosen; chosen=-1; unpack(rowArray_[1]); // compute dj in case linear { info->createGradient(this); double * gradient = info->gradient(); dualIn_=gradient[sequenceIn_]; int j; const double * element = rowArray_[1]->denseVector(); int numberNonZero=rowArray_[1]->getNumElements(); const int * whichRow = rowArray_[1]->getIndices(); bool packed = rowArray_[1]->packedMode(); if (packed) { for (j=0;j<numberNonZero; j++) { int iRow = whichRow[j]; dualIn_ -= element[j]*dual_[iRow]; } } else { for (j=0;j<numberNonZero; j++) { int iRow = whichRow[j]; dualIn_ -= element[iRow]*dual_[iRow]; } } } if ((!cleanupIteration&&sequenceIn_<numberXColumns)|| (cleanupIteration&&info->crucialSj()<numberColumns_&&info->crucialSj()>=numberRows_)) { int iSequence; if (!cleanupIteration) iSequence=sequenceIn_; else iSequence=info->crucialSj()-numberRows_; // may need to re-do row of basis int * impliedSj = info->impliedSj(); if (impliedSj[iSequence]>=0) { // mark as valid impliedSj[iSequence]=-1; ClpPackedMatrix* rowCopy = dynamic_cast< ClpPackedMatrix*>(rowCopy_); assert(rowCopy); const int * column = rowCopy->getIndices(); const CoinBigIndex * rowStart = rowCopy->getVectorStarts(); const double * rowElement = rowCopy->getElements(); int iRow=iSequence+numberXRows; int i; int * index = rowArray_[2]->getIndices(); double * element = rowArray_[2]->denseVector(); int n=0; int * basicRow = info->basicRow(); int * permute = factorization_->pivotColumn(); for (i=rowStart[iRow];i<rowStart[iRow+1];i++) { int iColumn = column[i]; iColumn = basicRow[iColumn]; if (iColumn>=0&&iColumn!=iRow) { index[n]=permute[iColumn]; element[n++]=rowElement[i]; } } factorization_->replaceRow(permute[iRow],n,index,element); memset(element,0,n*sizeof(double)); } } // Take out elements in implied Sj rows if (1) { int n=rowArray_[1]->getNumElements(); int * index = rowArray_[1]->getIndices(); double * element = rowArray_[1]->denseVector(); int i; int * impliedSj = info->impliedSj(); int n2=0; for (i=0;i<n;i++) { int iRow=index[i]-numberXRows; if (iRow<0||impliedSj[iRow]<0) { index[n2++]=iRow+numberXRows; } else { element[iRow+numberXRows]=0.0; } } rowArray_[1]->setNumElements(n2); } factorization_->updateColumnFT(rowArray_[2],rowArray_[1]); if (cleanupIteration) { // move back to a version of primalColumn? valueIn_=solution_[sequenceIn_]; // should keep pivot row of crucialSj as well (speed) int iSjRow=-1; if (crucialSj>=0) { double * work=rowArray_[1]->denseVector(); int number=rowArray_[1]->getNumElements(); int * which=rowArray_[1]->getIndices(); double tj = 0.0; int iIndex; for (iIndex=0;iIndex<number;iIndex++) { int iRow = which[iIndex]; double alpha = work[iRow]; int iPivot=pivotVariable_[iRow]; if (iPivot==crucialSj) { tj = alpha; iSjRow = iRow; double d2 = solution_[crucialSj]/tj; if (fabs(solution_[crucialSj])<1.0e-8) printf("zero crucialSj - pivot out - get right way\n"); // see which way to go if (d2>0) dj_[sequenceIn_]= -1.0; else dj_[sequenceIn_]= 1.0; break; } } if (!tj) { printf("trouble\n"); assert (sequenceIn_>numberXColumns&&sequenceIn_<numberColumns_); dj_[sequenceIn_]=solution_[sequenceIn_]; //assert(tj); } } dualIn_=dj_[sequenceIn_]; lowerIn_=lower_[sequenceIn_]; upperIn_=upper_[sequenceIn_]; if (dualIn_>0.0) directionIn_ = -1; else directionIn_ = 1; } else { // not clean up lowerIn_=lower_[sequenceIn_]; upperIn_=upper_[sequenceIn_]; dualIn_=dj_[sequenceIn_]; valueIn_=solution_[sequenceIn_]; if (dualIn_>0.0) directionIn_ = -1; else directionIn_ = 1; if (sequenceIn_<numberColumns_) { // Set dj as value of slack crucialSj = sequenceIn_+ numberQuadraticRows+numberXColumns; // sj which should go to 0.0 } else { // Set dj as value of pi int iRow = sequenceIn_-numberColumns_; crucialSj = iRow+numberXColumns; // pi which should go to 0.0 } if (crucialSj>=0&&getColumnStatus(crucialSj)!=basic) crucialSj=-1; info->setCrucialSj(crucialSj); } double oldSj=1.0e30; if (info->crucialSj()>=0&&cleanupIteration) oldSj= solution_[info->crucialSj()]; // save reduced cost //double saveDj = dualIn_; // do ratio test and re-compute dj // Note second parameter long enough for columns int result=primalRow(rowArray_[1],rowArray_[3], rowArray_[2],rowArray_[0], info, cleanupIteration); if (pivotRow_==-1&&(phase==2||cleanupIteration)&&fabs(dualIn_)<1.0e-3) { // try other way dualIn_=-dualIn_; directionIn_=-directionIn_; if (info->crucialSj()>=0) setColumnStatus(info->crucialSj(),basic); result=primalRow(rowArray_[1],rowArray_[3], rowArray_[2],rowArray_[0], info, cleanupIteration); } saveObjective = objectiveValue_; if (pivotRow_>=0) { // If sj out AND not gone out since last invert then add back // if stable replace in basis int updateStatus = 0; if (result<20) { double saveCheck = factorization_->getAccuracyCheck(); if (cleanupIteration) factorization_->relaxAccuracyCheck(1.0e3*saveCheck); updateStatus=factorization_->replaceColumn(rowArray_[2], pivotRow_, alpha_); factorization_->relaxAccuracyCheck(saveCheck); } if (result>=10) { updateStatus = max(updateStatus,result/10); result = result%10; if (updateStatus>1) { alpha_=0.0; // force re-factorization info->setSequenceIn(sequenceIn_); } } // if no pivots, bad update but reasonable alpha - take and invert if (updateStatus==2&& lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5) updateStatus=4; if (updateStatus==1||updateStatus==4) { // slight error if (factorization_->pivots()>5||updateStatus==4) { returnCode=-3; } } else if (updateStatus==2) { // major error // Reset sequenceIn_ // sequenceIn_=saveSequenceIn; nextSequenceIn=saveSequenceInInfo; info->setCrucialSj(saveCrucialSjInfo); if (saveCrucialSjInfo<0&&!phase) nextSequenceIn=-1; pivotRow_=-1; // better to have small tolerance even if slower factorization_->zeroTolerance(1.0e-15); int maxFactor = factorization_->maximumPivots(); if (maxFactor>10) { if (forceFactorization_<0) forceFactorization_= maxFactor; forceFactorization_ = max (1,(forceFactorization_>>1)); } // later we may need to unwind more e.g. fake bounds if(lastGoodIteration_ != numberIterations_||factorization_->pivots()) { rowArray_[1]->clear(); pivotRow_=-1; returnCode=-4; // retry on return if (info->crucialSj()>=0) nextSequenceIn = sequenceIn_; break; } else { // need to reject something char x = isColumn(sequenceIn_) ? 'C' :'R'; handler_->message(CLP_SIMPLEX_FLAG,messages_) <<x<<sequenceWithin(sequenceIn_) <<CoinMessageEol; setFlagged(sequenceIn_); lastBadIteration_ = numberIterations_; // say be more cautious rowArray_[1]->clear(); pivotRow_=-1; returnCode = -5; break; } } else if (updateStatus==3) { // out of memory // increase space if not many iterations if (factorization_->pivots()< 0.5*factorization_->maximumPivots()&& factorization_->pivots()<200) factorization_->areaFactor( factorization_->areaFactor() * 1.1); returnCode =-2; // factorize now } // here do part of steepest - ready for next iteration primalColumnPivot_->updateWeights(rowArray_[1]); } else { if (pivotRow_==-1) { // no outgoing row is valid rowArray_[0]->clear(); if (!factorization_->pivots()) { returnCode = 2; //say looks unbounded // do ray primalRay(rowArray_[1]); } else { returnCode = 4; //say looks unbounded but has iterated } break; } else { // flipping from bound to bound } } // update primal solution double objectiveChange=0.0; // Cost on pivot row may change - may need to change dualIn double oldCost=0.0; if (pivotRow_>=0) oldCost = cost(pivotVariable_[pivotRow_]); // rowArray_[1] is not empty - used to update djs updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,1); if (pivotRow_>=0) dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_])); double oldValue = valueIn_; if (directionIn_==-1) { // as if from upper bound if (sequenceIn_!=sequenceOut_) { // variable becoming basic valueIn_ -= fabs(theta_); } else { valueIn_=lowerIn_; } } else { // as if from lower bound if (sequenceIn_!=sequenceOut_) { // variable becoming basic valueIn_ += fabs(theta_); } else { valueIn_=upperIn_; } } objectiveChange += dualIn_*(valueIn_-oldValue); // outgoing if (sequenceIn_!=sequenceOut_) { if (directionOut_>0) { valueOut_ = lowerOut_; } else { valueOut_ = upperOut_; } double lowerValue = lower_[sequenceOut_]; double upperValue = upper_[sequenceOut_]; if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) { // really free but going to zero lowerValue=0.0; upperValue=0.0; } assert(valueOut_>=lowerValue-primalTolerance_&& valueOut_<=upperValue+primalTolerance_); // may not be exactly at bound and bounds may have changed #if 1 // Make sure outgoing looks feasible double useValue=valueOut_; if (valueOut_<=lowerValue+primalTolerance_) { directionOut_=1; useValue= lower_[sequenceOut_]; } else if (valueOut_>=upperValue-primalTolerance_) { directionOut_=-1; useValue= upper_[sequenceOut_]; } else { printf("*** variable wandered off bound %g %g %g!\n", lowerValue,valueOut_,upperValue); if (valueOut_-lowerValue<=upperValue-valueOut_) { useValue= lower_[sequenceOut_]; valueOut_=lowerValue; directionOut_=1; } else { useValue= upper_[sequenceOut_]; valueOut_=upperValue; directionOut_=-1; } } solution_[sequenceOut_]=valueOut_; nonLinearCost_->setOne(sequenceOut_,useValue); #else directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_); solution_[sequenceOut_]=valueOut_; #endif } // change cost and bounds on incoming if primal if (sequenceIn_==sequenceOut_&&(lowerOut_!=lowerIn_||upperOut_!=upperIn_)) { // linear variable superbasic setStatus(sequenceOut_,superBasic); } nonLinearCost_->setOne(sequenceIn_,valueIn_); int whatNext=housekeeping(objectiveChange); // and again as housekeeping may have changed if (sequenceIn_==sequenceOut_&&(lowerOut_!=lowerIn_||upperOut_!=upperIn_)) { // linear variable superbasic setStatus(sequenceOut_,superBasic); } // And update backward pointers to basic variables if (sequenceIn_!=sequenceOut_) { int * basicRow = info->basicRow(); basicRow[sequenceOut_]=-1; basicRow[sequenceIn_]=pivotRow_; } if (info->crucialSj()>=0) { assert(fabs(oldSj)>= fabs(solution_[info->crucialSj()])); printf("%d Sj value went from %g to %g\n",crucialSj,oldSj,solution_[info->crucialSj()]); } checkComplementarity (info,rowArray_[2],rowArray_[3]); printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n", objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost()); if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) { // really free but going to zero setStatus(sequenceOut_,isFree); if (sequenceOut_==info->crucialSj()) info->setCrucialSj(-1); } else if (info->crucialSj()>=0) { // Just possible crucialSj still in but original out again int crucialSj=info->crucialSj(); if (crucialSj>=numberXColumns+numberQuadraticRows) { if (sequenceOut_==crucialSj-numberXColumns-numberQuadraticRows) info->setCrucialSj(-1); } else { if (sequenceOut_-numberColumns_==crucialSj-numberXColumns) info->setCrucialSj(-1); } } if (info->crucialSj()<0) result=0; if (whatNext==1) { returnCode =-2; // refactorize } else if (whatNext==2) { // maximum iterations or equivalent returnCode=3; } else if(numberIterations_ == lastGoodIteration_ + 2 * factorization_->maximumPivots()) { // done a lot of flips - be safe returnCode =-2; // refactorize } // may need to go round again cleanupIteration = 1; // may not be correct on second time if (result&&(returnCode==-1||returnCode==-2||returnCode==-3)) { int crucialSj=info->crucialSj(); if (sequenceOut_<numberXColumns) { chosen =sequenceOut_ + numberQuadraticRows + numberXColumns; // sj variable } else if (sequenceOut_>=numberColumns_) { // Does this mean we can change pi int iRow = sequenceOut_-numberColumns_; if (iRow<numberXRows) { chosen=iRow+numberXColumns; } else { printf("Row %d is in column part\n",iRow); abort(); } } else if (sequenceOut_<numberQuadraticRows+numberXColumns) { // pi out chosen = sequenceOut_-numberXColumns+numberColumns_; } else { // sj out chosen = sequenceOut_-numberQuadraticRows-numberXColumns; } // But double check as original incoming might have gone out if (chosen==crucialSj) { chosen=-1; break; // means original has gone out } rowArray_[0]->clear(); columnArray_[0]->clear(); if (returnCode==-2) { // factorization nextSequenceIn = chosen; break; } } else { break; } } if (returnCode<-1&&returnCode>-5) { problemStatus_=-2; // } else if (returnCode==-5) { // something flagged - continue; } else if (returnCode==2) { problemStatus_=-5; // looks unbounded } else if (returnCode==4) { problemStatus_=-2; // looks unbounded but has iterated } else if (returnCode!=-1) { assert(returnCode==3); problemStatus_=3; } } else { // no pivot column #ifdef CLP_DEBUG if (handler_->logLevel()&32) printf("** no column pivot\n"); #endif if (nonLinearCost_->numberInfeasibilities()) problemStatus_=-4; // might be infeasible returnCode=0; break; } } info->setSequenceIn(nextSequenceIn); return returnCode; }

---

The documentation for this class was generated from the following files:

• ClpSimplexPrimalQuadratic.hpp
• ClpSimplexPrimalQuadratic.cpp

---