Changeset 2385 for trunk/Clp/src/ClpDynamicMatrix.cpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (19 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpDynamicMatrix.cpp
r2300 r2385 3 3 // Corporation and others. All Rights Reserved. 4 4 // This code is licensed under the terms of the Eclipse Public License (EPL). 5 6 5 7 6 #include <cstdio> … … 27 26 // Default Constructor 28 27 // 29 ClpDynamicMatrix::ClpDynamicMatrix 30 : ClpPackedMatrix(),31 sumDualInfeasibilities_(0.0),32 sumPrimalInfeasibilities_(0.0),33 sumOfRelaxedDualInfeasibilities_(0.0),34 sumOfRelaxedPrimalInfeasibilities_(0.0),35 savedBestGubDual_(0.0),36 savedBestSet_(0),37 backToPivotRow_(NULL),38 keyVariable_(NULL),39 toIndex_(NULL),40 fromIndex_(NULL),41 numberSets_(0),42 numberActiveSets_(0),43 objectiveOffset_(0.0),44 lowerSet_(NULL),45 upperSet_(NULL),46 status_(NULL),47 model_(NULL),48 firstAvailable_(0),49 firstAvailableBefore_(0),50 firstDynamic_(0),51 lastDynamic_(0),52 numberStaticRows_(0),53 numberElements_(0),54 numberDualInfeasibilities_(0),55 numberPrimalInfeasibilities_(0),56 noCheck_(1),57 infeasibilityWeight_(0.0),58 numberGubColumns_(0),59 maximumGubColumns_(0),60 maximumElements_(0),61 startSet_(NULL),62 next_(NULL),63 startColumn_(NULL),64 row_(NULL),65 element_(NULL),66 cost_(NULL),67 id_(NULL),68 dynamicStatus_(NULL),69 columnLower_(NULL),70 28 ClpDynamicMatrix::ClpDynamicMatrix() 29 : ClpPackedMatrix() 30 , sumDualInfeasibilities_(0.0) 31 , sumPrimalInfeasibilities_(0.0) 32 , sumOfRelaxedDualInfeasibilities_(0.0) 33 , sumOfRelaxedPrimalInfeasibilities_(0.0) 34 , savedBestGubDual_(0.0) 35 , savedBestSet_(0) 36 , backToPivotRow_(NULL) 37 , keyVariable_(NULL) 38 , toIndex_(NULL) 39 , fromIndex_(NULL) 40 , numberSets_(0) 41 , numberActiveSets_(0) 42 , objectiveOffset_(0.0) 43 , lowerSet_(NULL) 44 , upperSet_(NULL) 45 , status_(NULL) 46 , model_(NULL) 47 , firstAvailable_(0) 48 , firstAvailableBefore_(0) 49 , firstDynamic_(0) 50 , lastDynamic_(0) 51 , numberStaticRows_(0) 52 , numberElements_(0) 53 , numberDualInfeasibilities_(0) 54 , numberPrimalInfeasibilities_(0) 55 , noCheck_(1) 56 , infeasibilityWeight_(0.0) 57 , numberGubColumns_(0) 58 , maximumGubColumns_(0) 59 , maximumElements_(0) 60 , startSet_(NULL) 61 , next_(NULL) 62 , startColumn_(NULL) 63 , row_(NULL) 64 , element_(NULL) 65 , cost_(NULL) 66 , id_(NULL) 67 , dynamicStatus_(NULL) 68 , columnLower_(NULL) 69 , columnUpper_(NULL) 71 70 { 72 71 setType(15); 73 72 } 74 73 … … 76 75 // Copy constructor 77 76 // 78 ClpDynamicMatrix::ClpDynamicMatrix (const ClpDynamicMatrix &rhs)79 77 ClpDynamicMatrix::ClpDynamicMatrix(const ClpDynamicMatrix &rhs) 78 : ClpPackedMatrix(rhs) 80 79 { 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 status_ = ClpCopyOfArray(rhs.status_, static_cast<int>(2*numberSets_+4*sizeof(int)));97 98 sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;99 100 101 102 103 104 105 106 107 108 109 110 111 112 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);113 114 115 116 117 118 119 120 121 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);80 objectiveOffset_ = rhs.objectiveOffset_; 81 numberSets_ = rhs.numberSets_; 82 numberActiveSets_ = rhs.numberActiveSets_; 83 firstAvailable_ = rhs.firstAvailable_; 84 firstAvailableBefore_ = rhs.firstAvailableBefore_; 85 firstDynamic_ = rhs.firstDynamic_; 86 lastDynamic_ = rhs.lastDynamic_; 87 numberStaticRows_ = rhs.numberStaticRows_; 88 numberElements_ = rhs.numberElements_; 89 backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_); 90 keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); 91 toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); 92 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1  numberStaticRows_); 93 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_); 94 upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_); 95 status_ = ClpCopyOfArray(rhs.status_, static_cast< int >(2 * numberSets_ + 4 * sizeof(int))); 96 model_ = rhs.model_; 97 sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; 98 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; 99 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; 100 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; 101 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; 102 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; 103 savedBestGubDual_ = rhs.savedBestGubDual_; 104 savedBestSet_ = rhs.savedBestSet_; 105 noCheck_ = rhs.noCheck_; 106 infeasibilityWeight_ = rhs.infeasibilityWeight_; 107 // Now secondary data 108 numberGubColumns_ = rhs.numberGubColumns_; 109 maximumGubColumns_ = rhs.maximumGubColumns_; 110 maximumElements_ = rhs.maximumElements_; 111 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_ + 1); 112 next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_); 113 startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1); 114 row_ = ClpCopyOfArray(rhs.row_, maximumElements_); 115 element_ = ClpCopyOfArray(rhs.element_, maximumElements_); 116 cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_); 117 id_ = ClpCopyOfArray(rhs.id_, lastDynamic_  firstDynamic_); 118 columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_); 119 columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_); 120 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2 * maximumGubColumns_); 122 121 } 123 122 124 123 /* This is the real constructor*/ 125 ClpDynamicMatrix::ClpDynamicMatrix(ClpSimplex * 126 int numberGubColumns, const int *starts,127 const double * lower, const double *upper,128 const CoinBigIndex * startColumn, const int *row,129 const double * element, const double *cost,130 const double * columnLower, const double *columnUpper,131 const unsigned char *status,132 const unsigned char *dynamicStatus)133 124 ClpDynamicMatrix::ClpDynamicMatrix(ClpSimplex *model, int numberSets, 125 int numberGubColumns, const int *starts, 126 const double *lower, const double *upper, 127 const CoinBigIndex *startColumn, const int *row, 128 const double *element, const double *cost, 129 const double *columnLower, const double *columnUpper, 130 const unsigned char *status, 131 const unsigned char *dynamicStatus) 132 : ClpPackedMatrix() 134 133 { 135 setType(15); 136 objectiveOffset_ = model>objectiveOffset(); 137 model_ = model; 138 numberSets_ = numberSets; 139 numberGubColumns_ = numberGubColumns; 140 maximumGubColumns_ = numberGubColumns_; 141 if (numberGubColumns_) 142 maximumElements_ = startColumn[numberGubColumns_]; 143 else 144 maximumElements_ = 0; 145 startSet_ = new int [numberSets_+1]; 146 next_ = new int [maximumGubColumns_]; 147 // fill in startSet and next 148 int iSet; 149 if (numberGubColumns_) { 150 for (iSet = 0; iSet < numberSets_; iSet++) { 151 int first = starts[iSet]; 152 int last = starts[iSet+1]  1; 153 startSet_[iSet] = first; 154 for (int i = first; i < last; i++) 155 next_[i] = i + 1; 156 next_[last] = iSet  1; 157 } 158 startSet_[numberSets_] = starts[numberSets_]; 159 } 160 int numberColumns = model>numberColumns(); 161 int numberRows = model>numberRows(); 162 numberStaticRows_ = numberRows; 163 savedBestGubDual_ = 0.0; 164 savedBestSet_ = 0; 165 // Number of columns needed 166 int frequency = model>factorizationFrequency(); 167 int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4; 168 // But we may have two per row + one for incoming (make it two) 169 numberGubInSmall = CoinMax(2*numberRows+2,numberGubInSmall); 170 // for small problems this could be too big 171 //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_); 172 int numberNeeded = numberGubInSmall + numberColumns; 173 firstAvailable_ = numberColumns; 174 firstAvailableBefore_ = firstAvailable_; 175 firstDynamic_ = numberColumns; 176 lastDynamic_ = numberNeeded; 177 startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1); 178 if (!numberGubColumns_) { 179 if (!startColumn_) 180 startColumn_ = new CoinBigIndex [1]; 181 startColumn_[0] = 0; 182 } 183 CoinBigIndex numberElements = startColumn_[numberGubColumns_]; 184 row_ = ClpCopyOfArray(row, numberElements); 185 element_ = new double[numberElements]; 186 CoinBigIndex i; 187 for (i = 0; i < numberElements; i++) 188 element_[i] = element[i]; 189 cost_ = new double[numberGubColumns_]; 190 for (i = 0; i < numberGubColumns_; i++) { 191 cost_[i] = cost[i]; 192 // I don't think I need sorted but ... 193 CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i+1], element_ + startColumn_[i]); 194 } 195 if (columnLower) { 196 columnLower_ = new double[numberGubColumns_]; 197 for (i = 0; i < numberGubColumns_; i++) 198 columnLower_[i] = columnLower[i]; 199 } else { 200 columnLower_ = NULL; 201 } 202 if (columnUpper) { 203 columnUpper_ = new double[numberGubColumns_]; 204 for (i = 0; i < numberGubColumns_; i++) 205 columnUpper_[i] = columnUpper[i]; 206 } else { 207 columnUpper_ = NULL; 208 } 209 lowerSet_ = new double[numberSets_]; 210 for (i = 0; i < numberSets_; i++) { 211 if (lower[i] > 1.0e20) 212 lowerSet_[i] = lower[i]; 213 else 214 lowerSet_[i] = 1.0e30; 215 } 216 upperSet_ = new double[numberSets_]; 217 for (i = 0; i < numberSets_; i++) { 218 if (upper[i] < 1.0e20) 219 upperSet_[i] = upper[i]; 220 else 221 upperSet_[i] = 1.0e30; 222 } 223 id_ = new int[numberGubInSmall]; 224 for (i = 0; i < numberGubInSmall; i++) 225 id_[i] = 1; 226 ClpPackedMatrix* originalMatrixA = 227 dynamic_cast< ClpPackedMatrix*>(model>clpMatrix()); 228 assert (originalMatrixA); 229 CoinPackedMatrix * originalMatrix = originalMatrixA>getPackedMatrix(); 230 originalMatrixA>setMatrixNull(); // so can be deleted safely 231 // guess how much space needed 232 double guess = numberElements; 233 guess /= static_cast<double> (numberColumns); 234 guess *= 2 * numberGubInSmall; 235 numberElements_ = static_cast<int> (guess); 236 numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix>getNumElements(); 237 matrix_ = originalMatrix; 238 //delete originalMatrixA; 239 flags_ &= ~1; 240 // resize model (matrix stays same) 241 // modify frequency 242 if (frequency>=50) 243 frequency = 50+(frequency50)/2; 244 int newRowSize = numberRows + CoinMin(numberSets_, frequency+numberRows) + 1; 245 model>resize(newRowSize, numberNeeded); 246 for (i = numberRows; i < newRowSize; i++) 247 model>setRowStatus(i, ClpSimplex::basic); 248 if (columnUpper_) { 249 // set all upper bounds so we have enough space 250 double * columnUpper = model>columnUpper(); 251 for(i = firstDynamic_; i < lastDynamic_; i++) 252 columnUpper[i] = 1.0e10; 253 } 254 // resize matrix 255 // extra 1 is so can keep number of elements handy 256 originalMatrix>reserve(numberNeeded, numberElements_, true); 257 originalMatrix>reserve(numberNeeded + 1, numberElements_, false); 258 originalMatrix>getMutableVectorStarts()[numberColumns] = originalMatrix>getNumElements(); 259 originalMatrix>setDimensions(newRowSize, 1); 260 numberActiveColumns_ = firstDynamic_; 261 // redo number of columns 262 numberColumns = matrix_>getNumCols(); 263 backToPivotRow_ = new int[numberNeeded]; 264 keyVariable_ = new int[numberSets_]; 265 if (status) { 266 status_ = ClpCopyOfArray(status, static_cast<int>(2*numberSets_+4*sizeof(int))); 267 assert (dynamicStatus); 268 dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_); 269 } else { 270 status_ = new unsigned char [2*numberSets_+4*sizeof(int)]; 271 memset(status_, 0, numberSets_); 272 int i; 273 for (i = 0; i < numberSets_; i++) { 274 // make slack key 275 setStatus(i, ClpSimplex::basic); 276 } 277 dynamicStatus_ = new unsigned char [2*numberGubColumns_]; 278 memset(dynamicStatus_, 0, numberGubColumns_); // for clarity 279 for (i = 0; i < numberGubColumns_; i++) 280 setDynamicStatus(i, atLowerBound); 281 } 282 toIndex_ = new int[numberSets_]; 283 for (iSet = 0; iSet < numberSets_; iSet++) 284 toIndex_[iSet] = 1; 285 fromIndex_ = new int [newRowSizenumberStaticRows_+1]; 286 numberActiveSets_ = 0; 287 rhsOffset_ = NULL; 288 if (numberGubColumns_) { 289 if (!status) { 290 gubCrash(); 291 } else { 292 initialProblem(); 293 } 294 } 295 noCheck_ = 1; 296 infeasibilityWeight_ = 0.0; 134 setType(15); 135 objectiveOffset_ = model>objectiveOffset(); 136 model_ = model; 137 numberSets_ = numberSets; 138 numberGubColumns_ = numberGubColumns; 139 maximumGubColumns_ = numberGubColumns_; 140 if (numberGubColumns_) 141 maximumElements_ = startColumn[numberGubColumns_]; 142 else 143 maximumElements_ = 0; 144 startSet_ = new int[numberSets_ + 1]; 145 next_ = new int[maximumGubColumns_]; 146 // fill in startSet and next 147 int iSet; 148 if (numberGubColumns_) { 149 for (iSet = 0; iSet < numberSets_; iSet++) { 150 int first = starts[iSet]; 151 int last = starts[iSet + 1]  1; 152 startSet_[iSet] = first; 153 for (int i = first; i < last; i++) 154 next_[i] = i + 1; 155 next_[last] = iSet  1; 156 } 157 startSet_[numberSets_] = starts[numberSets_]; 158 } 159 int numberColumns = model>numberColumns(); 160 int numberRows = model>numberRows(); 161 numberStaticRows_ = numberRows; 162 savedBestGubDual_ = 0.0; 163 savedBestSet_ = 0; 164 // Number of columns needed 165 int frequency = model>factorizationFrequency(); 166 int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4; 167 // But we may have two per row + one for incoming (make it two) 168 numberGubInSmall = CoinMax(2 * numberRows + 2, numberGubInSmall); 169 // for small problems this could be too big 170 //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_); 171 int numberNeeded = numberGubInSmall + numberColumns; 172 firstAvailable_ = numberColumns; 173 firstAvailableBefore_ = firstAvailable_; 174 firstDynamic_ = numberColumns; 175 lastDynamic_ = numberNeeded; 176 startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1); 177 if (!numberGubColumns_) { 178 if (!startColumn_) 179 startColumn_ = new CoinBigIndex[1]; 180 startColumn_[0] = 0; 181 } 182 CoinBigIndex numberElements = startColumn_[numberGubColumns_]; 183 row_ = ClpCopyOfArray(row, numberElements); 184 element_ = new double[numberElements]; 185 CoinBigIndex i; 186 for (i = 0; i < numberElements; i++) 187 element_[i] = element[i]; 188 cost_ = new double[numberGubColumns_]; 189 for (i = 0; i < numberGubColumns_; i++) { 190 cost_[i] = cost[i]; 191 // I don't think I need sorted but ... 192 CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i + 1], element_ + startColumn_[i]); 193 } 194 if (columnLower) { 195 columnLower_ = new double[numberGubColumns_]; 196 for (i = 0; i < numberGubColumns_; i++) 197 columnLower_[i] = columnLower[i]; 198 } else { 199 columnLower_ = NULL; 200 } 201 if (columnUpper) { 202 columnUpper_ = new double[numberGubColumns_]; 203 for (i = 0; i < numberGubColumns_; i++) 204 columnUpper_[i] = columnUpper[i]; 205 } else { 206 columnUpper_ = NULL; 207 } 208 lowerSet_ = new double[numberSets_]; 209 for (i = 0; i < numberSets_; i++) { 210 if (lower[i] > 1.0e20) 211 lowerSet_[i] = lower[i]; 212 else 213 lowerSet_[i] = 1.0e30; 214 } 215 upperSet_ = new double[numberSets_]; 216 for (i = 0; i < numberSets_; i++) { 217 if (upper[i] < 1.0e20) 218 upperSet_[i] = upper[i]; 219 else 220 upperSet_[i] = 1.0e30; 221 } 222 id_ = new int[numberGubInSmall]; 223 for (i = 0; i < numberGubInSmall; i++) 224 id_[i] = 1; 225 ClpPackedMatrix *originalMatrixA = dynamic_cast< ClpPackedMatrix * >(model>clpMatrix()); 226 assert(originalMatrixA); 227 CoinPackedMatrix *originalMatrix = originalMatrixA>getPackedMatrix(); 228 originalMatrixA>setMatrixNull(); // so can be deleted safely 229 // guess how much space needed 230 double guess = numberElements; 231 guess /= static_cast< double >(numberColumns); 232 guess *= 2 * numberGubInSmall; 233 numberElements_ = static_cast< int >(guess); 234 numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix>getNumElements(); 235 matrix_ = originalMatrix; 236 //delete originalMatrixA; 237 flags_ &= ~1; 238 // resize model (matrix stays same) 239 // modify frequency 240 if (frequency >= 50) 241 frequency = 50 + (frequency  50) / 2; 242 int newRowSize = numberRows + CoinMin(numberSets_, frequency + numberRows) + 1; 243 model>resize(newRowSize, numberNeeded); 244 for (i = numberRows; i < newRowSize; i++) 245 model>setRowStatus(i, ClpSimplex::basic); 246 if (columnUpper_) { 247 // set all upper bounds so we have enough space 248 double *columnUpper = model>columnUpper(); 249 for (i = firstDynamic_; i < lastDynamic_; i++) 250 columnUpper[i] = 1.0e10; 251 } 252 // resize matrix 253 // extra 1 is so can keep number of elements handy 254 originalMatrix>reserve(numberNeeded, numberElements_, true); 255 originalMatrix>reserve(numberNeeded + 1, numberElements_, false); 256 originalMatrix>getMutableVectorStarts()[numberColumns] = originalMatrix>getNumElements(); 257 originalMatrix>setDimensions(newRowSize, 1); 258 numberActiveColumns_ = firstDynamic_; 259 // redo number of columns 260 numberColumns = matrix_>getNumCols(); 261 backToPivotRow_ = new int[numberNeeded]; 262 keyVariable_ = new int[numberSets_]; 263 if (status) { 264 status_ = ClpCopyOfArray(status, static_cast< int >(2 * numberSets_ + 4 * sizeof(int))); 265 assert(dynamicStatus); 266 dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2 * numberGubColumns_); 267 } else { 268 status_ = new unsigned char[2 * numberSets_ + 4 * sizeof(int)]; 269 memset(status_, 0, numberSets_); 270 int i; 271 for (i = 0; i < numberSets_; i++) { 272 // make slack key 273 setStatus(i, ClpSimplex::basic); 274 } 275 dynamicStatus_ = new unsigned char[2 * numberGubColumns_]; 276 memset(dynamicStatus_, 0, numberGubColumns_); // for clarity 277 for (i = 0; i < numberGubColumns_; i++) 278 setDynamicStatus(i, atLowerBound); 279 } 280 toIndex_ = new int[numberSets_]; 281 for (iSet = 0; iSet < numberSets_; iSet++) 282 toIndex_[iSet] = 1; 283 fromIndex_ = new int[newRowSize  numberStaticRows_ + 1]; 284 numberActiveSets_ = 0; 285 rhsOffset_ = NULL; 286 if (numberGubColumns_) { 287 if (!status) { 288 gubCrash(); 289 } else { 290 initialProblem(); 291 } 292 } 293 noCheck_ = 1; 294 infeasibilityWeight_ = 0.0; 297 295 } 298 296 … … 300 298 // Destructor 301 299 // 302 ClpDynamicMatrix::~ClpDynamicMatrix 300 ClpDynamicMatrix::~ClpDynamicMatrix() 303 301 { 304 delete[] backToPivotRow_;305 delete[] keyVariable_;306 delete[] toIndex_;307 delete[] fromIndex_;308 delete[] lowerSet_;309 delete[] upperSet_;310 delete[] status_;311 delete[] startSet_;312 delete[] next_;313 delete[] startColumn_;314 delete[] row_;315 delete[] element_;316 delete[] cost_;317 delete[] id_;318 delete[] dynamicStatus_;319 delete[] columnLower_;320 delete[] columnUpper_;302 delete[] backToPivotRow_; 303 delete[] keyVariable_; 304 delete[] toIndex_; 305 delete[] fromIndex_; 306 delete[] lowerSet_; 307 delete[] upperSet_; 308 delete[] status_; 309 delete[] startSet_; 310 delete[] next_; 311 delete[] startColumn_; 312 delete[] row_; 313 delete[] element_; 314 delete[] cost_; 315 delete[] id_; 316 delete[] dynamicStatus_; 317 delete[] columnLower_; 318 delete[] columnUpper_; 321 319 } 322 320 … … 325 323 // 326 324 ClpDynamicMatrix & 327 ClpDynamicMatrix::operator=(const ClpDynamicMatrix &rhs)325 ClpDynamicMatrix::operator=(const ClpDynamicMatrix &rhs) 328 326 { 329 330 331 delete[] backToPivotRow_;332 delete[] keyVariable_;333 delete[] toIndex_;334 delete[] fromIndex_;335 delete[] lowerSet_;336 delete[] upperSet_;337 delete[] status_;338 delete[] startSet_;339 delete[] next_;340 delete[] startColumn_;341 delete[] row_;342 delete[] element_;343 delete[] cost_;344 delete[] id_;345 delete[] dynamicStatus_;346 delete[] columnLower_;347 delete[] columnUpper_;348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 status_ = ClpCopyOfArray(rhs.status_, static_cast<int>(2*numberSets_+4*sizeof(int)));364 365 sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;366 367 368 369 370 371 372 373 374 375 376 377 378 379 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);380 381 382 383 384 385 386 387 388 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);389 390 327 if (this != &rhs) { 328 ClpPackedMatrix::operator=(rhs); 329 delete[] backToPivotRow_; 330 delete[] keyVariable_; 331 delete[] toIndex_; 332 delete[] fromIndex_; 333 delete[] lowerSet_; 334 delete[] upperSet_; 335 delete[] status_; 336 delete[] startSet_; 337 delete[] next_; 338 delete[] startColumn_; 339 delete[] row_; 340 delete[] element_; 341 delete[] cost_; 342 delete[] id_; 343 delete[] dynamicStatus_; 344 delete[] columnLower_; 345 delete[] columnUpper_; 346 objectiveOffset_ = rhs.objectiveOffset_; 347 numberSets_ = rhs.numberSets_; 348 numberActiveSets_ = rhs.numberActiveSets_; 349 firstAvailable_ = rhs.firstAvailable_; 350 firstAvailableBefore_ = rhs.firstAvailableBefore_; 351 firstDynamic_ = rhs.firstDynamic_; 352 lastDynamic_ = rhs.lastDynamic_; 353 numberStaticRows_ = rhs.numberStaticRows_; 354 numberElements_ = rhs.numberElements_; 355 backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_); 356 keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); 357 toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); 358 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1  numberStaticRows_); 359 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_); 360 upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_); 361 status_ = ClpCopyOfArray(rhs.status_, static_cast< int >(2 * numberSets_ + 4 * sizeof(int))); 362 model_ = rhs.model_; 363 sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; 364 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; 365 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; 366 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; 367 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; 368 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; 369 savedBestGubDual_ = rhs.savedBestGubDual_; 370 savedBestSet_ = rhs.savedBestSet_; 371 noCheck_ = rhs.noCheck_; 372 infeasibilityWeight_ = rhs.infeasibilityWeight_; 373 // Now secondary data 374 numberGubColumns_ = rhs.numberGubColumns_; 375 maximumGubColumns_ = rhs.maximumGubColumns_; 376 maximumElements_ = rhs.maximumElements_; 377 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_ + 1); 378 next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_); 379 startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1); 380 row_ = ClpCopyOfArray(rhs.row_, maximumElements_); 381 element_ = ClpCopyOfArray(rhs.element_, maximumElements_); 382 cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_); 383 id_ = ClpCopyOfArray(rhs.id_, lastDynamic_  firstDynamic_); 384 columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_); 385 columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_); 386 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2 * maximumGubColumns_); 387 } 388 return *this; 391 389 } 392 390 // 393 391 // Clone 394 392 // 395 ClpMatrixBase * 393 ClpMatrixBase *ClpDynamicMatrix::clone() const 396 394 { 397 395 return new ClpDynamicMatrix(*this); 398 396 } 399 397 // Partial pricing 400 void 401 ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 402 int & bestSequence, int & numberWanted) 398 void ClpDynamicMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction, 399 int &bestSequence, int &numberWanted) 403 400 { 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 int startG2 = static_cast<int>(startFraction * numberSets_);419 int endG2 = static_cast<int>(endFraction * numberSets_ + 0.1);420 421 422 423 424 double *reducedCost = model>djRegion();425 const double *duals = model>dualRowSolution();426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 401 numberWanted = currentWanted_; 402 assert(!model>rowScale()); 403 if (numberSets_) { 404 // Do packed part before gub 405 // always??? 406 //printf("normal packed price  start %d end %d (passed end %d, first %d)\n", 407 ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted); 408 } else { 409 // no gub 410 ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted); 411 return; 412 } 413 if (numberWanted > 0) { 414 // and do some proportion of full set 415 int startG2 = static_cast< int >(startFraction * numberSets_); 416 int endG2 = static_cast< int >(endFraction * numberSets_ + 0.1); 417 endG2 = CoinMin(endG2, numberSets_); 418 //printf("gub price  set start %d end %d\n", 419 // startG2,endG2); 420 double tolerance = model>currentDualTolerance(); 421 double *reducedCost = model>djRegion(); 422 const double *duals = model>dualRowSolution(); 423 double bestDj; 424 int numberRows = model>numberRows(); 425 int slackOffset = lastDynamic_ + numberRows; 426 int structuralOffset = slackOffset + numberSets_; 427 // If nothing found yet can go all the way to end 428 int endAll = endG2; 429 if (bestSequence < 0 && !startG2) 430 endAll = numberSets_; 431 if (bestSequence >= 0) { 432 if (bestSequence != savedBestSequence_) 433 bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent 434 else 435 bestDj = savedBestDj_; 436 } else { 437 bestDj = tolerance; 438 } 439 int saveSequence = bestSequence; 440 double djMod = 0.0; 441 double bestDjMod = 0.0; 442 //printf("iteration %d start %d end %d  wanted %d\n",model>numberIterations(), 443 // startG2,endG2,numberWanted); 444 int bestSet = 1; 448 445 #if 0 449 446 // make sure first available is clean (in case last iteration rejected) … … 457 454 } 458 455 #endif 459 int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_; 460 int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_; 461 for (int iSet = startG2; iSet < endAll; iSet++) { 462 if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) { 463 // give up 464 numberWanted = 0; 465 break; 466 } else if (iSet == endG2 && bestSequence >= 0) { 467 break; 468 } 469 int gubRow = toIndex_[iSet]; 470 if (gubRow >= 0) { 471 djMod = duals[gubRow+numberStaticRows_]; // have I got sign right? 472 } else { 473 int iBasic = keyVariable_[iSet]; 474 if (iBasic >= maximumGubColumns_) { 475 djMod = 0.0; // set not in 476 } else { 477 // get dj without 478 djMod = 0.0; 479 for (CoinBigIndex j = startColumn_[iBasic]; 480 j < startColumn_[iBasic+1]; j++) { 481 int jRow = row_[j]; 482 djMod = duals[jRow] * element_[j]; 483 } 484 djMod += cost_[iBasic]; 485 // See if gub slack possible  dj is djMod 486 if (getStatus(iSet) == ClpSimplex::atLowerBound) { 487 double value = djMod; 488 if (value > tolerance) { 489 numberWanted; 490 if (value > bestDj) { 491 // check flagged variable and correct dj 492 if (!flagged(iSet)) { 493 bestDj = value; 494 bestSequence = slackOffset + iSet; 495 bestDjMod = djMod; 496 bestSet = iSet; 497 } else { 498 // just to make sure we don't exit before got something 499 numberWanted++; 500 abort(); 501 } 502 } 503 } 504 } else if (getStatus(iSet) == ClpSimplex::atUpperBound) { 505 double value = djMod; 506 if (value > tolerance) { 507 numberWanted; 508 if (value > bestDj) { 509 // check flagged variable and correct dj 510 if (!flagged(iSet)) { 511 bestDj = value; 512 bestSequence = slackOffset + iSet; 513 bestDjMod = djMod; 514 bestSet = iSet; 515 } else { 516 // just to make sure we don't exit before got something 517 numberWanted++; 518 abort(); 519 } 520 } 521 } 522 } 523 } 524 } 525 int iSequence = startSet_[iSet]; 526 while (iSequence >= 0) { 527 DynamicStatus status = getDynamicStatus(iSequence); 528 if (status == atLowerBound  status == atUpperBound) { 529 double value = cost_[iSequence]  djMod; 530 for (CoinBigIndex j = startColumn_[iSequence]; 531 j < startColumn_[iSequence+1]; j++) { 532 int jRow = row_[j]; 533 value = duals[jRow] * element_[j]; 534 } 535 // change sign if at lower bound 536 if (status == atLowerBound) 537 value = value; 538 if (value > tolerance) { 539 numberWanted; 540 if (value > bestDj) { 541 // check flagged variable and correct dj 542 if (!flagged(iSequence)) { 543 if (false/*status == atLowerBound 544 &&keyValue(iSet)<1.0e7*/) { 545 // can't come in because 546 // of ones at ub 547 numberWanted++; 548 } else { 549 550 bestDj = value; 551 bestSequence = structuralOffset + iSequence; 552 bestDjMod = djMod; 553 bestSet = iSet; 554 } 555 } else { 556 // just to make sure we don't exit before got something 557 numberWanted++; 558 } 559 } 560 } 561 } 562 iSequence = next_[iSequence]; //onto next in set 563 } 564 if (numberWanted <= 0) { 565 numberWanted = 0; 566 break; 567 } 568 } 569 if (bestSequence != saveSequence) { 570 savedBestGubDual_ = bestDjMod; 571 savedBestDj_ = bestDj; 572 savedBestSequence_ = bestSequence; 573 savedBestSet_ = bestSet; 574 } 575 // See if may be finished 576 if (!startG2 && bestSequence < 0) 577 infeasibilityWeight_ = model_>infeasibilityCost(); 578 else if (bestSequence >= 0) 579 infeasibilityWeight_ = 1.0; 580 } 581 currentWanted_ = numberWanted; 456 int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_; 457 int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_; 458 for (int iSet = startG2; iSet < endAll; iSet++) { 459 if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) { 460 // give up 461 numberWanted = 0; 462 break; 463 } else if (iSet == endG2 && bestSequence >= 0) { 464 break; 465 } 466 int gubRow = toIndex_[iSet]; 467 if (gubRow >= 0) { 468 djMod = duals[gubRow + numberStaticRows_]; // have I got sign right? 469 } else { 470 int iBasic = keyVariable_[iSet]; 471 if (iBasic >= maximumGubColumns_) { 472 djMod = 0.0; // set not in 473 } else { 474 // get dj without 475 djMod = 0.0; 476 for (CoinBigIndex j = startColumn_[iBasic]; 477 j < startColumn_[iBasic + 1]; j++) { 478 int jRow = row_[j]; 479 djMod = duals[jRow] * element_[j]; 480 } 481 djMod += cost_[iBasic]; 482 // See if gub slack possible  dj is djMod 483 if (getStatus(iSet) == ClpSimplex::atLowerBound) { 484 double value = djMod; 485 if (value > tolerance) { 486 numberWanted; 487 if (value > bestDj) { 488 // check flagged variable and correct dj 489 if (!flagged(iSet)) { 490 bestDj = value; 491 bestSequence = slackOffset + iSet; 492 bestDjMod = djMod; 493 bestSet = iSet; 494 } else { 495 // just to make sure we don't exit before got something 496 numberWanted++; 497 abort(); 498 } 499 } 500 } 501 } else if (getStatus(iSet) == ClpSimplex::atUpperBound) { 502 double value = djMod; 503 if (value > tolerance) { 504 numberWanted; 505 if (value > bestDj) { 506 // check flagged variable and correct dj 507 if (!flagged(iSet)) { 508 bestDj = value; 509 bestSequence = slackOffset + iSet; 510 bestDjMod = djMod; 511 bestSet = iSet; 512 } else { 513 // just to make sure we don't exit before got something 514 numberWanted++; 515 abort(); 516 } 517 } 518 } 519 } 520 } 521 } 522 int iSequence = startSet_[iSet]; 523 while (iSequence >= 0) { 524 DynamicStatus status = getDynamicStatus(iSequence); 525 if (status == atLowerBound  status == atUpperBound) { 526 double value = cost_[iSequence]  djMod; 527 for (CoinBigIndex j = startColumn_[iSequence]; 528 j < startColumn_[iSequence + 1]; j++) { 529 int jRow = row_[j]; 530 value = duals[jRow] * element_[j]; 531 } 532 // change sign if at lower bound 533 if (status == atLowerBound) 534 value = value; 535 if (value > tolerance) { 536 numberWanted; 537 if (value > bestDj) { 538 // check flagged variable and correct dj 539 if (!flagged(iSequence)) { 540 if (false /*status == atLowerBound 541 &&keyValue(iSet)<1.0e7*/ 542 ) { 543 // can't come in because 544 // of ones at ub 545 numberWanted++; 546 } else { 547 548 bestDj = value; 549 bestSequence = structuralOffset + iSequence; 550 bestDjMod = djMod; 551 bestSet = iSet; 552 } 553 } else { 554 // just to make sure we don't exit before got something 555 numberWanted++; 556 } 557 } 558 } 559 } 560 iSequence = next_[iSequence]; //onto next in set 561 } 562 if (numberWanted <= 0) { 563 numberWanted = 0; 564 break; 565 } 566 } 567 if (bestSequence != saveSequence) { 568 savedBestGubDual_ = bestDjMod; 569 savedBestDj_ = bestDj; 570 savedBestSequence_ = bestSequence; 571 savedBestSet_ = bestSet; 572 } 573 // See if may be finished 574 if (!startG2 && bestSequence < 0) 575 infeasibilityWeight_ = model_>infeasibilityCost(); 576 else if (bestSequence >= 0) 577 infeasibilityWeight_ = 1.0; 578 } 579 currentWanted_ = numberWanted; 582 580 } 583 581 /* Returns effective RHS if it is being used. This is used for long problems … … 585 583 expensive. This may recompute */ 586 584 double * 587 ClpDynamicMatrix::rhsOffset(ClpSimplex * 588 585 ClpDynamicMatrix::rhsOffset(ClpSimplex *model, bool forceRefresh, 586 bool 589 587 #ifdef CLP_DEBUG 590 588 check 591 589 #endif 592 590 ) 593 591 { 594 595 596 597 592 // forceRefresh=true;printf("take out forceRefresh\n"); 593 if (!model_>numberIterations()) 594 forceRefresh = true; 595 //check=false; 598 596 #ifdef CLP_DEBUG 599 double *saveE = NULL;600 601 602 603 597 double *saveE = NULL; 598 if (rhsOffset_ && check) { 599 int numberRows = model>numberRows(); 600 saveE = new double[numberRows]; 601 } 604 602 #endif 605 603 if (rhsOffset_) { 606 604 #ifdef CLP_DEBUG 607 608 609 610 double *rhs = new double[numberRows];611 612 613 614 615 const double *smallSolution = model>solutionRegion();616 const double *element = matrix_>getElements();617 const int *row = matrix_>getIndices();618 const CoinBigIndex *startColumn = matrix_>getVectorStarts();619 const int *length = matrix_>getVectorLengths();620 621 622 623 624 625 626 627 628 629 630 631 632 633 double * solution = new double[numberGubColumns_];634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 int jFull = id_[iColumnfirstDynamic_];657 658 659 660 661 662 663 664 665 666 667 668 for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {669 670 671 672 673 674 675 676 677 678 delete[] solution;679 680 681 682 683 684 685 686 687 688 689 690 691 assert(iStatus != ClpSimplex::basic);692 693 694 695 696 697 698 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 delete[] rhs;716 605 if (check) { 606 // no need  but check anyway 607 int numberRows = model>numberRows(); 608 double *rhs = new double[numberRows]; 609 int iRow; 610 int iSet; 611 CoinZeroN(rhs, numberRows); 612 // do ones at bounds before gub 613 const double *smallSolution = model>solutionRegion(); 614 const double *element = matrix_>getElements(); 615 const int *row = matrix_>getIndices(); 616 const CoinBigIndex *startColumn = matrix_>getVectorStarts(); 617 const int *length = matrix_>getVectorLengths(); 618 int iColumn; 619 double objectiveOffset = 0.0; 620 for (iColumn = 0; iColumn < firstDynamic_; iColumn++) { 621 if (model>getStatus(iColumn) != ClpSimplex::basic) { 622 double value = smallSolution[iColumn]; 623 for (CoinBigIndex j = startColumn[iColumn]; 624 j < startColumn[iColumn] + length[iColumn]; j++) { 625 int jRow = row[j]; 626 rhs[jRow] = value * element[j]; 627 } 628 } 629 } 630 if (columnLower_  columnUpper_) { 631 double *solution = new double[numberGubColumns_]; 632 for (iSet = 0; iSet < numberSets_; iSet++) { 633 int j = startSet_[iSet]; 634 while (j >= 0) { 635 double value = 0.0; 636 if (getDynamicStatus(j) != inSmall) { 637 if (getDynamicStatus(j) == atLowerBound) { 638 if (columnLower_) 639 value = columnLower_[j]; 640 } else if (getDynamicStatus(j) == atUpperBound) { 641 value = columnUpper_[j]; 642 } else if (getDynamicStatus(j) == soloKey) { 643 value = keyValue(iSet); 644 } 645 objectiveOffset += value * cost_[j]; 646 } 647 solution[j] = value; 648 j = next_[j]; //onto next in set 649 } 650 } 651 // ones in gub and in small problem 652 for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) { 653 if (model_>getStatus(iColumn) != ClpSimplex::basic) { 654 int jFull = id_[iColumn  firstDynamic_]; 655 solution[jFull] = smallSolution[iColumn]; 656 } 657 } 658 for (iSet = 0; iSet < numberSets_; iSet++) { 659 int kRow = toIndex_[iSet]; 660 if (kRow >= 0) 661 kRow += numberStaticRows_; 662 int j = startSet_[iSet]; 663 while (j >= 0) { 664 double value = solution[j]; 665 if (value) { 666 for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j + 1]; k++) { 667 int iRow = row_[k]; 668 rhs[iRow] = element_[k] * value; 669 } 670 if (kRow >= 0) 671 rhs[kRow] = value; 672 } 673 j = next_[j]; //onto next in set 674 } 675 } 676 delete[] solution; 677 } else { 678 // bounds 679 ClpSimplex::Status iStatus; 680 for (int iSet = 0; iSet < numberSets_; iSet++) { 681 int kRow = toIndex_[iSet]; 682 if (kRow < 0) { 683 int iColumn = keyVariable_[iSet]; 684 if (iColumn < maximumGubColumns_) { 685 // key is not treated as basic 686 double b = 0.0; 687 // key is structural  where is slack 688 iStatus = getStatus(iSet); 689 assert(iStatus != ClpSimplex::basic); 690 if (iStatus == ClpSimplex::atLowerBound) 691 b = lowerSet_[iSet]; 692 else 693 b = upperSet_[iSet]; 694 if (b) { 695 objectiveOffset += b * cost_[iColumn]; 696 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) { 697 int iRow = row_[j]; 698 rhs[iRow] = element_[j] * b; 699 } 700 } 701 } 702 } 703 } 704 } 705 if (fabs(model>objectiveOffset()  objectiveOffset_  objectiveOffset) > 1.0e1) 706 printf("old offset %g, true %g\n", model>objectiveOffset(), 707 objectiveOffset_  objectiveOffset); 708 for (iRow = 0; iRow < numberRows; iRow++) { 709 if (fabs(rhs[iRow]  rhsOffset_[iRow]) > 1.0e3) 710 printf("** bad effective %d  true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]); 711 } 712 CoinMemcpyN(rhs, numberRows, saveE); 713 delete[] rhs; 714 } 717 715 #endif 718 if (forceRefresh  (refreshFrequency_ && model>numberIterations() >= 719 lastRefresh_ + refreshFrequency_)) { 720 int numberRows = model>numberRows(); 721 int iSet; 722 CoinZeroN(rhsOffset_, numberRows); 723 // do ones at bounds before gub 724 const double * smallSolution = model>solutionRegion(); 725 const double * element = matrix_>getElements(); 726 const int * row = matrix_>getIndices(); 727 const CoinBigIndex * startColumn = matrix_>getVectorStarts(); 728 const int * length = matrix_>getVectorLengths(); 729 int iColumn; 730 double objectiveOffset = 0.0; 731 for (iColumn = 0; iColumn < firstDynamic_; iColumn++) { 732 if (model>getStatus(iColumn) != ClpSimplex::basic) { 733 double value = smallSolution[iColumn]; 734 for (CoinBigIndex j = startColumn[iColumn]; 735 j < startColumn[iColumn] + length[iColumn]; j++) { 736 int jRow = row[j]; 737 rhsOffset_[jRow] = value * element[j]; 738 } 739 } 740 } 741 if (columnLower_  columnUpper_) { 742 double * solution = new double [numberGubColumns_]; 743 for (iSet = 0; iSet < numberSets_; iSet++) { 744 int j = startSet_[iSet]; 745 while (j >= 0) { 746 double value = 0.0; 747 if (getDynamicStatus(j) != inSmall) { 748 if (getDynamicStatus(j) == atLowerBound) { 749 if (columnLower_) 750 value = columnLower_[j]; 751 } else if (getDynamicStatus(j) == atUpperBound) { 752 value = columnUpper_[j]; 753 assert (value<1.0e30); 754 } else if (getDynamicStatus(j) == soloKey) { 755 value = keyValue(iSet); 756 } 757 objectiveOffset += value * cost_[j]; 758 } 759 solution[j] = value; 760 j = next_[j]; //onto next in set 761 } 762 } 763 // ones in gub and in small problem 764 for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) { 765 if (model_>getStatus(iColumn) != ClpSimplex::basic) { 766 int jFull = id_[iColumnfirstDynamic_]; 767 solution[jFull] = smallSolution[iColumn]; 768 } 769 } 770 for (iSet = 0; iSet < numberSets_; iSet++) { 771 int kRow = toIndex_[iSet]; 772 if (kRow >= 0) 773 kRow += numberStaticRows_; 774 int j = startSet_[iSet]; 775 while (j >= 0) { 776 //? use DynamicStatus status = getDynamicStatus(j); 777 double value = solution[j]; 778 if (value) { 779 for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) { 780 int iRow = row_[k]; 781 rhsOffset_[iRow] = element_[k] * value; 782 } 783 if (kRow >= 0) 784 rhsOffset_[kRow] = value; 785 } 786 j = next_[j]; //onto next in set 787 } 788 } 789 delete [] solution; 790 } else { 791 // bounds 792 ClpSimplex::Status iStatus; 793 for (int iSet = 0; iSet < numberSets_; iSet++) { 794 int kRow = toIndex_[iSet]; 795 if (kRow < 0) { 796 int iColumn = keyVariable_[iSet]; 797 if (iColumn < maximumGubColumns_) { 798 // key is not treated as basic 799 double b = 0.0; 800 // key is structural  where is slack 801 iStatus = getStatus(iSet); 802 assert (iStatus != ClpSimplex::basic); 803 if (iStatus == ClpSimplex::atLowerBound) 804 b = lowerSet_[iSet]; 805 else 806 b = upperSet_[iSet]; 807 if (b) { 808 objectiveOffset += b * cost_[iColumn]; 809 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) { 810 int iRow = row_[j]; 811 rhsOffset_[iRow] = element_[j] * b; 812 } 813 } 814 } 815 } 816 } 817 } 818 model>setObjectiveOffset(objectiveOffset_  objectiveOffset); 716 if (forceRefresh  (refreshFrequency_ && model>numberIterations() >= lastRefresh_ + refreshFrequency_)) { 717 int numberRows = model>numberRows(); 718 int iSet; 719 CoinZeroN(rhsOffset_, numberRows); 720 // do ones at bounds before gub 721 const double *smallSolution = model>solutionRegion(); 722 const double *element = matrix_>getElements(); 723 const int *row = matrix_>getIndices(); 724 const CoinBigIndex *startColumn = matrix_>getVectorStarts(); 725 const int *length = matrix_>getVectorLengths(); 726 int iColumn; 727 double objectiveOffset = 0.0; 728 for (iColumn = 0; iColumn < firstDynamic_; iColumn++) { 729 if (model>getStatus(iColumn) != ClpSimplex::basic) { 730 double value = smallSolution[iColumn]; 731 for (CoinBigIndex j = startColumn[iColumn]; 732 j < startColumn[iColumn] + length[iColumn]; j++) { 733 int jRow = row[j]; 734 rhsOffset_[jRow] = value * element[j]; 735 } 736 } 737 } 738 if (columnLower_  columnUpper_) { 739 double *solution = new double[numberGubColumns_]; 740 for (iSet = 0; iSet < numberSets_; iSet++) { 741 int j = startSet_[iSet]; 742 while (j >= 0) { 743 double value = 0.0; 744 if (getDynamicStatus(j) != inSmall) { 745 if (getDynamicStatus(j) == atLowerBound) { 746 if (columnLower_) 747 value = columnLower_[j]; 748 } else if (getDynamicStatus(j) == atUpperBound) { 749 value = columnUpper_[j]; 750 assert(value < 1.0e30); 751 } else if (getDynamicStatus(j) == soloKey) { 752 value = keyValue(iSet); 753 } 754 objectiveOffset += value * cost_[j]; 755 } 756 solution[j] = value; 757 j = next_[j]; //onto next in set 758 } 759 } 760 // ones in gub and in small problem 761 for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) { 762 if (model_>getStatus(iColumn) != ClpSimplex::basic) { 763 int jFull = id_[iColumn  firstDynamic_]; 764 solution[jFull] = smallSolution[iColumn]; 765 } 766 } 767 for (iSet = 0; iSet < numberSets_; iSet++) { 768 int kRow = toIndex_[iSet]; 769 if (kRow >= 0) 770 kRow += numberStaticRows_; 771 int j = startSet_[iSet]; 772 while (j >= 0) { 773 //? use DynamicStatus status = getDynamicStatus(j); 774 double value = solution[j]; 775 if (value) { 776 for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j + 1]; k++) { 777 int iRow = row_[k]; 778 rhsOffset_[iRow] = element_[k] * value; 779 } 780 if (kRow >= 0) 781 rhsOffset_[kRow] = value; 782 } 783 j = next_[j]; //onto next in set 784 } 785 } 786 delete[] solution; 787 } else { 788 // bounds 789 ClpSimplex::Status iStatus; 790 for (int iSet = 0; iSet < numberSets_; iSet++) { 791 int kRow = toIndex_[iSet]; 792 if (kRow < 0) { 793 int iColumn = keyVariable_[iSet]; 794 if (iColumn < maximumGubColumns_) { 795 // key is not treated as basic 796 double b = 0.0; 797 // key is structural  where is slack 798 iStatus = getStatus(iSet); 799 assert(iStatus != ClpSimplex::basic); 800 if (iStatus == ClpSimplex::atLowerBound) 801 b = lowerSet_[iSet]; 802 else 803 b = upperSet_[iSet]; 804 if (b) { 805 objectiveOffset += b * cost_[iColumn]; 806 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn + 1]; j++) { 807 int iRow = row_[j]; 808 rhsOffset_[iRow] = element_[j] * b; 809 } 810 } 811 } 812 } 813 } 814 } 815 model>setObjectiveOffset(objectiveOffset_  objectiveOffset); 819 816 #ifdef CLP_DEBUG 820 821 822 823 824 825 826 delete[] saveE;827 817 if (saveE) { 818 int iRow; 819 for (iRow = 0; iRow < numberRows; iRow++) { 820 if (fabs(saveE[iRow]  rhsOffset_[iRow]) > 1.0e3) 821 printf("** %d  old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]); 822 } 823 delete[] saveE; 824 } 828 825 #endif 829 830 831 832 826 lastRefresh_ = model>numberIterations(); 827 } 828 } 829 return rhsOffset_; 833 830 } 834 831 /* 835 832 update information for a pivot (and effective rhs) 836 833 */ 837 int 838 ClpDynamicMatrix::updatePivot(ClpSimplex * model, double oldInValue, double oldOutValue) 834 int ClpDynamicMatrix::updatePivot(ClpSimplex *model, double oldInValue, double oldOutValue) 839 835 { 840 836 841 // now update working model 842 int sequenceIn = model>sequenceIn(); 843 int sequenceOut = model>sequenceOut(); 844 int numberColumns = model>numberColumns(); 845 if (sequenceIn != sequenceOut && sequenceIn < numberColumns) 846 backToPivotRow_[sequenceIn] = model>pivotRow(); 847 if (sequenceIn >= firstDynamic_ && sequenceIn < numberColumns) { 848 int bigSequence = id_[sequenceInfirstDynamic_]; 849 if (getDynamicStatus(bigSequence) != inSmall) { 850 firstAvailable_++; 851 setDynamicStatus(bigSequence, inSmall); 852 } 853 } 854 // make sure slack is synchronized 855 if (sequenceIn >= numberColumns + numberStaticRows_) { 856 int iDynamic = sequenceIn  numberColumns  numberStaticRows_; 857 int iSet = fromIndex_[iDynamic]; 858 setStatus(iSet, model>getStatus(sequenceIn)); 859 } 860 if (sequenceOut >= numberColumns + numberStaticRows_) { 861 int iDynamic = sequenceOut  numberColumns  numberStaticRows_; 862 int iSet = fromIndex_[iDynamic]; 863 // out may have gone through barrier  so check 864 double valueOut = model>lowerRegion()[sequenceOut]; 865 if (fabs(valueOut  static_cast<double> (lowerSet_[iSet])) < 866 fabs(valueOut  static_cast<double> (upperSet_[iSet]))) 867 setStatus(iSet, ClpSimplex::atLowerBound); 868 else 869 setStatus(iSet, ClpSimplex::atUpperBound); 870 if (lowerSet_[iSet] == upperSet_[iSet]) 871 setStatus(iSet, ClpSimplex::isFixed); 837 // now update working model 838 int sequenceIn = model>sequenceIn(); 839 int sequenceOut = model>sequenceOut(); 840 int numberColumns = model>numberColumns(); 841 if (sequenceIn != sequenceOut && sequenceIn < numberColumns) 842 backToPivotRow_[sequenceIn] = model>pivotRow(); 843 if (sequenceIn >= firstDynamic_ && sequenceIn < numberColumns) { 844 int bigSequence = id_[sequenceIn  firstDynamic_]; 845 if (getDynamicStatus(bigSequence) != inSmall) { 846 firstAvailable_++; 847 setDynamicStatus(bigSequence, inSmall); 848 } 849 } 850 // make sure slack is synchronized 851 if (sequenceIn >= numberColumns + numberStaticRows_) { 852 int iDynamic = sequenceIn  numberColumns  numberStaticRows_; 853 int iSet = fromIndex_[iDynamic]; 854 setStatus(iSet, model>getStatus(sequenceIn)); 855 } 856 if (sequenceOut >= numberColumns + numberStaticRows_) { 857 int iDynamic = sequenceOut  numberColumns  numberStaticRows_; 858 int iSet = fromIndex_[iDynamic]; 859 // out may have gone through barrier  so check 860 double valueOut = model>lowerRegion()[sequenceOut]; 861 if (fabs(valueOut  static_cast< double >(lowerSet_[iSet])) < fabs(valueOut  static_cast< double >(upperSet_[iSet]))) 862 setStatus(iSet, ClpSimplex::atLowerBound); 863 else 864 setStatus(iSet, ClpSimplex::atUpperBound); 865 if (lowerSet_[iSet] == upperSet_[iSet]) 866 setStatus(iSet, ClpSimplex::isFixed); 872 867 #if 0 873 868 if (getStatus(iSet) != model>getStatus(sequenceOut)) … … 875 870 getStatus(iSet), model>getStatus(sequenceOut)); 876 871 #endif 877 878 872 } 873 ClpMatrixBase::updatePivot(model, oldInValue, oldOutValue); 879 874 #ifdef CLP_DEBUG 880 char * inSmall = new char[numberGubColumns_];881 882 const double *solution = model>solutionRegion();883 884 885 886 887 int k = id_[ifirstDynamic_];888 889 890 891 892 893 assert(!inSmall[i]);894 delete[] inSmall;875 char *inSmall = new char[numberGubColumns_]; 876 memset(inSmall, 0, numberGubColumns_); 877 const double *solution = model>solutionRegion(); 878 for (int i = 0; i < numberGubColumns_; i++) 879 if (getDynamicStatus(i) == ClpDynamicMatrix::inSmall) 880 inSmall[i] = 1; 881 for (int i = firstDynamic_; i < firstAvailable_; i++) { 882 int k = id_[i  firstDynamic_]; 883 inSmall[k] = 0; 884 //if (k>=23289&&k<23357&&solution[i]) 885 //printf("var %d (in small %d) has value %g\n",k,i,solution[i]); 886 } 887 for (int i = 0; i < numberGubColumns_; i++) 888 assert(!inSmall[i]); 889 delete[] inSmall; 895 890 #ifndef NDEBUG 896 897 898 assert(toIndex_[iSet] == i);899 900 901 902 903 904 905 906 907 891 for (int i = 0; i < numberActiveSets_; i++) { 892 int iSet = fromIndex_[i]; 893 assert(toIndex_[iSet] == i); 894 //if (getStatus(iSet)!=model>getRowStatus(i+numberStaticRows_)) 895 //printf("*** set %d status %d, var status %d\n",iSet, 896 // getStatus(iSet),model>getRowStatus(i+numberStaticRows_)); 897 //assert (model>getRowStatus(i+numberStaticRows_)==getStatus(iSet)); 898 //if (iSet==1035) { 899 //printf("rhs for set %d (%d) is %g %g  cost %g\n",iSet,i,model>lowerRegion(0)[i+numberStaticRows_], 900 // model>upperRegion(0)[i+numberStaticRows_],model>costRegion(0)[i+numberStaticRows_]); 901 //} 902 } 908 903 #endif 909 904 #endif 910 if (numberStaticRows_+numberActiveSets_<model>numberRows())911 912 913 905 if (numberStaticRows_ + numberActiveSets_ < model>numberRows()) 906 return 0; 907 else 908 return 1; 914 909 } 915 910 /* … … 918 913 Remember to update here when settled down 919 914 */ 920 void 921 ClpDynamicMatrix::dualExpanded(ClpSimplex * model, 922 CoinIndexedVector * /*array*/, 923 double * /*other*/, int mode) 915 void ClpDynamicMatrix::dualExpanded(ClpSimplex *model, 916 CoinIndexedVector * /*array*/, 917 double * /*other*/, int mode) 924 918 { 925 switch (mode) { 926 // modify costs before transposeUpdate 927 case 0: 928 break; 929 // create duals for key variables (without check on dual infeasible) 930 case 1: 931 break; 932 // as 1 but check slacks and compute djs 933 case 2: { 934 // do pivots 935 int * pivotVariable = model>pivotVariable(); 936 int numberRows = numberStaticRows_ + numberActiveSets_; 937 int numberColumns = model>numberColumns(); 938 for (int iRow = 0; iRow < numberRows; iRow++) { 939 int iPivot = pivotVariable[iRow]; 940 if (iPivot < numberColumns) 941 backToPivotRow_[iPivot] = iRow; 942 } 943 if (noCheck_ >= 0) { 944 if (infeasibilityWeight_ != model_>infeasibilityCost()) { 945 // don't bother checking 946 sumDualInfeasibilities_ = 100.0; 947 numberDualInfeasibilities_ = 1; 948 sumOfRelaxedDualInfeasibilities_ = 100.0; 949 return; 950 } 951 } 952 // In theory we should subtract out ones we have done but .... 953 // If key slack then dual 0.0 954 // If not then slack could be dual infeasible 955 // dj for key is zero so that defines dual on set 956 int i; 957 double * dual = model>dualRowSolution(); 958 double dualTolerance = model>dualTolerance(); 959 double relaxedTolerance = dualTolerance; 960 // we can't really trust infeasibilities if there is dual error 961 double error = CoinMin(1.0e2, model>largestDualError()); 962 // allow tolerance at least slightly bigger than standard 963 relaxedTolerance = relaxedTolerance + error; 964 // but we will be using difference 965 relaxedTolerance = dualTolerance; 966 sumDualInfeasibilities_ = 0.0; 967 numberDualInfeasibilities_ = 0; 968 sumOfRelaxedDualInfeasibilities_ = 0.0; 969 for (i = 0; i < numberSets_; i++) { 970 double value = 0.0; 971 int gubRow = toIndex_[i]; 972 if (gubRow < 0) { 973 int kColumn = keyVariable_[i]; 974 if (kColumn < maximumGubColumns_) { 975 // dj without set 976 value = cost_[kColumn]; 977 for (CoinBigIndex j = startColumn_[kColumn]; 978 j < startColumn_[kColumn+1]; j++) { 979 int iRow = row_[j]; 980 value = dual[iRow] * element_[j]; 981 } 982 double infeasibility = 0.0; 983 if (getStatus(i) == ClpSimplex::atLowerBound) { 984 if (value > dualTolerance) 985 infeasibility = value  dualTolerance; 986 } else if (getStatus(i) == ClpSimplex::atUpperBound) { 987 if (value > dualTolerance) 988 infeasibility = value  dualTolerance; 989 } 990 if (infeasibility > 0.0) { 991 sumDualInfeasibilities_ += infeasibility; 992 if (infeasibility > relaxedTolerance) 993 sumOfRelaxedDualInfeasibilities_ += infeasibility; 994 numberDualInfeasibilities_ ++; 995 } 996 } 997 } else { 998 value = dual[gubRow+numberStaticRows_]; 999 } 1000 // Now subtract out from all 1001 int k = startSet_[i]; 1002 while (k >= 0) { 1003 if (getDynamicStatus(k) != inSmall) { 1004 double djValue = cost_[k]  value; 1005 for (CoinBigIndex j = startColumn_[k]; 1006 j < startColumn_[k+1]; j++) { 1007 int iRow = row_[j]; 1008 djValue = dual[iRow] * element_[j]; 1009 } 1010 double infeasibility = 0.0; 1011 if (getDynamicStatus(k) == atLowerBound) { 1012 if (djValue < dualTolerance) 1013 infeasibility = djValue  dualTolerance; 1014 } else if (getDynamicStatus(k) == atUpperBound) { 1015 // at upper bound 1016 if (djValue > dualTolerance) 1017 infeasibility = djValue  dualTolerance; 1018 } 1019 if (infeasibility > 0.0) { 1020 sumDualInfeasibilities_ += infeasibility; 1021 if (infeasibility > relaxedTolerance) 1022 sumOfRelaxedDualInfeasibilities_ += infeasibility; 1023 numberDualInfeasibilities_ ++; 1024 } 1025 } 1026 k = next_[k]; //onto next in set 1027 } 1028 } 1029 } 1030 infeasibilityWeight_ = 1.0; 1031 break; 1032 // Report on infeasibilities of key variables 1033 case 3: { 1034 model>setSumDualInfeasibilities(model>sumDualInfeasibilities() + 1035 sumDualInfeasibilities_); 1036 model>setNumberDualInfeasibilities(model>numberDualInfeasibilities() + 1037 numberDualInfeasibilities_); 1038 model>setSumOfRelaxedDualInfeasibilities(model>sumOfRelaxedDualInfeasibilities() + 1039 sumOfRelaxedDualInfeasibilities_); 1040 } 1041 break; 1042 // modify costs before transposeUpdate for partial pricing 1043 case 4: 1044 break; 1045 } 919 switch (mode) { 920 // modify costs before transposeUpdate 921 case 0: 922 break; 923 // create duals for key variables (without check on dual infeasible) 924 case 1: 925 break; 926 // as 1 but check slacks and compute djs 927 case 2: { 928 // do pivots 929 int *pivotVariable = model>pivotVariable(); 930 int numberRows = numberStaticRows_ + numberActiveSets_; 931 int numberColumns = model>numberColumns(); 932 for (int iRow = 0; iRow < numberRows; iRow++) { 933 int iPivot = pivotVariable[iRow]; 934 if (iPivot < numberColumns) 935 backToPivotRow_[iPivot] = iRow; 936 } 937 if (noCheck_ >= 0) { 938 if (infeasibilityWeight_ != model_>infeasibilityCost()) { 939 // don't bother checking 940 sumDualInfeasibilities_ = 100.0; 941 numberDualInfeasibilities_ = 1; 942 sumOfRelaxedDualInfeasibilities_ = 100.0; 943 return; 944 } 945 } 946 // In theory we should subtract out ones we have done but .... 947 // If key slack then dual 0.0 948 // If not then slack could be dual infeasible 949 // dj for key is zero so that defines dual on set 950 int i; 951 double *dual = model>dualRowSolution(); 952 double dualTolerance = model>dualTolerance(); 953 double relaxedTolerance = dualTolerance; 954 // we can't really trust infeasibilities if there is dual error 955 double error = CoinMin(1.0e2, model>largestDualError()); 956 // allow tolerance at least slightly bigger than standard 957 relaxedTolerance = relaxedTolerance + error; 958 // but we will be using difference 959 relaxedTolerance = dualTolerance; 960 sumDualInfeasibilities_ = 0.0; 961 numberDualInfeasibilities_ = 0; 962 sumOfRelaxedDualInfeasibilities_ = 0.0; 963 for (i = 0; i < numberSets_; i++) { 964 double value = 0.0; 965 int gubRow = toIndex_[i]; 966 if (gubRow < 0) { 967 int kColumn = keyVariable_[i]; 968 if (kColumn < maximumGubColumns_) { 969 // dj without set 970 value = cost_[kColumn]; 971 for (CoinBigIndex j = startColumn_[kColumn]; 972 j < startColumn_[kColumn + 1]; j++) { 973 int iRow = row_[j]; 974 value = dual[iRow] * element_[j]; 975 } 976 double infeasibility = 0.0; 977 if (getStatus(i) == ClpSimplex::atLowerBound) { 978 if (value > dualTolerance) 979 infeasibility = value  dualTolerance; 980 } else if (getStatus(i) == ClpSimplex::atUpperBound) { 981 if (value > dualTolerance) 982 infeasibility = value  dualTolerance; 983 } 984 if (infeasibility > 0.0) { 985 sumDualInfeasibilities_ += infeasibility; 986 if (infeasibility > relaxedTolerance) 987 sumOfRelaxedDualInfeasibilities_ += infeasibility; 988 numberDualInfeasibilities_++; 989 } 990 } 991 } else { 992 value = dual[gubRow + numberStaticRows_]; 993 } 994 // Now subtract out from all 995 int k = startSet_[i]; 996 while (k >= 0) { 997 if (getDynamicStatus(k) != inSmall) { 998 double djValue = cost_[k]  value; 999 for (CoinBigIndex j = startColumn_[k]; 1000 j < startColumn_[k + 1]; j++) { 1001 int iRow = row_[j]; 1002 djValue = dual[iRow] * element_[j]; 1003 } 1004 double infeasibility = 0.0; 1005 if (getDynamicStatus(k) == atLowerBound) { 1006 if (djValue < dualTolerance) 1007 infeasibility = djValue  dualTolerance; 1008 } else if (getDynamicStatus(k) == atUpperBound) { 1009 // at upper bound 1010 if (djValue > dualTolerance) 1011 infeasibility = djValue  dualTolerance; 1012 } 1013 if (infeasibility > 0.0) { 1014 sumDualInfeasibilities_ += infeasibility; 1015 if (infeasibility > relaxedTolerance) 1016 sumOfRelaxedDualInfeasibilities_ += infeasibility; 1017 numberDualInfeasibilities_++; 1018 } 1019 } 1020 k = next_[k]; //onto next in set 1021 } 1022 } 1023 } 1024 infeasibilityWeight_ = 1.0; 1025 break; 1026 // Report on infeasibilities of key variables 1027 case 3: { 1028 model>setSumDualInfeasibilities(model>sumDualInfeasibilities() + sumDualInfeasibilities_); 1029 model>setNumberDualInfeasibilities(model>numberDualInfeasibilities() + numberDualInfeasibilities_); 1030 model>setSumOfRelaxedDualInfeasibilities(model>sumOfRelaxedDualInfeasibilities() + sumOfRelaxedDualInfeasibilities_); 1031 } break; 1032 // modify costs before transposeUpdate for partial pricing 1033 case 4: 1034 break; 1035 } 1046 1036 } 1047 1037 /* … … 1050 1040 Remember to update here when settled down 1051 1041 */ 1052 int 1053 ClpDynamicMatrix::generalExpanded(ClpSimplex * model, int mode, int &number) 1042 int ClpDynamicMatrix::generalExpanded(ClpSimplex *model, int mode, int &number) 1054 1043 { 1055 1044 int returnCode = 0; 1056 1045 #if 0 //ndef NDEBUG 1057 1046 { … … 1066 1055 } 1067 1056 #endif 1068 switch (mode) { 1069 // Fill in pivotVariable 1070 case 0: { 1071 // If no effective rhs  form it 1072 if (!rhsOffset_) { 1073 rhsOffset_ = new double[model>numberRows()]; 1074 rhsOffset(model, true); 1075 } 1076 int i; 1077 int numberBasic = number; 1078 int numberColumns = model>numberColumns(); 1079 // Use different array so can build from true pivotVariable_ 1080 //int * pivotVariable = model>pivotVariable(); 1081 int * pivotVariable = model>rowArray(0)>getIndices(); 1082 for (i = 0; i < numberColumns; i++) { 1083 if (model>getColumnStatus(i) == ClpSimplex::basic) 1084 pivotVariable[numberBasic++] = i; 1085 } 1086 number = numberBasic; 1087 } 1088 break; 1089 // Do initial extra rows + maximum basic 1090 case 2: { 1091 number = model>numberRows(); 1092 } 1093 break; 1094 // Before normal replaceColumn 1095 case 3: { 1096 if (numberActiveSets_ + numberStaticRows_ == model_>numberRows()) { 1097 // no space  refactorize 1098 returnCode = 4; 1099 number = 1; // say no need for normal replaceColumn 1100 } 1101 } 1102 break; 1103 // To see if can dual or primal 1104 case 4: { 1105 returnCode = 1; 1106 } 1107 break; 1108 // save status 1109 case 5: { 1110 memcpy(status_+numberSets_,status_,numberSets_); 1111 memcpy(status_+2*numberSets_,&numberActiveSets_,sizeof(int)); 1112 memcpy(dynamicStatus_+maximumGubColumns_, 1113 dynamicStatus_,maximumGubColumns_); 1114 } 1115 break; 1116 // restore status 1117 case 6: { 1118 memcpy(status_,status_+numberSets_,numberSets_); 1119 memcpy(&numberActiveSets_,status_+2*numberSets_,sizeof(int)); 1120 memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_, 1121 maximumGubColumns_); 1122 initialProblem(); 1123 } 1124 break; 1125 // unflag all variables 1126 case 8: { 1127 for (int i = 0; i < numberGubColumns_; i++) { 1128 if (flagged(i)) { 1129 unsetFlagged(i); 1130 returnCode++; 1131 } 1132 } 1133 } 1134 break; 1135 // redo costs in primal 1136 case 9: { 1137 double * cost = model>costRegion(); 1138 double * solution = model>solutionRegion(); 1139 double * columnLower = model>lowerRegion(); 1140 double * columnUpper = model>upperRegion(); 1141 int i; 1142 bool doCosts = (number & 4) != 0; 1143 bool doBounds = (number & 1) != 0; 1144 for ( i = firstDynamic_; i < firstAvailable_; i++) { 1145 int jColumn = id_[ifirstDynamic_]; 1146 if (doBounds) { 1147 if (!columnLower_ && !columnUpper_) { 1148 columnLower[i] = 0.0; 1149 columnUpper[i] = COIN_DBL_MAX; 1150 } else { 1151 if (columnLower_) 1152 columnLower[i] = columnLower_[jColumn]; 1153 else 1154 columnLower[i] = 0.0; 1155 if (columnUpper_) 1156 columnUpper[i] = columnUpper_[jColumn]; 1157 else 1158 columnUpper[i] = COIN_DBL_MAX; 1159 } 1160 } 1161 if (doCosts) { 1162 cost[i] = cost_[jColumn]; 1163 // Original bounds 1164 if (model>nonLinearCost()) 1165 model>nonLinearCost()>setOne(i, solution[i], 1166 this>columnLower(jColumn), 1167 this>columnUpper(jColumn), cost_[jColumn]); 1168 } 1169 } 1170 // and active sets 1171 for (i = 0; i < numberActiveSets_; i++ ) { 1172 int iSet = fromIndex_[i]; 1173 int iSequence = lastDynamic_ + numberStaticRows_ + i; 1174 if (doBounds) { 1175 if (lowerSet_[iSet] > 1.0e20) 1176 columnLower[iSequence] = lowerSet_[iSet]; 1177 else 1178 columnLower[iSequence] = COIN_DBL_MAX; 1179 if (upperSet_[iSet] < 1.0e20) 1180 columnUpper[iSequence] = upperSet_[iSet]; 1181 else 1182 columnUpper[iSequence] = COIN_DBL_MAX; 1183 } 1184 if (doCosts) { 1185 if (model>nonLinearCost()) { 1186 double trueLower; 1187 if (lowerSet_[iSet] > 1.0e20) 1188 trueLower = lowerSet_[iSet]; 1189 else 1190 trueLower = COIN_DBL_MAX; 1191 double trueUpper; 1192 if (upperSet_[iSet] < 1.0e20) 1193 trueUpper = upperSet_[iSet]; 1194 else 1195 trueUpper = COIN_DBL_MAX; 1196 model>nonLinearCost()>setOne(iSequence, solution[iSequence], 1197 trueLower, trueUpper, 0.0); 1198 } 1199 } 1200 } 1201 } 1202 break; 1203 // return 1 if there may be changing bounds on variable (column generation) 1204 case 10: { 1205 // return 1 as bounds on rhs will change 1206 returnCode = 1; 1207 } 1208 break; 1209 // make sure set is clean 1210 case 7: { 1211 // first flag 1212 if (number >= firstDynamic_ && number < lastDynamic_) { 1213 int sequence = id_[numberfirstDynamic_]; 1214 setFlagged(sequence); 1215 //model>clearFlagged(number); 1216 } else if (number>=model_>numberColumns()+numberStaticRows_) { 1217 // slack 1218 int iSet = fromIndex_[numbermodel_>numberColumns() 1219 numberStaticRows_]; 1220 setFlaggedSlack(iSet); 1221 //model>clearFlagged(number); 1222 } 1223 } 1224 case 11: { 1225 //int sequenceIn = model>sequenceIn(); 1226 if (number >= firstDynamic_ && number < lastDynamic_) { 1227 //assert (number == model>sequenceIn()); 1228 // take out variable (but leave key) 1229 double * cost = model>costRegion(); 1230 double * columnLower = model>lowerRegion(); 1231 double * columnUpper = model>upperRegion(); 1232 double * solution = model>solutionRegion(); 1233 int * length = matrix_>getMutableVectorLengths(); 1057 switch (mode) { 1058 // Fill in pivotVariable 1059 case 0: { 1060 // If no effective rhs  form it 1061 if (!rhsOffset_) { 1062 rhsOffset_ = new double[model>numberRows()]; 1063 rhsOffset(model, true); 1064 } 1065 int i; 1066 int numberBasic = number; 1067 int numberColumns = model>numberColumns(); 1068 // Use different array so can build from true pivotVariable_ 1069 //int * pivotVariable = model>pivotVariable(); 1070 int *pivotVariable = model>rowArray(0)>getIndices(); 1071 for (i = 0; i < numberColumns; i++) { 1072 if (model>getColumnStatus(i) == ClpSimplex::basic) 1073 pivotVariable[numberBasic++] = i; 1074 } 1075 number = numberBasic; 1076 } break; 1077 // Do initial extra rows + maximum basic 1078 case 2: { 1079 number = model>numberRows(); 1080 } break; 1081 // Before normal replaceColumn 1082 case 3: { 1083 if (numberActiveSets_ + numberStaticRows_ == model_>numberRows()) { 1084 // no space  refactorize 1085 returnCode = 4; 1086 number = 1; // say no need for normal replaceColumn 1087 } 1088 } break; 1089 // To see if can dual or primal 1090 case 4: { 1091 returnCode = 1; 1092 } break; 1093 // save status 1094 case 5: { 1095 memcpy(status_ + numberSets_, status_, numberSets_); 1096 memcpy(status_ + 2 * numberSets_, &numberActiveSets_, sizeof(int)); 1097 memcpy(dynamicStatus_ + maximumGubColumns_, 1098 dynamicStatus_, maximumGubColumns_); 1099 } break; 1100 // restore status 1101 case 6: { 1102 memcpy(status_, status_ + numberSets_, numberSets_); 1103 memcpy(&numberActiveSets_, status_ + 2 * numberSets_, sizeof(int)); 1104 memcpy(dynamicStatus_, dynamicStatus_ + maximumGubColumns_, 1105 maximumGubColumns_); 1106 initialProblem(); 1107 } break; 1108 // unflag all variables 1109 case 8: { 1110 for (int i = 0; i < numberGubColumns_; i++) { 1111 if (flagged(i)) { 1112 unsetFlagged(i); 1113 returnCode++; 1114 } 1115 } 1116 } break; 1117 // redo costs in primal 1118 case 9: { 1119 double *cost = model>costRegion(); 1120 double *solution = model>solutionRegion(); 1121 double *columnLower = model>lowerRegion(); 1122 double *columnUpper = model>upperRegion(); 1123 int i; 1124 bool doCosts = (number & 4) != 0; 1125 bool doBounds = (number & 1) != 0; 1126 for (i = firstDynamic_; i < firstAvailable_; i++) { 1127 int jColumn = id_[i  firstDynamic_]; 1128 if (doBounds) { 1129 if (!columnLower_ && !columnUpper_) { 1130 columnLower[i] = 0.0; 1131 columnUpper[i] = COIN_DBL_MAX; 1132 } else { 1133 if (columnLower_) 1134 columnLower[i] = columnLower_[jColumn]; 1135 else 1136 columnLower[i] = 0.0; 1137 if (columnUpper_) 1138 columnUpper[i] = columnUpper_[jColumn]; 1139 else 1140 columnUpper[i] = COIN_DBL_MAX; 1141 } 1142 } 1143 if (doCosts) { 1144 cost[i] = cost_[jColumn]; 1145 // Original bounds 1146 if (model>nonLinearCost()) 1147 model>nonLinearCost()>setOne(i, solution[i], 1148 this>columnLower(jColumn), 1149 this>columnUpper(jColumn), cost_[jColumn]); 1150 } 1151 } 1152 // and active sets 1153 for (i = 0; i < numberActiveSets_; i++) { 1154 int iSet = fromIndex_[i]; 1155 int iSequence = lastDynamic_ + numberStaticRows_ + i; 1156 if (doBounds) { 1157 if (lowerSet_[iSet] > 1.0e20) 1158 columnLower[iSequence] = lowerSet_[iSet]; 1159 else 1160 columnLower[iSequence] = COIN_DBL_MAX; 1161 if (upperSet_[iSet] < 1.0e20) 1162 columnUpper[iSequence] = upperSet_[iSet]; 1163 else 1164 columnUpper[iSequence] = COIN_DBL_MAX; 1165 } 1166 if (doCosts) { 1167 if (model>nonLinearCost()) { 1168 double trueLower; 1169 if (lowerSet_[iSet] > 1.0e20) 1170 trueLower = lowerSet_[iSet]; 1171 else 1172 trueLower = COIN_DBL_MAX; 1173 double trueUpper; 1174 if (upperSet_[iSet] < 1.0e20) 1175 trueUpper = upperSet_[iSet]; 1176 else 1177 trueUpper = COIN_DBL_MAX; 1178 model>nonLinearCost()>setOne(iSequence, solution[iSequence], 1179 trueLower, trueUpper, 0.0); 1180 } 1181 } 1182 } 1183 } break; 1184 // return 1 if there may be changing bounds on variable (column generation) 1185 case 10: { 1186 // return 1 as bounds on rhs will change 1187 returnCode = 1; 1188 } break; 1189 // make sure set is clean 1190 case 7: { 1191 // first flag 1192 if (number >= firstDynamic_ && number < lastDynamic_) { 1193 int sequence = id_[number  firstDynamic_]; 1194 setFlagged(sequence); 1195 //model>clearFlagged(number); 1196 } else if (number >= model_>numberColumns() + numberStaticRows_) { 1197 // slack 1198 int iSet = fromIndex_[number  model_>numberColumns()  numberStaticRows_]; 1199 setFlaggedSlack(iSet); 1200 //model>clearFlagged(number); 1201 } 1202 } 1203 case 11: { 1204 //int sequenceIn = model>sequenceIn(); 1205 if (number >= firstDynamic_ && number < lastDynamic_) { 1206 //assert (number == model>sequenceIn()); 1207 // take out variable (but leave key) 1208 double *cost = model>costRegion(); 1209 double *columnLower = model>lowerRegion(); 1210 double *columnUpper = model>upperRegion(); 1211 double *solution = model>solutionRegion(); 1212 int *length = matrix_>getMutableVectorLengths(); 1234 1213 #ifndef NDEBUG 1235 1236 int *row = matrix_>getMutableIndices();1237 CoinBigIndex *startColumn = matrix_>getMutableVectorStarts();1238 1239 1240 assert(which >= 0);1241 1242 assert(toIndex_[iSet] == which);1243 1214 if (length[number]) { 1215 int *row = matrix_>getMutableIndices(); 1216 CoinBigIndex *startColumn = matrix_>getMutableVectorStarts(); 1217 int iRow = row[startColumn[number] + length[number]  1]; 1218 int which = iRow  numberStaticRows_; 1219 assert(which >= 0); 1220 int iSet = fromIndex_[which]; 1221 assert(toIndex_[iSet] == which); 1222 } 1244 1223 #endif 1245 1246 1247 1248 1249 1250 1251 1252 1224 // no need firstAvailable_; 1225 solution[firstAvailable_] = 0.0; 1226 cost[firstAvailable_] = 0.0; 1227 length[firstAvailable_] = 0; 1228 model>nonLinearCost()>setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0); 1229 model>setStatus(firstAvailable_, ClpSimplex::atLowerBound); 1230 columnLower[firstAvailable_] = 0.0; 1231 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1253 1232 1254 // not really in small problem 1255 int iBig = id_[numberfirstDynamic_]; 1256 if (model>getStatus(number) == ClpSimplex::atLowerBound) { 1257 setDynamicStatus(iBig, atLowerBound); 1258 if (columnLower_) 1259 modifyOffset(number, columnLower_[iBig]); 1260 } else { 1261 setDynamicStatus(iBig, atUpperBound); 1262 modifyOffset(number, columnUpper_[iBig]); 1263 } 1264 } else if (number>=model_>numberColumns()+numberStaticRows_) { 1265 // slack 1266 int iSet = fromIndex_[numbermodel_>numberColumns() 1267 numberStaticRows_]; 1268 printf("what now  set %d\n",iSet); 1269 } 1270 } 1271 break; 1272 default: 1273 break; 1274 } 1275 return returnCode; 1233 // not really in small problem 1234 int iBig = id_[number  firstDynamic_]; 1235 if (model>getStatus(number) == ClpSimplex::atLowerBound) { 1236 setDynamicStatus(iBig, atLowerBound); 1237 if (columnLower_) 1238 modifyOffset(number, columnLower_[iBig]); 1239 } else { 1240 setDynamicStatus(iBig, atUpperBound); 1241 modifyOffset(number, columnUpper_[iBig]); 1242 } 1243 } else if (number >= model_>numberColumns() + numberStaticRows_) { 1244 // slack 1245 int iSet = fromIndex_[number  model_>numberColumns()  numberStaticRows_]; 1246 printf("what now  set %d\n", iSet); 1247 } 1248 } break; 1249 default: 1250 break; 1251 } 1252 return returnCode; 1276 1253 } 1277 1254 /* Purely for column generation and similar ideas. Allows … … 1279 1256 Returns nonzero if any changes. 1280 1257 */ 1281 int 1282 ClpDynamicMatrix::refresh(ClpSimplex * model) 1258 int ClpDynamicMatrix::refresh(ClpSimplex *model) 1283 1259 { 1284 1285 1286 1287 1288 1289 1290 1291 1260 // If at beginning then create initial problem 1261 if (firstDynamic_ == firstAvailable_) { 1262 initialProblem(); 1263 return 1; 1264 } else if (!model>nonLinearCost()) { 1265 // will be same as last time 1266 return 1; 1267 } 1292 1268 #ifndef NDEBUG 1293 1294 1295 1296 int *pivotVariable = model>pivotVariable();1297 for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {1298 assert (pivotVariable[i]==i+numberColumns);1299 1300 1269 { 1270 int numberColumns = model>numberColumns(); 1271 int numberRows = model>numberRows(); 1272 int *pivotVariable = model>pivotVariable(); 1273 for (int i = numberStaticRows_ + numberActiveSets_; i < numberRows; i++) { 1274 assert(pivotVariable[i] == i + numberColumns); 1275 } 1276 } 1301 1277 #endif 1302 1303 int * active = new int[numberActiveSets_];1304 1305 1306 1307 double * element =matrix_>getMutableElements();1308 int *row = matrix_>getMutableIndices();1309 CoinBigIndex *startColumn = matrix_>getMutableVectorStarts();1310 int *length = matrix_>getMutableVectorLengths();1311 double *cost = model>costRegion();1312 double *columnLower = model>lowerRegion();1313 double *columnUpper = model>upperRegion();1314 1315 1316 1317 1318 1319 double *solution = model>solutionRegion();1320 1321 1322 1323 1324 assert(which >= 0);1325 1326 assert(toIndex_[iSet] == which);1327 1328 1329 1330 keyVariable_[iSet] = id_[iColumnfirstDynamic_];1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1278 // lookup array 1279 int *active = new int[numberActiveSets_]; 1280 CoinZeroN(active, numberActiveSets_); 1281 int iColumn; 1282 int numberColumns = model>numberColumns(); 1283 double *element = matrix_>getMutableElements(); 1284 int *row = matrix_>getMutableIndices(); 1285 CoinBigIndex *startColumn = matrix_>getMutableVectorStarts(); 1286 int *length = matrix_>getMutableVectorLengths(); 1287 double *cost = model>costRegion(); 1288 double *columnLower = model>lowerRegion(); 1289 double *columnUpper = model>upperRegion(); 1290 CoinBigIndex numberElements = startColumn[firstDynamic_]; 1291 // first just do lookup and basic stuff 1292 int currentNumber = firstAvailable_; 1293 firstAvailable_ = firstDynamic_; 1294 double objectiveChange = 0.0; 1295 double *solution = model>solutionRegion(); 1296 int currentNumberActiveSets = numberActiveSets_; 1297 for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) { 1298 int iRow = row[startColumn[iColumn] + length[iColumn]  1]; 1299 int which = iRow  numberStaticRows_; 1300 assert(which >= 0); 1301 int iSet = fromIndex_[which]; 1302 assert(toIndex_[iSet] == which); 1303 if (model>getStatus(iColumn) == ClpSimplex::basic) { 1304 active[which]++; 1305 // may as well make key 1306 keyVariable_[iSet] = id_[iColumn  firstDynamic_]; 1307 } 1308 } 1309 int i; 1310 numberActiveSets_ = 0; 1311 int numberDeleted = 0; 1312 for (i = 0; i < currentNumberActiveSets; i++) { 1313 int iRow = i + numberStaticRows_; 1314 int numberActive = active[i]; 1315 int iSet = fromIndex_[i]; 1316 if (model>getRowStatus(iRow) == ClpSimplex::basic) { 1317 numberActive++; 1318 // may as well make key 1319 keyVariable_[iSet] = maximumGubColumns_ + iSet; 1320 } 1321 if (numberActive > 1) { 1322 // keep 1323 active[i] = numberActiveSets_; 1324 numberActiveSets_++; 1325 } else { 1326 active[i] = 1; 1327 } 1328 } 1353 1329 1354 1355 1356 1357 int jColumn = id_[iColumnfirstDynamic_];1358 if (model>getStatus(iColumn) == ClpSimplex::atLowerBound 1359 model>getStatus(iColumn) == ClpSimplex::atUpperBound) {1360 double value = solution[iColumn];1361 bool toLowerBound = true;1362 assert (jColumn>=0);assert (iColumn>=0);1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 id_[firstAvailable_firstDynamic_] = jColumn;1394 int numberThis = startColumn_[jColumn+1]  startColumn_[jColumn] + 1;1395 1396 1397 1398 1399 1400 1401 1402 1403 row[numberElements] = row_[base+j];1404 element[numberElements++] = element_[base+j];1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 for (iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 assert(toIndex_[iSet] == i);1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 rhsOffset_[numberActiveSets_+numberStaticRows_] = rhsOffset_[i+numberStaticRows_];1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 delete[] active;1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 rhsOffset_[i+numberStaticRows_] = 0.0;1483 1484 1485 1486 int *pivotVariable = model>pivotVariable();1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 if (iPut<numberStaticRows_+numberActiveSets_) {1516 1517 numberStaticRows_+numberActiveSets_iPut);1518 iPut = numberStaticRows_+numberActiveSets_;1519 1520 1521 1522 1523 1524 1330 for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) { 1331 int iRow = row[startColumn[iColumn] + length[iColumn]  1]; 1332 int which = iRow  numberStaticRows_; 1333 int jColumn = id_[iColumn  firstDynamic_]; 1334 if (model>getStatus(iColumn) == ClpSimplex::atLowerBound  model>getStatus(iColumn) == ClpSimplex::atUpperBound) { 1335 double value = solution[iColumn]; 1336 bool toLowerBound = true; 1337 assert(jColumn >= 0); 1338 assert(iColumn >= 0); 1339 if (columnUpper_) { 1340 if (!columnLower_) { 1341 if (fabs(value  columnUpper_[jColumn]) < fabs(value)) 1342 toLowerBound = false; 1343 } else if (fabs(value  columnUpper_[jColumn]) < fabs(value  columnLower_[jColumn])) { 1344 toLowerBound = false; 1345 } 1346 } 1347 if (toLowerBound) { 1348 // throw out to lower bound 1349 if (columnLower_) { 1350 setDynamicStatus(jColumn, atLowerBound); 1351 // treat solution as if exactly at a bound 1352 double value = columnLower[iColumn]; 1353 objectiveChange += cost[iColumn] * value; 1354 } else { 1355 setDynamicStatus(jColumn, atLowerBound); 1356 } 1357 } else { 1358 // throw out to upper bound 1359 setDynamicStatus(jColumn, atUpperBound); 1360 double value = columnUpper[iColumn]; 1361 objectiveChange += cost[iColumn] * value; 1362 } 1363 numberDeleted++; 1364 } else { 1365 assert(model>getStatus(iColumn) != ClpSimplex::superBasic); // deal with later 1366 int iPut = active[which]; 1367 if (iPut >= 0) { 1368 // move 1369 id_[firstAvailable_  firstDynamic_] = jColumn; 1370 int numberThis = startColumn_[jColumn + 1]  startColumn_[jColumn] + 1; 1371 length[firstAvailable_] = numberThis; 1372 cost[firstAvailable_] = cost[iColumn]; 1373 columnLower[firstAvailable_] = columnLower[iColumn]; 1374 columnUpper[firstAvailable_] = columnUpper[iColumn]; 1375 model>nonLinearCost()>setOne(firstAvailable_, solution[iColumn], 0.0, COIN_DBL_MAX, 1376 cost_[jColumn]); 1377 CoinBigIndex base = startColumn_[jColumn]; 1378 for (int j = 0; j < numberThis  1; j++) { 1379 row[numberElements] = row_[base + j]; 1380 element[numberElements++] = element_[base + j]; 1381 } 1382 row[numberElements] = iPut + numberStaticRows_; 1383 element[numberElements++] = 1.0; 1384 model>setStatus(firstAvailable_, model>getStatus(iColumn)); 1385 solution[firstAvailable_] = solution[iColumn]; 1386 firstAvailable_++; 1387 startColumn[firstAvailable_] = numberElements; 1388 } 1389 } 1390 } 1391 // zero solution 1392 CoinZeroN(solution + firstAvailable_, currentNumber  firstAvailable_); 1393 // zero cost 1394 CoinZeroN(cost + firstAvailable_, currentNumber  firstAvailable_); 1395 // zero lengths 1396 CoinZeroN(length + firstAvailable_, currentNumber  firstAvailable_); 1397 for (iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) { 1398 model>nonLinearCost()>setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0); 1399 model>setStatus(iColumn, ClpSimplex::atLowerBound); 1400 columnLower[iColumn] = 0.0; 1401 columnUpper[iColumn] = COIN_DBL_MAX; 1402 } 1403 // move rhs and set rest to infinite 1404 numberActiveSets_ = 0; 1405 for (i = 0; i < currentNumberActiveSets; i++) { 1406 int iSet = fromIndex_[i]; 1407 assert(toIndex_[iSet] == i); 1408 int iRow = i + numberStaticRows_; 1409 int iPut = active[i]; 1410 if (iPut >= 0) { 1411 int in = iRow + numberColumns; 1412 int out = iPut + numberColumns + numberStaticRows_; 1413 solution[out] = solution[in]; 1414 columnLower[out] = columnLower[in]; 1415 columnUpper[out] = columnUpper[in]; 1416 cost[out] = cost[in]; 1417 double trueLower; 1418 if (lowerSet_[iSet] > 1.0e20) 1419 trueLower = lowerSet_[iSet]; 1420 else 1421 trueLower = COIN_DBL_MAX; 1422 double trueUpper; 1423 if (upperSet_[iSet] < 1.0e20) 1424 trueUpper = upperSet_[iSet]; 1425 else 1426 trueUpper = COIN_DBL_MAX; 1427 model>nonLinearCost()>setOne(out, solution[out], 1428 trueLower, trueUpper, 0.0); 1429 model>setStatus(out, model>getStatus(in)); 1430 toIndex_[iSet] = numberActiveSets_; 1431 rhsOffset_[numberActiveSets_ + numberStaticRows_] = rhsOffset_[i + numberStaticRows_]; 1432 fromIndex_[numberActiveSets_++] = fromIndex_[i]; 1433 } else { 1434 // adjust offset 1435 // put out as key 1436 int jColumn = keyVariable_[iSet]; 1437 toIndex_[iSet] = 1; 1438 if (jColumn < maximumGubColumns_) { 1439 setDynamicStatus(jColumn, soloKey); 1440 double value = keyValue(iSet); 1441 objectiveChange += cost_[jColumn] * value; 1442 modifyOffset(jColumn, value); 1443 } 1444 } 1445 } 1446 model>setObjectiveOffset(model>objectiveOffset()  objectiveChange); 1447 delete[] active; 1448 for (i = numberActiveSets_; i < currentNumberActiveSets; i++) { 1449 int iSequence = i + numberStaticRows_ + numberColumns; 1450 solution[iSequence] = 0.0; 1451 columnLower[iSequence] = COIN_DBL_MAX; 1452 columnUpper[iSequence] = COIN_DBL_MAX; 1453 cost[iSequence] = 0.0; 1454 model>nonLinearCost()>setOne(iSequence, solution[iSequence], 1455 columnLower[iSequence], 1456 columnUpper[iSequence], 0.0); 1457 model>setStatus(iSequence, ClpSimplex::basic); 1458 rhsOffset_[i + numberStaticRows_] = 0.0; 1459 } 1460 if (currentNumberActiveSets != numberActiveSets_  numberDeleted) { 1461 // deal with pivotVariable 1462 int *pivotVariable = model>pivotVariable(); 1463 int sequence = firstDynamic_; 1464 int iRow = 0; 1465 int base1 = firstDynamic_; 1466 int base2 = lastDynamic_; 1467 int base3 = numberColumns + numberStaticRows_; 1468 int numberRows = numberStaticRows_ + currentNumberActiveSets; 1469 while (sequence < firstAvailable_) { 1470 int iPivot = pivotVariable[iRow]; 1471 while (iPivot < base1  (iPivot >= base2 && iPivot < base3)) { 1472 iPivot = pivotVariable[++iRow]; 1473 } 1474 pivotVariable[iRow++] = sequence; 1475 sequence++; 1476 } 1477 // move normal basic ones down 1478 int iPut = iRow; 1479 for (; iRow < numberRows; iRow++) { 1480 int iPivot = pivotVariable[iRow]; 1481 if (iPivot < base1  (iPivot >= base2 && iPivot < base3)) { 1482 pivotVariable[iPut++] = iPivot; 1483 } 1484 } 1485 // and basic others 1486 for (i = 0; i < numberActiveSets_; i++) { 1487 if (model>getRowStatus(i + numberStaticRows_) == ClpSimplex::basic) { 1488 pivotVariable[iPut++] = i + base3; 1489 } 1490 } 1491 if (iPut < numberStaticRows_ + numberActiveSets_) { 1492 printf("lost %d sets\n", 1493 numberStaticRows_ + numberActiveSets_  iPut); 1494 iPut = numberStaticRows_ + numberActiveSets_; 1495 } 1496 for (i = numberActiveSets_; i < currentNumberActiveSets; i++) { 1497 pivotVariable[iPut++] = i + base3; 1498 } 1499 //assert (iPut == numberRows); 1500 } 1525 1501 #ifdef CLP_DEBUG 1526 1502 #if 0 … … 1535 1511 #endif 1536 1512 #ifdef CLP_DEBUG 1537 1538 1539 1540 1541 1542 int * back = new int[numberGubColumns_];1543 1544 1545 1546 back[id_[ifirstDynamic_]] = i;1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 delete[] back;1566 1567 1513 { 1514 // debug coding to analyze set 1515 int i; 1516 int kSet = 6; 1517 if (kSet >= 0) { 1518 int *back = new int[numberGubColumns_]; 1519 for (i = 0; i < numberGubColumns_; i++) 1520 back[i] = 1; 1521 for (i = firstDynamic_; i < firstAvailable_; i++) 1522 back[id_[i  firstDynamic_]] = i; 1523 int sequence = startSet_[kSet]; 1524 if (toIndex_[kSet] < 0) { 1525 printf("Not in  Set %d  slack status %d, key %d\n", kSet, status_[kSet], keyVariable_[kSet]); 1526 while (sequence >= 0) { 1527 printf("( %d status %d ) ", sequence, getDynamicStatus(sequence)); 1528 sequence = next_[sequence]; 1529 } 1530 } else { 1531 int iRow = numberStaticRows_ + toIndex_[kSet]; 1532 printf("In  Set %d  slack status %d, key %d offset %g slack %g\n", kSet, status_[kSet], keyVariable_[kSet], 1533 rhsOffset_[iRow], model>solutionRegion(0)[iRow]); 1534 while (sequence >= 0) { 1535 int iBack = back[sequence]; 1536 printf("( %d status %d value %g) ", sequence, getDynamicStatus(sequence), model>solutionRegion()[iBack]); 1537 sequence = next_[sequence]; 1538 } 1539 } 1540 printf("\n"); 1541 delete[] back; 1542 } 1543 } 1568 1544 #endif 1569 1570 1571 1572 1573 1574 1575 int k=0;1576 for (int j=startSet_[i];j<startSet_[i+1];j++) {1577 if (getDynamicStatus(j)==soloKey)1578 1579 1580 assert (k<2);1581 1582 assert(n == numberSets_);1545 int n = numberActiveSets_; 1546 for (i = 0; i < numberSets_; i++) { 1547 if (toIndex_[i] < 0) { 1548 //assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]); 1549 n++; 1550 } 1551 int k = 0; 1552 for (int j = startSet_[i]; j < startSet_[i + 1]; j++) { 1553 if (getDynamicStatus(j) == soloKey) 1554 k++; 1555 } 1556 assert(k < 2); 1557 } 1558 assert(n == numberSets_); 1583 1559 #endif 1584 1560 return 1; 1585 1561 } 1586 void 1587 ClpDynamicMatrix::times(double scalar, 1588 const double * x, double * y) const 1562 void ClpDynamicMatrix::times(double scalar, 1563 const double *x, double *y) const 1589 1564 { 1590 1591 1592 1593 1594 const double * element =matrix_>getElements();1595 const int *row = matrix_>getIndices();1596 const CoinBigIndex *startColumn = matrix_>getVectorStarts();1597 const int *length = matrix_>getVectorLengths();1598 int *pivotVariable = model_>pivotVariable();1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1565 if (model_>specialOptions() != 16) { 1566 ClpPackedMatrix::times(scalar, x, y); 1567 } else { 1568 int iRow; 1569 const double *element = matrix_>getElements(); 1570 const int *row = matrix_>getIndices(); 1571 const CoinBigIndex *startColumn = matrix_>getVectorStarts(); 1572 const int *length = matrix_>getVectorLengths(); 1573 int *pivotVariable = model_>pivotVariable(); 1574 for (iRow = 0; iRow < numberStaticRows_ + numberActiveSets_; iRow++) { 1575 y[iRow] = scalar * rhsOffset_[iRow]; 1576 int iColumn = pivotVariable[iRow]; 1577 if (iColumn < lastDynamic_) { 1578 CoinBigIndex j; 1579 double value = scalar * x[iColumn]; 1580 if (value) { 1581 for (j = startColumn[iColumn]; 1582 j < startColumn[iColumn] + length[iColumn]; j++) { 1583 int jRow = row[j]; 1584 y[jRow] += value * element[j]; 1585 } 1586 } 1587 } 1588 } 1589 } 1615 1590 } 1616 1591 // Modifies rhs offset 1617 void 1618 ClpDynamicMatrix::modifyOffset(int sequence, double amount) 1592 void ClpDynamicMatrix::modifyOffset(int sequence, double amount) 1619 1593 { 1620 1621 assert(rhsOffset_);1622 1623 for (j = startColumn_[sequence]; j < startColumn_[sequence+1]; j++) {1624 1625 1626 1627 1594 if (amount) { 1595 assert(rhsOffset_); 1596 CoinBigIndex j; 1597 for (j = startColumn_[sequence]; j < startColumn_[sequence + 1]; j++) { 1598 int iRow = row_[j]; 1599 rhsOffset_[iRow] += amount * element_[j]; 1600 } 1601 } 1628 1602 } 1629 1603 // Gets key value when none in small … … 1631 1605 ClpDynamicMatrix::keyValue(int iSet) const 1632 1606 { 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 assert(status != inSmall);1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 assert(numberKey == 1);1656 1657 1658 1659 1660 assert(status != inSmall);1661 assert(status != soloKey);1662 1663 1664 1665 1666 1667 1668 1607 double value = 0.0; 1608 if (toIndex_[iSet] < 0) { 1609 int key = keyVariable_[iSet]; 1610 if (key < maximumGubColumns_) { 1611 if (getStatus(iSet) == ClpSimplex::atLowerBound) 1612 value = lowerSet_[iSet]; 1613 else 1614 value = upperSet_[iSet]; 1615 int numberKey = 0; 1616 int j = startSet_[iSet]; 1617 while (j >= 0) { 1618 DynamicStatus status = getDynamicStatus(j); 1619 assert(status != inSmall); 1620 if (status == soloKey) { 1621 numberKey++; 1622 } else if (status == atUpperBound) { 1623 value = columnUpper_[j]; 1624 } else if (columnLower_) { 1625 value = columnLower_[j]; 1626 } 1627 j = next_[j]; //onto next in set 1628 } 1629 assert(numberKey == 1); 1630 } else { 1631 int j = startSet_[iSet]; 1632 while (j >= 0) { 1633 DynamicStatus status = getDynamicStatus(j); 1634 assert(status != inSmall); 1635 assert(status != soloKey); 1636 if (status == atUpperBound) { 1637 value += columnUpper_[j]; 1638 } else if (columnLower_) { 1639 value += columnLower_[j]; 1640 } 1641 j = next_[j]; //onto next in set 1642 } 1669 1643 #if 0 1670 1644 // slack must be feasible … … 1676 1650 value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]); 1677 1651 #endif 1678 1679 1680 1652 } 1653 } 1654 return value; 1681 1655 } 1682 1656 // Switches off dj checking each factorization (for BIG models) 1683 void 1684 ClpDynamicMatrix::switchOffCheck() 1657 void ClpDynamicMatrix::switchOffCheck() 1685 1658 { 1686 1687 1659 noCheck_ = 0; 1660 infeasibilityWeight_ = 0.0; 1688 1661 } 1689 1662 /* Creates a variable. This is called after partial pricing and may modify matrix. 1690 1663 May update bestSequence. 1691 1664 */ 1692 void 1693 ClpDynamicMatrix::createVariable(ClpSimplex * model, int & bestSequence) 1665 void ClpDynamicMatrix::createVariable(ClpSimplex *model, int &bestSequence) 1694 1666 { 1695 1696 1697 1698 1699 1700 double *columnLower = model>lowerRegion();1701 double *columnUpper = model>upperRegion();1702 double *solution = model>solutionRegion();1703 double *reducedCost = model>djRegion();1704 const double *duals = model>dualRowSolution();1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1667 int numberRows = model>numberRows(); 1668 int slackOffset = lastDynamic_ + numberRows; 1669 int structuralOffset = slackOffset + numberSets_; 1670 int bestSequence2 = savedBestSequence_  structuralOffset; 1671 if (bestSequence >= slackOffset) { 1672 double *columnLower = model>lowerRegion(); 1673 double *columnUpper = model>upperRegion(); 1674 double *solution = model>solutionRegion(); 1675 double *reducedCost = model>djRegion(); 1676 const double *duals = model>dualRowSolution(); 1677 if (toIndex_[savedBestSet_] < 0) { 1678 // need to put key into basis 1679 int newRow = numberActiveSets_ + numberStaticRows_; 1680 model>dualRowSolution()[newRow] = savedBestGubDual_; 1681 double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set 1682 toIndex_[savedBestSet_] = numberActiveSets_; 1683 fromIndex_[numberActiveSets_++] = savedBestSet_; 1684 int iSequence = lastDynamic_ + newRow; 1685 // we need to get lower and upper correct 1686 double shift = 0.0; 1687 int j = startSet_[savedBestSet_]; 1688 while (j >= 0) { 1689 if (getDynamicStatus(j) == atUpperBound) 1690 shift += columnUpper_[j]; 1691 else if (getDynamicStatus(j) == atLowerBound && columnLower_) 1692 shift += columnLower_[j]; 1693 j = next_[j]; //onto next in set 1694 } 1695 if (lowerSet_[savedBestSet_] > 1.0e20) 1696 columnLower[iSequence] = lowerSet_[savedBestSet_]; 1697 else 1698 columnLower[iSequence] = COIN_DBL_MAX; 1699 if (upperSet_[savedBestSet_] < 1.0e20) 1700 columnUpper[iSequence] = upperSet_[savedBestSet_]; 1701 else 1702 columnUpper[iSequence] = COIN_DBL_MAX; 1731 1703 #ifdef CLP_DEBUG 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1704 if (model_>logLevel() == 63) { 1705 printf("%d in in set %d, key is %d rhs %g %g  keyvalue %g\n", 1706 bestSequence2, savedBestSet_, keyVariable_[savedBestSet_], 1707 columnLower[iSequence], columnUpper[iSequence], valueOfKey); 1708 int j = startSet_[savedBestSet_]; 1709 while (j >= 0) { 1710 if (getDynamicStatus(j) == atUpperBound) 1711 printf("%d atup ", j); 1712 else if (getDynamicStatus(j) == atLowerBound) 1713 printf("%d atlo ", j); 1714 else if (getDynamicStatus(j) == soloKey) 1715 printf("%d solo ", j); 1716 else 1717 abort(); 1718 j = next_[j]; //onto next in set 1719 } 1720 printf("\n"); 1721 } 1750 1722 #endif 1751 if (keyVariable_[savedBestSet_] < maximumGubColumns_) { 1752 // slack not key 1753 model_>pivotVariable()[newRow] = firstAvailable_; 1754 backToPivotRow_[firstAvailable_] = newRow; 1755 model>setStatus(iSequence, getStatus(savedBestSet_)); 1756 model>djRegion()[iSequence] = savedBestGubDual_; 1757 solution[iSequence] = valueOfKey; 1758 // create variable and pivot in 1759 int key = keyVariable_[savedBestSet_]; 1760 setDynamicStatus(key, inSmall); 1761 double * element = matrix_>getMutableElements(); 1762 int * row = matrix_>getMutableIndices(); 1763 CoinBigIndex * startColumn = matrix_>getMutableVectorStarts(); 1764 int * length = matrix_>getMutableVectorLengths(); 1765 CoinBigIndex numberElements = startColumn[firstAvailable_]; 1766 int numberThis = startColumn_[key+1]  startColumn_[key] + 1; 1767 if (numberElements + numberThis > numberElements_) { 1768 // need to redo 1769 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); 1770 matrix_>reserve(lastDynamic_, numberElements_); 1771 element = matrix_>getMutableElements(); 1772 row = matrix_>getMutableIndices(); 1773 // these probably okay but be safe 1774 startColumn = matrix_>getMutableVectorStarts(); 1775 length = matrix_>getMutableVectorLengths(); 1776 } 1777 // already set startColumn[firstAvailable_]=numberElements; 1778 length[firstAvailable_] = numberThis; 1779 model>costRegion()[firstAvailable_] = cost_[key]; 1780 CoinBigIndex base = startColumn_[key]; 1781 for (int j = 0; j < numberThis  1; j++) { 1782 row[numberElements] = row_[base+j]; 1783 element[numberElements++] = element_[base+j]; 1784 } 1785 row[numberElements] = newRow; 1786 element[numberElements++] = 1.0; 1787 id_[firstAvailable_firstDynamic_] = key; 1788 model>setObjectiveOffset(model>objectiveOffset() + cost_[key]* 1789 valueOfKey); 1790 model>solutionRegion()[firstAvailable_] = valueOfKey; 1791 model>setStatus(firstAvailable_, ClpSimplex::basic); 1792 // ***** need to adjust effective rhs 1793 if (!columnLower_ && !columnUpper_) { 1794 columnLower[firstAvailable_] = 0.0; 1795 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1796 } else { 1797 if (columnLower_) 1798 columnLower[firstAvailable_] = columnLower_[key]; 1799 else 1800 columnLower[firstAvailable_] = 0.0; 1801 if (columnUpper_) 1802 columnUpper[firstAvailable_] = columnUpper_[key]; 1803 else 1804 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1805 } 1806 model>nonLinearCost()>setOne(firstAvailable_, solution[firstAvailable_], 1807 columnLower[firstAvailable_], 1808 columnUpper[firstAvailable_], cost_[key]); 1809 startColumn[firstAvailable_+1] = numberElements; 1810 reducedCost[firstAvailable_] = 0.0; 1811 modifyOffset(key, valueOfKey); 1812 rhsOffset_[newRow] = shift; // sign? 1723 if (keyVariable_[savedBestSet_] < maximumGubColumns_) { 1724 // slack not key 1725 model_>pivotVariable()[newRow] = firstAvailable_; 1726 backToPivotRow_[firstAvailable_] = newRow; 1727 model>setStatus(iSequence, getStatus(savedBestSet_)); 1728 model>djRegion()[iSequence] = savedBestGubDual_; 1729 solution[iSequence] = valueOfKey; 1730 // create variable and pivot in 1731 int key = keyVariable_[savedBestSet_]; 1732 setDynamicStatus(key, inSmall); 1733 double *element = matrix_>getMutableElements(); 1734 int *row = matrix_>getMutableIndices(); 1735 CoinBigIndex *startColumn = matrix_>getMutableVectorStarts(); 1736 int *length = matrix_>getMutableVectorLengths(); 1737 CoinBigIndex numberElements = startColumn[firstAvailable_]; 1738 int numberThis = startColumn_[key + 1]  startColumn_[key] + 1; 1739 if (numberElements + numberThis > numberElements_) { 1740 // need to redo 1741 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); 1742 matrix_>reserve(lastDynamic_, numberElements_); 1743 element = matrix_>getMutableElements(); 1744 row = matrix_>getMutableIndices(); 1745 // these probably okay but be safe 1746 startColumn = matrix_>getMutableVectorStarts(); 1747 length = matrix_>getMutableVectorLengths(); 1748 } 1749 // already set startColumn[firstAvailable_]=numberElements; 1750 length[firstAvailable_] = numberThis; 1751 model>costRegion()[firstAvailable_] = cost_[key]; 1752 CoinBigIndex base = startColumn_[key]; 1753 for (int j = 0; j < numberThis  1; j++) { 1754 row[numberElements] = row_[base + j]; 1755 element[numberElements++] = element_[base + j]; 1756 } 1757 row[numberElements] = newRow; 1758 element[numberElements++] = 1.0; 1759 id_[firstAvailable_  firstDynamic_] = key; 1760 model>setObjectiveOffset(model>objectiveOffset() + cost_[key] * valueOfKey); 1761 model>solutionRegion()[firstAvailable_] = valueOfKey; 1762 model>setStatus(firstAvailable_, ClpSimplex::basic); 1763 // ***** need to adjust effective rhs 1764 if (!columnLower_ && !columnUpper_) { 1765 columnLower[firstAvailable_] = 0.0; 1766 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1767 } else { 1768 if (columnLower_) 1769 columnLower[firstAvailable_] = columnLower_[key]; 1770 else 1771 columnLower[firstAvailable_] = 0.0; 1772 if (columnUpper_) 1773 columnUpper[firstAvailable_] = columnUpper_[key]; 1774 else 1775 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1776 } 1777 model>nonLinearCost()>setOne(firstAvailable_, solution[firstAvailable_], 1778 columnLower[firstAvailable_], 1779 columnUpper[firstAvailable_], cost_[key]); 1780 startColumn[firstAvailable_ + 1] = numberElements; 1781 reducedCost[firstAvailable_] = 0.0; 1782 modifyOffset(key, valueOfKey); 1783 rhsOffset_[newRow] = shift; // sign? 1813 1784 #ifdef CLP_DEBUG 1814 1785 model>rowArray(1)>checkClear(); 1815 1786 #endif 1816 // now pivot in 1817 unpack(model, model>rowArray(1), firstAvailable_); 1818 model>factorization()>updateColumnFT(model>rowArray(2), model>rowArray(1)); 1819 double alpha = model>rowArray(1)>denseVector()[newRow]; 1820 int updateStatus = 1821 model>factorization()>replaceColumn(model, 1822 model>rowArray(2), 1823 model>rowArray(1), 1824 newRow, alpha); 1825 model>rowArray(1)>clear(); 1826 if (updateStatus) { 1827 if (updateStatus == 3) { 1828 // out of memory 1829 // increase space if not many iterations 1830 if (model>factorization()>pivots() < 1831 0.5 * model>factorization()>maximumPivots() && 1832 model>factorization()>pivots() < 400) 1833 model>factorization()>areaFactor( 1834 model>factorization()>areaFactor() * 1.1); 1835 } else { 1836 printf("Bad returncode %d from replaceColumn\n", updateStatus); 1837 } 1838 bestSequence = 1; 1839 return; 1840 } 1841 // firstAvailable_ only finally updated if good pivot (in updatePivot) 1842 // otherwise it reverts to firstAvailableBefore_ 1843 firstAvailable_++; 1844 } else { 1845 // slack key 1846 model>setStatus(iSequence, ClpSimplex::basic); 1847 model>djRegion()[iSequence] = 0.0; 1848 solution[iSequence] = valueOfKey+shift; 1849 rhsOffset_[newRow] = shift; // sign? 1850 } 1851 // correct slack 1852 model>costRegion()[iSequence] = 0.0; 1853 model>nonLinearCost()>setOne(iSequence, solution[iSequence], columnLower[iSequence], 1854 columnUpper[iSequence], 0.0); 1855 } 1856 if (savedBestSequence_ >= structuralOffset) { 1857 // recompute dj and create 1858 double value = cost_[bestSequence2]  savedBestGubDual_; 1859 for (CoinBigIndex jBigIndex = startColumn_[bestSequence2]; 1860 jBigIndex < startColumn_[bestSequence2+1]; jBigIndex++) { 1861 int jRow = row_[jBigIndex]; 1862 value = duals[jRow] * element_[jBigIndex]; 1863 } 1864 int gubRow = toIndex_[savedBestSet_] + numberStaticRows_; 1865 double * element = matrix_>getMutableElements(); 1866 int * row = matrix_>getMutableIndices(); 1867 CoinBigIndex * startColumn = matrix_>getMutableVectorStarts(); 1868 int * length = matrix_>getMutableVectorLengths(); 1869 CoinBigIndex numberElements = startColumn[firstAvailable_]; 1870 int numberThis = startColumn_[bestSequence2+1]  startColumn_[bestSequence2] + 1; 1871 if (numberElements + numberThis > numberElements_) { 1872 // need to redo 1873 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); 1874 matrix_>reserve(lastDynamic_, numberElements_); 1875 element = matrix_>getMutableElements(); 1876 row = matrix_>getMutableIndices(); 1877 // these probably okay but be safe 1878 startColumn = matrix_>getMutableVectorStarts(); 1879 length = matrix_>getMutableVectorLengths(); 1880 } 1881 // already set startColumn[firstAvailable_]=numberElements; 1882 length[firstAvailable_] = numberThis; 1883 model>costRegion()[firstAvailable_] = cost_[bestSequence2]; 1884 CoinBigIndex base = startColumn_[bestSequence2]; 1885 for (int j = 0; j < numberThis  1; j++) { 1886 row[numberElements] = row_[base+j]; 1887 element[numberElements++] = element_[base+j]; 1888 } 1889 row[numberElements] = gubRow; 1890 element[numberElements++] = 1.0; 1891 id_[firstAvailable_firstDynamic_] = bestSequence2; 1892 //printf("best %d\n",bestSequence2); 1893 model>solutionRegion()[firstAvailable_] = 0.0; 1894 model>clearFlagged(firstAvailable_); 1895 if (!columnLower_ && !columnUpper_) { 1896 model>setStatus(firstAvailable_, ClpSimplex::atLowerBound); 1897 columnLower[firstAvailable_] = 0.0; 1898 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1899 } else { 1900 DynamicStatus status = getDynamicStatus(bestSequence2); 1901 if (columnLower_) 1902 columnLower[firstAvailable_] = columnLower_[bestSequence2]; 1903 else 1904 columnLower[firstAvailable_] = 0.0; 1905 if (columnUpper_) 1906 columnUpper[firstAvailable_] = columnUpper_[bestSequence2]; 1907 else 1908 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1909 if (status == atLowerBound) { 1910 solution[firstAvailable_] = columnLower[firstAvailable_]; 1911 model>setStatus(firstAvailable_, ClpSimplex::atLowerBound); 1912 } else { 1913 solution[firstAvailable_] = columnUpper[firstAvailable_]; 1914 model>setStatus(firstAvailable_, ClpSimplex::atUpperBound); 1915 } 1916 } 1917 model>setObjectiveOffset(model>objectiveOffset() + cost_[bestSequence2]* 1918 solution[firstAvailable_]); 1919 model>nonLinearCost()>setOne(firstAvailable_, solution[firstAvailable_], 1920 columnLower[firstAvailable_], 1921 columnUpper[firstAvailable_], cost_[bestSequence2]); 1922 bestSequence = firstAvailable_; 1923 // firstAvailable_ only updated if good pivot (in updatePivot) 1924 startColumn[firstAvailable_+1] = numberElements; 1925 //printf("price struct %d  dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_); 1926 reducedCost[bestSequence] = value; 1787 // now pivot in 1788 unpack(model, model>rowArray(1), firstAvailable_); 1789 model>factorization()>updateColumnFT(model>rowArray(2), model>rowArray(1)); 1790 double alpha = model>rowArray(1)>denseVector()[newRow]; 1791 int updateStatus = model>factorization()>replaceColumn(model, 1792 model>rowArray(2), 1793 model>rowArray(1), 1794 newRow, alpha); 1795 model>rowArray(1)>clear(); 1796 if (updateStatus) { 1797 if (updateStatus == 3) { 1798 // out of memory 1799 // increase space if not many iterations 1800 if (model>factorization()>pivots() < 0.5 * model>factorization()>maximumPivots() && model>factorization()>pivots() < 400) 1801 model>factorization()>areaFactor( 1802 model>factorization()>areaFactor() * 1.1); 1927 1803 } else { 1928 // slack  row must just have been created 1929 assert (toIndex_[savedBestSet_] == numberActiveSets_  1); 1930 int newRow = numberStaticRows_ + numberActiveSets_  1; 1931 bestSequence = lastDynamic_ + newRow; 1932 reducedCost[bestSequence] = savedBestGubDual_; 1933 } 1934 } 1935 // clear for next iteration 1936 savedBestSequence_ = 1; 1804 printf("Bad returncode %d from replaceColumn\n", updateStatus); 1805 } 1806 bestSequence = 1; 1807 return; 1808 } 1809 // firstAvailable_ only finally updated if good pivot (in updatePivot) 1810 // otherwise it reverts to firstAvailableBefore_ 1811 firstAvailable_++; 1812 } else { 1813 // slack key 1814 model>setStatus(iSequence, ClpSimplex::basic); 1815 model>djRegion()[iSequence] = 0.0; 1816 solution[iSequence] = valueOfKey + shift; 1817 rhsOffset_[newRow] = shift; // sign? 1818 } 1819 // correct slack 1820 model>costRegion()[iSequence] = 0.0; 1821 model>nonLinearCost()>setOne(iSequence, solution[iSequence], columnLower[iSequence], 1822 columnUpper[iSequence], 0.0); 1823 } 1824 if (savedBestSequence_ >= structuralOffset) { 1825 // recompute dj and create 1826 double value = cost_[bestSequence2]  savedBestGubDual_; 1827 for (CoinBigIndex jBigIndex = startColumn_[bestSequence2]; 1828 jBigIndex < startColumn_[bestSequence2 + 1]; jBigIndex++) { 1829 int jRow = row_[jBigIndex]; 1830 value = duals[jRow] * element_[jBigIndex]; 1831 } 1832 int gubRow = toIndex_[savedBestSet_] + numberStaticRows_; 1833 double *element = matrix_>getMutableElements(); 1834 int *row = matrix_>getMutableIndices(); 1835 CoinBigIndex *startColumn = matrix_>getMutableVectorStarts(); 1836 int *length = matrix_>getMutableVectorLengths(); 1837 CoinBigIndex numberElements = startColumn[firstAvailable_]; 1838 int numberThis = startColumn_[bestSequence2 + 1]  startColumn_[bestSequence2] + 1; 1839 if (numberElements + numberThis > numberElements_) { 1840 // need to redo 1841 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); 1842 matrix_>reserve(lastDynamic_, numberElements_); 1843 element = matrix_>getMutableElements(); 1844 row = matrix_>getMutableIndices(); 1845 // these probably okay but be safe 1846 startColumn = matrix_>getMutableVectorStarts(); 1847 length = matrix_>getMutableVectorLengths(); 1848 } 1849 // already set startColumn[firstAvailable_]=numberElements; 1850 length[firstAvailable_] = numberThis; 1851 model>costRegion()[firstAvailable_] = cost_[bestSequence2]; 1852 CoinBigIndex base = startColumn_[bestSequence2]; 1853 for (int j = 0; j < numberThis  1; j++) { 1854 row[numberElements] = row_[base + j]; 1855 element[numberElements++] = element_[base + j]; 1856 } 1857 row[numberElements] = gubRow; 1858 element[numberElements++] = 1.0; 1859 id_[firstAvailable_  firstDynamic_] = bestSequence2; 1860 //printf("best %d\n",bestSequence2); 1861 model>solutionRegion()[firstAvailable_] = 0.0; 1862 model>clearFlagged(firstAvailable_); 1863 if (!columnLower_ && !columnUpper_) { 1864 model>setStatus(firstAvailable_, ClpSimplex::atLowerBound); 1865 columnLower[firstAvailable_] = 0.0; 1866 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1867 } else { 1868 DynamicStatus status = getDynamicStatus(bestSequence2); 1869 if (columnLower_) 1870 columnLower[firstAvailable_] = columnLower_[bestSequence2]; 1871 else 1872 columnLower[firstAvailable_] = 0.0; 1873 if (columnUpper_) 1874 columnUpper[firstAvailable_] = columnUpper_[bestSequence2]; 1875 else 1876 columnUpper[firstAvailable_] = COIN_DBL_MAX; 1877 if (status == atLowerBound) { 1878 solution[firstAvailable_] = columnLower[firstAvailable_]; 1879 model>setStatus(firstAvailable_, ClpSimplex::atLowerBound); 1880 } else { 1881 solution[firstAvailable_] = columnUpper[firstAvailable_]; 1882 model>setStatus(firstAvailable_, ClpSimplex::atUpperBound); 1883 } 1884 } 1885 model>setObjectiveOffset(model>objectiveOffset() + cost_[bestSequence2] * solution[firstAvailable_]); 1886 model>nonLinearCost()>setOne(firstAvailable_, solution[firstAvailable_], 1887 columnLower[firstAvailable_], 1888 columnUpper[firstAvailable_], cost_[bestSequence2]); 1889 bestSequence = firstAvailable_; 1890 // firstAvailable_ only updated if good pivot (in updatePivot) 1891 startColumn[firstAvailable_ + 1] = numberElements; 1892 //printf("price struct %d  dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_); 1893 reducedCost[bestSequence] = value; 1894 } else { 1895 // slack  row must just have been created 1896 assert(toIndex_[savedBestSet_] == numberActiveSets_  1); 1897 int newRow = numberStaticRows_ + numberActiveSets_  1; 1898 bestSequence = lastDynamic_ + newRow; 1899 reducedCost[bestSequence] = savedBestGubDual_; 1900 } 1901 } 1902 // clear for next iteration 1903 savedBestSequence_ = 1; 1937 1904 } 1938 1905 // Returns reduced cost of a variable 1939 1906 double 1940 ClpDynamicMatrix::reducedCost(ClpSimplex * 1907 ClpDynamicMatrix::reducedCost(ClpSimplex *model, int sequence) const 1941 1908 { 1942 1943 1944 1945 1946 1947 1909 int numberRows = model>numberRows(); 1910 int slackOffset = lastDynamic_ + numberRows; 1911 if (sequence < slackOffset) 1912 return model>djRegion()[sequence]; 1913 else 1914 return savedBestDj_; 1948 1915 } 1949 1916 // Does gub crash 1950 void 1951 ClpDynamicMatrix::gubCrash() 1917 void ClpDynamicMatrix::gubCrash() 1952 1918 { 1953 // Do basis  cheapest or slack if feasible 1954 int longestSet = 0; 1955 int iSet; 1956 for (iSet = 0; iSet < numberSets_; iSet++) { 1957 int n = 0; 1958 int j = startSet_[iSet]; 1959 while (j >= 0) { 1960 n++; 1961 j = next_[j]; 1962 } 1963 longestSet = CoinMax(longestSet, n); 1964 } 1965 double * upper = new double[longestSet+1]; 1966 double * cost = new double[longestSet+1]; 1967 double * lower = new double[longestSet+1]; 1968 double * solution = new double[longestSet+1]; 1969 int * back = new int[longestSet+1]; 1970 double tolerance = model_>primalTolerance(); 1971 double objectiveOffset = 0.0; 1972 for (iSet = 0; iSet < numberSets_; iSet++) { 1973 int iBasic = 1; 1974 double value = 0.0; 1975 // find cheapest 1976 int numberInSet = 0; 1977 int j = startSet_[iSet]; 1978 while (j >= 0) { 1979 if (!columnLower_) 1980 lower[numberInSet] = 0.0; 1981 else 1982 lower[numberInSet] = columnLower_[j]; 1983 if (!columnUpper_) 1984 upper[numberInSet] = COIN_DBL_MAX; 1985 else 1986 upper[numberInSet] = columnUpper_[j]; 1987 back[numberInSet++] = j; 1988 j = next_[j]; 1989 } 1990 CoinFillN(solution, numberInSet, 0.0); 1991 // and slack 1992 iBasic = numberInSet; 1993 solution[iBasic] = value; 1994 lower[iBasic] = upperSet_[iSet]; 1995 upper[iBasic] = lowerSet_[iSet]; 1996 int kphase; 1997 if (value < lowerSet_[iSet]  tolerance  value > upperSet_[iSet] + tolerance) { 1998 // infeasible 1999 kphase = 0; 2000 // remember bounds are flipped so opposite to natural 2001 if (value < lowerSet_[iSet]  tolerance) 2002 cost[iBasic] = 1.0; 2003 else 2004 cost[iBasic] = 1.0; 2005 CoinZeroN(cost, numberInSet); 2006 double dualTolerance = model_>dualTolerance(); 2007 for (int iphase = kphase; iphase < 2; iphase++) { 2008 if (iphase) { 2009 cost[numberInSet] = 0.0; 2010 for (int j = 0; j < numberInSet; j++) 2011 cost[j] = cost_[back[j]]; 2012 } 2013 // now do one row lp 2014 bool improve = true; 2015 while (improve) { 2016 improve = false; 2017 double dual = cost[iBasic]; 2018 int chosen = 1; 2019 double best = dualTolerance; 2020 int way = 0; 2021 for (int i = 0; i <= numberInSet; i++) { 2022 double dj = cost[i]  dual; 2023 double improvement = 0.0; 2024 if (iphase  i < numberInSet) 2025 assert (solution[i] >= lower[i] && solution[i] <= upper[i]); 2026 if (dj > dualTolerance) 2027 improvement = dj * (solution[i]  lower[i]); 2028 else if (dj < dualTolerance) 2029 improvement = dj * (solution[i]  upper[i]); 2030 if (improvement > best) { 2031 best = improvement; 2032 chosen = i; 2033 if (dj < 0.0) { 2034 way = 1; 2035 } else { 2036 way = 1; 2037 } 2038 } 2039 } 2040 if (chosen >= 0) { 2041 improve = true; 2042 // now see how far 2043 if (way > 0) { 2044 // incoming increasing so basic decreasing 2045 // if phase 0 then go to nearest bound 2046 double distance = upper[chosen]  solution[chosen]; 2047 double basicDistance; 2048 if (!iphase) { 2049 assert (iBasic == numberInSet); 2050 assert (solution[iBasic] > upper[iBasic]); 2051 basicDistance = solution[iBasic]  upper[iBasic]; 2052 } else { 2053 basicDistance = solution[iBasic]  lower[iBasic]; 2054 } 2055 // need extra coding for unbounded 2056 assert (CoinMin(distance, basicDistance) < 1.0e20); 2057 if (distance > basicDistance) { 2058 // incoming becomes basic 2059 solution[chosen] += basicDistance; 2060 if (!iphase) 2061 solution[iBasic] = upper[iBasic]; 2062 else 2063 solution[iBasic] = lower[iBasic]; 2064 iBasic = chosen; 2065 } else { 2066 // flip 2067 solution[chosen] = upper[chosen]; 2068 solution[iBasic] = distance; 2069 } 2070 } else { 2071 // incoming decreasing so basic increasing 2072 // if phase 0 then go to nearest bound 2073 double distance = solution[chosen]  lower[chosen]; 2074 double basicDistance; 2075 if (!iphase) { 2076 assert (iBasic == numberInSet); 2077 assert (solution[iBasic] < lower[iBasic]); 2078 basicDistance = lower[iBasic]  solution[iBasic]; 2079 } else { 2080 basicDistance = upper[iBasic]  solution[iBasic]; 2081 } 2082 // need extra coding for unbounded  for now just exit 2083 if (CoinMin(distance, basicDistance) > 1.0e20) { 2084 printf("unbounded on set %d\n", iSet); 2085 iphase = 1; 2086 iBasic = numberInSet; 2087 break; 2088 } 2089 if (distance > basicDistance) { 2090 // incoming becomes basic 2091 solution[chosen] = basicDistance; 2092 if (!iphase) 2093 solution[iBasic] = lower[iBasic]; 2094 else 2095 solution[iBasic] = upper[iBasic]; 2096 iBasic = chosen; 2097 } else { 2098 // flip 2099 solution[chosen] = lower[chosen]; 2100 solution[iBasic] += distance; 2101 } 2102 } 2103 if (!iphase) { 2104 if(iBasic < numberInSet) 2105 break; // feasible 2106 else if (solution[iBasic] >= lower[iBasic] && 2107 solution[iBasic] <= upper[iBasic]) 2108 break; // feasible (on flip) 2109 } 2110 } 2111 } 2112 } 2113 } 2114 // do solution i.e. bounds 2115 if (columnLower_  columnUpper_) { 2116 for (int j = 0; j < numberInSet; j++) { 2117 if (j != iBasic) { 2118 objectiveOffset += solution[j] * cost[j]; 2119 if (columnLower_ && columnUpper_) { 2120 if (fabs(solution[j]  columnLower_[back[j]]) > 2121 fabs(solution[j]  columnUpper_[back[j]])) 2122 setDynamicStatus(back[j], atUpperBound); 2123 } else if (columnUpper_ && solution[j] > 0.0) { 2124 setDynamicStatus(back[j], atUpperBound); 2125 } else { 2126 setDynamicStatus(back[j], atLowerBound); 2127 assert(!solution[j]); 2128 } 2129 } 2130 } 2131 } 2132 // convert iBasic back and do bounds 2133 if (iBasic == numberInSet) { 2134 // slack basic 2135 setStatus(iSet, ClpSimplex::basic); 2136 iBasic = iSet + maximumGubColumns_; 1919 // Do basis  cheapest or slack if feasible 1920 int longestSet = 0; 1921 int iSet; 1922 for (iSet = 0; iSet < numberSets_; iSet++) { 1923 int n = 0; 1924 int j = startSet_[iSet]; 1925 while (j >= 0) { 1926 n++; 1927 j = next_[j]; 1928 } 1929 longestSet = CoinMax(longestSet, n); 1930 } 1931 double *upper = new double[longestSet + 1]; 1932 double *cost = new double[longestSet + 1]; 1933 double *lower = new double[longestSet + 1]; 1934 double *solution = new double[longestSet + 1]; 1935 int *back = new int[longestSet + 1]; 1936 double tolerance = model_>primalTolerance(); 1937 double objectiveOffset = 0.0; 1938 for (iSet = 0; iSet < numberSets_; iSet++) { 1939 int iBasic = 1; 1940 double value = 0.0; 1941 // find cheapest 1942 int numberInSet = 0; 1943 int j = startSet_[iSet]; 1944 while (j >= 0) { 1945 if (!columnLower_) 1946 lower[numberInSet] = 0.0; 1947 else 1948 lower[numberInSet] = columnLower_[j]; 1949 if (!columnUpper_) 1950 upper[numberInSet] = COIN_DBL_MAX; 1951 else 1952 upper[numberInSet] = columnUpper_[j]; 1953 back[numberInSet++] = j; 1954 j = next_[j]; 1955 } 1956 CoinFillN(solution, numberInSet, 0.0); 1957 // and slack 1958 iBasic = numberInSet; 1959 solution[iBasic] = value; 1960 lower[iBasic] = upperSet_[iSet]; 1961 upper[iBasic] = lowerSet_[iSet]; 1962 int kphase; 1963 if (value < lowerSet_[iSet]  tolerance  value > upperSet_[iSet] + tolerance) { 1964 // infeasible 1965 kphase = 0; 1966 // remember bounds are flipped so opposite to natural 1967 if (value < lowerSet_[iSet]  tolerance) 1968 cost[iBasic] = 1.0; 1969 else 1970 cost[iBasic] = 1.0; 1971 CoinZeroN(cost, numberInSet); 1972 double dualTolerance = model_>dualTolerance(); 1973 for (int iphase = kphase; iphase < 2; iphase++) { 1974 if (iphase) { 1975 cost[numberInSet] = 0.0; 1976 for (int j = 0; j < numberInSet; j++) 1977 cost[j] = cost_[back[j]]; 1978 } 1979 // now do one row lp 1980 bool improve = true; 1981 while (improve) { 1982 improve = false; 1983 double dual = cost[iBasic]; 1984 int chosen = 1; 1985 double best = dualTolerance; 1986 int way = 0; 1987 for (int i = 0; i <= numberInSet; i++) { 1988 double dj = cost[i]  dual; 1989 double improvement = 0.0; 1990 if (iphase  i < numberInSet) 1991 assert(solution[i] >= lower[i] && solution[i] <= upper[i]); 1992 if (dj > dualTolerance) 1993 improvement = dj * (solution[i]  lower[i]); 1994 else if (dj < dualTolerance) 1995 improvement = dj * (solution[i]  upper[i]); 1996 if (improvement > best) { 1997 best = improvement; 1998 chosen = i; 1999 if (dj < 0.0) { 2000 way = 1; 2001 } else { 2002 way = 1; 2003 } 2004 } 2005 } 2006 if (chosen >= 0) { 2007 improve = true; 2008 // now see how far 2009 if (way > 0) { 2010 // incoming increasing so basic decreasing 2011 // if phase 0 then go to nearest bound 2012 double distance = upper[chosen]  solution[chosen]; 2013 double basicDistance; 2014 if (!iphase) { 2015 assert(iBasic == numberInSet); 2016 assert(solution[iBasic] > upper[iBasic]); 2017 basicDistance = solution[iBasic]  upper[iBasic]; 2018 } else { 2019 basicDistance = solution[iBasic]  lower[iBasic]; 2020 } 2021 // need extra coding for unbounded 2022 assert(CoinMin(distance, basicDistance) < 1.0e20); 2023 if (distance > basicDistance) { 2024 // incoming becomes basic 2025 solution[chosen] += basicDistance; 2026 if (!iphase) 2027 solution[iBasic] = upper[iBasic]; 2028 else 2029 solution[iBasic] = lower[iBasic]; 2030 iBasic = chosen; 2031 } else { 2032 // flip 2033 solution[chosen] = upper[chosen]; 2034 solution[iBasic] = distance; 2035 } 2036 } else { 2037 // incoming decreasing so basic increasing 2038 // if phase 0 then go to nearest bound 2039 double distance = solution[chosen]  lower[chosen]; 2040 double basicDistance; 2041 if (!iphase) { 2042 assert(iBasic == numberInSet); 2043 assert(solution[iBasic] < lower[iBasic]); 2044 basicDistance = lower[iBasic]  solution[iBasic]; 2045 } else { 2046 basicDistance = upper[iBasic]  solution[iBasic]; 2047 } 2048 // need extra coding for unbounded  for now just exit 2049 if (CoinMin(distance, basicDistance) > 1.0e20) { 2050 printf("unbounded on set %d\n", iSet); 2051 iphase = 1; 2052 iBasic = numberInSet; 2053 break; 2054 } 2055 if (distance > basicDistance) { 2056 // incoming becomes basic 2057 solution[chosen] = basicDistance; 2058 if (!iphase) 2059 solution[iBasic] = lower[iBasic]; 2060 else 2061 solution[iBasic] = upper[iBasic]; 2062 iBasic = chosen; 2063 } else { 2064 // flip 2065 solution[chosen] = lower[chosen]; 2066 solution[iBasic] += distance; 2067 } 2068 } 2069 if (!iphase) { 2070 if (iBasic < numberInSet) 2071 break; // feasible 2072 else if (solution[iBasic] >= lower[iBasic] && solution[iBasic] <= upper[iBasic]) 2073 break; // feasible (on flip) 2074 } 2075 } 2076 } 2077 } 2078 } 2079 // do solution i.e. bounds 2080 if (columnLower_  columnUpper_) { 2081 for (int j = 0; j < numberInSet; j++) { 2082 if (j != iBasic) { 2083 objectiveOffset += solution[j] * cost[j]; 2084 if (columnLower_ && columnUpper_) { 2085 if (fabs(solution[j]  columnLower_[back[j]]) > fabs(solution[j]  columnUpper_[back[j]])) 2086 setDynamicStatus(back[j], atUpperBound); 2087 } else if (columnUpper_ && solution[j] > 0.0) { 2088 setDynamicStatus(back[j], atUpperBound); 2137 2089 } else { 2138 iBasic = back[iBasic]; 2139 setDynamicStatus(iBasic, soloKey); 2140 // remember bounds flipped 2141 if (upper[numberInSet] == lower[numberInSet]) 2142 setStatus(iSet, ClpSimplex::isFixed); 2143 else if (solution[numberInSet] == upper[numberInSet]) 2144 setStatus(iSet, ClpSimplex::atLowerBound); 2145 else if (solution[numberInSet] == lower[numberInSet]) 2146 setStatus(iSet, ClpSimplex::atUpperBound); 2147 else 2148 abort(); 2149 } 2150 keyVariable_[iSet] = iBasic; 2151 } 2152 model_>setObjectiveOffset(objectiveOffset_  objectiveOffset); 2153 delete [] lower; 2154 delete [] solution; 2155 delete [] upper; 2156 delete [] cost; 2157 delete [] back; 2158 // make sure matrix is in good shape 2159 matrix_>orderMatrix(); 2090 setDynamicStatus(back[j], atLowerBound); 2091 assert(!solution[j]); 2092 } 2093 } 2094 } 2095 } 2096 // convert iBasic back and do bounds 2097 if (iBasic == numberInSet) { 2098 // slack basic 2099 setStatus(iSet, ClpSimplex::basic); 2100 iBasic = iSet + maximumGubColumns_; 2101 } else { 2102 iBasic = back[iBasic]; 2103 setDynamicStatus(iBasic, soloKey); 2104 // remember bounds flipped 2105 if (upper[numberInSet] == lower[numberInSet]) 2106 setStatus(iSet, ClpSimplex::isFixed); 2107 else if (solution[numberInSet] == upper[numberInSet]) 2108 setStatus(iSet, ClpSimplex::atLowerBound); 2109 else if (solution[numberInSet] == lower[numberInSet]) 2110 setStatus(iSet, ClpSimplex::atUpperBound); 2111 else 2112 abort(); 2113 } 2114 keyVariable_[iSet] = iBasic; 2115 } 2116 model_>setObjectiveOffset(objectiveOffset_  objectiveOffset); 2117 delete[] lower; 2118 delete[] solution; 2119 delete[] upper; 2120 delete[] cost; 2121 delete[] back; 2122 // make sure matrix is in good shape 2123 matrix_>orderMatrix(); 2160 2124 } 2161 2125 // Populates initial matrix from dynamic status 2162 void 2163 ClpDynamicMatrix::initialProblem() 2126 void ClpDynamicMatrix::initialProblem() 2164 2127 { 2165 2166 double * element =matrix_>getMutableElements();2167 int *row = matrix_>getMutableIndices();2168 CoinBigIndex *startColumn = matrix_>getMutableVectorStarts();2169 int *length = matrix_>getMutableVectorLengths();2170 double *cost = model_>objective();2171 double *solution = model_>primalColumnSolution();2172 double *columnLower = model_>columnLower();2173 double *columnUpper = model_>columnUpper();2174 double *rowSolution = model_>primalRowSolution();2175 double *rowLower = model_>rowLower();2176 double *rowUpper = model_>rowUpper();2177 2128 int iSet; 2129 double *element = matrix_>getMutableElements(); 2130 int *row = matrix_>getMutableIndices(); 2131 CoinBigIndex *startColumn = matrix_>getMutableVectorStarts(); 2132 int *length = matrix_>getMutableVectorLengths(); 2133 double *cost = model_>objective(); 2134 double *solution = model_>primalColumnSolution(); 2135 double *columnLower = model_>columnLower(); 2136 double *columnUpper = model_>columnUpper(); 2137 double *rowSolution = model_>primalRowSolution(); 2138 double *rowLower = model_>rowLower(); 2139 double *rowUpper = model_>rowUpper(); 2140 CoinBigIndex numberElements = startColumn[firstDynamic_]; 2178 2141 2179 firstAvailable_ = firstDynamic_; 2180 numberActiveSets_ = 0; 2181 for (iSet = 0; iSet < numberSets_; iSet++) { 2182 toIndex_[iSet] = 1; 2183 int numberActive = 0; 2184 int whichKey = 1; 2185 if (getStatus(iSet) == ClpSimplex::basic) { 2186 whichKey = maximumGubColumns_ + iSet; 2187 numberActive = 1; 2142 firstAvailable_ = firstDynamic_; 2143 numberActiveSets_ = 0; 2144 for (iSet = 0; iSet < numberSets_; iSet++) { 2145 toIndex_[iSet] = 1; 2146 int numberActive = 0; 2147 int whichKey = 1; 2148 if (getStatus(iSet) == ClpSimplex::basic) { 2149 whichKey = maximumGubColumns_ + iSet; 2150 numberActive = 1; 2151 } else { 2152 whichKey = 1; 2153 } 2154 int j = startSet_[iSet]; 2155 while (j >= 0) { 2156 assert(getDynamicStatus(j) != soloKey  whichKey == 1); 2157 if (getDynamicStatus(j) == inSmall) { 2158 numberActive++; 2159 } else if (getDynamicStatus(j) == soloKey) { 2160 whichKey = j; 2161 numberActive++; 2162 } 2163 j = next_[j]; //onto next in set 2164 } 2165 if (numberActive > 1) { 2166 int iRow = numberActiveSets_ + numberStaticRows_; 2167 rowSolution[iRow] = 0.0; 2168 double lowerValue; 2169 if (lowerSet_[iSet] > 1.0e20) 2170 lowerValue = lowerSet_[iSet]; 2171 else 2172 lowerValue = COIN_DBL_MAX; 2173 double upperValue; 2174 if (upperSet_[iSet] < 1.0e20) 2175 upperValue = upperSet_[iSet]; 2176 else 2177 upperValue = COIN_DBL_MAX; 2178 rowLower[iRow] = lowerValue; 2179 rowUpper[iRow] = upperValue; 2180 if (getStatus(iSet) == ClpSimplex::basic) { 2181 model_>setRowStatus(iRow, ClpSimplex::basic); 2182 rowSolution[iRow] = 0.0; 2183 } else if (getStatus(iSet) == ClpSimplex::atLowerBound) { 2184 model_>setRowStatus(iRow, ClpSimplex::atLowerBound); 2185 rowSolution[iRow] = lowerValue; 2186 } else { 2187 model_>setRowStatus(iRow, ClpSimplex::atUpperBound); 2188 rowSolution[iRow] = upperValue; 2189 } 2190 j = startSet_[iSet]; 2191 while (j >= 0) { 2192 DynamicStatus status = getDynamicStatus(j); 2193 if (status == inSmall) { 2194 int numberThis = startColumn_[j + 1]  startColumn_[j] + 1; 2195 if (numberElements + numberThis > numberElements_) { 2196 // need to redo 2197 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); 2198 matrix_>reserve(lastDynamic_, numberElements_); 2199 element = matrix_>getMutableElements(); 2200 row = matrix_>getMutableIndices(); 2201 // these probably okay but be safe 2202 startColumn = matrix_>getMutableVectorStarts(); 2203 length = matrix_>getMutableVectorLengths(); 2204 } 2205 length[firstAvailable_] = numberThis; 2206 cost[firstAvailable_] = cost_[j]; 2207 CoinBigIndex base = startColumn_[j]; 2208 for (int k = 0; k < numberThis  1; k++) { 2209 row[numberElements] = row_[base + k]; 2210 element[numberElements++] = element_[base + k]; 2211 } 2212 row[numberElements] = iRow; 2213 element[numberElements++] = 1.0; 2214 id_[firstAvailable_  firstDynamic_] = j; 2215 solution[firstAvailable_] = 0.0; 2216 model_>setStatus(firstAvailable_, ClpSimplex::basic); 2217 if (!columnLower_ && !columnUpper_) { 2218 columnLower[firstAvailable_] = 0.0; 2219 columnUpper[firstAvailable_] = COIN_DBL_MAX; 2188 2220 } else { 2189 whichKey = 1; 2190 } 2191 int j = startSet_[iSet]; 2192 while (j >= 0) { 2193 assert (getDynamicStatus(j) != soloKey  whichKey == 1); 2194 if (getDynamicStatus(j) == inSmall) { 2195 numberActive++; 2196 } else if (getDynamicStatus(j) == soloKey) { 2197 whichKey = j; 2198 numberActive++; 2199 } 2200 j = next_[j]; //onto next in set 2201 } 2202 if (numberActive > 1) { 2203 int iRow = numberActiveSets_ + numberStaticRows_; 2204 rowSolution[iRow] = 0.0; 2205 double lowerValue; 2206 if (lowerSet_[iSet] > 1.0e20) 2207 lowerValue = lowerSet_[iSet]; 2208 else 2209 lowerValue = COIN_DBL_MAX; 2210 double upperValue; 2211 if (upperSet_[iSet] < 1.0e20) 2212 upperValue = upperSet_[iSet]; 2213 else 2214 upperValue = COIN_DBL_MAX; 2215 rowLower[iRow] = lowerValue; 2216 rowUpper[iRow] = upperValue; 2217 if (getStatus(iSet) == ClpSimplex::basic) { 2218 model_>setRowStatus(iRow, ClpSimplex::basic); 2219 rowSolution[iRow] = 0.0; 2220 } else if (getStatus(iSet) == ClpSimplex::atLowerBound) { 2221 model_>setRowStatus(iRow, ClpSimplex::atLowerBound); 2222 rowSolution[iRow] = lowerValue; 2223 } else { 2224 model_>setRowStatus(iRow, ClpSimplex::atUpperBound); 2225 rowSolution[iRow] = upperValue; 2226 } 2227 j = startSet_[iSet]; 2228 while (j >= 0) { 2229 DynamicStatus status = getDynamicStatus(j); 2230 if (status == inSmall) { 2231 int numberThis = startColumn_[j+1]  startColumn_[j] + 1; 2232 if (numberElements + numberThis > numberElements_) { 2233 // need to redo 2234 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); 2235 matrix_>reserve(lastDynamic_, numberElements_); 2236 element = matrix_>getMutableElements(); 2237 row = matrix_>getMutableIndices(); 2238 // these probably okay but be safe 2239 startColumn = matrix_>getMutableVectorStarts(); 2240 length = matrix_>getMutableVectorLengths(); 2241 } 2242 length[firstAvailable_] = numberThis; 2243 cost[firstAvailable_] = cost_[j]; 2244 CoinBigIndex base = startColumn_[j]; 2245 for (int k = 0; k < numberThis  1; k++) { 2246 row[numberElements] = row_[base+k]; 2247 element[numberElements++] = element_[base+k]; 2248 } 2249 row[numberElements] = iRow; 2250 element[numberElements++] = 1.0; 2251 id_[firstAvailable_firstDynamic_] = j; 2252 solution[firstAvailable_] = 0.0; 2253 model_>setStatus(firstAvailable_, ClpSimplex::basic); 2254 if (!columnLower_ && !columnUpper_) { 2255 columnLower[firstAvailable_] = 0.0; 2256 columnUpper[firstAvailable_] = COIN_DBL_MAX; 2257 } else { 2258 if (columnLower_) 2259 columnLower[firstAvailable_] = columnLower_[j]; 2260 else 2261 columnLower[firstAvailable_] = 0.0; 2262 if (columnUpper_) 2263 columnUpper[firstAvailable_] = columnUpper_[j]; 2264 else 2265 columnUpper[firstAvailable_] = COIN_DBL_MAX; 2266 if (status != atUpperBound) { 2267 solution[firstAvailable_] = columnLower[firstAvailable_]; 2268 } else { 2269 solution[firstAvailable_] = columnUpper[firstAvailable_]; 2270 } 2271 } 2272 firstAvailable_++; 2273 startColumn[firstAvailable_] = numberElements; 2274 } 2275 j = next_[j]; //onto next in set 2276 } 2277 model_>setRowStatus(numberActiveSets_ + numberStaticRows_, getStatus(iSet)); 2278 toIndex_[iSet] = numberActiveSets_; 2279 fromIndex_[numberActiveSets_++] = iSet; 2280 } else { 2281 // solo key 2282 bool needKey=false; 2283 if (numberActive) { 2284 if (whichKey<maximumGubColumns_) { 2285 // structural  assume ok 2286 needKey = false; 2287 } else { 2288 // slack 2289 keyVariable_[iSet] = maximumGubColumns_ + iSet; 2290 double value = keyValue(iSet); 2291 if (value<lowerSet_[iSet]1.0e8 2292 value>upperSet_[iSet]+1.0e8) 2293 needKey=true; 2294 } 2295 } else { 2296 needKey = true; 2297 } 2298 if (needKey) { 2299 // all to lb then up some (slack/null if possible) 2300 int length=99999999; 2301 int which=1; 2302 double sum=0.0; 2303 for (int iColumn=startSet_[iSet];iColumn<startSet_[iSet+1];iColumn++) { 2304 setDynamicStatus(iColumn,atLowerBound); 2305 sum += columnLower_[iColumn]; 2306 if (length>startColumn_[iColumn+1]startColumn_[iColumn]) { 2307 which=iColumn; 2308 length=startColumn_[iColumn+1]startColumn_[iColumn]; 2309 } 2310 } 2311 if (sum>lowerSet_[iSet]1.0e8) { 2312 // slack can be basic 2313 setStatus(iSet,ClpSimplex::basic); 2314 keyVariable_[iSet] = maximumGubColumns_ + iSet; 2315 } else { 2316 // use shortest 2317 setDynamicStatus(which,soloKey); 2318 keyVariable_[iSet] = which; 2319 setStatus(iSet,ClpSimplex::atLowerBound); 2320 } 2321 } 2322 } 2323 assert (toIndex_[iSet] >= 0  whichKey >= 0); 2324 keyVariable_[iSet] = whichKey; 2325 } 2326 // clean up pivotVariable 2327 int numberColumns = model_>numberColumns(); 2328 int numberRows = model_>numberRows(); 2329 int * pivotVariable = model_>pivotVariable(); 2330 if (pivotVariable) { 2331 for (int i=0; i<numberStaticRows_+numberActiveSets_;i++) { 2332 if (model_>getRowStatus(i)!=ClpSimplex::basic) 2333 pivotVariable[i]=1; 2334 else 2335 pivotVariable[i]=numberColumns+i; 2336 } 2337 for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) { 2338 pivotVariable[i]=i+numberColumns; 2339 } 2340 int put=1; 2341 for (int i=0;i<numberColumns;i++) { 2342 if (model_>getColumnStatus(i)==ClpSimplex::basic) { 2343 while(put<numberRows) { 2344 put++; 2345 if (pivotVariable[put]==1) { 2346 pivotVariable[put]=i; 2347 break; 2348 } 2349 } 2350 } 2351 } 2352 for (int i=CoinMax(put,0);i<numberRows;i++) { 2353 if (pivotVariable[i]==1) 2354 pivotVariable[i]=i+numberColumns; 2355 } 2356 } 2357 if (rhsOffset_ && model_>costRegion()) { 2358 double * cost = model_>costRegion(); 2359 double * columnLower = model_>lowerRegion(); 2360 double * columnUpper = model_>upperRegion(); 2361 double * solution = model_>solutionRegion(); 2362 int numberRows = model_>numberRows(); 2363 for (int i = numberActiveSets_; i < numberRowsnumberStaticRows_; i++) { 2364 int iSequence = i + numberStaticRows_ + numberColumns; 2365 solution[iSequence] = 0.0; 2366 columnLower[iSequence] = COIN_DBL_MAX; 2367 columnUpper[iSequence] = COIN_DBL_MAX; 2368 cost[iSequence] = 0.0; 2369 model_>nonLinearCost()>setOne(iSequence, solution[iSequence], 2370 columnLower[iSequence], 2371 columnUpper[iSequence], 0.0); 2372 model_>setStatus(iSequence, ClpSimplex::basic); 2373 rhsOffset_[i+numberStaticRows_] = 0.0; 2374 } 2221 if (columnLower_) 2222 columnLower[firstAvailable_] = columnLower_[j]; 2223 else 2224 columnLower[firstAvailable_] = 0.0; 2225 if (columnUpper_) 2226 columnUpper[firstAvailable_] = columnUpper_[j]; 2227 else 2228 columnUpper[firstAvailable_] = COIN_DBL_MAX; 2229 if (status != atUpperBound) { 2230 solution[firstAvailable_] = columnLower[firstAvailable_]; 2231 } else { 2232 solution[firstAvailable_] = columnUpper[firstAvailable_]; 2233 } 2234 } 2235 firstAvailable_++; 2236 startColumn[firstAvailable_] = numberElements; 2237 } 2238 j = next_[j]; //onto next in set 2239 } 2240 model_>setRowStatus(numberActiveSets_ + numberStaticRows_, getStatus(iSet)); 2241 toIndex_[iSet] = numberActiveSets_; 2242 fromIndex_[numberActiveSets_++] = iSet; 2243 } else { 2244 // solo key 2245 bool needKey = false; 2246 if (numberActive) { 2247 if (whichKey < maximumGubColumns_) { 2248 // structural  assume ok 2249 needKey = false; 2250 } else { 2251 // slack 2252 keyVariable_[iSet] = maximumGubColumns_ + iSet; 2253 double value = keyValue(iSet); 2254 if (value < lowerSet_[iSet]  1.0e8  value > upperSet_[iSet] + 1.0e8) 2255 needKey = true; 2256 } 2257 } else { 2258 needKey = true; 2259 } 2260 if (needKey) { 2261 // all to lb then up some (slack/null if possible) 2262 int length = 99999999; 2263 int which = 1; 2264 double sum = 0.0; 2265 for (int iColumn = startSet_[iSet]; iColumn < startSet_[iSet + 1]; iColumn++) { 2266 setDynamicStatus(iColumn, atLowerBound); 2267 sum += columnLower_[iColumn]; 2268 if (length > startColumn_[iColumn + 1]  startColumn_[iColumn]) { 2269 which = iColumn; 2270 length = startColumn_[iColumn + 1]  startColumn_[iColumn]; 2271 } 2272 } 2273 if (sum > lowerSet_[iSet]  1.0e8) { 2274 // slack can be basic 2275 setStatus(iSet, ClpSimplex::basic); 2276 keyVariable_[iSet] = maximumGubColumns_ + iSet; 2277 } else { 2278 // use shortest 2279 setDynamicStatus(which, soloKey); 2280 keyVariable_[iSet] = which; 2281 setStatus(iSet, ClpSimplex::atLowerBound); 2282 } 2283 } 2284 } 2285 assert(toIndex_[iSet] >= 0  whichKey >= 0); 2286 keyVariable_[iSet] = whichKey; 2287 } 2288 // clean up pivotVariable 2289 int numberColumns = model_>numberColumns(); 2290 int numberRows = model_>numberRows(); 2291 int *pivotVariable = model_>pivotVariable(); 2292 if (pivotVariable) { 2293 for (int i = 0; i < numberStaticRows_ + numberActiveSets_; i++) { 2294 if (model_>getRowStatus(i) != ClpSimplex::basic) 2295 pivotVariable[i] = 1; 2296 else 2297 pivotVariable[i] = numberColumns + i; 2298 } 2299 for (int i = numberStaticRows_ + numberActiveSets_; i < numberRows; i++) { 2300 pivotVariable[i] = i + numberColumns; 2301 } 2302 int put = 1; 2303 for (int i = 0; i < numberColumns; i++) { 2304 if (model_>getColumnStatus(i) == ClpSimplex::basic) { 2305 while (put < numberRows) { 2306 put++; 2307 if (pivotVariable[put] == 1) { 2308 pivotVariable[put] = i; 2309 break; 2310 } 2311 } 2312 } 2313 } 2314 for (int i = CoinMax(put, 0); i < numberRows; i++) { 2315 if (pivotVariable[i] == 1) 2316 pivotVariable[i] = i + numberColumns; 2317 } 2318 } 2319 if (rhsOffset_ && model_>costRegion()) { 2320 double *cost = model_>costRegion(); 2321 double *columnLower = model_>lowerRegion(); 2322 double *columnUpper = model_>upperRegion(); 2323 double *solution = model_>solutionRegion(); 2324 int numberRows = model_>numberRows(); 2325 for (int i = numberActiveSets_; i < numberRows  numberStaticRows_; i++) { 2326 int iSequence = i + numberStaticRows_ + numberColumns; 2327 solution[iSequence] = 0.0; 2328 columnLower[iSequence] = COIN_DBL_MAX; 2329 columnUpper[iSequence] = COIN_DBL_MAX; 2330 cost[iSequence] = 0.0; 2331 model_>nonLinearCost()>setOne(iSequence, solution[iSequence], 2332 columnLower[iSequence], 2333 columnUpper[iSequence], 0.0); 2334 model_>setStatus(iSequence, ClpSimplex::basic); 2335 rhsOffset_[i + numberStaticRows_] = 0.0; 2336 } 2375 2337 #if 0 2376 2338 for (int i=0;i<numberStaticRows_;i++) … … 2378 2340 i,rhsOffset_[i]); 2379 2341 #endif 2380 2381 2342 } 2343 numberActiveColumns_ = firstAvailable_; 2382 2344 #if 0 2383 2345 for (iSet = 0; iSet < numberSets_; iSet++) { … … 2390 2352 } 2391 2353 #endif 2392 2354 return; 2393 2355 } 2394 2356 // Writes out model (without names) 2395 void 2396 ClpDynamicMatrix::writeMps(const char * name) 2357 void ClpDynamicMatrix::writeMps(const char *name) 2397 2358 { 2398 int numberTotalRows = numberStaticRows_ +numberSets_;2399 int numberTotalColumns = firstDynamic_ +numberGubColumns_;2359 int numberTotalRows = numberStaticRows_ + numberSets_; 2360 int numberTotalColumns = firstDynamic_ + numberGubColumns_; 2400 2361 // over estimate 2401 int numberElements = getNumElements() +startColumn_[numberGubColumns_]2362 int numberElements = getNumElements() + startColumn_[numberGubColumns_] 2402 2363 + numberGubColumns_; 2403 double * columnLower = new double[numberTotalColumns];2404 double * columnUpper = new double[numberTotalColumns];2405 double * cost = new double[numberTotalColumns];2406 double * rowLower = new double[numberTotalRows];2407 double * rowUpper = new double[numberTotalRows];2408 CoinBigIndex * start = new CoinBigIndex[numberTotalColumns+1];2409 int * row = new int[numberElements];2410 double * element = new double[numberElements];2364 double *columnLower = new double[numberTotalColumns]; 2365 double *columnUpper = new double[numberTotalColumns]; 2366 double *cost = new double[numberTotalColumns]; 2367 double *rowLower = new double[numberTotalRows]; 2368 double *rowUpper = new double[numberTotalRows]; 2369 CoinBigIndex *start = new CoinBigIndex[numberTotalColumns + 1]; 2370 int *row = new int[numberElements]; 2371 double *element = new double[numberElements]; 2411 2372 // Fill in 2412 const CoinBigIndex * 2413 const int * 2414 const int * 2415 const double * 2416 const double * 2417 const double * 2418 const double * 2419 const double * 2420 const double * 2421 start[0] =0;2422 numberElements =0;2423 for (int i =0;i<firstDynamic_;i++) {2373 const CoinBigIndex *startA = getVectorStarts(); 2374 const int *lengthA = getVectorLengths(); 2375 const int *rowA = getIndices(); 2376 const double *elementA = getElements(); 2377 const double *columnLowerA = model_>columnLower(); 2378 const double *columnUpperA = model_>columnUpper(); 2379 const double *costA = model_>objective(); 2380 const double *rowLowerA = model_>rowLower(); 2381 const double *rowUpperA = model_>rowUpper(); 2382 start[0] = 0; 2383 numberElements = 0; 2384 for (int i = 0; i < firstDynamic_; i++) { 2424 2385 columnLower[i] = columnLowerA[i]; 2425 2386 columnUpper[i] = columnUpperA[i]; 2426 2387 cost[i] = costA[i]; 2427 for (CoinBigIndex j = startA[i]; j<startA[i]+lengthA[i];j++) {2388 for (CoinBigIndex j = startA[i]; j < startA[i] + lengthA[i]; j++) { 2428 2389 row[numberElements] = rowA[j]; 2429 element[numberElements++] =elementA[j];2430 } 2431 start[i +1]=numberElements;2432 } 2433 for (int i =0;i<numberStaticRows_;i++) {2390 element[numberElements++] = elementA[j]; 2391 } 2392 start[i + 1] = numberElements; 2393 } 2394 for (int i = 0; i < numberStaticRows_; i++) { 2434 2395 rowLower[i] = rowLowerA[i]; 2435 2396 rowUpper[i] = rowUpperA[i]; 2436 2397 } 2437 int putC =firstDynamic_;2438 int putR =numberStaticRows_;2439 for (int i =0;i<numberSets_;i++) {2440 rowLower[putR] =lowerSet_[i];2441 rowUpper[putR] =upperSet_[i];2442 for (CoinBigIndex k =startSet_[i];k<startSet_[i+1];k++) {2443 columnLower[putC] =columnLower_[k];2444 columnUpper[putC] =columnUpper_[k];2445 cost[putC] =cost_[k];2398 int putC = firstDynamic_; 2399 int putR = numberStaticRows_; 2400 for (int i = 0; i < numberSets_; i++) { 2401 rowLower[putR] = lowerSet_[i]; 2402 rowUpper[putR] = upperSet_[i]; 2403 for (CoinBigIndex k = startSet_[i]; k < startSet_[i + 1]; k++) { 2404 columnLower[putC] = columnLower_[k]; 2405 columnUpper[putC] = columnUpper_[k]; 2406 cost[putC] = cost_[k]; 2446 2407 putC++; 2447 for (CoinBigIndex j = startColumn_[k]; j<startColumn_[k+1];j++) {2448 2449 element[numberElements++]=element_[j];2408 for (CoinBigIndex j = startColumn_[k]; j < startColumn_[k + 1]; j++) { 2409 row[numberElements] = row_[j]; 2410 element[numberElements++] = element_[j]; 2450 2411 } 2451 2412 row[numberElements] = putR; 2452 element[numberElements++] =1.0;2453 start[putC] =numberElements;2413 element[numberElements++] = 1.0; 2414 start[putC] = numberElements; 2454 2415 } 2455 2416 putR++; 2456 2417 } 2457 2418 2458 assert (putR==numberTotalRows);2459 assert (putC==numberTotalColumns);2419 assert(putR == numberTotalRows); 2420 assert(putC == numberTotalColumns); 2460 2421 ClpSimplex modelOut; 2461 modelOut.loadProblem(numberTotalColumns, numberTotalRows,2462 start,row,element,2463 columnLower,columnUpper,cost,2464 rowLower,rowUpper);2422 modelOut.loadProblem(numberTotalColumns, numberTotalRows, 2423 start, row, element, 2424 columnLower, columnUpper, cost, 2425 rowLower, rowUpper); 2465 2426 modelOut.writeMps(name); 2466 delete 2467 delete 2468 delete 2469 delete 2470 delete 2471 delete 2472 delete [] row;2473 delete 2427 delete[] columnLower; 2428 delete[] columnUpper; 2429 delete[] cost; 2430 delete[] rowLower; 2431 delete[] rowUpper; 2432 delete[] start; 2433 delete[] row; 2434 delete[] element; 2474 2435 } 2475 2436 // Adds in a column to gub structure (called from descendant) 2476 int 2477 ClpDynamicMatrix::addColumn(CoinBigIndex numberEntries, const int * row, const double * element, 2478 double cost, double lower, double upper, int iSet, 2479 DynamicStatus status) 2437 int ClpDynamicMatrix::addColumn(CoinBigIndex numberEntries, const int *row, const double *element, 2438 double cost, double lower, double upper, int iSet, 2439 DynamicStatus status) 2480 2440 { 2481 2482 2483 2484 if (startColumn_[j+1]  startColumn_[j] == numberEntries) {2485 const int *row2 = row_ + startColumn_[j];2486 const double *element2 = element_ + startColumn_[j];2487 2488 2489 2490 2491 2492 2493 2494 2495 2496 2497 2498 2499 2500 2501 2502 2503 2504 2505 2506 2507 2508 2509 2510 2511 2512 2513 2514 2441 // check if already in 2442 int j = startSet_[iSet]; 2443 while (j >= 0) { 2444 if (startColumn_[j + 1]  startColumn_[j] == numberEntries) { 2445 const int *row2 = row_ + startColumn_[j]; 2446 const double *element2 = element_ + startColumn_[j]; 2447 bool same = true; 2448 for (int k = 0; k < numberEntries; k++) { 2449 if (row[k] != row2[k]  element[k] != element2[k]) { 2450 same = false; 2451 break; 2452 } 2453 } 2454 if (same) { 2455 bool odd = false; 2456 if (cost != cost_[j]) 2457 odd = true; 2458 if (columnLower_ && lower != columnLower_[j]) 2459 odd = true; 2460 if (columnUpper_ && upper != columnUpper_[j]) 2461 odd = true; 2462 if (odd) { 2463 printf("seems odd  same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n", 2464 cost, lower, upper, cost_[j], 2465 columnLower_ ? columnLower_[j] : 0.0, 2466 columnUpper_ ? columnUpper_[j] : 1.0e100); 2467 } else { 2468 setDynamicStatus(j, status); 2469 return j; 2470 } 2471 } 2472 } 2473 j = next_[j]; 2474 } 2515 2475 2516 if (numberGubColumns_ == maximumGubColumns_  2517 startColumn_[numberGubColumns_] + numberEntries > maximumElements_) { 2518 CoinBigIndex j; 2519 int i; 2520 int put = 0; 2521 int numberElements = 0; 2522 CoinBigIndex start = 0; 2523 // compress  leave ones at ub and basic 2524 int * which = new int [numberGubColumns_]; 2525 for (i = 0; i < numberGubColumns_; i++) { 2526 CoinBigIndex end = startColumn_[i+1]; 2527 // what about ubs if column generation? 2528 if (getDynamicStatus(i) != atLowerBound) { 2529 // keep in 2530 for (j = start; j < end; j++) { 2531 row_[numberElements] = row_[j]; 2532 element_[numberElements++] = element_[j]; 2533 } 2534 startColumn_[put+1] = numberElements; 2535 cost_[put] = cost_[i]; 2536 if (columnLower_) 2537 columnLower_[put] = columnLower_[i]; 2538 if (columnUpper_) 2539 columnUpper_[put] = columnUpper_[i]; 2540 dynamicStatus_[put] = dynamicStatus_[i]; 2541 id_[put] = id_[i]; 2542 which[i] = put; 2543 put++; 2544 } else { 2545 // out 2546 which[i] = 1; 2547 } 2548 start = end; 2549 } 2550 // now redo startSet_ and next_ 2551 int * newNext = new int [maximumGubColumns_]; 2552 for (int jSet = 0; jSet < numberSets_; jSet++) { 2553 int sequence = startSet_[jSet]; 2554 while (which[sequence] < 0) { 2555 // out 2556 assert (next_[sequence] >= 0); 2557 sequence = next_[sequence]; 2558 } 2559 startSet_[jSet] = which[sequence]; 2560 int last = which[sequence]; 2561 while (next_[sequence] >= 0) { 2562 sequence = next_[sequence]; 2563 if(which[sequence] >= 0) { 2564 // keep 2565 newNext[last] = which[sequence]; 2566 last = which[sequence]; 2567 } 2568 } 2569 newNext[last] = jSet  1; 2570 } 2571 delete [] next_; 2572 next_ = newNext; 2573 delete [] which; 2574 abort(); 2575 } 2576 CoinBigIndex start = startColumn_[numberGubColumns_]; 2577 CoinMemcpyN(row, numberEntries, row_ + start); 2578 CoinMemcpyN(element, numberEntries, element_ + start); 2579 startColumn_[numberGubColumns_+1] = start + numberEntries; 2580 cost_[numberGubColumns_] = cost; 2581 if (columnLower_) 2582 columnLower_[numberGubColumns_] = lower; 2583 else 2584 assert (!lower); 2585 if (columnUpper_) 2586 columnUpper_[numberGubColumns_] = upper; 2587 else 2588 assert (upper > 1.0e20); 2589 setDynamicStatus(numberGubColumns_, status); 2590 // Do next_ 2591 j = startSet_[iSet]; 2592 startSet_[iSet] = numberGubColumns_; 2593 next_[numberGubColumns_] = j; 2594 numberGubColumns_++; 2595 return numberGubColumns_  1; 2476 if (numberGubColumns_ == maximumGubColumns_  startColumn_[numberGubColumns_] + numberEntries > maximumElements_) { 2477 CoinBigIndex j; 2478 int i; 2479 int put = 0; 2480 int numberElements = 0; 2481 CoinBigIndex start = 0; 2482 // compress  leave ones at ub and basic 2483 int *which = new int[numberGubColumns_]; 2484 for (i = 0; i < numberGubColumns_; i++) { 2485 CoinBigIndex end = startColumn_[i + 1]; 2486 // what about ubs if column generation? 2487 if (getDynamicStatus(i) != atLowerBound) { 2488 // keep in 2489 for (j = start; j < end; j++) { 2490 row_[numberElements] = row_[j]; 2491 element_[numberElements++] = element_[j]; 2492 } 2493 startColumn_[put + 1] = numberElements; 2494 cost_[put] = cost_[i]; 2495 if (columnLower_) 2496 columnLower_[put] = columnLower_[i]; 2497 if (columnUpper_) 2498 columnUpper_[put] = columnUpper_[i]; 2499 dynamicStatus_[put] = dynamicStatus_[i]; 2500 id_[put] = id_[i]; 2501 which[i] = put; 2502 put++; 2503 } else { 2504 // out 2505 which[i] = 1; 2506 } 2507 start = end; 2508 } 2509 // now redo startSet_ and next_ 2510 int *newNext = new int[maximumGubColumns_]; 2511 for (int jSet = 0; jSet < numberSets_; jSet++) { 2512 int sequence = startSet_[jSet]; 2513 while (which[sequence] < 0) { 2514 // out 2515 assert(next_[sequence] >= 0); 2516 sequence = next_[sequence]; 2517 } 2518 startSet_[jSet] = which[sequence]; 2519 int last = which[sequence]; 2520 while (next_[sequence] >= 0) { 2521 sequence = next_[sequence]; 2522 if (which[sequence] >= 0) { 2523 // keep 2524 newNext[last] = which[sequence]; 2525 last = which[sequence]; 2526 } 2527 } 2528 newNext[last] = jSet  1; 2529 } 2530 delete[] next_; 2531 next_ = newNext; 2532 delete[] which; 2533 abort(); 2534 } 2535 CoinBigIndex start = startColumn_[numberGubColumns_]; 2536 CoinMemcpyN(row, numberEntries, row_ + start); 2537 CoinMemcpyN(element, numberEntries, element_ + start); 2538 startColumn_[numberGubColumns_ + 1] = start + numberEntries; 2539 cost_[numberGubColumns_] = cost; 2540 if (columnLower_) 2541 columnLower_[numberGubColumns_] = lower; 2542 else 2543 assert(!lower); 2544 if (columnUpper_) 2545 columnUpper_[numberGubColumns_] = upper; 2546 else 2547 assert(upper > 1.0e20); 2548 setDynamicStatus(numberGubColumns_, status); 2549 // Do next_ 2550 j = startSet_[iSet]; 2551 startSet_[iSet] = numberGubColumns_; 2552 next_[numberGubColumns_] = j; 2553 numberGubColumns_++; 2554 return numberGubColumns_  1; 2596 2555 } 2597 2556 // Returns which set a variable is in 2598 int 2599 ClpDynamicMatrix::whichSet (int sequence) const 2557 int ClpDynamicMatrix::whichSet(int sequence) const 2600 2558 { 2601 2602 2603 int iSet = next_[sequence]  1;2604 2559 while (next_[sequence] >= 0) 2560 sequence = next_[sequence]; 2561 int iSet = next_[sequence]  1; 2562 return iSet; 2605 2563 } 2564 2565 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 2566 */
Note: See TracChangeset
for help on using the changeset viewer.