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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpDualRowSteepest.cpp
r2070 r2385 18 18 // Default Constructor 19 19 // 20 ClpDualRowSteepest::ClpDualRowSteepest 21 : ClpDualRowPivot(),22 state_(1),23 mode_(mode),24 persistence_(normal),25 weights_(NULL),26 infeasible_(NULL),27 alternateWeights_(NULL),28 savedWeights_(NULL),29 30 { 31 20 ClpDualRowSteepest::ClpDualRowSteepest(int mode) 21 : ClpDualRowPivot() 22 , state_(1) 23 , mode_(mode) 24 , persistence_(normal) 25 , weights_(NULL) 26 , infeasible_(NULL) 27 , alternateWeights_(NULL) 28 , savedWeights_(NULL) 29 , dubiousWeights_(NULL) 30 { 31 type_ = 2 + 64 * mode; 32 32 } 33 33 … … 35 35 // Copy constructor 36 36 // 37 ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest &rhs)38 39 { 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 37 ClpDualRowSteepest::ClpDualRowSteepest(const ClpDualRowSteepest &rhs) 38 : ClpDualRowPivot(rhs) 39 { 40 state_ = rhs.state_; 41 mode_ = rhs.mode_; 42 persistence_ = rhs.persistence_; 43 model_ = rhs.model_; 44 if ((model_ && model_>whatsChanged() & 1) != 0) { 45 int number = model_>numberRows(); 46 if (rhs.savedWeights_) 47 number = CoinMin(number, rhs.savedWeights_>capacity()); 48 if (rhs.infeasible_) { 49 infeasible_ = new CoinIndexedVector(rhs.infeasible_); 50 } else { 51 infeasible_ = NULL; 52 } 53 if (rhs.weights_) { 54 weights_ = new double[number]; 55 ClpDisjointCopyN(rhs.weights_, number, weights_); 56 } else { 57 weights_ = NULL; 58 } 59 if (rhs.alternateWeights_) { 60 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_); 61 } else { 62 alternateWeights_ = NULL; 63 } 64 if (rhs.savedWeights_) { 65 savedWeights_ = new CoinIndexedVector(rhs.savedWeights_); 66 } else { 67 savedWeights_ = NULL; 68 } 69 if (rhs.dubiousWeights_) { 70 assert(model_); 71 int number = model_>numberRows(); 72 dubiousWeights_ = new int[number]; 73 ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_); 74 } else { 75 dubiousWeights_ = NULL; 76 } 77 } else { 78 infeasible_ = NULL; 79 weights_ = NULL; 80 alternateWeights_ = NULL; 81 savedWeights_ = NULL; 82 dubiousWeights_ = NULL; 83 } 84 84 } 85 85 … … 87 87 // Destructor 88 88 // 89 ClpDualRowSteepest::~ClpDualRowSteepest 90 { 91 delete[] weights_;92 delete[] dubiousWeights_;93 94 95 89 ClpDualRowSteepest::~ClpDualRowSteepest() 90 { 91 delete[] weights_; 92 delete[] dubiousWeights_; 93 delete infeasible_; 94 delete alternateWeights_; 95 delete savedWeights_; 96 96 } 97 97 … … 100 100 // 101 101 ClpDualRowSteepest & 102 ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs) 103 { 104 if (this != &rhs) { 105 ClpDualRowPivot::operator=(rhs); 106 state_ = rhs.state_; 107 mode_ = rhs.mode_; 108 persistence_ = rhs.persistence_; 109 model_ = rhs.model_; 110 delete [] weights_; 111 delete [] dubiousWeights_; 112 delete infeasible_; 113 delete alternateWeights_; 114 delete savedWeights_; 115 assert(model_); 116 int number = model_>numberRows(); 117 if (rhs.savedWeights_) 118 number = CoinMin(number, rhs.savedWeights_>capacity()); 119 if (rhs.infeasible_ != NULL) { 120 infeasible_ = new CoinIndexedVector(rhs.infeasible_); 102 ClpDualRowSteepest::operator=(const ClpDualRowSteepest &rhs) 103 { 104 if (this != &rhs) { 105 ClpDualRowPivot::operator=(rhs); 106 state_ = rhs.state_; 107 mode_ = rhs.mode_; 108 persistence_ = rhs.persistence_; 109 model_ = rhs.model_; 110 delete[] weights_; 111 delete[] dubiousWeights_; 112 delete infeasible_; 113 delete alternateWeights_; 114 delete savedWeights_; 115 assert(model_); 116 int number = model_>numberRows(); 117 if (rhs.savedWeights_) 118 number = CoinMin(number, rhs.savedWeights_>capacity()); 119 if (rhs.infeasible_ != NULL) { 120 infeasible_ = new CoinIndexedVector(rhs.infeasible_); 121 } else { 122 infeasible_ = NULL; 123 } 124 if (rhs.weights_ != NULL) { 125 weights_ = new double[number]; 126 ClpDisjointCopyN(rhs.weights_, number, weights_); 127 } else { 128 weights_ = NULL; 129 } 130 if (rhs.alternateWeights_ != NULL) { 131 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_); 132 } else { 133 alternateWeights_ = NULL; 134 } 135 if (rhs.savedWeights_ != NULL) { 136 savedWeights_ = new CoinIndexedVector(rhs.savedWeights_); 137 } else { 138 savedWeights_ = NULL; 139 } 140 if (rhs.dubiousWeights_) { 141 assert(model_); 142 int number = model_>numberRows(); 143 dubiousWeights_ = new int[number]; 144 ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_); 145 } else { 146 dubiousWeights_ = NULL; 147 } 148 } 149 return *this; 150 } 151 // Fill most values 152 void ClpDualRowSteepest::fill(const ClpDualRowSteepest &rhs) 153 { 154 state_ = rhs.state_; 155 mode_ = rhs.mode_; 156 persistence_ = rhs.persistence_; 157 assert(model_>numberRows() == rhs.model_>numberRows()); 158 model_ = rhs.model_; 159 assert(model_); 160 int number = model_>numberRows(); 161 if (rhs.savedWeights_) 162 number = CoinMin(number, rhs.savedWeights_>capacity()); 163 if (rhs.infeasible_ != NULL) { 164 if (!infeasible_) 165 infeasible_ = new CoinIndexedVector(rhs.infeasible_); 166 else 167 *infeasible_ = *rhs.infeasible_; 168 } else { 169 delete infeasible_; 170 infeasible_ = NULL; 171 } 172 if (rhs.weights_ != NULL) { 173 if (!weights_) 174 weights_ = new double[number]; 175 ClpDisjointCopyN(rhs.weights_, number, weights_); 176 } else { 177 delete[] weights_; 178 weights_ = NULL; 179 } 180 if (rhs.alternateWeights_ != NULL) { 181 if (!alternateWeights_) 182 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_); 183 else 184 *alternateWeights_ = *rhs.alternateWeights_; 185 } else { 186 delete alternateWeights_; 187 alternateWeights_ = NULL; 188 } 189 if (rhs.savedWeights_ != NULL) { 190 if (!savedWeights_) 191 savedWeights_ = new CoinIndexedVector(rhs.savedWeights_); 192 else 193 *savedWeights_ = *rhs.savedWeights_; 194 } else { 195 delete savedWeights_; 196 savedWeights_ = NULL; 197 } 198 if (rhs.dubiousWeights_) { 199 assert(model_); 200 int number = model_>numberRows(); 201 if (!dubiousWeights_) 202 dubiousWeights_ = new int[number]; 203 ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_); 204 } else { 205 delete[] dubiousWeights_; 206 dubiousWeights_ = NULL; 207 } 208 } 209 // Returns pivot row, 1 if none 210 int ClpDualRowSteepest::pivotRow() 211 { 212 assert(model_); 213 int i, iRow; 214 double *infeas = infeasible_>denseVector(); 215 double largest = 0.0; 216 int *index = infeasible_>getIndices(); 217 int number = infeasible_>getNumElements(); 218 const int *pivotVariable = model_>pivotVariable(); 219 int chosenRow = 1; 220 int lastPivotRow = model_>pivotRow(); 221 assert(lastPivotRow < model_>numberRows()); 222 double tolerance = model_>currentPrimalTolerance(); 223 // we can't really trust infeasibilities if there is primal error 224 // this coding has to mimic coding in checkPrimalSolution 225 double error = CoinMin(1.0e2, model_>largestPrimalError()); 226 // allow tolerance at least slightly bigger than standard 227 tolerance = tolerance + error; 228 // But cap 229 tolerance = CoinMin(1000.0, tolerance); 230 tolerance *= tolerance; // as we are using squares 231 bool toleranceChanged = false; 232 double *solution = model_>solutionRegion(); 233 double *lower = model_>lowerRegion(); 234 double *upper = model_>upperRegion(); 235 // do last pivot row one here 236 //#define CLP_DUAL_FIXED_COLUMN_MULTIPLIER 10.0 237 if (lastPivotRow >= 0 && lastPivotRow < model_>numberRows()) { 238 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 239 int numberColumns = model_>numberColumns(); 240 #endif 241 int iPivot = pivotVariable[lastPivotRow]; 242 double value = solution[iPivot]; 243 double lower = model_>lower(iPivot); 244 double upper = model_>upper(iPivot); 245 if (value > upper + tolerance) { 246 value = upper; 247 value *= value; 248 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 249 if (iPivot < numberColumns) 250 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 251 #endif 252 // store square in list 253 if (infeas[lastPivotRow]) 254 infeas[lastPivotRow] = value; // already there 255 else 256 infeasible_>quickAdd(lastPivotRow, value); 257 } else if (value < lower  tolerance) { 258 value = lower; 259 value *= value; 260 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 261 if (iPivot < numberColumns) 262 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 263 #endif 264 // store square in list 265 if (infeas[lastPivotRow]) 266 infeas[lastPivotRow] = value; // already there 267 else 268 infeasible_>add(lastPivotRow, value); 269 } else { 270 // feasible  was it infeasible  if so set tiny 271 if (infeas[lastPivotRow]) 272 infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 273 } 274 number = infeasible_>getNumElements(); 275 } 276 if (model_>numberIterations() < model_>lastBadIteration() + 200) { 277 // we can't really trust infeasibilities if there is dual error 278 if (model_>largestDualError() > model_>largestPrimalError()) { 279 tolerance *= CoinMin(model_>largestDualError() / model_>largestPrimalError(), 1000.0); 280 toleranceChanged = true; 281 } 282 } 283 int numberWanted; 284 if (mode_ < 2) { 285 numberWanted = number + 1; 286 } else if (mode_ == 2) { 287 numberWanted = CoinMax(2000, number / 8); 288 } else { 289 int numberElements = model_>factorization()>numberElements(); 290 double ratio = static_cast< double >(numberElements) / static_cast< double >(model_>numberRows()); 291 numberWanted = CoinMax(2000, number / 8); 292 if (ratio < 1.0) { 293 numberWanted = CoinMax(2000, number / 20); 294 } else if (ratio > 10.0) { 295 ratio = number * (ratio / 80.0); 296 if (ratio > number) 297 numberWanted = number + 1; 298 else 299 numberWanted = CoinMax(2000, static_cast< int >(ratio)); 300 } 301 } 302 if (model_>largestPrimalError() > 1.0e3) 303 numberWanted = number + 1; // be safe 304 int iPass; 305 // Setup two passes 306 int start[4]; 307 start[1] = number; 308 start[2] = 0; 309 double dstart = static_cast< double >(number) * model_>randomNumberGenerator()>randomDouble(); 310 start[0] = static_cast< int >(dstart); 311 start[3] = start[0]; 312 //double largestWeight=0.0; 313 //double smallestWeight=1.0e100; 314 for (iPass = 0; iPass < 2; iPass++) { 315 int end = start[2 * iPass + 1]; 316 for (i = start[2 * iPass]; i < end; i++) { 317 iRow = index[i]; 318 double value = infeas[iRow]; 319 if (value > tolerance) { 320 //#define OUT_EQ 321 #ifdef OUT_EQ 322 { 323 int iSequence = pivotVariable[iRow]; 324 if (upper[iSequence] == lower[iSequence]) 325 value *= 2.0; 326 } 327 #endif 328 double weight = CoinMin(weights_[iRow], 1.0e50); 329 //largestWeight = CoinMax(largestWeight,weight); 330 //smallestWeight = CoinMin(smallestWeight,weight); 331 //double dubious = dubiousWeights_[iRow]; 332 //weight *= dubious; 333 //if (value>2.0*largest*weight(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) { 334 if (value > largest * weight) { 335 // make last pivot row last resort choice 336 if (iRow == lastPivotRow) { 337 if (value * 1.0e10 < largest * weight) 338 continue; 339 else 340 value *= 1.0e10; 341 } 342 int iSequence = pivotVariable[iRow]; 343 if (!model_>flagged(iSequence)) { 344 //#define CLP_DEBUG 3 345 #ifdef CLP_DEBUG 346 double value2 = 0.0; 347 if (solution[iSequence] > upper[iSequence] + tolerance) 348 value2 = solution[iSequence]  upper[iSequence]; 349 else if (solution[iSequence] < lower[iSequence]  tolerance) 350 value2 = solution[iSequence]  lower[iSequence]; 351 assert(fabs(value2 * value2  infeas[iRow]) < 1.0e8 * CoinMin(value2 * value2, infeas[iRow])); 352 #endif 353 if (solution[iSequence] > upper[iSequence] + tolerance  solution[iSequence] < lower[iSequence]  tolerance) { 354 chosenRow = iRow; 355 largest = value / weight; 356 //largestWeight = dubious; 357 } 121 358 } else { 122 infeasible_ = NULL; 359 // just to make sure we don't exit before got something 360 numberWanted++; 123 361 } 124 if (rhs.weights_ != NULL) { 125 weights_ = new double[number]; 126 ClpDisjointCopyN(rhs.weights_, number, weights_); 127 } else { 128 weights_ = NULL; 129 } 130 if (rhs.alternateWeights_ != NULL) { 131 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_); 132 } else { 133 alternateWeights_ = NULL; 134 } 135 if (rhs.savedWeights_ != NULL) { 136 savedWeights_ = new CoinIndexedVector(rhs.savedWeights_); 137 } else { 138 savedWeights_ = NULL; 139 } 140 if (rhs.dubiousWeights_) { 141 assert(model_); 142 int number = model_>numberRows(); 143 dubiousWeights_ = new int[number]; 144 ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_); 145 } else { 146 dubiousWeights_ = NULL; 147 } 148 } 149 return *this; 150 } 151 // Fill most values 152 void 153 ClpDualRowSteepest::fill(const ClpDualRowSteepest& rhs) 154 { 155 state_ = rhs.state_; 156 mode_ = rhs.mode_; 157 persistence_ = rhs.persistence_; 158 assert (model_>numberRows() == rhs.model_>numberRows()); 159 model_ = rhs.model_; 160 assert(model_); 161 int number = model_>numberRows(); 162 if (rhs.savedWeights_) 163 number = CoinMin(number, rhs.savedWeights_>capacity()); 164 if (rhs.infeasible_ != NULL) { 165 if (!infeasible_) 166 infeasible_ = new CoinIndexedVector(rhs.infeasible_); 167 else 168 *infeasible_ = *rhs.infeasible_; 169 } else { 170 delete infeasible_; 171 infeasible_ = NULL; 172 } 173 if (rhs.weights_ != NULL) { 174 if (!weights_) 175 weights_ = new double[number]; 176 ClpDisjointCopyN(rhs.weights_, number, weights_); 177 } else { 178 delete [] weights_; 179 weights_ = NULL; 180 } 181 if (rhs.alternateWeights_ != NULL) { 182 if (!alternateWeights_) 183 alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_); 184 else 185 *alternateWeights_ = *rhs.alternateWeights_; 186 } else { 187 delete alternateWeights_; 188 alternateWeights_ = NULL; 189 } 190 if (rhs.savedWeights_ != NULL) { 191 if (!savedWeights_) 192 savedWeights_ = new CoinIndexedVector(rhs.savedWeights_); 193 else 194 *savedWeights_ = *rhs.savedWeights_; 195 } else { 196 delete savedWeights_; 197 savedWeights_ = NULL; 198 } 199 if (rhs.dubiousWeights_) { 200 assert(model_); 201 int number = model_>numberRows(); 202 if (!dubiousWeights_) 203 dubiousWeights_ = new int[number]; 204 ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_); 205 } else { 206 delete [] dubiousWeights_; 207 dubiousWeights_ = NULL; 208 } 209 } 210 // Returns pivot row, 1 if none 211 int 212 ClpDualRowSteepest::pivotRow() 213 { 214 assert(model_); 215 int i, iRow; 216 double * infeas = infeasible_>denseVector(); 217 double largest = 0.0; 218 int * index = infeasible_>getIndices(); 219 int number = infeasible_>getNumElements(); 220 const int * pivotVariable = model_>pivotVariable(); 221 int chosenRow = 1; 222 int lastPivotRow = model_>pivotRow(); 223 assert (lastPivotRow < model_>numberRows()); 224 double tolerance = model_>currentPrimalTolerance(); 225 // we can't really trust infeasibilities if there is primal error 226 // this coding has to mimic coding in checkPrimalSolution 227 double error = CoinMin(1.0e2, model_>largestPrimalError()); 228 // allow tolerance at least slightly bigger than standard 229 tolerance = tolerance + error; 230 // But cap 231 tolerance = CoinMin(1000.0, tolerance); 232 tolerance *= tolerance; // as we are using squares 233 bool toleranceChanged = false; 234 double * solution = model_>solutionRegion(); 235 double * lower = model_>lowerRegion(); 236 double * upper = model_>upperRegion(); 237 // do last pivot row one here 238 //#define CLP_DUAL_FIXED_COLUMN_MULTIPLIER 10.0 239 if (lastPivotRow >= 0 && lastPivotRow < model_>numberRows()) { 240 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 241 int numberColumns = model_>numberColumns(); 242 #endif 243 int iPivot = pivotVariable[lastPivotRow]; 244 double value = solution[iPivot]; 245 double lower = model_>lower(iPivot); 246 double upper = model_>upper(iPivot); 247 if (value > upper + tolerance) { 248 value = upper; 249 value *= value; 250 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 251 if (iPivot < numberColumns) 252 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 253 #endif 254 // store square in list 255 if (infeas[lastPivotRow]) 256 infeas[lastPivotRow] = value; // already there 257 else 258 infeasible_>quickAdd(lastPivotRow, value); 259 } else if (value < lower  tolerance) { 260 value = lower; 261 value *= value; 262 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 263 if (iPivot < numberColumns) 264 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 265 #endif 266 // store square in list 267 if (infeas[lastPivotRow]) 268 infeas[lastPivotRow] = value; // already there 269 else 270 infeasible_>add(lastPivotRow, value); 271 } else { 272 // feasible  was it infeasible  if so set tiny 273 if (infeas[lastPivotRow]) 274 infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 275 } 276 number = infeasible_>getNumElements(); 277 } 278 if(model_>numberIterations() < model_>lastBadIteration() + 200) { 279 // we can't really trust infeasibilities if there is dual error 280 if (model_>largestDualError() > model_>largestPrimalError()) { 281 tolerance *= CoinMin(model_>largestDualError() / model_>largestPrimalError(), 1000.0); 282 toleranceChanged = true; 283 } 284 } 285 int numberWanted; 286 if (mode_ < 2 ) { 287 numberWanted = number + 1; 288 } else if (mode_ == 2) { 289 numberWanted = CoinMax(2000, number / 8); 290 } else { 291 int numberElements = model_>factorization()>numberElements(); 292 double ratio = static_cast<double> (numberElements) / 293 static_cast<double> (model_>numberRows()); 294 numberWanted = CoinMax(2000, number / 8); 295 if (ratio < 1.0) { 296 numberWanted = CoinMax(2000, number / 20); 297 } else if (ratio > 10.0) { 298 ratio = number * (ratio / 80.0); 299 if (ratio > number) 300 numberWanted = number + 1; 301 else 302 numberWanted = CoinMax(2000, static_cast<int> (ratio)); 303 } 304 } 305 if (model_>largestPrimalError() > 1.0e3) 306 numberWanted = number + 1; // be safe 307 int iPass; 308 // Setup two passes 309 int start[4]; 310 start[1] = number; 311 start[2] = 0; 312 double dstart = static_cast<double> (number) * model_>randomNumberGenerator()>randomDouble(); 313 start[0] = static_cast<int> (dstart); 314 start[3] = start[0]; 315 //double largestWeight=0.0; 316 //double smallestWeight=1.0e100; 317 for (iPass = 0; iPass < 2; iPass++) { 318 int end = start[2*iPass+1]; 319 for (i = start[2*iPass]; i < end; i++) { 320 iRow = index[i]; 321 double value = infeas[iRow]; 322 if (value > tolerance) { 323 //#define OUT_EQ 324 #ifdef OUT_EQ 325 { 326 int iSequence = pivotVariable[iRow]; 327 if (upper[iSequence] == lower[iSequence]) 328 value *= 2.0; 329 } 330 #endif 331 double weight = CoinMin(weights_[iRow], 1.0e50); 332 //largestWeight = CoinMax(largestWeight,weight); 333 //smallestWeight = CoinMin(smallestWeight,weight); 334 //double dubious = dubiousWeights_[iRow]; 335 //weight *= dubious; 336 //if (value>2.0*largest*weight(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) { 337 if (value > largest * weight) { 338 // make last pivot row last resort choice 339 if (iRow == lastPivotRow) { 340 if (value * 1.0e10 < largest * weight) 341 continue; 342 else 343 value *= 1.0e10; 344 } 345 int iSequence = pivotVariable[iRow]; 346 if (!model_>flagged(iSequence)) { 347 //#define CLP_DEBUG 3 348 #ifdef CLP_DEBUG 349 double value2 = 0.0; 350 if (solution[iSequence] > upper[iSequence] + tolerance) 351 value2 = solution[iSequence]  upper[iSequence]; 352 else if (solution[iSequence] < lower[iSequence]  tolerance) 353 value2 = solution[iSequence]  lower[iSequence]; 354 assert(fabs(value2 * value2  infeas[iRow]) < 1.0e8 * CoinMin(value2 * value2, infeas[iRow])); 355 #endif 356 if (solution[iSequence] > upper[iSequence] + tolerance  357 solution[iSequence] < lower[iSequence]  tolerance) { 358 chosenRow = iRow; 359 largest = value / weight; 360 //largestWeight = dubious; 361 } 362 } else { 363 // just to make sure we don't exit before got something 364 numberWanted++; 365 } 366 } 367 numberWanted; 368 if (!numberWanted) 369 break; 370 } 371 } 372 if (!numberWanted) 373 break; 374 } 375 //printf("smallest %g largest %g\n",smallestWeight,largestWeight); 376 if (chosenRow < 0 && toleranceChanged) { 377 // won't line up with checkPrimalSolution  do again 378 double saveError = model_>largestDualError(); 379 model_>setLargestDualError(0.0); 380 // can't loop 381 chosenRow = pivotRow(); 382 model_>setLargestDualError(saveError); 383 } 384 if (chosenRow < 0 && lastPivotRow < 0) { 385 int nLeft = 0; 386 for (int i = 0; i < number; i++) { 387 int iRow = index[i]; 388 if (fabs(infeas[iRow]) > 1.0e50) { 389 index[nLeft++] = iRow; 390 } else { 391 infeas[iRow] = 0.0; 392 } 393 } 394 infeasible_>setNumElements(nLeft); 395 model_>setNumberPrimalInfeasibilities(nLeft); 396 } 397 return chosenRow; 362 } 363 numberWanted; 364 if (!numberWanted) 365 break; 366 } 367 } 368 if (!numberWanted) 369 break; 370 } 371 //printf("smallest %g largest %g\n",smallestWeight,largestWeight); 372 if (chosenRow < 0 && toleranceChanged) { 373 // won't line up with checkPrimalSolution  do again 374 double saveError = model_>largestDualError(); 375 model_>setLargestDualError(0.0); 376 // can't loop 377 chosenRow = pivotRow(); 378 model_>setLargestDualError(saveError); 379 } 380 if (chosenRow < 0 && lastPivotRow < 0) { 381 int nLeft = 0; 382 for (int i = 0; i < number; i++) { 383 int iRow = index[i]; 384 if (fabs(infeas[iRow]) > 1.0e50) { 385 index[nLeft++] = iRow; 386 } else { 387 infeas[iRow] = 0.0; 388 } 389 } 390 infeasible_>setNumElements(nLeft); 391 model_>setNumberPrimalInfeasibilities(nLeft); 392 } 393 return chosenRow; 398 394 } 399 395 #if 0 … … 407 403 Also does FT update */ 408 404 double 409 ClpDualRowSteepest::updateWeights(CoinIndexedVector * input, 410 CoinIndexedVector * spare, 411 CoinIndexedVector * spare2, 412 CoinIndexedVector * updatedColumn) 413 { 414 //#define CLP_DEBUG 3 415 #if CLP_DEBUG>2 416 // Very expensive debug 417 { 418 int numberRows = model_>numberRows(); 419 CoinIndexedVector * temp = new CoinIndexedVector(); 420 temp>reserve(numberRows + 421 model_>factorization()>maximumPivots()); 422 double * array = alternateWeights_>denseVector(); 423 int * which = alternateWeights_>getIndices(); 424 alternateWeights_>clear(); 425 int i; 426 for (i = 0; i < numberRows; i++) { 427 double value = 0.0; 428 array[i] = 1.0; 429 which[0] = i; 430 alternateWeights_>setNumElements(1); 431 model_>factorization()>updateColumnTranspose(temp, 432 alternateWeights_); 433 int number = alternateWeights_>getNumElements(); 434 int j; 435 for (j = 0; j < number; j++) { 436 int iRow = which[j]; 437 value += array[iRow] * array[iRow]; 438 array[iRow] = 0.0; 439 } 440 alternateWeights_>setNumElements(0); 441 double w = CoinMax(weights_[i], value) * .1; 442 if (fabs(weights_[i]  value) > w) { 443 printf("%d old %g, true %g\n", i, weights_[i], value); 444 weights_[i] = value; // to reduce printout 445 } 446 //else 447 //printf("%d matches %g\n",i,value); 448 } 449 delete temp; 450 } 451 #endif 452 assert (input>packedMode()); 453 if (!updatedColumn>packedMode()) { 454 // I think this means empty 405 ClpDualRowSteepest::updateWeights(CoinIndexedVector *input, 406 CoinIndexedVector *spare, 407 CoinIndexedVector *spare2, 408 CoinIndexedVector *updatedColumn) 409 { 410 //#define CLP_DEBUG 3 411 #if CLP_DEBUG > 2 412 // Very expensive debug 413 { 414 int numberRows = model_>numberRows(); 415 CoinIndexedVector *temp = new CoinIndexedVector(); 416 temp>reserve(numberRows + model_>factorization()>maximumPivots()); 417 double *array = alternateWeights_>denseVector(); 418 int *which = alternateWeights_>getIndices(); 419 alternateWeights_>clear(); 420 int i; 421 for (i = 0; i < numberRows; i++) { 422 double value = 0.0; 423 array[i] = 1.0; 424 which[0] = i; 425 alternateWeights_>setNumElements(1); 426 model_>factorization()>updateColumnTranspose(temp, 427 alternateWeights_); 428 int number = alternateWeights_>getNumElements(); 429 int j; 430 for (j = 0; j < number; j++) { 431 int iRow = which[j]; 432 value += array[iRow] * array[iRow]; 433 array[iRow] = 0.0; 434 } 435 alternateWeights_>setNumElements(0); 436 double w = CoinMax(weights_[i], value) * .1; 437 if (fabs(weights_[i]  value) > w) { 438 printf("%d old %g, true %g\n", i, weights_[i], value); 439 weights_[i] = value; // to reduce printout 440 } 441 //else 442 //printf("%d matches %g\n",i,value); 443 } 444 delete temp; 445 } 446 #endif 447 assert(input>packedMode()); 448 if (!updatedColumn>packedMode()) { 449 // I think this means empty 455 450 #ifdef COIN_DEVELOP 456 457 #endif 458 459 460 461 462 463 464 465 466 double *work = input>denseVector();467 468 int *which = input>getIndices();469 double *work2 = spare>denseVector();470 int *which2 = spare>getIndices();471 472 473 474 475 476 477 478 for ( i = 0; i < numberNonZero; i ++) {479 480 481 482 483 484 485 486 487 for ( i = 0; i < numberNonZero; i ++) {488 489 490 491 492 493 494 495 496 spare>setNumElements ( numberNonZero);497 451 printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n"); 452 #endif 453 return 0.0; 454 } 455 double alpha = 0.0; 456 if (!model_>factorization()>networkBasis()) { 457 // clear other region 458 alternateWeights_>clear(); 459 double norm = 0.0; 460 int i; 461 double *work = input>denseVector(); 462 int numberNonZero = input>getNumElements(); 463 int *which = input>getIndices(); 464 double *work2 = spare>denseVector(); 465 int *which2 = spare>getIndices(); 466 // ftran 467 //permute and move indices into index array 468 //also compute norm 469 //int *regionIndex = alternateWeights_>getIndices ( ); 470 const int *permute = model_>factorization()>permute(); 471 //double * region = alternateWeights_>denseVector(); 472 if (permute) { 473 for (i = 0; i < numberNonZero; i++) { 474 int iRow = which[i]; 475 double value = work[i]; 476 norm += value * value; 477 iRow = permute[iRow]; 478 work2[iRow] = value; 479 which2[i] = iRow; 480 } 481 } else { 482 for (i = 0; i < numberNonZero; i++) { 483 int iRow = which[i]; 484 double value = work[i]; 485 norm += value * value; 486 //iRow = permute[iRow]; 487 work2[iRow] = value; 488 which2[i] = iRow; 489 } 490 } 491 spare>setNumElements(numberNonZero); 492 // Do FT update 498 493 #if 0 499 494 ft_count_in += updatedColumn>getNumElements(); 500 495 up_count_in += spare>getNumElements(); 501 496 #endif 502 503 #if CLP_DEBUG >2504 505 506 #endif 507 508 509 #if CLP_DEBUG >2510 511 512 #endif 513 514 515 516 517 497 if (permute  true) { 498 #if CLP_DEBUG > 2 499 printf("REGION before %d els\n", spare>getNumElements()); 500 spare>print(); 501 #endif 502 model_>factorization()>updateTwoColumnsFT(spare2, updatedColumn, 503 spare, permute != NULL); 504 #if CLP_DEBUG > 2 505 printf("REGION after %d els\n", spare>getNumElements()); 506 spare>print(); 507 #endif 508 } else { 509 // Leave as old way 510 model_>factorization()>updateColumnFT(spare2, updatedColumn); 511 model_>factorization()>updateColumn(spare2, spare, false); 512 } 518 513 #undef CLP_DEBUG 519 514 #if 0 … … 525 520 ft_count_in, up_count_in); 526 521 #endif 527 528 529 522 numberNonZero = spare>getNumElements(); 523 // alternateWeights_ should still be empty 524 int pivotRow = model_>pivotRow(); 530 525 #ifdef CLP_DEBUG 531 if ( model_>logLevel ( ) > 4 && 532 fabs(norm  weights_[pivotRow]) > 1.0e3 * (1.0 + norm)) 533 printf("on row %d, true weight %g, old %g\n", 534 pivotRow, sqrt(norm), sqrt(weights_[pivotRow])); 535 #endif 536 // could reinitialize here (could be expensive) 537 norm /= model_>alpha() * model_>alpha(); 538 assert(model_>alpha()); 539 assert(norm); 540 // pivot element 541 alpha = 0.0; 542 double multiplier = 2.0 / model_>alpha(); 543 // look at updated column 544 work = updatedColumn>denseVector(); 545 numberNonZero = updatedColumn>getNumElements(); 546 which = updatedColumn>getIndices(); 526 if (model_>logLevel() > 4 && fabs(norm  weights_[pivotRow]) > 1.0e3 * (1.0 + norm)) 527 printf("on row %d, true weight %g, old %g\n", 528 pivotRow, sqrt(norm), sqrt(weights_[pivotRow])); 529 #endif 530 // could reinitialize here (could be expensive) 531 norm /= model_>alpha() * model_>alpha(); 532 assert(model_>alpha()); 533 assert(norm); 534 // pivot element 535 alpha = 0.0; 536 double multiplier = 2.0 / model_>alpha(); 537 // look at updated column 538 work = updatedColumn>denseVector(); 539 numberNonZero = updatedColumn>getNumElements(); 540 which = updatedColumn>getIndices(); 547 541 548 549 double *work3 = alternateWeights_>denseVector();550 int *which3 = alternateWeights_>getIndices();551 const int *pivotColumn = model_>factorization()>pivotColumn();552 553 554 555 556 557 558 559 560 561 562 563 devex +=theta * (theta * norm + value * multiplier);564 565 566 567 568 569 570 571 572 573 574 575 542 int nSave = 0; 543 double *work3 = alternateWeights_>denseVector(); 544 int *which3 = alternateWeights_>getIndices(); 545 const int *pivotColumn = model_>factorization()>pivotColumn(); 546 for (i = 0; i < numberNonZero; i++) { 547 int iRow = which[i]; 548 double theta = work[i]; 549 if (iRow == pivotRow) 550 alpha = theta; 551 double devex = weights_[iRow]; 552 work3[nSave] = devex; // save old 553 which3[nSave++] = iRow; 554 // transform to match spare 555 int jRow = permute ? pivotColumn[iRow] : iRow; 556 double value = work2[jRow]; 557 devex += theta * (theta * norm + value * multiplier); 558 if (devex < DEVEX_TRY_NORM) 559 devex = DEVEX_TRY_NORM; 560 weights_[iRow] = devex; 561 } 562 alternateWeights_>setPackedMode(true); 563 alternateWeights_>setNumElements(nSave); 564 if (norm < DEVEX_TRY_NORM) 565 norm = DEVEX_TRY_NORM; 566 // Try this to make less likely will happen again and stop cycling 567 //norm *= 1.02; 568 weights_[pivotRow] = norm; 569 spare>clear(); 576 570 #ifdef CLP_DEBUG 577 578 #endif 579 580 581 582 583 584 585 586 double *work = input>denseVector();587 588 int *which = input>getIndices();589 double *work2 = spare>denseVector();590 int *which2 = spare>getIndices();591 592 593 594 595 596 597 598 599 571 spare>checkClear(); 572 #endif 573 } else { 574 // Do FT update 575 model_>factorization()>updateColumnFT(spare, updatedColumn); 576 // clear other region 577 alternateWeights_>clear(); 578 double norm = 0.0; 579 int i; 580 double *work = input>denseVector(); 581 int number = input>getNumElements(); 582 int *which = input>getIndices(); 583 double *work2 = spare>denseVector(); 584 int *which2 = spare>getIndices(); 585 for (i = 0; i < number; i++) { 586 int iRow = which[i]; 587 double value = work[i]; 588 norm += value * value; 589 work2[iRow] = value; 590 which2[i] = iRow; 591 } 592 spare>setNumElements(number); 593 // ftran 600 594 #ifndef NDEBUG 601 602 #endif 603 604 595 alternateWeights_>checkClear(); 596 #endif 597 model_>factorization()>updateColumn(alternateWeights_, spare); 598 // alternateWeights_ should still be empty 605 599 #ifndef NDEBUG 606 607 #endif 608 600 alternateWeights_>checkClear(); 601 #endif 602 int pivotRow = model_>pivotRow(); 609 603 #ifdef CLP_DEBUG 610 if ( model_>logLevel ( ) > 4 && 611 fabs(norm  weights_[pivotRow]) > 1.0e3 * (1.0 + norm)) 612 printf("on row %d, true weight %g, old %g\n", 613 pivotRow, sqrt(norm), sqrt(weights_[pivotRow])); 614 #endif 615 // could reinitialize here (could be expensive) 616 norm /= model_>alpha() * model_>alpha(); 604 if (model_>logLevel() > 4 && fabs(norm  weights_[pivotRow]) > 1.0e3 * (1.0 + norm)) 605 printf("on row %d, true weight %g, old %g\n", 606 pivotRow, sqrt(norm), sqrt(weights_[pivotRow])); 607 #endif 608 // could reinitialize here (could be expensive) 609 norm /= model_>alpha() * model_>alpha(); 617 610 618 619 620 621 622 623 624 625 626 627 611 assert(norm); 612 //if (norm < DEVEX_TRY_NORM) 613 //norm = DEVEX_TRY_NORM; 614 // pivot element 615 alpha = 0.0; 616 double multiplier = 2.0 / model_>alpha(); 617 // look at updated column 618 work = updatedColumn>denseVector(); 619 number = updatedColumn>getNumElements(); 620 which = updatedColumn>getIndices(); 628 621 629 630 double *work3 = alternateWeights_>denseVector();631 int *which3 = alternateWeights_>getIndices();632 633 634 635 636 637 638 639 640 641 devex +=theta * (theta * norm + value * multiplier);642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 622 int nSave = 0; 623 double *work3 = alternateWeights_>denseVector(); 624 int *which3 = alternateWeights_>getIndices(); 625 for (i = 0; i < number; i++) { 626 int iRow = which[i]; 627 double theta = work[i]; 628 if (iRow == pivotRow) 629 alpha = theta; 630 double devex = weights_[iRow]; 631 work3[nSave] = devex; // save old 632 which3[nSave++] = iRow; 633 double value = work2[iRow]; 634 devex += theta * (theta * norm + value * multiplier); 635 if (devex < DEVEX_TRY_NORM) 636 devex = DEVEX_TRY_NORM; 637 weights_[iRow] = devex; 638 } 639 if (!alpha) { 640 // error  but carry on 641 alpha = 1.0e50; 642 } 643 alternateWeights_>setPackedMode(true); 644 alternateWeights_>setNumElements(nSave); 645 if (norm < DEVEX_TRY_NORM) 646 norm = DEVEX_TRY_NORM; 647 weights_[pivotRow] = norm; 648 spare>clear(); 649 } 657 650 #ifdef CLP_DEBUG 658 659 #endif 660 651 spare>checkClear(); 652 #endif 653 return alpha; 661 654 } 662 655 … … 665 658 Computes change in objective function 666 659 */ 667 void 668 ClpDualRowSteepest::updatePrimalSolution( 669 CoinIndexedVector * primalUpdate, 670 double primalRatio, 671 double & objectiveChange) 672 { 673 double * COIN_RESTRICT work = primalUpdate>denseVector(); 674 int number = primalUpdate>getNumElements(); 675 int * COIN_RESTRICT which = primalUpdate>getIndices(); 676 int i; 677 double changeObj = 0.0; 678 double tolerance = model_>currentPrimalTolerance(); 679 const int * COIN_RESTRICT pivotVariable = model_>pivotVariable(); 680 double * COIN_RESTRICT infeas = infeasible_>denseVector(); 681 double * COIN_RESTRICT solution = model_>solutionRegion(); 682 const double * COIN_RESTRICT costModel = model_>costRegion(); 683 const double * COIN_RESTRICT lowerModel = model_>lowerRegion(); 684 const double * COIN_RESTRICT upperModel = model_>upperRegion(); 660 void ClpDualRowSteepest::updatePrimalSolution( 661 CoinIndexedVector *primalUpdate, 662 double primalRatio, 663 double &objectiveChange) 664 { 665 double *COIN_RESTRICT work = primalUpdate>denseVector(); 666 int number = primalUpdate>getNumElements(); 667 int *COIN_RESTRICT which = primalUpdate>getIndices(); 668 int i; 669 double changeObj = 0.0; 670 double tolerance = model_>currentPrimalTolerance(); 671 const int *COIN_RESTRICT pivotVariable = model_>pivotVariable(); 672 double *COIN_RESTRICT infeas = infeasible_>denseVector(); 673 double *COIN_RESTRICT solution = model_>solutionRegion(); 674 const double *COIN_RESTRICT costModel = model_>costRegion(); 675 const double *COIN_RESTRICT lowerModel = model_>lowerRegion(); 676 const double *COIN_RESTRICT upperModel = model_>upperRegion(); 685 677 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 686 687 #endif 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 678 int numberColumns = model_>numberColumns(); 679 #endif 680 if (primalUpdate>packedMode()) { 681 for (i = 0; i < number; i++) { 682 int iRow = which[i]; 683 int iPivot = pivotVariable[iRow]; 684 double value = solution[iPivot]; 685 double cost = costModel[iPivot]; 686 double change = primalRatio * work[i]; 687 work[i] = 0.0; 688 value = change; 689 changeObj = change * cost; 690 double lower = lowerModel[iPivot]; 691 double upper = upperModel[iPivot]; 692 solution[iPivot] = value; 693 if (value < lower  tolerance) { 694 value = lower; 695 value *= value; 704 696 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 705 706 697 if (iPivot < numberColumns) 698 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 707 699 #endif 708 700 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 709 710 711 #endif 712 713 714 715 716 717 718 719 701 if (lower == upper) 702 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 703 #endif 704 // store square in list 705 if (infeas[iRow]) 706 infeas[iRow] = value; // already there 707 else 708 infeasible_>quickAdd(iRow, value); 709 } else if (value > upper + tolerance) { 710 value = upper; 711 value *= value; 720 712 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 721 722 713 if (iPivot < numberColumns) 714 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 723 715 #endif 724 716 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 725 726 727 #endif 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 717 if (lower == upper) 718 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 719 #endif 720 // store square in list 721 if (infeas[iRow]) 722 infeas[iRow] = value; // already there 723 else 724 infeasible_>quickAdd(iRow, value); 725 } else { 726 // feasible  was it infeasible  if so set tiny 727 if (infeas[iRow]) 728 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 729 } 730 } 731 } else { 732 for (i = 0; i < number; i++) { 733 int iRow = which[i]; 734 int iPivot = pivotVariable[iRow]; 735 double value = solution[iPivot]; 736 double cost = costModel[iPivot]; 737 double change = primalRatio * work[iRow]; 738 value = change; 739 changeObj = change * cost; 740 double lower = lowerModel[iPivot]; 741 double upper = upperModel[iPivot]; 742 solution[iPivot] = value; 743 if (value < lower  tolerance) { 744 value = lower; 745 value *= value; 754 746 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 755 756 747 if (iPivot < numberColumns) 748 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 757 749 #endif 758 750 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 759 760 761 #endif 762 763 764 765 766 767 768 769 751 if (lower == upper) 752 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 753 #endif 754 // store square in list 755 if (infeas[iRow]) 756 infeas[iRow] = value; // already there 757 else 758 infeasible_>quickAdd(iRow, value); 759 } else if (value > upper + tolerance) { 760 value = upper; 761 value *= value; 770 762 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 771 772 763 if (iPivot < numberColumns) 764 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 773 765 #endif 774 766 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 775 776 777 #endif 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 767 if (lower == upper) 768 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 769 #endif 770 // store square in list 771 if (infeas[iRow]) 772 infeas[iRow] = value; // already there 773 else 774 infeasible_>quickAdd(iRow, value); 775 } else { 776 // feasible  was it infeasible  if so set tiny 777 if (infeas[iRow]) 778 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 779 } 780 work[iRow] = 0.0; 781 } 782 } 783 // Do pivot row 784 { 785 int iRow = model_>pivotRow(); 786 // feasible  was it infeasible  if so set tiny 787 //assert (infeas[iRow]); 788 if (infeas[iRow]) 789 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 790 } 791 primalUpdate>setNumElements(0); 792 objectiveChange += changeObj; 801 793 } 802 794 /* Saves any weights round factorization as pivot rows may change … … 809 801 7) for strong branching  initialize full weights , infeasibilities 810 802 */ 811 void 812 ClpDualRowSteepest::saveWeights(ClpSimplex * model, int mode) 813 { 814 // alternateWeights_ is defined as indexed but is treated oddly 815 model_ = model; 816 int numberRows = model_>numberRows(); 817 int numberColumns = model_>numberColumns(); 818 const int * pivotVariable = model_>pivotVariable(); 819 int i; 820 if (mode == 1) { 821 if(weights_) { 822 // Check if size has changed 823 if (infeasible_>capacity() == numberRows) { 824 alternateWeights_>clear(); 825 // change from row numbers to sequence numbers 826 int * which = alternateWeights_>getIndices(); 827 for (i = 0; i < numberRows; i++) { 828 int iPivot = pivotVariable[i]; 829 which[i] = iPivot; 830 } 831 state_ = 1; 832 // clear marker on savedWeights_ 833 assert (savedWeights_); 834 if (savedWeights_>packedMode()) { 835 savedWeights_>setPackedMode(false); 836 savedWeights_>setNumElements(0); 837 } 838 } else { 839 // size has changed  clear everything 840 delete [] weights_; 841 weights_ = NULL; 842 delete [] dubiousWeights_; 843 dubiousWeights_ = NULL; 844 delete infeasible_; 845 infeasible_ = NULL; 846 delete alternateWeights_; 847 alternateWeights_ = NULL; 848 delete savedWeights_; 849 savedWeights_ = NULL; 850 state_ = 1; 851 } 803 void ClpDualRowSteepest::saveWeights(ClpSimplex *model, int mode) 804 { 805 // alternateWeights_ is defined as indexed but is treated oddly 806 model_ = model; 807 int numberRows = model_>numberRows(); 808 int numberColumns = model_>numberColumns(); 809 const int *pivotVariable = model_>pivotVariable(); 810 int i; 811 if (mode == 1) { 812 if (weights_) { 813 // Check if size has changed 814 if (infeasible_>capacity() == numberRows) { 815 alternateWeights_>clear(); 816 // change from row numbers to sequence numbers 817 int *which = alternateWeights_>getIndices(); 818 for (i = 0; i < numberRows; i++) { 819 int iPivot = pivotVariable[i]; 820 which[i] = iPivot; 821 } 822 state_ = 1; 823 // clear marker on savedWeights_ 824 assert(savedWeights_); 825 if (savedWeights_>packedMode()) { 826 savedWeights_>setPackedMode(false); 827 savedWeights_>setNumElements(0); 828 } 829 } else { 830 // size has changed  clear everything 831 delete[] weights_; 832 weights_ = NULL; 833 delete[] dubiousWeights_; 834 dubiousWeights_ = NULL; 835 delete infeasible_; 836 infeasible_ = NULL; 837 delete alternateWeights_; 838 alternateWeights_ = NULL; 839 delete savedWeights_; 840 savedWeights_ = NULL; 841 state_ = 1; 842 } 843 } 844 } else if (mode == 2  mode == 4  mode >= 5) { 845 // restore 846 if (!weights_  state_ == 1  mode == 5  mode == 7) { 847 // initialize weights 848 delete[] weights_; 849 delete alternateWeights_; 850 weights_ = new double[numberRows]; 851 // initialize to 1.0 (can we do better?) 852 for (i = 0; i < numberRows; i++) { 853 weights_[i] = 1.0; 854 } 855 alternateWeights_ = new CoinIndexedVector(); 856 // enough space so can use it for factorization 857 alternateWeights_>reserve(numberRows + model_>factorization()>maximumPivots()); 858 if (mode_ != 1  mode == 5) { 859 } else { 860 CoinIndexedVector *temp = new CoinIndexedVector(); 861 temp>reserve(numberRows + model_>factorization()>maximumPivots()); 862 double *array = alternateWeights_>denseVector(); 863 int *which = alternateWeights_>getIndices(); 864 int firstRow = 0; 865 int lastRow = numberRows; 866 if (mode == 7) { 867 // use info passed in 868 firstRow = model>spareIntArray_[0]; 869 lastRow = model>spareIntArray_[1]; 870 } 871 for (i = firstRow; i < lastRow; i++) { 872 double value = 0.0; 873 array[0] = 1.0; 874 which[0] = i; 875 alternateWeights_>setNumElements(1); 876 alternateWeights_>setPackedMode(true); 877 model_>factorization()>updateColumnTranspose(temp, 878 alternateWeights_); 879 int number = alternateWeights_>getNumElements(); 880 int j; 881 for (j = 0; j < number; j++) { 882 value += array[j] * array[j]; 883 array[j] = 0.0; 852 884 } 853 } else if (mode == 2  mode == 4  mode >= 5) { 854 // restore 855 if (!weights_  state_ == 1  mode == 5  mode == 7) { 856 // initialize weights 857 delete [] weights_; 858 delete alternateWeights_; 859 weights_ = new double[numberRows]; 860 // initialize to 1.0 (can we do better?) 861 for (i = 0; i < numberRows; i++) { 862 weights_[i] = 1.0; 863 } 864 alternateWeights_ = new CoinIndexedVector(); 865 // enough space so can use it for factorization 866 alternateWeights_>reserve(numberRows + 867 model_>factorization()>maximumPivots()); 868 if (mode_ != 1  mode == 5) { 869 } else { 870 CoinIndexedVector * temp = new CoinIndexedVector(); 871 temp>reserve(numberRows + 872 model_>factorization()>maximumPivots()); 873 double * array = alternateWeights_>denseVector(); 874 int * which = alternateWeights_>getIndices(); 875 int firstRow=0; 876 int lastRow=numberRows; 877 if (mode==7) { 878 // use info passed in 879 firstRow=model>spareIntArray_[0]; 880 lastRow=model>spareIntArray_[1]; 881 } 882 for (i = firstRow; i < lastRow; i++) { 883 double value = 0.0; 884 array[0] = 1.0; 885 which[0] = i; 886 alternateWeights_>setNumElements(1); 887 alternateWeights_>setPackedMode(true); 888 model_>factorization()>updateColumnTranspose(temp, 889 alternateWeights_); 890 int number = alternateWeights_>getNumElements(); 891 int j; 892 for (j = 0; j < number; j++) { 893 value += array[j] * array[j]; 894 array[j] = 0.0; 895 } 896 alternateWeights_>setNumElements(0); 897 weights_[i] = value; 898 } 899 delete temp; 900 } 901 // create saved weights (not really indexedvector) 902 savedWeights_ = new CoinIndexedVector(); 903 savedWeights_>reserve(numberRows); 904 for (int i=0;i<model_>numberRows();i++) 905 savedWeights_>denseVector()[i]=1.0; 885 alternateWeights_>setNumElements(0); 886 weights_[i] = value; 887 } 888 delete temp; 889 } 890 // create saved weights (not really indexedvector) 891 savedWeights_ = new CoinIndexedVector(); 892 savedWeights_>reserve(numberRows); 893 for (int i = 0; i < model_>numberRows(); i++) 894 savedWeights_>denseVector()[i] = 1.0; 906 895 907 double *array = savedWeights_>denseVector();908 int *which = savedWeights_>getIndices();909 910 911 912 913 if (mode==7) {914 915 916 917 918 int *which = alternateWeights_>getIndices();919 CoinIndexedVector *rowArray3 = model_>rowArray(3);920 921 int *back = rowArray3>getIndices();922 923 924 925 926 927 CoinMemcpyN(which,numberRows, savedWeights_>getIndices());928 CoinMemcpyN(weights_,numberRows, savedWeights_>denseVector());929 930 931 932 933 934 935 936 937 938 double *array = savedWeights_>denseVector();939 940 941 942 943 int firstRow=0;944 int lastRow=numberRows;945 if (mode==7) {946 947 firstRow=model>spareIntArray_[0];948 lastRow=model>spareIntArray_[1];949 950 951 952 953 954 955 956 957 958 959 960 961 962 896 double *array = savedWeights_>denseVector(); 897 int *which = savedWeights_>getIndices(); 898 for (int i = 0; i < numberRows; i++) { 899 array[i] = weights_[i]; 900 which[i] = pivotVariable[i]; 901 } 902 if (mode == 7) { 903 savedWeights_>setNumElements(numberRows); //flag as special 904 savedWeights_>setPackedMode(true); 905 } 906 } else if (mode != 6) { 907 int *which = alternateWeights_>getIndices(); 908 CoinIndexedVector *rowArray3 = model_>rowArray(3); 909 rowArray3>clear(); 910 int *back = rowArray3>getIndices(); 911 // In case something went wrong 912 for (i = 0; i < numberRows + numberColumns; i++) 913 back[i] = 1; 914 if (mode != 4) { 915 // save 916 CoinMemcpyN(which, numberRows, savedWeights_>getIndices()); 917 CoinMemcpyN(weights_, numberRows, savedWeights_>denseVector()); 918 } else { 919 // restore 920 //memcpy(which,savedWeights_>getIndices(), 921 // numberRows*sizeof(int)); 922 //memcpy(weights_,savedWeights_>denseVector(), 923 // numberRows*sizeof(double)); 924 which = savedWeights_>getIndices(); 925 } 926 // restore (a bit slow  but only every refactorization) 927 double *array = savedWeights_>denseVector(); 928 for (i = 0; i < numberRows; i++) { 929 int iSeq = which[i]; 930 back[iSeq] = i; 931 } 932 int firstRow = 0; 933 int lastRow = numberRows; 934 if (mode == 7) { 935 // use info passed in 936 firstRow = model>spareIntArray_[0]; 937 lastRow = model>spareIntArray_[1]; 938 } 939 for (i = firstRow; i < lastRow; i++) { 940 int iPivot = pivotVariable[i]; 941 iPivot = back[iPivot]; 942 if (iPivot >= 0) { 943 weights_[i] = array[iPivot]; 944 if (weights_[i] < DEVEX_TRY_NORM) 945 weights_[i] = DEVEX_TRY_NORM; // may need to check more 946 } else { 947 // odd 948 weights_[i] = 1.0; 949 //printf("bad pivot row %d (%d>%d) iPivot %d\n",i,firstRow,lastRow,iPivot); 950 } 951 } 963 952 #if 0 964 953 printf("mode %d mode_ %d state_ %d\n",mode,mode_,state_); … … 969 958 } 970 959 #endif 971 972 973 974 double allowed;975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 const int *pivotVariable = model_>pivotVariable();1008 1009 1010 1011 1012 1013 1014 1015 1016 960 } else { 961 // mode 6  scale back weights as primal errors 962 double primalError = model_>largestPrimalError(); 963 double allowed; 964 if (primalError > 1.0e3) 965 allowed = 10.0; 966 else if (primalError > 1.0e2) 967 allowed = 50.0; 968 else if (primalError > 1.0e1) 969 allowed = 100.0; 970 else 971 allowed = 1000.0; 972 double allowedInv = 1.0 / allowed; 973 for (i = 0; i < numberRows; i++) { 974 double value = weights_[i]; 975 if (value < allowedInv) 976 value = allowedInv; 977 else if (value > allowed) 978 value = allowed; 979 weights_[i] = allowed; 980 } 981 } 982 state_ = 0; 983 // set up infeasibilities 984 if (!infeasible_) { 985 infeasible_ = new CoinIndexedVector(); 986 infeasible_>reserve(numberRows); 987 } 988 } 989 if (mode >= 2) { 990 // Get dubious weights 991 //if (!dubiousWeights_) 992 //dubiousWeights_=new int[numberRows]; 993 //model_>factorization()>getWeights(dubiousWeights_); 994 infeasible_>clear(); 995 int iRow; 996 const int *pivotVariable = model_>pivotVariable(); 997 double tolerance = model_>currentPrimalTolerance(); 998 for (iRow = 0; iRow < numberRows; iRow++) { 999 int iPivot = pivotVariable[iRow]; 1000 double value = model_>solution(iPivot); 1001 double lower = model_>lower(iPivot); 1002 double upper = model_>upper(iPivot); 1003 if (value < lower  tolerance) { 1004 value = lower; 1005 value *= value; 1017 1006 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 1018 1019 1007 if (iPivot < numberColumns) 1008 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 1020 1009 #endif 1021 1010 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1022 1023 1024 #endif 1025 1026 1027 1028 1029 1011 if (lower == upper) 1012 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1013 #endif 1014 // store square in list 1015 infeasible_>quickAdd(iRow, value); 1016 } else if (value > upper + tolerance) { 1017 value = upper; 1018 value *= value; 1030 1019 #ifdef CLP_DUAL_COLUMN_MULTIPLIER 1031 1032 1020 if (iPivot < numberColumns) 1021 value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns 1033 1022 #endif 1034 1023 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1035 1036 1037 #endif 1038 1039 1040 1041 1042 1043 1044 if (mode==2&&!model>numberIterations()) {1045 int options=model>specialOptions();1046 if ((options&16384)!=0&&true) {1047 1048 1049 if ((options&524288)!=0&&false) {1050 // fathom 1051 for (int i = 0; i < numberRows; i++) 1052 1053 1054 1055 for (int i = 0; i < numberRows; i++) 1056 1057 1058 1059 1024 if (lower == upper) 1025 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1026 #endif 1027 // store square in list 1028 infeasible_>quickAdd(iRow, value); 1029 } 1030 } 1031 } 1032 // see where coming from 1033 if (mode == 2 && !model>numberIterations()) { 1034 int options = model>specialOptions(); 1035 if ((options & 16384) != 0 && true) { 1036 // fast of some sort  could be clever??? 1037 // for now initialize 1038 if ((options & 524288) != 0 && false) { 1039 // fathom 1040 for (int i = 0; i < numberRows; i++) 1041 weights_[i] = 1.0; 1042 } else if (true) { 1043 // strong 1044 for (int i = 0; i < numberRows; i++) 1045 weights_[i] = 1.0; 1046 } 1047 } 1048 } 1060 1049 } 1061 1050 // Pass in saved weights 1062 void 1063 ClpDualRowSteepest::passInSavedWeights(const CoinIndexedVector * saved) 1051 void ClpDualRowSteepest::passInSavedWeights(const CoinIndexedVector *saved) 1064 1052 { 1065 1053 delete savedWeights_; 1066 savedWeights_ =new CoinIndexedVector(*saved);1054 savedWeights_ = new CoinIndexedVector(*saved); 1067 1055 } 1068 1056 // Gets rid of last update 1069 void 1070 ClpDualRowSteepest::unrollWeights() 1071 { 1072 double * saved = alternateWeights_>denseVector(); 1073 int number = alternateWeights_>getNumElements(); 1074 int * which = alternateWeights_>getIndices(); 1075 int i; 1076 if (alternateWeights_>packedMode()) { 1077 for (i = 0; i < number; i++) { 1078 int iRow = which[i]; 1079 weights_[iRow] = saved[i]; 1080 saved[i] = 0.0; 1081 } 1082 } else { 1083 for (i = 0; i < number; i++) { 1084 int iRow = which[i]; 1085 weights_[iRow] = saved[iRow]; 1086 saved[iRow] = 0.0; 1087 } 1088 } 1089 alternateWeights_>setNumElements(0); 1057 void ClpDualRowSteepest::unrollWeights() 1058 { 1059 double *saved = alternateWeights_>denseVector(); 1060 int number = alternateWeights_>getNumElements(); 1061 int *which = alternateWeights_>getIndices(); 1062 int i; 1063 if (alternateWeights_>packedMode()) { 1064 for (i = 0; i < number; i++) { 1065 int iRow = which[i]; 1066 weights_[iRow] = saved[i]; 1067 saved[i] = 0.0; 1068 } 1069 } else { 1070 for (i = 0; i < number; i++) { 1071 int iRow = which[i]; 1072 weights_[iRow] = saved[iRow]; 1073 saved[iRow] = 0.0; 1074 } 1075 } 1076 alternateWeights_>setNumElements(0); 1090 1077 } 1091 1078 // 1092 1079 // Clone 1093 1080 // 1094 ClpDualRowPivot * 1095 { 1096 1097 1098 1099 1100 1081 ClpDualRowPivot *ClpDualRowSteepest::clone(bool CopyData) const 1082 { 1083 if (CopyData) { 1084 return new ClpDualRowSteepest(*this); 1085 } else { 1086 return new ClpDualRowSteepest(); 1087 } 1101 1088 } 1102 1089 // Gets rid of all arrays 1103 void 1104 ClpDualRowSteepest::clearArrays() 1105 { 1106 if (persistence_ == normal) { 1107 delete [] weights_; 1108 weights_ = NULL; 1109 delete [] dubiousWeights_; 1110 dubiousWeights_ = NULL; 1111 delete infeasible_; 1112 infeasible_ = NULL; 1113 delete alternateWeights_; 1114 alternateWeights_ = NULL; 1115 delete savedWeights_; 1116 savedWeights_ = NULL; 1117 } 1118 state_ = 1; 1090 void ClpDualRowSteepest::clearArrays() 1091 { 1092 if (persistence_ == normal) { 1093 delete[] weights_; 1094 weights_ = NULL; 1095 delete[] dubiousWeights_; 1096 dubiousWeights_ = NULL; 1097 delete infeasible_; 1098 infeasible_ = NULL; 1099 delete alternateWeights_; 1100 alternateWeights_ = NULL; 1101 delete savedWeights_; 1102 savedWeights_ = NULL; 1103 } 1104 state_ = 1; 1119 1105 } 1120 1106 // Returns true if would not find any row 1121 bool 1122 ClpDualRowSteepest::looksOptimal() const 1123 { 1124 int iRow; 1125 const int * pivotVariable = model_>pivotVariable(); 1126 double tolerance = model_>currentPrimalTolerance(); 1127 // we can't really trust infeasibilities if there is primal error 1128 // this coding has to mimic coding in checkPrimalSolution 1129 double error = CoinMin(1.0e2, model_>largestPrimalError()); 1130 // allow tolerance at least slightly bigger than standard 1131 tolerance = tolerance + error; 1132 // But cap 1133 tolerance = CoinMin(1000.0, tolerance); 1134 int numberRows = model_>numberRows(); 1135 int numberInfeasible = 0; 1136 for (iRow = 0; iRow < numberRows; iRow++) { 1137 int iPivot = pivotVariable[iRow]; 1138 double value = model_>solution(iPivot); 1139 double lower = model_>lower(iPivot); 1140 double upper = model_>upper(iPivot); 1141 if (value < lower  tolerance) { 1142 numberInfeasible++; 1143 } else if (value > upper + tolerance) { 1144 numberInfeasible++; 1145 } 1146 } 1147 return (numberInfeasible == 0); 1107 bool ClpDualRowSteepest::looksOptimal() const 1108 { 1109 int iRow; 1110 const int *pivotVariable = model_>pivotVariable(); 1111 double tolerance = model_>currentPrimalTolerance(); 1112 // we can't really trust infeasibilities if there is primal error 1113 // this coding has to mimic coding in checkPrimalSolution 1114 double error = CoinMin(1.0e2, model_>largestPrimalError()); 1115 // allow tolerance at least slightly bigger than standard 1116 tolerance = tolerance + error; 1117 // But cap 1118 tolerance = CoinMin(1000.0, tolerance); 1119 int numberRows = model_>numberRows(); 1120 int numberInfeasible = 0; 1121 for (iRow = 0; iRow < numberRows; iRow++) { 1122 int iPivot = pivotVariable[iRow]; 1123 double value = model_>solution(iPivot); 1124 double lower = model_>lower(iPivot); 1125 double upper = model_>upper(iPivot); 1126 if (value < lower  tolerance) { 1127 numberInfeasible++; 1128 } else if (value > upper + tolerance) { 1129 numberInfeasible++; 1130 } 1131 } 1132 return (numberInfeasible == 0); 1148 1133 } 1149 1134 // Called when maximum pivots changes 1150 void 1151 ClpDualRowSteepest::maximumPivotsChanged() 1152 { 1153 if (alternateWeights_ && 1154 alternateWeights_>capacity() != model_>numberRows() + 1155 model_>factorization()>maximumPivots()) { 1156 delete alternateWeights_; 1157 alternateWeights_ = new CoinIndexedVector(); 1158 // enough space so can use it for factorization 1159 alternateWeights_>reserve(model_>numberRows() + 1160 model_>factorization()>maximumPivots()); 1161 } 1135 void ClpDualRowSteepest::maximumPivotsChanged() 1136 { 1137 if (alternateWeights_ && alternateWeights_>capacity() != model_>numberRows() + model_>factorization()>maximumPivots()) { 1138 delete alternateWeights_; 1139 alternateWeights_ = new CoinIndexedVector(); 1140 // enough space so can use it for factorization 1141 alternateWeights_>reserve(model_>numberRows() + model_>factorization()>maximumPivots()); 1142 } 1162 1143 } 1163 1144 1145 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 1146 */
Note: See TracChangeset
for help on using the changeset viewer.