/* $Id: ClpFactorization.cpp 1903 2013-01-03 22:49:30Z stefan $ */ // Copyright (C) 2002, International Business Machines // Corporation and others. All Rights Reserved. // This code is licensed under the terms of the Eclipse Public License (EPL). #include "CoinPragma.hpp" #include "ClpFactorization.hpp" #ifndef SLIM_CLP #include "ClpQuadraticObjective.hpp" #endif #include "CoinHelperFunctions.hpp" #include "CoinIndexedVector.hpp" #include "ClpSimplex.hpp" #include "ClpSimplexDual.hpp" #include "ClpMatrixBase.hpp" #ifndef SLIM_CLP #include "ClpNetworkBasis.hpp" #include "ClpNetworkMatrix.hpp" //#define CHECK_NETWORK #ifdef CHECK_NETWORK const static bool doCheck = true; #else const static bool doCheck = false; #endif #endif //#define CLP_FACTORIZATION_INSTRUMENT #ifdef CLP_FACTORIZATION_INSTRUMENT #include "CoinTime.hpp" double factorization_instrument(int type) { static int times[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; static double startTime = 0.0; static double totalTimes [10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; if (type < 0) { assert (!startTime); startTime = CoinCpuTime(); return 0.0; } else if (type > 0) { times[type]++; double difference = CoinCpuTime() - startTime; totalTimes[type] += difference; startTime = 0.0; return difference; } else { // report const char * types[10] = { "", "fac=rhs_etc", "factorize", "replace", "update_FT", "update", "update_transpose", "gosparse", "getWeights!", "update2_FT" }; double total = 0.0; for (int i = 1; i < 10; i++) { if (times[i]) { printf("%s was called %d times taking %g seconds\n", types[i], times[i], totalTimes[i]); total += totalTimes[i]; times[i] = 0; totalTimes[i] = 0.0; } } return total; } } #endif //############################################################################# // Constructors / Destructor / Assignment //############################################################################# #ifndef CLP_MULTIPLE_FACTORIZATIONS //------------------------------------------------------------------- // Default Constructor //------------------------------------------------------------------- ClpFactorization::ClpFactorization () : CoinFactorization() { #ifndef SLIM_CLP networkBasis_ = NULL; #endif } //------------------------------------------------------------------- // Copy constructor //------------------------------------------------------------------- ClpFactorization::ClpFactorization (const ClpFactorization & rhs, int dummyDenseIfSmaller) : CoinFactorization(rhs) { #ifndef SLIM_CLP if (rhs.networkBasis_) networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_)); else networkBasis_ = NULL; #endif } ClpFactorization::ClpFactorization (const CoinFactorization & rhs) : CoinFactorization(rhs) { #ifndef SLIM_CLP networkBasis_ = NULL; #endif } //------------------------------------------------------------------- // Destructor //------------------------------------------------------------------- ClpFactorization::~ClpFactorization () { #ifndef SLIM_CLP delete networkBasis_; #endif } //---------------------------------------------------------------- // Assignment operator //------------------------------------------------------------------- ClpFactorization & ClpFactorization::operator=(const ClpFactorization& rhs) { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif if (this != &rhs) { CoinFactorization::operator=(rhs); #ifndef SLIM_CLP delete networkBasis_; if (rhs.networkBasis_) networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_)); else networkBasis_ = NULL; #endif } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(1); #endif return *this; } #if 0 static unsigned int saveList[10000]; int numberSave = -1; inline bool isDense(int i) { return ((saveList[i>>5] >> (i & 31)) & 1) != 0; } inline void setDense(int i) { unsigned int & value = saveList[i>>5]; int bit = i & 31; value |= (1 << bit); } #endif int ClpFactorization::factorize ( ClpSimplex * model, int solveType, bool valuesPass) { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif ClpMatrixBase * matrix = model->clpMatrix(); int numberRows = model->numberRows(); int numberColumns = model->numberColumns(); if (!numberRows) return 0; // If too many compressions increase area if (numberPivots_ > 1 && numberCompressions_ * 10 > numberPivots_ + 10) { areaFactor_ *= 1.1; } //int numberPivots=numberPivots_; #if 0 if (model->algorithm() > 0) numberSave = -1; #endif #ifndef SLIM_CLP if (!networkBasis_ || doCheck) { #endif status_ = -99; int * pivotVariable = model->pivotVariable(); //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */ while (status_ < -98) { int i; int numberBasic = 0; int numberRowBasic; // Move pivot variables across if they look good int * pivotTemp = model->rowArray(0)->getIndices(); assert (!model->rowArray(0)->getNumElements()); if (!matrix->rhsOffset(model)) { #if 0 if (numberSave > 0) { int nStill = 0; int nAtBound = 0; int nZeroDual = 0; CoinIndexedVector * array = model->rowArray(3); CoinIndexedVector * objArray = model->columnArray(1); array->clear(); objArray->clear(); double * cost = model->costRegion(); double tolerance = model->primalTolerance(); double offset = 0.0; for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot < numberColumns && isDense(iPivot)) { if (model->getColumnStatus(iPivot) == ClpSimplex::basic) { nStill++; double value = model->solutionRegion()[iPivot]; double dual = model->dualRowSolution()[i]; double lower = model->lowerRegion()[iPivot]; double upper = model->upperRegion()[iPivot]; ClpSimplex::Status status; if (fabs(value - lower) < tolerance) { status = ClpSimplex::atLowerBound; nAtBound++; } else if (fabs(value - upper) < tolerance) { nAtBound++; status = ClpSimplex::atUpperBound; } else if (value > lower && value < upper) { status = ClpSimplex::superBasic; } else { status = ClpSimplex::basic; } if (status != ClpSimplex::basic) { if (model->getRowStatus(i) != ClpSimplex::basic) { model->setColumnStatus(iPivot, ClpSimplex::atLowerBound); model->setRowStatus(i, ClpSimplex::basic); pivotVariable[i] = i + numberColumns; model->dualRowSolution()[i] = 0.0; model->djRegion(0)[i] = 0.0; array->add(i, dual); offset += dual * model->solutionRegion(0)[i]; } } if (fabs(dual) < 1.0e-5) nZeroDual++; } } } printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n", numberSave, nStill, nAtBound, nZeroDual, offset); if (array->getNumElements()) { // modify costs model->clpMatrix()->transposeTimes(model, 1.0, array, model->columnArray(0), objArray); array->clear(); int n = objArray->getNumElements(); int * indices = objArray->getIndices(); double * elements = objArray->denseVector(); for (i = 0; i < n; i++) { int iColumn = indices[i]; cost[iColumn] -= elements[iColumn]; elements[iColumn] = 0.0; } objArray->setNumElements(0); } } #endif // Seems to prefer things in order so quickest // way is to go though like this for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) pivotTemp[numberBasic++] = i; } numberRowBasic = numberBasic; /* Put column basic variables into pivotVariable This is done by ClpMatrixBase to allow override for gub */ matrix->generalExpanded(model, 0, numberBasic); } else { // Long matrix - do a different way bool fullSearch = false; for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot >= numberColumns) { pivotTemp[numberBasic++] = iPivot - numberColumns; } } numberRowBasic = numberBasic; for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot < numberColumns) { if (iPivot >= 0) { pivotTemp[numberBasic++] = iPivot; } else { // not full basis fullSearch = true; break; } } } if (fullSearch) { // do slow way numberBasic = 0; for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) pivotTemp[numberBasic++] = i; } numberRowBasic = numberBasic; /* Put column basic variables into pivotVariable This is done by ClpMatrixBase to allow override for gub */ matrix->generalExpanded(model, 0, numberBasic); } } if (numberBasic > model->maximumBasic()) { #if 0 // ndef NDEBUG printf("%d basic - should only be %d\n", numberBasic, numberRows); #endif // Take out some numberBasic = numberRowBasic; for (int i = 0; i < numberColumns; i++) { if (model->getColumnStatus(i) == ClpSimplex::basic) { if (numberBasic < numberRows) numberBasic++; else model->setColumnStatus(i, ClpSimplex::superBasic); } } numberBasic = numberRowBasic; matrix->generalExpanded(model, 0, numberBasic); } #ifndef SLIM_CLP // see if matrix a network #ifndef NO_RTTI ClpNetworkMatrix* networkMatrix = dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix()); #else ClpNetworkMatrix* networkMatrix = NULL; if (model->clpMatrix()->type() == 11) networkMatrix = static_cast< ClpNetworkMatrix*>(model->clpMatrix()); #endif // If network - still allow ordinary factorization first time for laziness if (networkMatrix) biasLU_ = 0; // All to U if network //int saveMaximumPivots = maximumPivots(); delete networkBasis_; networkBasis_ = NULL; if (networkMatrix && !doCheck) maximumPivots(1); #endif //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR()); while (status_ == -99) { // maybe for speed will be better to leave as many regions as possible gutsOfDestructor(); gutsOfInitialize(2); CoinBigIndex numberElements = numberRowBasic; // compute how much in basis int i; // can change for gub int numberColumnBasic = numberBasic - numberRowBasic; numberElements += matrix->countBasis(model, pivotTemp + numberRowBasic, numberRowBasic, numberColumnBasic); // and recompute as network side say different if (model->numberIterations()) numberRowBasic = numberBasic - numberColumnBasic; numberElements = 3 * numberBasic + 3 * numberElements + 20000; #if 0 // If iteration not zero then may be compressed getAreas ( !model->numberIterations() ? numberRows : numberBasic, numberRowBasic + numberColumnBasic, numberElements, 2 * numberElements ); #else getAreas ( numberRows, numberRowBasic + numberColumnBasic, numberElements, 2 * numberElements ); #endif //fill // Fill in counts so we can skip part of preProcess int * numberInRow = numberInRow_.array(); int * numberInColumn = numberInColumn_.array(); CoinZeroN ( numberInRow, numberRows_ + 1 ); CoinZeroN ( numberInColumn, maximumColumnsExtra_ + 1 ); double * elementU = elementU_.array(); int * indexRowU = indexRowU_.array(); CoinBigIndex * startColumnU = startColumnU_.array(); for (i = 0; i < numberRowBasic; i++) { int iRow = pivotTemp[i]; indexRowU[i] = iRow; startColumnU[i] = i; elementU[i] = slackValue_; numberInRow[iRow] = 1; numberInColumn[i] = 1; } startColumnU[numberRowBasic] = numberRowBasic; // can change for gub so redo numberColumnBasic = numberBasic - numberRowBasic; matrix->fillBasis(model, pivotTemp + numberRowBasic, numberColumnBasic, indexRowU, startColumnU + numberRowBasic, numberInRow, numberInColumn + numberRowBasic, elementU); #if 0 { printf("%d row basic, %d column basic\n", numberRowBasic, numberColumnBasic); for (int i = 0; i < numberElements; i++) printf("row %d col %d value %g\n", indexRowU_.array()[i], indexColumnU_[i], elementU_.array()[i]); } #endif // recompute number basic numberBasic = numberRowBasic + numberColumnBasic; if (numberBasic) numberElements = startColumnU[numberBasic-1] + numberInColumn[numberBasic-1]; else numberElements = 0; lengthU_ = numberElements; //saveFactorization("dump.d"); if (biasLU_ >= 3 || numberRows_ != numberColumns_) preProcess ( 2 ); else preProcess ( 3 ); // no row copy factor ( ); if (status_ == -99) { // get more memory areaFactor(2.0 * areaFactor()); } else if (status_ == -1 && model->numberIterations() == 0 && denseThreshold_) { // Round again without dense denseThreshold_ = 0; status_ = -99; } } // If we get here status is 0 or -1 if (status_ == 0) { // We may need to tamper with order and redo - e.g. network with side int useNumberRows = numberRows; // **** we will also need to add test in dual steepest to do // as we do for network matrix->generalExpanded(model, 12, useNumberRows); const int * permuteBack = permuteBack_.array(); const int * back = pivotColumnBack_.array(); //int * pivotTemp = pivotColumn_.array(); //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp ); // Redo pivot order for (i = 0; i < numberRowBasic; i++) { int k = pivotTemp[i]; // so rowIsBasic[k] would be permuteBack[back[i]] pivotVariable[permuteBack[back[i]]] = k + numberColumns; } for (; i < useNumberRows; i++) { int k = pivotTemp[i]; // so rowIsBasic[k] would be permuteBack[back[i]] pivotVariable[permuteBack[back[i]]] = k; } #if 0 if (numberSave >= 0) { numberSave = numberDense_; memset(saveList, 0, ((numberRows_ + 31) >> 5)*sizeof(int)); for (i = numberRows_ - numberSave; i < numberRows_; i++) { int k = pivotTemp[pivotColumn_.array()[i]]; setDense(k); } } #endif // Set up permutation vector // these arrays start off as copies of permute // (and we could use permute_ instead of pivotColumn (not back though)) ClpDisjointCopyN ( permute_.array(), useNumberRows , pivotColumn_.array() ); ClpDisjointCopyN ( permuteBack_.array(), useNumberRows , pivotColumnBack_.array() ); #ifndef SLIM_CLP if (networkMatrix) { maximumPivots(CoinMax(2000, maximumPivots())); // redo arrays for (int iRow = 0; iRow < 4; iRow++) { int length = model->numberRows() + maximumPivots(); if (iRow == 3 || model->objectiveAsObject()->type() > 1) length += model->numberColumns(); model->rowArray(iRow)->reserve(length); } // create network factorization if (doCheck) delete networkBasis_; // temp networkBasis_ = new ClpNetworkBasis(model, numberRows_, pivotRegion_.array(), permuteBack_.array(), startColumnU_.array(), numberInColumn_.array(), indexRowU_.array(), elementU_.array()); // kill off arrays in ordinary factorization if (!doCheck) { gutsOfDestructor(); // but make sure numberRows_ set numberRows_ = model->numberRows(); status_ = 0; #if 0 // but put back permute arrays so odd things will work int numberRows = model->numberRows(); pivotColumnBack_ = new int [numberRows]; permute_ = new int [numberRows]; int i; for (i = 0; i < numberRows; i++) { pivotColumnBack_[i] = i; permute_[i] = i; } #endif } } else { #endif // See if worth going sparse and when if (numberFtranCounts_ > 100) { ftranCountInput_ = CoinMax(ftranCountInput_, 1.0); ftranAverageAfterL_ = CoinMax(ftranCountAfterL_ / ftranCountInput_, 1.0); ftranAverageAfterR_ = CoinMax(ftranCountAfterR_ / ftranCountAfterL_, 1.0); ftranAverageAfterU_ = CoinMax(ftranCountAfterU_ / ftranCountAfterR_, 1.0); if (btranCountInput_ && btranCountAfterU_ && btranCountAfterR_) { btranAverageAfterU_ = CoinMax(btranCountAfterU_ / btranCountInput_, 1.0); btranAverageAfterR_ = CoinMax(btranCountAfterR_ / btranCountAfterU_, 1.0); btranAverageAfterL_ = CoinMax(btranCountAfterL_ / btranCountAfterR_, 1.0); } else { // we have not done any useful btrans (values pass?) btranAverageAfterU_ = 1.0; btranAverageAfterR_ = 1.0; btranAverageAfterL_ = 1.0; } } // scale back ftranCountInput_ *= 0.8; ftranCountAfterL_ *= 0.8; ftranCountAfterR_ *= 0.8; ftranCountAfterU_ *= 0.8; btranCountInput_ *= 0.8; btranCountAfterU_ *= 0.8; btranCountAfterR_ *= 0.8; btranCountAfterL_ *= 0.8; #ifndef SLIM_CLP } #endif } else if (status_ == -1 && (solveType == 0 || solveType == 2)) { // This needs redoing as it was merged coding - does not need array int numberTotal = numberRows + numberColumns; int * isBasic = new int [numberTotal]; int * rowIsBasic = isBasic + numberColumns; int * columnIsBasic = isBasic; for (i = 0; i < numberTotal; i++) isBasic[i] = -1; for (i = 0; i < numberRowBasic; i++) { int iRow = pivotTemp[i]; rowIsBasic[iRow] = 1; } for (; i < numberBasic; i++) { int iColumn = pivotTemp[i]; columnIsBasic[iColumn] = 1; } numberBasic = 0; for (i = 0; i < numberRows; i++) pivotVariable[i] = -1; // mark as basic or non basic const int * pivotColumn = pivotColumn_.array(); for (i = 0; i < numberRows; i++) { if (rowIsBasic[i] >= 0) { if (pivotColumn[numberBasic] >= 0) { rowIsBasic[i] = pivotColumn[numberBasic]; } else { rowIsBasic[i] = -1; model->setRowStatus(i, ClpSimplex::superBasic); } numberBasic++; } } for (i = 0; i < numberColumns; i++) { if (columnIsBasic[i] >= 0) { if (pivotColumn[numberBasic] >= 0) columnIsBasic[i] = pivotColumn[numberBasic]; else columnIsBasic[i] = -1; numberBasic++; } } // leave pivotVariable in useful form for cleaning basis int * pivotVariable = model->pivotVariable(); for (i = 0; i < numberRows; i++) { pivotVariable[i] = -1; } for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) { int iPivot = rowIsBasic[i]; if (iPivot >= 0) pivotVariable[iPivot] = i + numberColumns; } } for (i = 0; i < numberColumns; i++) { if (model->getColumnStatus(i) == ClpSimplex::basic) { int iPivot = columnIsBasic[i]; if (iPivot >= 0) pivotVariable[iPivot] = i; } } delete [] isBasic; double * columnLower = model->lowerRegion(); double * columnUpper = model->upperRegion(); double * columnActivity = model->solutionRegion(); double * rowLower = model->lowerRegion(0); double * rowUpper = model->upperRegion(0); double * rowActivity = model->solutionRegion(0); //redo basis - first take ALL columns out int iColumn; double largeValue = model->largeValue(); for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getColumnStatus(iColumn) == ClpSimplex::basic) { // take out if (!valuesPass) { double lower = columnLower[iColumn]; double upper = columnUpper[iColumn]; double value = columnActivity[iColumn]; if (lower > -largeValue || upper < largeValue) { if (fabs(value - lower) < fabs(value - upper)) { model->setColumnStatus(iColumn, ClpSimplex::atLowerBound); columnActivity[iColumn] = lower; } else { model->setColumnStatus(iColumn, ClpSimplex::atUpperBound); columnActivity[iColumn] = upper; } } else { model->setColumnStatus(iColumn, ClpSimplex::isFree); } } else { model->setColumnStatus(iColumn, ClpSimplex::superBasic); } } } int iRow; for (iRow = 0; iRow < numberRows; iRow++) { int iSequence = pivotVariable[iRow]; if (iSequence >= 0) { // basic if (iSequence >= numberColumns) { // slack in - leave //assert (iSequence-numberColumns==iRow); } else { assert(model->getRowStatus(iRow) != ClpSimplex::basic); // put back structural model->setColumnStatus(iSequence, ClpSimplex::basic); } } else { // put in slack model->setRowStatus(iRow, ClpSimplex::basic); } } // Put back any key variables for gub (status_ not touched) matrix->generalExpanded(model, 1, status_); // signal repeat status_ = -99; // set fixed if they are for (iRow = 0; iRow < numberRows; iRow++) { if (model->getRowStatus(iRow) != ClpSimplex::basic ) { if (rowLower[iRow] == rowUpper[iRow]) { rowActivity[iRow] = rowLower[iRow]; model->setRowStatus(iRow, ClpSimplex::isFixed); } } } for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getColumnStatus(iColumn) != ClpSimplex::basic ) { if (columnLower[iColumn] == columnUpper[iColumn]) { columnActivity[iColumn] = columnLower[iColumn]; model->setColumnStatus(iColumn, ClpSimplex::isFixed); } } } } } #ifndef SLIM_CLP } else { // network - fake factorization - do nothing status_ = 0; numberPivots_ = 0; } #endif #ifndef SLIM_CLP if (!status_) { // take out part if quadratic if (model->algorithm() == 2) { ClpObjective * obj = model->objectiveAsObject(); #ifndef NDEBUG ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj)); assert (quadraticObj); #else ClpQuadraticObjective * quadraticObj = (static_cast< ClpQuadraticObjective*>(obj)); #endif CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective(); int numberXColumns = quadratic->getNumCols(); assert (numberXColumns < numberColumns); int base = numberColumns - numberXColumns; int * which = new int [numberXColumns]; int * pivotVariable = model->pivotVariable(); int * permute = pivotColumn(); int i; int n = 0; for (i = 0; i < numberRows; i++) { int iSj = pivotVariable[i] - base; if (iSj >= 0 && iSj < numberXColumns) which[n++] = permute[i]; } if (n) emptyRows(n, which); delete [] which; } } #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(2); #endif return status_; } /* Replaces one Column to basis, returns 0=OK, 1=Probably OK, 2=singular, 3=no room If checkBeforeModifying is true will do all accuracy checks before modifying factorization. Whether to set this depends on speed considerations. You could just do this on first iteration after factorization and thereafter re-factorize partial update already in U */ int ClpFactorization::replaceColumn ( const ClpSimplex * model, CoinIndexedVector * regionSparse, CoinIndexedVector * tableauColumn, int pivotRow, double pivotCheck , bool checkBeforeModifying, double acceptablePivot) { int returnCode; #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif // see if FT if (doForrestTomlin_) { returnCode = CoinFactorization::replaceColumn(regionSparse, pivotRow, pivotCheck, checkBeforeModifying, acceptablePivot); } else { returnCode = CoinFactorization::replaceColumnPFI(tableauColumn, pivotRow, pivotCheck); // Note array } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(3); #endif #ifndef SLIM_CLP } else { if (doCheck) { returnCode = CoinFactorization::replaceColumn(regionSparse, pivotRow, pivotCheck, checkBeforeModifying, acceptablePivot); networkBasis_->replaceColumn(regionSparse, pivotRow); } else { // increase number of pivots numberPivots_++; returnCode = networkBasis_->replaceColumn(regionSparse, pivotRow); } } #endif return returnCode; } /* Updates one column (FTRAN) from region2 number returned is negative if no room region1 starts as zero and is zero at end */ int ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse, CoinIndexedVector * regionSparse2) { #ifdef CLP_DEBUG regionSparse->checkClear(); #endif if (!numberRows_) return 0; #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif collectStatistics_ = true; int returnCode = CoinFactorization::updateColumnFT(regionSparse, regionSparse2); collectStatistics_ = false; #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(4); #endif return returnCode; #ifndef SLIM_CLP } else { #ifdef CHECK_NETWORK CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2); double * check = new double[numberRows_]; int returnCode = CoinFactorization::updateColumnFT(regionSparse, regionSparse2); networkBasis_->updateColumn(regionSparse, save, -1); int i; double * array = regionSparse2->denseVector(); int * indices = regionSparse2->getIndices(); int n = regionSparse2->getNumElements(); memset(check, 0, numberRows_ * sizeof(double)); double * array2 = save->denseVector(); int * indices2 = save->getIndices(); int n2 = save->getNumElements(); assert (n == n2); if (save->packedMode()) { for (i = 0; i < n; i++) { check[indices[i]] = array[i]; } for (i = 0; i < n; i++) { double value2 = array2[i]; assert (check[indices2[i]] == value2); } } else { for (i = 0; i < numberRows_; i++) { double value1 = array[i]; double value2 = array2[i]; assert (value1 == value2); } } delete save; delete [] check; return returnCode; #else networkBasis_->updateColumn(regionSparse, regionSparse2, -1); return 1; #endif } #endif } /* Updates one column (FTRAN) from region2 number returned is negative if no room region1 starts as zero and is zero at end */ int ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse, CoinIndexedVector * regionSparse2, bool noPermute) const { #ifdef CLP_DEBUG if (!noPermute) regionSparse->checkClear(); #endif if (!numberRows_) return 0; #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif collectStatistics_ = true; int returnCode = CoinFactorization::updateColumn(regionSparse, regionSparse2, noPermute); collectStatistics_ = false; #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(5); #endif return returnCode; #ifndef SLIM_CLP } else { #ifdef CHECK_NETWORK CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2); double * check = new double[numberRows_]; int returnCode = CoinFactorization::updateColumn(regionSparse, regionSparse2, noPermute); networkBasis_->updateColumn(regionSparse, save, -1); int i; double * array = regionSparse2->denseVector(); int * indices = regionSparse2->getIndices(); int n = regionSparse2->getNumElements(); memset(check, 0, numberRows_ * sizeof(double)); double * array2 = save->denseVector(); int * indices2 = save->getIndices(); int n2 = save->getNumElements(); assert (n == n2); if (save->packedMode()) { for (i = 0; i < n; i++) { check[indices[i]] = array[i]; } for (i = 0; i < n; i++) { double value2 = array2[i]; assert (check[indices2[i]] == value2); } } else { for (i = 0; i < numberRows_; i++) { double value1 = array[i]; double value2 = array2[i]; assert (value1 == value2); } } delete save; delete [] check; return returnCode; #else networkBasis_->updateColumn(regionSparse, regionSparse2, -1); return 1; #endif } #endif } /* Updates one column (FTRAN) from region2 Tries to do FT update number returned is negative if no room. Also updates region3 region1 starts as zero and is zero at end */ int ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1, CoinIndexedVector * regionSparse2, CoinIndexedVector * regionSparse3, bool noPermuteRegion3) { int returnCode = updateColumnFT(regionSparse1, regionSparse2); updateColumn(regionSparse1, regionSparse3, noPermuteRegion3); return returnCode; } /* Updates one column (FTRAN) from region2 number returned is negative if no room region1 starts as zero and is zero at end */ int ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse, CoinIndexedVector * regionSparse2, bool noPermute) const { if (!noPermute) regionSparse->checkClear(); if (!numberRows_) return 0; collectStatistics_ = false; int returnCode = CoinFactorization::updateColumn(regionSparse, regionSparse2, noPermute); return returnCode; } /* Updates one column (BTRAN) from region2 region1 starts as zero and is zero at end */ int ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse, CoinIndexedVector * regionSparse2) const { if (!numberRows_) return 0; #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif collectStatistics_ = true; int returnCode = CoinFactorization::updateColumnTranspose(regionSparse, regionSparse2); collectStatistics_ = false; #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(6); #endif return returnCode; #ifndef SLIM_CLP } else { #ifdef CHECK_NETWORK CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2); double * check = new double[numberRows_]; int returnCode = CoinFactorization::updateColumnTranspose(regionSparse, regionSparse2); networkBasis_->updateColumnTranspose(regionSparse, save); int i; double * array = regionSparse2->denseVector(); int * indices = regionSparse2->getIndices(); int n = regionSparse2->getNumElements(); memset(check, 0, numberRows_ * sizeof(double)); double * array2 = save->denseVector(); int * indices2 = save->getIndices(); int n2 = save->getNumElements(); assert (n == n2); if (save->packedMode()) { for (i = 0; i < n; i++) { check[indices[i]] = array[i]; } for (i = 0; i < n; i++) { double value2 = array2[i]; assert (check[indices2[i]] == value2); } } else { for (i = 0; i < numberRows_; i++) { double value1 = array[i]; double value2 = array2[i]; assert (value1 == value2); } } delete save; delete [] check; return returnCode; #else return networkBasis_->updateColumnTranspose(regionSparse, regionSparse2); #endif } #endif } /* makes a row copy of L for speed and to allow very sparse problems */ void ClpFactorization::goSparse() { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif #ifndef SLIM_CLP if (!networkBasis_) #endif CoinFactorization::goSparse(); #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(7); #endif } // Cleans up i.e. gets rid of network basis void ClpFactorization::cleanUp() { #ifndef SLIM_CLP delete networkBasis_; networkBasis_ = NULL; #endif resetStatistics(); } /// Says whether to redo pivot order bool ClpFactorization::needToReorder() const { #ifdef CHECK_NETWORK return true; #endif #ifndef SLIM_CLP if (!networkBasis_) #endif return true; #ifndef SLIM_CLP else return false; #endif } // Get weighted row list void ClpFactorization::getWeights(int * weights) const { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif #ifndef SLIM_CLP if (networkBasis_) { // Network - just unit for (int i = 0; i < numberRows_; i++) weights[i] = 1; return; } #endif int * numberInRow = numberInRow_.array(); int * numberInColumn = numberInColumn_.array(); int * permuteBack = pivotColumnBack_.array(); int * indexRowU = indexRowU_.array(); const CoinBigIndex * startColumnU = startColumnU_.array(); const CoinBigIndex * startRowL = startRowL_.array(); if (!startRowL || !numberInRow_.array()) { int * temp = new int[numberRows_]; memset(temp, 0, numberRows_ * sizeof(int)); int i; for (i = 0; i < numberRows_; i++) { // one for pivot temp[i]++; CoinBigIndex j; for (j = startColumnU[i]; j < startColumnU[i] + numberInColumn[i]; j++) { int iRow = indexRowU[j]; temp[iRow]++; } } CoinBigIndex * startColumnL = startColumnL_.array(); int * indexRowL = indexRowL_.array(); for (i = baseL_; i < baseL_ + numberL_; i++) { CoinBigIndex j; for (j = startColumnL[i]; j < startColumnL[i+1]; j++) { int iRow = indexRowL[j]; temp[iRow]++; } } for (i = 0; i < numberRows_; i++) { int number = temp[i]; int iPermute = permuteBack[i]; weights[iPermute] = number; } delete [] temp; } else { int i; for (i = 0; i < numberRows_; i++) { int number = startRowL[i+1] - startRowL[i] + numberInRow[i] + 1; //number = startRowL[i+1]-startRowL[i]+1; //number = numberInRow[i]+1; int iPermute = permuteBack[i]; weights[iPermute] = number; } } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(8); #endif } #else // This one allows multiple factorizations #if CLP_MULTIPLE_FACTORIZATIONS == 1 typedef CoinDenseFactorization CoinOtherFactorization; typedef CoinOslFactorization CoinOtherFactorization; #elif CLP_MULTIPLE_FACTORIZATIONS == 2 #include "CoinSimpFactorization.hpp" typedef CoinSimpFactorization CoinOtherFactorization; typedef CoinOslFactorization CoinOtherFactorization; #elif CLP_MULTIPLE_FACTORIZATIONS == 3 #include "CoinSimpFactorization.hpp" #define CoinOslFactorization CoinDenseFactorization #elif CLP_MULTIPLE_FACTORIZATIONS == 4 #include "CoinSimpFactorization.hpp" //#define CoinOslFactorization CoinDenseFactorization #include "CoinOslFactorization.hpp" #endif //------------------------------------------------------------------- // Default Constructor //------------------------------------------------------------------- ClpFactorization::ClpFactorization () { #ifndef SLIM_CLP networkBasis_ = NULL; #endif //coinFactorizationA_ = NULL; coinFactorizationA_ = new CoinFactorization() ; coinFactorizationB_ = NULL; //coinFactorizationB_ = new CoinOtherFactorization(); forceB_ = 0; goOslThreshold_ = -1; goDenseThreshold_ = -1; goSmallThreshold_ = -1; } //------------------------------------------------------------------- // Copy constructor //------------------------------------------------------------------- ClpFactorization::ClpFactorization (const ClpFactorization & rhs, int denseIfSmaller) { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif #ifndef SLIM_CLP if (rhs.networkBasis_) networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_)); else networkBasis_ = NULL; #endif forceB_ = rhs.forceB_; goOslThreshold_ = rhs.goOslThreshold_; goDenseThreshold_ = rhs.goDenseThreshold_; goSmallThreshold_ = rhs.goSmallThreshold_; int goDense = 0; #ifdef CLP_REUSE_ETAS model_=rhs.model_; #endif if (denseIfSmaller > 0 && denseIfSmaller <= goDenseThreshold_) { CoinDenseFactorization * denseR = dynamic_cast(rhs.coinFactorizationB_); if (!denseR) goDense = 1; } if (denseIfSmaller > 0 && !rhs.coinFactorizationB_) { if (denseIfSmaller <= goDenseThreshold_) goDense = 1; else if (denseIfSmaller <= goSmallThreshold_) goDense = 2; else if (denseIfSmaller <= goOslThreshold_) goDense = 3; } else if (denseIfSmaller < 0) { if (-denseIfSmaller <= goDenseThreshold_) goDense = 1; else if (-denseIfSmaller <= goSmallThreshold_) goDense = 2; else if (-denseIfSmaller <= goOslThreshold_) goDense = 3; } if (rhs.coinFactorizationA_ && !goDense) coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_)); else coinFactorizationA_ = NULL; if (rhs.coinFactorizationB_ && (denseIfSmaller >= 0 || !goDense)) coinFactorizationB_ = rhs.coinFactorizationB_->clone(); else coinFactorizationB_ = NULL; if (goDense) { delete coinFactorizationB_; if (goDense == 1) coinFactorizationB_ = new CoinDenseFactorization(); else if (goDense == 2) coinFactorizationB_ = new CoinSimpFactorization(); else coinFactorizationB_ = new CoinOslFactorization(); if (rhs.coinFactorizationA_) { coinFactorizationB_->maximumPivots(rhs.coinFactorizationA_->maximumPivots()); coinFactorizationB_->pivotTolerance(rhs.coinFactorizationA_->pivotTolerance()); coinFactorizationB_->zeroTolerance(rhs.coinFactorizationA_->zeroTolerance()); } else { assert (coinFactorizationB_); coinFactorizationB_->maximumPivots(rhs.coinFactorizationB_->maximumPivots()); coinFactorizationB_->pivotTolerance(rhs.coinFactorizationB_->pivotTolerance()); coinFactorizationB_->zeroTolerance(rhs.coinFactorizationB_->zeroTolerance()); } } assert (!coinFactorizationA_ || !coinFactorizationB_); #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(1); #endif } ClpFactorization::ClpFactorization (const CoinFactorization & rhs) { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif #ifndef SLIM_CLP networkBasis_ = NULL; #endif coinFactorizationA_ = new CoinFactorization(rhs); coinFactorizationB_ = NULL; #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(1); #endif forceB_ = 0; goOslThreshold_ = -1; goDenseThreshold_ = -1; goSmallThreshold_ = -1; assert (!coinFactorizationA_ || !coinFactorizationB_); } ClpFactorization::ClpFactorization (const CoinOtherFactorization & rhs) { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif #ifndef SLIM_CLP networkBasis_ = NULL; #endif coinFactorizationA_ = NULL; coinFactorizationB_ = rhs.clone(); //coinFactorizationB_ = new CoinOtherFactorization(rhs); forceB_ = 0; goOslThreshold_ = -1; goDenseThreshold_ = -1; goSmallThreshold_ = -1; #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(1); #endif assert (!coinFactorizationA_ || !coinFactorizationB_); } //------------------------------------------------------------------- // Destructor //------------------------------------------------------------------- ClpFactorization::~ClpFactorization () { #ifndef SLIM_CLP delete networkBasis_; #endif delete coinFactorizationA_; delete coinFactorizationB_; } //---------------------------------------------------------------- // Assignment operator //------------------------------------------------------------------- ClpFactorization & ClpFactorization::operator=(const ClpFactorization& rhs) { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif if (this != &rhs) { #ifndef SLIM_CLP delete networkBasis_; if (rhs.networkBasis_) networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_)); else networkBasis_ = NULL; #endif forceB_ = rhs.forceB_; #ifdef CLP_REUSE_ETAS model_=rhs.model_; #endif goOslThreshold_ = rhs.goOslThreshold_; goDenseThreshold_ = rhs.goDenseThreshold_; goSmallThreshold_ = rhs.goSmallThreshold_; if (rhs.coinFactorizationA_) { if (coinFactorizationA_) *coinFactorizationA_ = *(rhs.coinFactorizationA_); else coinFactorizationA_ = new CoinFactorization(*rhs.coinFactorizationA_); } else { delete coinFactorizationA_; coinFactorizationA_ = NULL; } if (rhs.coinFactorizationB_) { if (coinFactorizationB_) { CoinDenseFactorization * denseR = dynamic_cast(rhs.coinFactorizationB_); CoinDenseFactorization * dense = dynamic_cast(coinFactorizationB_); CoinOslFactorization * oslR = dynamic_cast(rhs.coinFactorizationB_); CoinOslFactorization * osl = dynamic_cast(coinFactorizationB_); CoinSimpFactorization * simpR = dynamic_cast(rhs.coinFactorizationB_); CoinSimpFactorization * simp = dynamic_cast(coinFactorizationB_); if (dense && denseR) { *dense = *denseR; } else if (osl && oslR) { *osl = *oslR; } else if (simp && simpR) { *simp = *simpR; } else { delete coinFactorizationB_; coinFactorizationB_ = rhs.coinFactorizationB_->clone(); } } else { coinFactorizationB_ = rhs.coinFactorizationB_->clone(); } } else { delete coinFactorizationB_; coinFactorizationB_ = NULL; } } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(1); #endif assert (!coinFactorizationA_ || !coinFactorizationB_); return *this; } // Go over to dense code void ClpFactorization::goDenseOrSmall(int numberRows) { if (!forceB_) { if (numberRows <= goDenseThreshold_) { delete coinFactorizationA_; delete coinFactorizationB_; coinFactorizationA_ = NULL; coinFactorizationB_ = new CoinDenseFactorization(); //printf("going dense\n"); } else if (numberRows <= goSmallThreshold_) { delete coinFactorizationA_; delete coinFactorizationB_; coinFactorizationA_ = NULL; coinFactorizationB_ = new CoinSimpFactorization(); //printf("going small\n"); } else if (numberRows <= goOslThreshold_) { delete coinFactorizationA_; delete coinFactorizationB_; coinFactorizationA_ = NULL; coinFactorizationB_ = new CoinOslFactorization(); //printf("going small\n"); } } assert (!coinFactorizationA_ || !coinFactorizationB_); } // If nonzero force use of 1,dense 2,small 3,osl void ClpFactorization::forceOtherFactorization(int which) { delete coinFactorizationB_; forceB_ = 0; coinFactorizationB_ = NULL; if (which > 0 && which < 4) { delete coinFactorizationA_; coinFactorizationA_ = NULL; forceB_ = which; switch (which) { case 1: coinFactorizationB_ = new CoinDenseFactorization(); goDenseThreshold_ = COIN_INT_MAX; break; case 2: coinFactorizationB_ = new CoinSimpFactorization(); goSmallThreshold_ = COIN_INT_MAX; break; case 3: coinFactorizationB_ = new CoinOslFactorization(); goOslThreshold_ = COIN_INT_MAX; break; } } else if (!coinFactorizationA_) { coinFactorizationA_ = new CoinFactorization(); goOslThreshold_ = -1; goDenseThreshold_ = -1; goSmallThreshold_ = -1; } } int ClpFactorization::factorize ( ClpSimplex * model, int solveType, bool valuesPass) { #ifdef CLP_REUSE_ETAS model_= model; #endif //if ((model->specialOptions()&16384)) //printf("factor at %d iterations\n",model->numberIterations()); ClpMatrixBase * matrix = model->clpMatrix(); int numberRows = model->numberRows(); int numberColumns = model->numberColumns(); if (!numberRows) return 0; #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif bool anyChanged = false; if (coinFactorizationB_) { coinFactorizationB_->setStatus(-99); int * pivotVariable = model->pivotVariable(); //returns 0 -okay, -1 singular, -2 too many in basis */ // allow dense int solveMode = 2; if (model->numberIterations()) solveMode += 8; if (valuesPass) solveMode += 4; coinFactorizationB_->setSolveMode(solveMode); while (status() < -98) { int i; int numberBasic = 0; int numberRowBasic; // Move pivot variables across if they look good int * pivotTemp = model->rowArray(0)->getIndices(); assert (!model->rowArray(0)->getNumElements()); if (!matrix->rhsOffset(model)) { // Seems to prefer things in order so quickest // way is to go though like this for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) pivotTemp[numberBasic++] = i; } numberRowBasic = numberBasic; /* Put column basic variables into pivotVariable This is done by ClpMatrixBase to allow override for gub */ matrix->generalExpanded(model, 0, numberBasic); } else { // Long matrix - do a different way bool fullSearch = false; for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot >= numberColumns) { pivotTemp[numberBasic++] = iPivot - numberColumns; } } numberRowBasic = numberBasic; for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot < numberColumns) { if (iPivot >= 0) { pivotTemp[numberBasic++] = iPivot; } else { // not full basis fullSearch = true; break; } } } if (fullSearch) { // do slow way numberBasic = 0; for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) pivotTemp[numberBasic++] = i; } numberRowBasic = numberBasic; /* Put column basic variables into pivotVariable This is done by ClpMatrixBase to allow override for gub */ matrix->generalExpanded(model, 0, numberBasic); } } if (numberBasic > model->maximumBasic()) { // Take out some numberBasic = numberRowBasic; for (int i = 0; i < numberColumns; i++) { if (model->getColumnStatus(i) == ClpSimplex::basic) { if (numberBasic < numberRows) numberBasic++; else model->setColumnStatus(i, ClpSimplex::superBasic); } } numberBasic = numberRowBasic; matrix->generalExpanded(model, 0, numberBasic); } else if (numberBasic < numberRows) { // add in some int needed = numberRows - numberBasic; // move up columns for (i = numberBasic - 1; i >= numberRowBasic; i--) pivotTemp[i+needed] = pivotTemp[i]; numberRowBasic = 0; numberBasic = numberRows; for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) { pivotTemp[numberRowBasic++] = i; } else if (needed) { needed--; model->setRowStatus(i, ClpSimplex::basic); pivotTemp[numberRowBasic++] = i; } } } CoinBigIndex numberElements = numberRowBasic; // compute how much in basis // can change for gub int numberColumnBasic = numberBasic - numberRowBasic; numberElements += matrix->countBasis(pivotTemp + numberRowBasic, numberColumnBasic); #ifndef NDEBUG //#define CHECK_CLEAN_BASIS #ifdef CHECK_CLEAN_BASIS int saveNumberElements = numberElements; #endif #endif // Not needed for dense numberElements = 3 * numberBasic + 3 * numberElements + 20000; int numberIterations = model->numberIterations(); coinFactorizationB_->setUsefulInformation(&numberIterations, 0); coinFactorizationB_->getAreas ( numberRows, numberRowBasic + numberColumnBasic, numberElements, 2 * numberElements ); // Fill in counts so we can skip part of preProcess // This is NOT needed for dense but would be needed for later versions CoinFactorizationDouble * elementU; int * indexRowU; CoinBigIndex * startColumnU; int * numberInRow; int * numberInColumn; elementU = coinFactorizationB_->elements(); indexRowU = coinFactorizationB_->indices(); startColumnU = coinFactorizationB_->starts(); #ifdef CHECK_CLEAN_BASIS for (int i=0;islackValue(); #else #define slackValue -1.0 #endif numberInRow = coinFactorizationB_->numberInRow(); numberInColumn = coinFactorizationB_->numberInColumn(); CoinZeroN ( numberInRow, numberRows ); CoinZeroN ( numberInColumn, numberRows ); for (i = 0; i < numberRowBasic; i++) { int iRow = pivotTemp[i]; // Change pivotTemp to correct sequence pivotTemp[i] = iRow + numberColumns; indexRowU[i] = iRow; startColumnU[i] = i; elementU[i] = slackValue; numberInRow[iRow] = 1; numberInColumn[i] = 1; } startColumnU[numberRowBasic] = numberRowBasic; // can change for gub so redo numberColumnBasic = numberBasic - numberRowBasic; matrix->fillBasis(model, pivotTemp + numberRowBasic, numberColumnBasic, indexRowU, startColumnU + numberRowBasic, numberInRow, numberInColumn + numberRowBasic, elementU); // recompute number basic numberBasic = numberRowBasic + numberColumnBasic; for (i = numberBasic; i < numberRows; i++) pivotTemp[i] = -1; // mark not there if (numberBasic) numberElements = startColumnU[numberBasic-1] + numberInColumn[numberBasic-1]; else numberElements = 0; #ifdef CHECK_CLEAN_BASIS assert (!startColumnU[0]); int lastStart=0; for (int i=0;ilastStart); lastStart=startColumnU[i+1]; } assert (lastStart==saveNumberElements); for (int i=0;i=0&&indexRowU[i]preProcess ( ); coinFactorizationB_->factor ( ); if (coinFactorizationB_->status() == -1 && (coinFactorizationB_->solveMode() % 3) != 0) { int solveMode = coinFactorizationB_->solveMode(); solveMode -= solveMode % 3; // so bottom will be 0 coinFactorizationB_->setSolveMode(solveMode); coinFactorizationB_->setStatus(-99); } if (coinFactorizationB_->status() == -99) continue; // If we get here status is 0 or -1 if (coinFactorizationB_->status() == 0 && numberBasic == numberRows) { coinFactorizationB_->postProcess(pivotTemp, pivotVariable); } else if (solveType == 0 || solveType == 2/*||solveType==1*/) { // Change pivotTemp to be correct list anyChanged = true; coinFactorizationB_->makeNonSingular(pivotTemp, numberColumns); double * columnLower = model->lowerRegion(); double * columnUpper = model->upperRegion(); double * columnActivity = model->solutionRegion(); double * rowLower = model->lowerRegion(0); double * rowUpper = model->upperRegion(0); double * rowActivity = model->solutionRegion(0); //redo basis - first take ALL out int iColumn; double largeValue = model->largeValue(); for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getColumnStatus(iColumn) == ClpSimplex::basic) { // take out if (!valuesPass) { double lower = columnLower[iColumn]; double upper = columnUpper[iColumn]; double value = columnActivity[iColumn]; if (lower > -largeValue || upper < largeValue) { if (fabs(value - lower) < fabs(value - upper)) { model->setColumnStatus(iColumn, ClpSimplex::atLowerBound); columnActivity[iColumn] = lower; } else { model->setColumnStatus(iColumn, ClpSimplex::atUpperBound); columnActivity[iColumn] = upper; } } else { model->setColumnStatus(iColumn, ClpSimplex::isFree); } } else { model->setColumnStatus(iColumn, ClpSimplex::superBasic); } } } int iRow; for (iRow = 0; iRow < numberRows; iRow++) { if (model->getRowStatus(iRow) == ClpSimplex::basic) { // take out if (!valuesPass) { double lower = columnLower[iRow]; double upper = columnUpper[iRow]; double value = columnActivity[iRow]; if (lower > -largeValue || upper < largeValue) { if (fabs(value - lower) < fabs(value - upper)) { model->setRowStatus(iRow, ClpSimplex::atLowerBound); columnActivity[iRow] = lower; } else { model->setRowStatus(iRow, ClpSimplex::atUpperBound); columnActivity[iRow] = upper; } } else { model->setRowStatus(iRow, ClpSimplex::isFree); } } else { model->setRowStatus(iRow, ClpSimplex::superBasic); } } } for (iRow = 0; iRow < numberRows; iRow++) { int iSequence = pivotTemp[iRow]; assert (iSequence >= 0); // basic model->setColumnStatus(iSequence, ClpSimplex::basic); } // signal repeat coinFactorizationB_->setStatus(-99); // set fixed if they are for (iRow = 0; iRow < numberRows; iRow++) { if (model->getRowStatus(iRow) != ClpSimplex::basic ) { if (rowLower[iRow] == rowUpper[iRow]) { rowActivity[iRow] = rowLower[iRow]; model->setRowStatus(iRow, ClpSimplex::isFixed); } } } for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getColumnStatus(iColumn) != ClpSimplex::basic ) { if (columnLower[iColumn] == columnUpper[iColumn]) { columnActivity[iColumn] = columnLower[iColumn]; model->setColumnStatus(iColumn, ClpSimplex::isFixed); } } } } } #ifdef CLP_DEBUG // check basic CoinIndexedVector region1(2 * numberRows); CoinIndexedVector region2B(2 * numberRows); int iPivot; double * arrayB = region2B.denseVector(); int i; for (iPivot = 0; iPivot < numberRows; iPivot++) { int iSequence = pivotVariable[iPivot]; model->unpack(®ion2B, iSequence); coinFactorizationB_->updateColumn(®ion1, ®ion2B); if (fabs(arrayB[iPivot] - 1.0) < 1.0e-4) { // OK? arrayB[iPivot] = 0.0; } else { assert (fabs(arrayB[iPivot]) < 1.0e-4); for (i = 0; i < numberRows; i++) { if (fabs(arrayB[i] - 1.0) < 1.0e-4) break; } assert (i < numberRows); printf("variable on row %d landed up on row %d\n", iPivot, i); arrayB[i] = 0.0; } for (i = 0; i < numberRows; i++) assert (fabs(arrayB[i]) < 1.0e-4); region2B.clear(); } #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(2); #endif if ( anyChanged && model->algorithm() < 0 && solveType > 0) { double dummyCost; static_cast (model)->changeBounds(3, NULL, dummyCost); } return coinFactorizationB_->status(); } // If too many compressions increase area if (coinFactorizationA_->pivots() > 1 && coinFactorizationA_->numberCompressions() * 10 > coinFactorizationA_->pivots() + 10) { coinFactorizationA_->areaFactor( coinFactorizationA_->areaFactor() * 1.1); } //int numberPivots=coinFactorizationA_->pivots(); #if 0 if (model->algorithm() > 0) numberSave = -1; #endif #ifndef SLIM_CLP if (!networkBasis_ || doCheck) { #endif coinFactorizationA_->setStatus(-99); int * pivotVariable = model->pivotVariable(); int nTimesRound = 0; //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */ while (coinFactorizationA_->status() < -98) { nTimesRound++; int i; int numberBasic = 0; int numberRowBasic; // Move pivot variables across if they look good int * pivotTemp = model->rowArray(0)->getIndices(); assert (!model->rowArray(0)->getNumElements()); if (!matrix->rhsOffset(model)) { #if 0 if (numberSave > 0) { int nStill = 0; int nAtBound = 0; int nZeroDual = 0; CoinIndexedVector * array = model->rowArray(3); CoinIndexedVector * objArray = model->columnArray(1); array->clear(); objArray->clear(); double * cost = model->costRegion(); double tolerance = model->primalTolerance(); double offset = 0.0; for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot < numberColumns && isDense(iPivot)) { if (model->getColumnStatus(iPivot) == ClpSimplex::basic) { nStill++; double value = model->solutionRegion()[iPivot]; double dual = model->dualRowSolution()[i]; double lower = model->lowerRegion()[iPivot]; double upper = model->upperRegion()[iPivot]; ClpSimplex::Status status; if (fabs(value - lower) < tolerance) { status = ClpSimplex::atLowerBound; nAtBound++; } else if (fabs(value - upper) < tolerance) { nAtBound++; status = ClpSimplex::atUpperBound; } else if (value > lower && value < upper) { status = ClpSimplex::superBasic; } else { status = ClpSimplex::basic; } if (status != ClpSimplex::basic) { if (model->getRowStatus(i) != ClpSimplex::basic) { model->setColumnStatus(iPivot, ClpSimplex::atLowerBound); model->setRowStatus(i, ClpSimplex::basic); pivotVariable[i] = i + numberColumns; model->dualRowSolution()[i] = 0.0; model->djRegion(0)[i] = 0.0; array->add(i, dual); offset += dual * model->solutionRegion(0)[i]; } } if (fabs(dual) < 1.0e-5) nZeroDual++; } } } printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n", numberSave, nStill, nAtBound, nZeroDual, offset); if (array->getNumElements()) { // modify costs model->clpMatrix()->transposeTimes(model, 1.0, array, model->columnArray(0), objArray); array->clear(); int n = objArray->getNumElements(); int * indices = objArray->getIndices(); double * elements = objArray->denseVector(); for (i = 0; i < n; i++) { int iColumn = indices[i]; cost[iColumn] -= elements[iColumn]; elements[iColumn] = 0.0; } objArray->setNumElements(0); } } #endif // Seems to prefer things in order so quickest // way is to go though like this for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) pivotTemp[numberBasic++] = i; } numberRowBasic = numberBasic; /* Put column basic variables into pivotVariable This is done by ClpMatrixBase to allow override for gub */ matrix->generalExpanded(model, 0, numberBasic); } else { // Long matrix - do a different way bool fullSearch = false; for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot >= numberColumns) { pivotTemp[numberBasic++] = iPivot - numberColumns; } } numberRowBasic = numberBasic; for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot < numberColumns) { if (iPivot >= 0) { pivotTemp[numberBasic++] = iPivot; } else { // not full basis fullSearch = true; break; } } } if (fullSearch) { // do slow way numberBasic = 0; for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) pivotTemp[numberBasic++] = i; } numberRowBasic = numberBasic; /* Put column basic variables into pivotVariable This is done by ClpMatrixBase to allow override for gub */ matrix->generalExpanded(model, 0, numberBasic); } } if (numberBasic > model->maximumBasic()) { #if 0 // ndef NDEBUG printf("%d basic - should only be %d\n", numberBasic, numberRows); #endif // Take out some numberBasic = numberRowBasic; for (int i = 0; i < numberColumns; i++) { if (model->getColumnStatus(i) == ClpSimplex::basic) { if (numberBasic < numberRows) numberBasic++; else model->setColumnStatus(i, ClpSimplex::superBasic); } } numberBasic = numberRowBasic; matrix->generalExpanded(model, 0, numberBasic); } #ifndef SLIM_CLP // see if matrix a network #ifndef NO_RTTI ClpNetworkMatrix* networkMatrix = dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix()); #else ClpNetworkMatrix* networkMatrix = NULL; if (model->clpMatrix()->type() == 11) networkMatrix = static_cast< ClpNetworkMatrix*>(model->clpMatrix()); #endif // If network - still allow ordinary factorization first time for laziness if (networkMatrix) coinFactorizationA_->setBiasLU(0); // All to U if network //int saveMaximumPivots = maximumPivots(); delete networkBasis_; networkBasis_ = NULL; if (networkMatrix && !doCheck) maximumPivots(1); #endif //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR()); while (coinFactorizationA_->status() == -99) { // maybe for speed will be better to leave as many regions as possible coinFactorizationA_->gutsOfDestructor(); coinFactorizationA_->gutsOfInitialize(2); CoinBigIndex numberElements = numberRowBasic; // compute how much in basis int i; // can change for gub int numberColumnBasic = numberBasic - numberRowBasic; numberElements += matrix->countBasis( pivotTemp + numberRowBasic, numberColumnBasic); // and recompute as network side say different if (model->numberIterations()) numberRowBasic = numberBasic - numberColumnBasic; numberElements = 3 * numberBasic + 3 * numberElements + 20000; coinFactorizationA_->getAreas ( numberRows, numberRowBasic + numberColumnBasic, numberElements, 2 * numberElements ); //fill // Fill in counts so we can skip part of preProcess int * numberInRow = coinFactorizationA_->numberInRow(); int * numberInColumn = coinFactorizationA_->numberInColumn(); CoinZeroN ( numberInRow, coinFactorizationA_->numberRows() + 1 ); CoinZeroN ( numberInColumn, coinFactorizationA_->maximumColumnsExtra() + 1 ); CoinFactorizationDouble * elementU = coinFactorizationA_->elementU(); int * indexRowU = coinFactorizationA_->indexRowU(); CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU(); #ifndef COIN_FAST_CODE double slackValue = coinFactorizationA_->slackValue(); #endif for (i = 0; i < numberRowBasic; i++) { int iRow = pivotTemp[i]; indexRowU[i] = iRow; startColumnU[i] = i; elementU[i] = slackValue; numberInRow[iRow] = 1; numberInColumn[i] = 1; } startColumnU[numberRowBasic] = numberRowBasic; // can change for gub so redo numberColumnBasic = numberBasic - numberRowBasic; matrix->fillBasis(model, pivotTemp + numberRowBasic, numberColumnBasic, indexRowU, startColumnU + numberRowBasic, numberInRow, numberInColumn + numberRowBasic, elementU); #if 0 { printf("%d row basic, %d column basic\n", numberRowBasic, numberColumnBasic); for (int i = 0; i < numberElements; i++) printf("row %d col %d value %g\n", indexRowU[i], indexColumnU_[i], elementU[i]); } #endif // recompute number basic numberBasic = numberRowBasic + numberColumnBasic; if (numberBasic) numberElements = startColumnU[numberBasic-1] + numberInColumn[numberBasic-1]; else numberElements = 0; coinFactorizationA_->setNumberElementsU(numberElements); //saveFactorization("dump.d"); if (coinFactorizationA_->biasLU() >= 3 || coinFactorizationA_->numberRows() != coinFactorizationA_->numberColumns()) coinFactorizationA_->preProcess ( 2 ); else coinFactorizationA_->preProcess ( 3 ); // no row copy coinFactorizationA_->factor ( ); if (coinFactorizationA_->status() == -99) { // get more memory coinFactorizationA_->areaFactor(2.0 * coinFactorizationA_->areaFactor()); } else if (coinFactorizationA_->status() == -1 && (model->numberIterations() == 0 || nTimesRound > 2) && coinFactorizationA_->denseThreshold()) { // Round again without dense coinFactorizationA_->setDenseThreshold(0); coinFactorizationA_->setStatus(-99); } } // If we get here status is 0 or -1 if (coinFactorizationA_->status() == 0) { // We may need to tamper with order and redo - e.g. network with side int useNumberRows = numberRows; // **** we will also need to add test in dual steepest to do // as we do for network matrix->generalExpanded(model, 12, useNumberRows); const int * permuteBack = coinFactorizationA_->permuteBack(); const int * back = coinFactorizationA_->pivotColumnBack(); //int * pivotTemp = pivotColumn_.array(); //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp ); #ifndef NDEBUG CoinFillN(pivotVariable, numberRows, -1); #endif // Redo pivot order for (i = 0; i < numberRowBasic; i++) { int k = pivotTemp[i]; // so rowIsBasic[k] would be permuteBack[back[i]] int j = permuteBack[back[i]]; assert (pivotVariable[j] == -1); pivotVariable[j] = k + numberColumns; } for (; i < useNumberRows; i++) { int k = pivotTemp[i]; // so rowIsBasic[k] would be permuteBack[back[i]] int j = permuteBack[back[i]]; assert (pivotVariable[j] == -1); pivotVariable[j] = k; } #if 0 if (numberSave >= 0) { numberSave = numberDense_; memset(saveList, 0, ((coinFactorizationA_->numberRows() + 31) >> 5)*sizeof(int)); for (i = coinFactorizationA_->numberRows() - numberSave; i < coinFactorizationA_->numberRows(); i++) { int k = pivotTemp[pivotColumn_.array()[i]]; setDense(k); } } #endif // Set up permutation vector // these arrays start off as copies of permute // (and we could use permute_ instead of pivotColumn (not back though)) ClpDisjointCopyN ( coinFactorizationA_->permute(), useNumberRows , coinFactorizationA_->pivotColumn() ); ClpDisjointCopyN ( coinFactorizationA_->permuteBack(), useNumberRows , coinFactorizationA_->pivotColumnBack() ); #ifndef SLIM_CLP if (networkMatrix) { maximumPivots(CoinMax(2000, maximumPivots())); // redo arrays for (int iRow = 0; iRow < 4; iRow++) { int length = model->numberRows() + maximumPivots(); if (iRow == 3 || model->objectiveAsObject()->type() > 1) length += model->numberColumns(); model->rowArray(iRow)->reserve(length); } // create network factorization if (doCheck) delete networkBasis_; // temp networkBasis_ = new ClpNetworkBasis(model, coinFactorizationA_->numberRows(), coinFactorizationA_->pivotRegion(), coinFactorizationA_->permuteBack(), coinFactorizationA_->startColumnU(), coinFactorizationA_->numberInColumn(), coinFactorizationA_->indexRowU(), coinFactorizationA_->elementU()); // kill off arrays in ordinary factorization if (!doCheck) { coinFactorizationA_->gutsOfDestructor(); // but make sure coinFactorizationA_->numberRows() set coinFactorizationA_->setNumberRows(model->numberRows()); coinFactorizationA_->setStatus(0); #if 0 // but put back permute arrays so odd things will work int numberRows = model->numberRows(); pivotColumnBack_ = new int [numberRows]; permute_ = new int [numberRows]; int i; for (i = 0; i < numberRows; i++) { pivotColumnBack_[i] = i; permute_[i] = i; } #endif } } else { #endif // See if worth going sparse and when coinFactorizationA_->checkSparse(); #ifndef SLIM_CLP } #endif } else if (coinFactorizationA_->status() == -1 && (solveType == 0 || solveType == 2)) { // This needs redoing as it was merged coding - does not need array int numberTotal = numberRows + numberColumns; int * isBasic = new int [numberTotal]; int * rowIsBasic = isBasic + numberColumns; int * columnIsBasic = isBasic; for (i = 0; i < numberTotal; i++) isBasic[i] = -1; for (i = 0; i < numberRowBasic; i++) { int iRow = pivotTemp[i]; rowIsBasic[iRow] = 1; } for (; i < numberBasic; i++) { int iColumn = pivotTemp[i]; columnIsBasic[iColumn] = 1; } numberBasic = 0; for (i = 0; i < numberRows; i++) pivotVariable[i] = -1; // mark as basic or non basic const int * pivotColumn = coinFactorizationA_->pivotColumn(); for (i = 0; i < numberRows; i++) { if (rowIsBasic[i] >= 0) { if (pivotColumn[numberBasic] >= 0) { rowIsBasic[i] = pivotColumn[numberBasic]; } else { rowIsBasic[i] = -1; model->setRowStatus(i, ClpSimplex::superBasic); } numberBasic++; } } for (i = 0; i < numberColumns; i++) { if (columnIsBasic[i] >= 0) { if (pivotColumn[numberBasic] >= 0) columnIsBasic[i] = pivotColumn[numberBasic]; else columnIsBasic[i] = -1; numberBasic++; } } // leave pivotVariable in useful form for cleaning basis int * pivotVariable = model->pivotVariable(); for (i = 0; i < numberRows; i++) { pivotVariable[i] = -1; } for (i = 0; i < numberRows; i++) { if (model->getRowStatus(i) == ClpSimplex::basic) { int iPivot = rowIsBasic[i]; if (iPivot >= 0) pivotVariable[iPivot] = i + numberColumns; } } for (i = 0; i < numberColumns; i++) { if (model->getColumnStatus(i) == ClpSimplex::basic) { int iPivot = columnIsBasic[i]; if (iPivot >= 0) pivotVariable[iPivot] = i; } } delete [] isBasic; double * columnLower = model->lowerRegion(); double * columnUpper = model->upperRegion(); double * columnActivity = model->solutionRegion(); double * rowLower = model->lowerRegion(0); double * rowUpper = model->upperRegion(0); double * rowActivity = model->solutionRegion(0); //redo basis - first take ALL columns out int iColumn; double largeValue = model->largeValue(); for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getColumnStatus(iColumn) == ClpSimplex::basic) { // take out if (!valuesPass) { double lower = columnLower[iColumn]; double upper = columnUpper[iColumn]; double value = columnActivity[iColumn]; if (lower > -largeValue || upper < largeValue) { if (fabs(value - lower) < fabs(value - upper)) { model->setColumnStatus(iColumn, ClpSimplex::atLowerBound); columnActivity[iColumn] = lower; } else { model->setColumnStatus(iColumn, ClpSimplex::atUpperBound); columnActivity[iColumn] = upper; } } else { model->setColumnStatus(iColumn, ClpSimplex::isFree); } } else { model->setColumnStatus(iColumn, ClpSimplex::superBasic); } } } int iRow; for (iRow = 0; iRow < numberRows; iRow++) { int iSequence = pivotVariable[iRow]; if (iSequence >= 0) { // basic if (iSequence >= numberColumns) { // slack in - leave //assert (iSequence-numberColumns==iRow); } else { assert(model->getRowStatus(iRow) != ClpSimplex::basic); // put back structural model->setColumnStatus(iSequence, ClpSimplex::basic); } } else { // put in slack model->setRowStatus(iRow, ClpSimplex::basic); } } // Put back any key variables for gub int dummy; matrix->generalExpanded(model, 1, dummy); // signal repeat coinFactorizationA_->setStatus(-99); // set fixed if they are for (iRow = 0; iRow < numberRows; iRow++) { if (model->getRowStatus(iRow) != ClpSimplex::basic ) { if (rowLower[iRow] == rowUpper[iRow]) { rowActivity[iRow] = rowLower[iRow]; model->setRowStatus(iRow, ClpSimplex::isFixed); } } } for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getColumnStatus(iColumn) != ClpSimplex::basic ) { if (columnLower[iColumn] == columnUpper[iColumn]) { columnActivity[iColumn] = columnLower[iColumn]; model->setColumnStatus(iColumn, ClpSimplex::isFixed); } } } } } #ifndef SLIM_CLP } else { // network - fake factorization - do nothing coinFactorizationA_->setStatus(0); coinFactorizationA_->setPivots(0); } #endif #ifndef SLIM_CLP if (!coinFactorizationA_->status()) { // take out part if quadratic if (model->algorithm() == 2) { ClpObjective * obj = model->objectiveAsObject(); #ifndef NDEBUG ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj)); assert (quadraticObj); #else ClpQuadraticObjective * quadraticObj = (static_cast< ClpQuadraticObjective*>(obj)); #endif CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective(); int numberXColumns = quadratic->getNumCols(); assert (numberXColumns < numberColumns); int base = numberColumns - numberXColumns; int * which = new int [numberXColumns]; int * pivotVariable = model->pivotVariable(); int * permute = pivotColumn(); int i; int n = 0; for (i = 0; i < numberRows; i++) { int iSj = pivotVariable[i] - base; if (iSj >= 0 && iSj < numberXColumns) which[n++] = permute[i]; } if (n) coinFactorizationA_->emptyRows(n, which); delete [] which; } } #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(2); #endif return coinFactorizationA_->status(); } /* Replaces one Column in basis, returns 0=OK, 1=Probably OK, 2=singular, 3=no room If checkBeforeModifying is true will do all accuracy checks before modifying factorization. Whether to set this depends on speed considerations. You could just do this on first iteration after factorization and thereafter re-factorize partial update already in U */ int ClpFactorization::replaceColumn ( const ClpSimplex * model, CoinIndexedVector * regionSparse, CoinIndexedVector * tableauColumn, int pivotRow, double pivotCheck , bool checkBeforeModifying, double acceptablePivot) { #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif int returnCode; // see if FT if (!coinFactorizationA_ || coinFactorizationA_->forrestTomlin()) { if (coinFactorizationA_) { returnCode = coinFactorizationA_->replaceColumn(regionSparse, pivotRow, pivotCheck, checkBeforeModifying, acceptablePivot); } else { bool tab = coinFactorizationB_->wantsTableauColumn(); #ifdef CLP_REUSE_ETAS int tempInfo[2]; tempInfo[1] = model_->sequenceOut(); #else int tempInfo[1]; #endif tempInfo[0] = model->numberIterations(); coinFactorizationB_->setUsefulInformation(tempInfo, 1); returnCode = coinFactorizationB_->replaceColumn(tab ? tableauColumn : regionSparse, pivotRow, pivotCheck, checkBeforeModifying, acceptablePivot); #ifdef CLP_DEBUG // check basic int numberRows = coinFactorizationB_->numberRows(); CoinIndexedVector region1(2 * numberRows); CoinIndexedVector region2A(2 * numberRows); CoinIndexedVector region2B(2 * numberRows); int iPivot; double * arrayB = region2B.denseVector(); int * pivotVariable = model->pivotVariable(); int i; for (iPivot = 0; iPivot < numberRows; iPivot++) { int iSequence = pivotVariable[iPivot]; if (iPivot == pivotRow) iSequence = model->sequenceIn(); model->unpack(®ion2B, iSequence); coinFactorizationB_->updateColumn(®ion1, ®ion2B); assert (fabs(arrayB[iPivot] - 1.0) < 1.0e-4); arrayB[iPivot] = 0.0; for (i = 0; i < numberRows; i++) assert (fabs(arrayB[i]) < 1.0e-4); region2B.clear(); } #endif } } else { returnCode = coinFactorizationA_->replaceColumnPFI(tableauColumn, pivotRow, pivotCheck); // Note array } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(3); #endif return returnCode; #ifndef SLIM_CLP } else { if (doCheck) { int returnCode = coinFactorizationA_->replaceColumn(regionSparse, pivotRow, pivotCheck, checkBeforeModifying, acceptablePivot); networkBasis_->replaceColumn(regionSparse, pivotRow); return returnCode; } else { // increase number of pivots coinFactorizationA_->setPivots(coinFactorizationA_->pivots() + 1); return networkBasis_->replaceColumn(regionSparse, pivotRow); } } #endif } /* Updates one column (FTRAN) from region2 number returned is negative if no room region1 starts as zero and is zero at end */ int ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse, CoinIndexedVector * regionSparse2) { #ifdef CLP_DEBUG regionSparse->checkClear(); #endif if (!numberRows()) return 0; #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif int returnCode; if (coinFactorizationA_) { coinFactorizationA_->setCollectStatistics(true); returnCode = coinFactorizationA_->updateColumnFT(regionSparse, regionSparse2); coinFactorizationA_->setCollectStatistics(false); } else { #ifdef CLP_REUSE_ETAS int tempInfo[2]; tempInfo[0] = model_->numberIterations(); tempInfo[1] = model_->sequenceIn(); coinFactorizationB_->setUsefulInformation(tempInfo, 2); #endif returnCode = coinFactorizationB_->updateColumnFT(regionSparse, regionSparse2); } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(4); #endif return returnCode; #ifndef SLIM_CLP } else { #ifdef CHECK_NETWORK CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2); double * check = new double[coinFactorizationA_->numberRows()]; int returnCode = coinFactorizationA_->updateColumnFT(regionSparse, regionSparse2); networkBasis_->updateColumn(regionSparse, save, -1); int i; double * array = regionSparse2->denseVector(); int * indices = regionSparse2->getIndices(); int n = regionSparse2->getNumElements(); memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double)); double * array2 = save->denseVector(); int * indices2 = save->getIndices(); int n2 = save->getNumElements(); assert (n == n2); if (save->packedMode()) { for (i = 0; i < n; i++) { check[indices[i]] = array[i]; } for (i = 0; i < n; i++) { double value2 = array2[i]; assert (check[indices2[i]] == value2); } } else { int numberRows = coinFactorizationA_->numberRows(); for (i = 0; i < numberRows; i++) { double value1 = array[i]; double value2 = array2[i]; assert (value1 == value2); } } delete save; delete [] check; return returnCode; #else networkBasis_->updateColumn(regionSparse, regionSparse2, -1); return 1; #endif } #endif } /* Updates one column (FTRAN) from region2 number returned is negative if no room region1 starts as zero and is zero at end */ int ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse, CoinIndexedVector * regionSparse2, bool noPermute) const { #ifdef CLP_DEBUG if (!noPermute) regionSparse->checkClear(); #endif if (!numberRows()) return 0; #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif int returnCode; if (coinFactorizationA_) { coinFactorizationA_->setCollectStatistics(true); returnCode = coinFactorizationA_->updateColumn(regionSparse, regionSparse2, noPermute); coinFactorizationA_->setCollectStatistics(false); } else { returnCode = coinFactorizationB_->updateColumn(regionSparse, regionSparse2, noPermute); } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(5); #endif //#define PRINT_VECTOR #ifdef PRINT_VECTOR printf("Update\n"); regionSparse2->print(); #endif return returnCode; #ifndef SLIM_CLP } else { #ifdef CHECK_NETWORK CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2); double * check = new double[coinFactorizationA_->numberRows()]; int returnCode = coinFactorizationA_->updateColumn(regionSparse, regionSparse2, noPermute); networkBasis_->updateColumn(regionSparse, save, -1); int i; double * array = regionSparse2->denseVector(); int * indices = regionSparse2->getIndices(); int n = regionSparse2->getNumElements(); memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double)); double * array2 = save->denseVector(); int * indices2 = save->getIndices(); int n2 = save->getNumElements(); assert (n == n2); if (save->packedMode()) { for (i = 0; i < n; i++) { check[indices[i]] = array[i]; } for (i = 0; i < n; i++) { double value2 = array2[i]; assert (check[indices2[i]] == value2); } } else { int numberRows = coinFactorizationA_->numberRows(); for (i = 0; i < numberRows; i++) { double value1 = array[i]; double value2 = array2[i]; assert (value1 == value2); } } delete save; delete [] check; return returnCode; #else networkBasis_->updateColumn(regionSparse, regionSparse2, -1); return 1; #endif } #endif } /* Updates one column (FTRAN) from region2 Tries to do FT update number returned is negative if no room. Also updates region3 region1 starts as zero and is zero at end */ int ClpFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1, CoinIndexedVector * regionSparse2, CoinIndexedVector * regionSparse3, bool noPermuteRegion3) { #ifdef CLP_DEBUG regionSparse1->checkClear(); #endif if (!numberRows()) return 0; int returnCode = 0; #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif if (coinFactorizationA_) { coinFactorizationA_->setCollectStatistics(true); if (coinFactorizationA_->spaceForForrestTomlin()) { assert (regionSparse2->packedMode()); assert (!regionSparse3->packedMode()); returnCode = coinFactorizationA_->updateTwoColumnsFT(regionSparse1, regionSparse2, regionSparse3, noPermuteRegion3); } else { returnCode = coinFactorizationA_->updateColumnFT(regionSparse1, regionSparse2); coinFactorizationA_->updateColumn(regionSparse1, regionSparse3, noPermuteRegion3); } coinFactorizationA_->setCollectStatistics(false); } else { #if 0 CoinSimpFactorization * fact = dynamic_cast< CoinSimpFactorization*>(coinFactorizationB_); if (!fact) { returnCode = coinFactorizationB_->updateColumnFT(regionSparse1, regionSparse2); coinFactorizationB_->updateColumn(regionSparse1, regionSparse3, noPermuteRegion3); } else { returnCode = fact->updateTwoColumnsFT(regionSparse1, regionSparse2, regionSparse3, noPermuteRegion3); } #else #ifdef CLP_REUSE_ETAS int tempInfo[2]; tempInfo[0] = model_->numberIterations(); tempInfo[1] = model_->sequenceIn(); coinFactorizationB_->setUsefulInformation(tempInfo, 3); #endif returnCode = coinFactorizationB_->updateTwoColumnsFT( regionSparse1, regionSparse2, regionSparse3, noPermuteRegion3); #endif } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(9); #endif #ifdef PRINT_VECTOR printf("UpdateTwoFT\n"); regionSparse2->print(); regionSparse3->print(); #endif return returnCode; #ifndef SLIM_CLP } else { returnCode = updateColumnFT(regionSparse1, regionSparse2); updateColumn(regionSparse1, regionSparse3, noPermuteRegion3); } #endif return returnCode; } /* Updates one column (FTRAN) from region2 number returned is negative if no room region1 starts as zero and is zero at end */ int ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse, CoinIndexedVector * regionSparse2, bool noPermute) const { if (!noPermute) regionSparse->checkClear(); if (!coinFactorizationA_->numberRows()) return 0; coinFactorizationA_->setCollectStatistics(false); int returnCode = coinFactorizationA_->updateColumn(regionSparse, regionSparse2, noPermute); return returnCode; } /* Updates one column (BTRAN) from region2 region1 starts as zero and is zero at end */ int ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse, CoinIndexedVector * regionSparse2) const { if (!numberRows()) return 0; #ifndef SLIM_CLP if (!networkBasis_) { #endif #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif int returnCode; if (coinFactorizationA_) { coinFactorizationA_->setCollectStatistics(true); returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse, regionSparse2); coinFactorizationA_->setCollectStatistics(false); } else { returnCode = coinFactorizationB_->updateColumnTranspose(regionSparse, regionSparse2); } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(6); #endif #ifdef PRINT_VECTOR printf("UpdateTranspose\n"); regionSparse2->print(); #endif return returnCode; #ifndef SLIM_CLP } else { #ifdef CHECK_NETWORK CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2); double * check = new double[coinFactorizationA_->numberRows()]; int returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse, regionSparse2); networkBasis_->updateColumnTranspose(regionSparse, save); int i; double * array = regionSparse2->denseVector(); int * indices = regionSparse2->getIndices(); int n = regionSparse2->getNumElements(); memset(check, 0, coinFactorizationA_->numberRows()*sizeof(double)); double * array2 = save->denseVector(); int * indices2 = save->getIndices(); int n2 = save->getNumElements(); assert (n == n2); if (save->packedMode()) { for (i = 0; i < n; i++) { check[indices[i]] = array[i]; } for (i = 0; i < n; i++) { double value2 = array2[i]; assert (check[indices2[i]] == value2); } } else { int numberRows = coinFactorizationA_->numberRows(); for (i = 0; i < numberRows; i++) { double value1 = array[i]; double value2 = array2[i]; assert (value1 == value2); } } delete save; delete [] check; return returnCode; #else return networkBasis_->updateColumnTranspose(regionSparse, regionSparse2); #endif } #endif } /* makes a row copy of L for speed and to allow very sparse problems */ void ClpFactorization::goSparse() { #ifndef SLIM_CLP if (!networkBasis_) { #endif if (coinFactorizationA_) { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif coinFactorizationA_->goSparse(); #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(7); #endif } } } // Cleans up i.e. gets rid of network basis void ClpFactorization::cleanUp() { #ifndef SLIM_CLP delete networkBasis_; networkBasis_ = NULL; #endif if (coinFactorizationA_) coinFactorizationA_->resetStatistics(); } /// Says whether to redo pivot order bool ClpFactorization::needToReorder() const { #ifdef CHECK_NETWORK return true; #endif #ifndef SLIM_CLP if (!networkBasis_) #endif return true; #ifndef SLIM_CLP else return false; #endif } // Get weighted row list void ClpFactorization::getWeights(int * weights) const { #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(-1); #endif #ifndef SLIM_CLP if (networkBasis_) { // Network - just unit int numberRows = coinFactorizationA_->numberRows(); for (int i = 0; i < numberRows; i++) weights[i] = 1; return; } #endif int * numberInRow = coinFactorizationA_->numberInRow(); int * numberInColumn = coinFactorizationA_->numberInColumn(); int * permuteBack = coinFactorizationA_->pivotColumnBack(); int * indexRowU = coinFactorizationA_->indexRowU(); const CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU(); const CoinBigIndex * startRowL = coinFactorizationA_->startRowL(); int numberRows = coinFactorizationA_->numberRows(); if (!startRowL || !coinFactorizationA_->numberInRow()) { int * temp = new int[numberRows]; memset(temp, 0, numberRows * sizeof(int)); int i; for (i = 0; i < numberRows; i++) { // one for pivot temp[i]++; CoinBigIndex j; for (j = startColumnU[i]; j < startColumnU[i] + numberInColumn[i]; j++) { int iRow = indexRowU[j]; temp[iRow]++; } } CoinBigIndex * startColumnL = coinFactorizationA_->startColumnL(); int * indexRowL = coinFactorizationA_->indexRowL(); int numberL = coinFactorizationA_->numberL(); CoinBigIndex baseL = coinFactorizationA_->baseL(); for (i = baseL; i < baseL + numberL; i++) { CoinBigIndex j; for (j = startColumnL[i]; j < startColumnL[i+1]; j++) { int iRow = indexRowL[j]; temp[iRow]++; } } for (i = 0; i < numberRows; i++) { int number = temp[i]; int iPermute = permuteBack[i]; weights[iPermute] = number; } delete [] temp; } else { int i; for (i = 0; i < numberRows; i++) { int number = startRowL[i+1] - startRowL[i] + numberInRow[i] + 1; //number = startRowL[i+1]-startRowL[i]+1; //number = numberInRow[i]+1; int iPermute = permuteBack[i]; weights[iPermute] = number; } } #ifdef CLP_FACTORIZATION_INSTRUMENT factorization_instrument(8); #endif } // Set tolerances to safer of existing and given void ClpFactorization::saferTolerances ( double zeroValue, double pivotValue) { double newValue; // better to have small tolerance even if slower if (zeroValue > 0.0) newValue = zeroValue; else newValue = -zeroTolerance() * zeroValue; zeroTolerance(CoinMin(zeroTolerance(), zeroValue)); // better to have large tolerance even if slower if (pivotValue > 0.0) newValue = pivotValue; else newValue = -pivotTolerance() * pivotValue; pivotTolerance(CoinMin(CoinMax(pivotTolerance(), newValue), 0.999)); } // Sets factorization void ClpFactorization::setFactorization(ClpFactorization & rhs) { ClpFactorization::operator=(rhs); } #endif