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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/AbcDualRowSteepest.cpp
r1910 r2385 20 20 // Default Constructor 21 21 // 22 AbcDualRowSteepest::AbcDualRowSteepest 23 : AbcDualRowPivot() ,24 norm_(0.0),25 factorizationRatio_(10.0),26 state_(1),27 mode_(mode),28 persistence_(normal),29 weights_(NULL),30 infeasible_(NULL),31 22 AbcDualRowSteepest::AbcDualRowSteepest(int mode) 23 : AbcDualRowPivot() 24 , norm_(0.0) 25 , factorizationRatio_(10.0) 26 , state_(1) 27 , mode_(mode) 28 , persistence_(normal) 29 , weights_(NULL) 30 , infeasible_(NULL) 31 , savedWeights_(NULL) 32 32 { 33 33 type_ = 2 + 64 * mode; … … 38 38 // Copy constructor 39 39 // 40 AbcDualRowSteepest::AbcDualRowSteepest (const AbcDualRowSteepest &rhs)40 AbcDualRowSteepest::AbcDualRowSteepest(const AbcDualRowSteepest &rhs) 41 41 : AbcDualRowPivot(rhs) 42 42 { … … 77 77 // Destructor 78 78 // 79 AbcDualRowSteepest::~AbcDualRowSteepest 79 AbcDualRowSteepest::~AbcDualRowSteepest() 80 80 { 81 81 //printf("XX %x AbcDualRowSteepest destructor\n",this); … … 89 89 // 90 90 AbcDualRowSteepest & 91 AbcDualRowSteepest::operator=(const AbcDualRowSteepest &rhs)91 AbcDualRowSteepest::operator=(const AbcDualRowSteepest &rhs) 92 92 { 93 93 if (this != &rhs) { … … 125 125 } 126 126 // Fill most values 127 void 128 AbcDualRowSteepest::fill(const AbcDualRowSteepest& rhs) 127 void AbcDualRowSteepest::fill(const AbcDualRowSteepest &rhs) 129 128 { 130 129 norm_ = rhs.norm_; … … 133 132 mode_ = rhs.mode_; 134 133 persistence_ = rhs.persistence_; 135 assert 134 assert(model_>numberRows() == rhs.model_>numberRows()); 136 135 model_ = rhs.model_; 137 136 assert(model_); … … 166 165 } 167 166 #ifdef CHECK_NUMBER_WANTED 168 static int xTimes =0;169 static int xWanted =0;170 #endif 171 #if ABC_PARALLEL ==2167 static int xTimes = 0; 168 static int xWanted = 0; 169 #endif 170 #if ABC_PARALLEL == 2 172 171 #define DO_REDUCE 2 173 172 #ifdef DO_REDUCE 174 #if DO_REDUCE ==1173 #if DO_REDUCE == 1 175 174 #include <cilk/reducer_max.h> 176 static void choose(int & chosenRow,double & largest, int n, 177 const int * index, const double * infeas, 178 const double * weights,double tolerance) 179 { 180 cilk::reducer_max_index<int,double> maximumIndex(chosenRow,largest); 181 #pragma cilk_grainsize=128 182 cilk_for (int i = 0; i < n; i++) { 175 static void choose(int &chosenRow, double &largest, int n, 176 const int *index, const double *infeas, 177 const double *weights, double tolerance) 178 { 179 cilk::reducer_max_index< int, double > maximumIndex(chosenRow, largest); 180 #pragma cilk_grainsize = 128 181 cilk_for(int i = 0; i < n; i++) 182 { 183 183 int iRow = index[i]; 184 184 double value = infeas[iRow]; 185 185 if (value > tolerance) { 186 186 double thisWeight = CoinMin(weights[iRow], 1.0e50); 187 maximumIndex.calc_max(iRow, value/thisWeight);188 } 189 } 190 chosenRow =maximumIndex.get_index();191 largest =maximumIndex.get_value();187 maximumIndex.calc_max(iRow, value / thisWeight); 188 } 189 } 190 chosenRow = maximumIndex.get_index(); 191 largest = maximumIndex.get_value(); 192 192 } 193 193 #else 194 static void choose(AbcDualRowSteepest * 195 int & chosenRowSave,double &largestSave, int first, int last,196 197 { 198 if (last first>256) {199 int mid =(last+first)>>1;200 int chosenRow2 =chosenRowSave;201 double largest2 =largestSave;202 cilk_spawn choose(steepest, chosenRow2,largest2, first, mid,203 204 choose(steepest, chosenRowSave,largestSave, mid, last,205 194 static void choose(AbcDualRowSteepest *steepest, 195 int &chosenRowSave, double &largestSave, int first, int last, 196 double tolerance) 197 { 198 if (last  first > 256) { 199 int mid = (last + first) >> 1; 200 int chosenRow2 = chosenRowSave; 201 double largest2 = largestSave; 202 cilk_spawn choose(steepest, chosenRow2, largest2, first, mid, 203 tolerance); 204 choose(steepest, chosenRowSave, largestSave, mid, last, 205 tolerance); 206 206 cilk_sync; 207 if (largest2 >largestSave) {208 largestSave =largest2;209 chosenRowSave =chosenRow2;207 if (largest2 > largestSave) { 208 largestSave = largest2; 209 chosenRowSave = chosenRow2; 210 210 } 211 211 } else { 212 const int * index=steepest>infeasible()>getIndices();213 const double * infeas=steepest>infeasible()>denseVector();214 const double * weights=steepest>weights()>denseVector();215 double largest =largestSave;216 int chosenRow =chosenRowSave;217 if ((steepest>model()>stateOfProblem() &VALUES_PASS2)==0) {212 const int *index = steepest>infeasible()>getIndices(); 213 const double *infeas = steepest>infeasible()>denseVector(); 214 const double *weights = steepest>weights()>denseVector(); 215 double largest = largestSave; 216 int chosenRow = chosenRowSave; 217 if ((steepest>model()>stateOfProblem() & VALUES_PASS2) == 0) { 218 218 for (int i = first; i < last; i++) { 219 220 221 222 223 if (value>largest*thisWeight) {224 largest=value/thisWeight;225 chosenRow=iRow;226 227 219 int iRow = index[i]; 220 double value = infeas[iRow]; 221 if (value > tolerance) { 222 double thisWeight = CoinMin(weights[iRow], 1.0e50); 223 if (value > largest * thisWeight) { 224 largest = value / thisWeight; 225 chosenRow = iRow; 226 } 227 } 228 228 } 229 229 } else { 230 const double * 231 const int * 230 const double *fakeDjs = steepest>model()>fakeDjs(); 231 const int *pivotVariable = steepest>model()>pivotVariable(); 232 232 for (int i = first; i < last; i++) { 233 234 235 236 int iSequence=pivotVariable[iRow];237 value *=(fabs(fakeDjs[iSequence])+1.0e6);238 239 233 int iRow = index[i]; 234 double value = infeas[iRow]; 235 if (value > tolerance) { 236 int iSequence = pivotVariable[iRow]; 237 value *= (fabs(fakeDjs[iSequence]) + 1.0e6); 238 double thisWeight = CoinMin(weights[iRow], 1.0e50); 239 /* 240 240 Ideas 241 241 always use fake … … 244 244 how to switch off if cilking  easiest at saveWeights 245 245 */ 246 if (value>largest*thisWeight) {247 largest=value/thisWeight;248 chosenRow=iRow;249 250 251 } 252 } 253 chosenRowSave =chosenRow;254 largestSave =largest;255 } 256 } 257 static void choose2(AbcDualRowSteepest * 258 int & chosenRowSave,double &largestSave, int first, int last,259 260 { 261 if (last first>256) {262 int mid =(last+first)>>1;263 int chosenRow2 =chosenRowSave;264 double largest2 =largestSave;265 cilk_spawn choose2(steepest, chosenRow2,largest2, first, mid,266 267 choose2(steepest, chosenRowSave,largestSave, mid, last,268 246 if (value > largest * thisWeight) { 247 largest = value / thisWeight; 248 chosenRow = iRow; 249 } 250 } 251 } 252 } 253 chosenRowSave = chosenRow; 254 largestSave = largest; 255 } 256 } 257 static void choose2(AbcDualRowSteepest *steepest, 258 int &chosenRowSave, double &largestSave, int first, int last, 259 double tolerance) 260 { 261 if (last  first > 256) { 262 int mid = (last + first) >> 1; 263 int chosenRow2 = chosenRowSave; 264 double largest2 = largestSave; 265 cilk_spawn choose2(steepest, chosenRow2, largest2, first, mid, 266 tolerance); 267 choose2(steepest, chosenRowSave, largestSave, mid, last, 268 tolerance); 269 269 cilk_sync; 270 if (largest2 >largestSave) {271 largestSave =largest2;272 chosenRowSave =chosenRow2;270 if (largest2 > largestSave) { 271 largestSave = largest2; 272 chosenRowSave = chosenRow2; 273 273 } 274 274 } else { 275 const int * index=steepest>infeasible()>getIndices();276 const double * infeas=steepest>infeasible()>denseVector();277 double largest =largestSave;278 int chosenRow =chosenRowSave;275 const int *index = steepest>infeasible()>getIndices(); 276 const double *infeas = steepest>infeasible()>denseVector(); 277 double largest = largestSave; 278 int chosenRow = chosenRowSave; 279 279 for (int i = first; i < last; i++) { 280 280 int iRow = index[i]; 281 281 double value = infeas[iRow]; 282 282 if (value > largest) { 283 largest=value;284 chosenRow=iRow;285 } 286 } 287 chosenRowSave =chosenRow;288 largestSave =largest;283 largest = value; 284 chosenRow = iRow; 285 } 286 } 287 chosenRowSave = chosenRow; 288 largestSave = largest; 289 289 } 290 290 } … … 293 293 #endif 294 294 // Returns pivot row, 1 if none 295 int 296 AbcDualRowSteepest::pivotRow() 295 int AbcDualRowSteepest::pivotRow() 297 296 { 298 297 assert(model_); 299 double * 300 int * 298 double *COIN_RESTRICT infeas = infeasible_>denseVector(); 299 int *COIN_RESTRICT index = infeasible_>getIndices(); 301 300 int number = infeasible_>getNumElements(); 302 301 int lastPivotRow = model_>lastPivotRow(); 303 assert 302 assert(lastPivotRow < model_>numberRows()); 304 303 double tolerance = model_>currentPrimalTolerance(); 305 304 // we can't really trust infeasibilities if there is primal error … … 307 306 double error = CoinMin(1.0e2, model_>largestPrimalError()); 308 307 // allow tolerance at least slightly bigger than standard 309 tolerance = tolerance + 308 tolerance = tolerance + error; 310 309 // But cap 311 310 tolerance = CoinMin(1000.0, tolerance); 312 311 tolerance *= tolerance; // as we are using squares 313 312 bool toleranceChanged = false; 314 const double * 315 const double * 316 const double * 313 const double *COIN_RESTRICT solutionBasic = model_>solutionBasic(); 314 const double *COIN_RESTRICT lowerBasic = model_>lowerBasic(); 315 const double *COIN_RESTRICT upperBasic = model_>upperBasic(); 317 316 // do last pivot row one here 318 317 if (lastPivotRow >= 0 && lastPivotRow < model_>numberRows()) { … … 325 324 // store square in list 326 325 if (infeas[lastPivotRow]) 327 326 infeas[lastPivotRow] = value; // already there 328 327 else 329 328 infeasible_>quickAdd(lastPivotRow, value); 330 329 } else if (value < lower  tolerance) { 331 330 value = lower; … … 333 332 // store square in list 334 333 if (infeas[lastPivotRow]) 335 334 infeas[lastPivotRow] = value; // already there 336 335 else 337 336 infeasible_>add(lastPivotRow, value); 338 337 } else { 339 338 // feasible  was it infeasible  if so set tiny 340 339 if (infeas[lastPivotRow]) 341 340 infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 342 341 } 343 342 number = infeasible_>getNumElements(); 344 343 } 345 if (model_>numberIterations() < model_>lastBadIteration() + 200) {344 if (model_>numberIterations() < model_>lastBadIteration() + 200) { 346 345 // we can't really trust infeasibilities if there is dual error 347 346 if (model_>largestDualError() > model_>largestPrimalError()) { … … 351 350 } 352 351 int numberWanted; 353 if (mode_ < 2 352 if (mode_ < 2) { 354 353 numberWanted = number + 1; 355 354 } else { 356 if (model_>factorization()>pivots() ==0) {355 if (model_>factorization()>pivots() == 0) { 357 356 int numberElements = model_>factorization()>numberElements(); 358 357 int numberRows = model_>numberRows(); … … 362 361 static_cast<double> (numberRows+1numberSlacks); 363 362 #else 364 factorizationRatio_ = static_cast<double> (numberElements) / 365 static_cast<double> (numberRows); 363 factorizationRatio_ = static_cast< double >(numberElements) / static_cast< double >(numberRows); 366 364 #endif 367 365 //printf("%d iterations, numberElements %d ratio %g\n", 368 366 // model_>numberIterations(),numberElements,factorizationRatio_); 369 367 // Bias so less likely to go to steepest if large 370 double weights[] ={1.0,2.0,5.0};368 double weights[] = { 1.0, 2.0, 5.0 }; 371 369 double modifier; 372 if (numberRows <100000)373 modifier=weights[0];374 else if (numberRows <200000)375 modifier=weights[1];370 if (numberRows < 100000) 371 modifier = weights[0]; 372 else if (numberRows < 200000) 373 modifier = weights[1]; 376 374 else 377 modifier=weights[2];378 if (mode_ ==2&&factorizationRatio_>modifier) {379 380 mode_=3;381 saveWeights(model_,5);375 modifier = weights[2]; 376 if (mode_ == 2 && factorizationRatio_ > modifier) { 377 // switch from dantzig to steepest 378 mode_ = 3; 379 saveWeights(model_, 5); 382 380 } 383 381 } … … 407 405 } 408 406 #endif 409 numberWanted =CoinMin(numberWanted,number+1);407 numberWanted = CoinMin(numberWanted, number + 1); 410 408 } 411 409 if (model_>largestPrimalError() > 1.0e3) 412 410 numberWanted = number + 1; // be safe 413 411 #ifdef CHECK_NUMBER_WANTED 414 xWanted +=numberWanted;412 xWanted += numberWanted; 415 413 xTimes++; 416 if ((xTimes %1000)==0)417 printf("time %d average wanted %g\n", xTimes,static_cast<double>(xWanted)/xTimes);414 if ((xTimes % 1000) == 0) 415 printf("time %d average wanted %g\n", xTimes, static_cast< double >(xWanted) / xTimes); 418 416 #endif 419 417 int iPass; … … 422 420 start[1] = number; 423 421 start[2] = 0; 424 double dstart = static_cast< double>(number) * model_>randomNumberGenerator()>randomDouble();425 start[0] = static_cast< int>(dstart);422 double dstart = static_cast< double >(number) * model_>randomNumberGenerator()>randomDouble(); 423 start[0] = static_cast< int >(dstart); 426 424 start[3] = start[0]; 427 425 //double largestWeight=0.0; 428 426 //double smallestWeight=1.0e100; 429 const double * 430 const int * 431 double savePivotWeight =weights_>denseVector()[lastPivotRow];427 const double *weights = weights_>denseVector(); 428 const int *pivotVariable = model_>pivotVariable(); 429 double savePivotWeight = weights_>denseVector()[lastPivotRow]; 432 430 weights_>denseVector()[lastPivotRow] *= 1.0e10; 433 431 double largest = 0.0; 434 432 int chosenRow = 1; 435 int saveNumberWanted =numberWanted;433 int saveNumberWanted = numberWanted; 436 434 #ifdef DO_REDUCE 437 bool doReduce =true;438 int lastChosen =1;439 double lastLargest =0.0;435 bool doReduce = true; 436 int lastChosen = 1; 437 double lastLargest = 0.0; 440 438 #endif 441 439 for (iPass = 0; iPass < 2; iPass++) { 442 int endThis = start[2 *iPass+1];443 int startThis =start[2*iPass];444 while (startThis <endThis) {445 int end = CoinMin(endThis, startThis+numberWanted);446 #ifdef DO_REDUCE 440 int endThis = start[2 * iPass + 1]; 441 int startThis = start[2 * iPass]; 442 while (startThis < endThis) { 443 int end = CoinMin(endThis, startThis + numberWanted); 444 #ifdef DO_REDUCE 447 445 if (doReduce) { 448 #if DO_REDUCE ==1449 choose(chosenRow,largest,endstartThis,450 index+startThis,infeas,weights,tolerance);446 #if DO_REDUCE == 1 447 choose(chosenRow, largest, end  startThis, 448 index + startThis, infeas, weights, tolerance); 451 449 #else 452 if (mode_!=2) 453 choose(this,chosenRow,largest,startThis,end,tolerance); 454 else 455 choose2(this,chosenRow,largest,startThis,end,tolerance); 456 #endif 457 if (chosenRow!=lastChosen) { 458 assert (chosenRow>=0); 459 if (model_>flagged(pivotVariable[chosenRow]) 460 (solutionBasic[chosenRow] <= upperBasic[chosenRow] + tolerance && 461 solutionBasic[chosenRow] >= lowerBasic[chosenRow]  tolerance)) { 462 doReduce=false; 463 chosenRow=lastChosen; 464 largest=lastLargest; 465 } else { 466 lastChosen=chosenRow; 467 lastLargest=largest; 468 } 469 } 450 if (mode_ != 2) 451 choose(this, chosenRow, largest, startThis, end, tolerance); 452 else 453 choose2(this, chosenRow, largest, startThis, end, tolerance); 454 #endif 455 if (chosenRow != lastChosen) { 456 assert(chosenRow >= 0); 457 if (model_>flagged(pivotVariable[chosenRow])  (solutionBasic[chosenRow] <= upperBasic[chosenRow] + tolerance && solutionBasic[chosenRow] >= lowerBasic[chosenRow]  tolerance)) { 458 doReduce = false; 459 chosenRow = lastChosen; 460 largest = lastLargest; 461 } else { 462 lastChosen = chosenRow; 463 lastLargest = largest; 464 } 465 } 470 466 } 471 467 if (!doReduce) { 472 468 #endif 473 for (int i = startThis; i < end; i++) {474 475 476 477 469 for (int i = startThis; i < end; i++) { 470 int iRow = index[i]; 471 double value = infeas[iRow]; 472 if (value > tolerance) { 473 //#define OUT_EQ 478 474 #ifdef OUT_EQ 479 480 481 #endif 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 475 if (upper[iRow] == lower[iRow]) 476 value *= 2.0; 477 #endif 478 double thisWeight = CoinMin(weights[iRow], 1.0e50); 479 //largestWeight = CoinMax(largestWeight,weight); 480 //smallestWeight = CoinMin(smallestWeight,weight); 481 //double dubious = dubiousWeights_[iRow]; 482 //weight *= dubious; 483 //if (value>2.0*largest*weight(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) { 484 if (value > largest * thisWeight) { 485 // make last pivot row last resort choice 486 //if (iRow == lastPivotRow) { 487 //if (value * 1.0e10 < largest * thisWeight) 488 // continue; 489 //else 490 // value *= 1.0e10; 491 //} 492 if (!model_>flagged(pivotVariable[iRow])) { 493 //#define CLP_DEBUG 3 498 494 #ifdef CLP_DEBUG 499 double value2 = 0.0; 500 if (solutionBasic[iRow] > upperBasic[iRow] + tolerance) 501 value2 = solutionBasic[iRow]  upperBasic[iRow]; 502 else if (solutionBasic[iRow] < lowerBasic[iRow]  tolerance) 503 value2 = solutionBasic[iRow]  lowerBasic[iRow]; 504 assert(fabs(value2 * value2  infeas[iRow]) < 1.0e8 * CoinMin(value2 * value2, infeas[iRow])); 505 #endif 506 if (solutionBasic[iRow] > upperBasic[iRow] + tolerance  507 solutionBasic[iRow] < lowerBasic[iRow]  tolerance) { 508 chosenRow = iRow; 509 largest = value / thisWeight; 510 //largestWeight = dubious; 511 } 512 } 513 } 514 } 515 } 495 double value2 = 0.0; 496 if (solutionBasic[iRow] > upperBasic[iRow] + tolerance) 497 value2 = solutionBasic[iRow]  upperBasic[iRow]; 498 else if (solutionBasic[iRow] < lowerBasic[iRow]  tolerance) 499 value2 = solutionBasic[iRow]  lowerBasic[iRow]; 500 assert(fabs(value2 * value2  infeas[iRow]) < 1.0e8 * CoinMin(value2 * value2, infeas[iRow])); 501 #endif 502 if (solutionBasic[iRow] > upperBasic[iRow] + tolerance  solutionBasic[iRow] < lowerBasic[iRow]  tolerance) { 503 chosenRow = iRow; 504 largest = value / thisWeight; 505 //largestWeight = dubious; 506 } 507 } 508 } 509 } 510 } 516 511 #ifdef DO_REDUCE 517 512 } 518 513 #endif 519 numberWanted =(endstartThis);514 numberWanted = (end  startThis); 520 515 if (!numberWanted) { 521 if(chosenRow>=0)522 523 524 numberWanted=(saveNumberWanted+1)>>1;525 } 526 startThis =end;516 if (chosenRow >= 0) 517 break; 518 else 519 numberWanted = (saveNumberWanted + 1) >> 1; 520 } 521 startThis = end; 527 522 } 528 523 if (!numberWanted) { 529 if (chosenRow>=0)530 524 if (chosenRow >= 0) 525 break; 531 526 else 532 numberWanted=(saveNumberWanted+1)>>1;533 } 534 } 535 weights_>denseVector()[lastPivotRow] =savePivotWeight;527 numberWanted = (saveNumberWanted + 1) >> 1; 528 } 529 } 530 weights_>denseVector()[lastPivotRow] = savePivotWeight; 536 531 //printf("smallest %g largest %g\n",smallestWeight,largestWeight); 537 532 if (chosenRow < 0 && toleranceChanged) { … … 548 543 int iRow = index[i]; 549 544 if (fabs(infeas[iRow]) > 1.0e50) { 550 545 index[nLeft++] = iRow; 551 546 } else { 552 547 infeas[iRow] = 0.0; 553 548 } 554 549 } … … 556 551 model_>setNumberPrimalInfeasibilities(nLeft); 557 552 } 558 if (true &&(model_>stateOfProblem()&VALUES_PASS2)!=0) {559 if (chosenRow >=0) {560 double * 553 if (true && (model_>stateOfProblem() & VALUES_PASS2) != 0) { 554 if (chosenRow >= 0) { 555 double *fakeDjs = model_>fakeDjs(); 561 556 //const int * pivotVariable = model_>pivotVariable(); 562 fakeDjs[pivotVariable[chosenRow]] =0.0;557 fakeDjs[pivotVariable[chosenRow]] = 0.0; 563 558 } 564 559 } … … 567 562 //#define PRINT_WEIGHTS_FREQUENCY 568 563 #ifdef PRINT_WEIGHTS_FREQUENCY 569 int xweights1 =0;570 int xweights2 =0;571 int xweights3 =0;564 int xweights1 = 0; 565 int xweights2 = 0; 566 int xweights3 = 0; 572 567 #endif 573 568 /* Updates weights and returns pivot alpha. 574 569 Also does FT update */ 575 570 double 576 AbcDualRowSteepest::updateWeights(CoinIndexedVector & 577 CoinIndexedVector &updateColumn)578 { 579 double *COIN_RESTRICT weights = weights_>denseVector();580 #if CLP_DEBUG >2571 AbcDualRowSteepest::updateWeights(CoinIndexedVector &input, 572 CoinIndexedVector &updateColumn) 573 { 574 double *COIN_RESTRICT weights = weights_>denseVector(); 575 #if CLP_DEBUG > 2 581 576 // Very expensive debug 582 577 { … … 584 579 CoinIndexedVector temp; 585 580 temp.reserve(numberRows); 586 double * 587 int * 581 double *array = temp.denseVector(); 582 int *which = temp.getIndices(); 588 583 for (int i = 0; i < numberRows; i++) { 589 584 double value = 0.0; … … 594 589 int number = temp.getNumElements(); 595 590 for (int j = 0; j < number; j++) { 596 597 598 591 int iRow = which[j]; 592 value += array[iRow] * array[iRow]; 593 array[iRow] = 0.0; 599 594 } 600 595 temp.setNumElements(0); 601 596 double w = CoinMax(weights[i], value) * .1; 602 597 if (fabs(weights[i]  value) > w) { 603 604 598 printf("%d old %g, true %g\n", i, weights[i], value); 599 weights[i] = value; // to reduce printout 605 600 } 606 601 } … … 610 605 double norm = 0.0; 611 606 int i; 612 double * 607 double *COIN_RESTRICT work2 = input.denseVector(); 613 608 int numberNonZero = input.getNumElements(); 614 int * 615 609 int *COIN_RESTRICT which = input.getIndices(); 610 616 611 //compute norm 617 for (int i = 0; i < numberNonZero; i ++) {612 for (int i = 0; i < numberNonZero; i++) { 618 613 int iRow = which[i]; 619 614 double value = work2[iRow]; … … 622 617 // Do FT update 623 618 if (true) { 624 model_>factorization()>updateTwoColumnsFT(updateColumn, input);625 #if CLP_DEBUG >3619 model_>factorization()>updateTwoColumnsFT(updateColumn, input); 620 #if CLP_DEBUG > 3 626 621 printf("REGION after %d els\n", input.getNumElements()); 627 622 input.print(); … … 635 630 numberNonZero = input.getNumElements(); 636 631 int pivotRow = model_>lastPivotRow(); 637 assert (model_>pivotRow()==model_>lastPivotRow());632 assert(model_>pivotRow() == model_>lastPivotRow()); 638 633 #ifdef CLP_DEBUG 639 if ( model_>logLevel ( ) > 4 && 640 fabs(norm  weights[pivotRow]) > 1.0e2 * (1.0 + norm)) 634 if (model_>logLevel() > 4 && fabs(norm  weights[pivotRow]) > 1.0e2 * (1.0 + norm)) 641 635 printf("on row %d, true weight %g, old %g\n", 642 636 pivotRow, sqrt(norm), sqrt(weights[pivotRow])); 643 637 #endif 644 638 // could reinitialize here (could be expensive) … … 649 643 #ifdef PRINT_WEIGHTS_FREQUENCY 650 644 xweights1++; 651 if ((xweights1 +xweights2+xweights3)%100==0) {645 if ((xweights1 + xweights2 + xweights3) % 100 == 0) { 652 646 printf("ZZ iteration %d sum weights %d  %d %d %d\n", 653 model_>numberIterations(),xweights1+xweights2+xweights3,xweights1,xweights2,xweights3);654 assert(abs(model_>numberIterations() (xweights1+xweights2+xweights3))<=1);647 model_>numberIterations(), xweights1 + xweights2 + xweights3, xweights1, xweights2, xweights3); 648 assert(abs(model_>numberIterations()  (xweights1 + xweights2 + xweights3)) <= 1); 655 649 } 656 650 #endif 657 651 // look at updated column 658 double * 652 double *COIN_RESTRICT work = updateColumn.denseVector(); 659 653 numberNonZero = updateColumn.getNumElements(); 660 654 which = updateColumn.getIndices(); 661 655 662 656 // pivot element 663 alpha =work[pivotRow];657 alpha = work[pivotRow]; 664 658 for (i = 0; i < numberNonZero; i++) { 665 659 int iRow = which[i]; … … 669 663 //which3[nSave++] = iRow; 670 664 double value = work2[iRow]; 671 devex += 665 devex += theta * (theta * norm + value * multiplier); 672 666 if (devex < DEVEX_TRY_NORM) 673 667 devex = DEVEX_TRY_NORM; … … 688 682 return alpha; 689 683 } 690 double 691 AbcDualRowSteepest::updateWeights1(CoinIndexedVector & input,CoinIndexedVector &updateColumn)692 { 693 if (mode_ ==2) {684 double 685 AbcDualRowSteepest::updateWeights1(CoinIndexedVector &input, CoinIndexedVector &updateColumn) 686 { 687 if (mode_ == 2) { 694 688 abort(); 695 689 // no update … … 699 693 double alpha = 0.0; 700 694 // look at updated column 701 double * 695 double *work = updateColumn.denseVector(); 702 696 int pivotRow = model_>lastPivotRow(); 703 assert 704 705 assert 697 assert(pivotRow == model_>pivotRow()); 698 699 assert(!updateColumn.packedMode()); 706 700 alpha = work[pivotRow]; 707 701 return alpha; 708 702 } 709 #if CLP_DEBUG >2703 #if CLP_DEBUG > 2 710 704 // Very expensive debug 711 705 { 712 double * 706 double *weights = weights_>denseVector(); 713 707 int numberRows = model_>numberRows(); 714 const int * 708 const int *pivotVariable = model_>pivotVariable(); 715 709 CoinIndexedVector temp; 716 710 temp.reserve(numberRows); 717 double * 718 int * 711 double *array = temp.denseVector(); 712 int *which = temp.getIndices(); 719 713 for (int i = 0; i < numberRows; i++) { 720 714 double value = 0.0; … … 725 719 int number = temp.getNumElements(); 726 720 for (int j = 0; j < number; j++) { 727 728 729 721 int iRow = which[j]; 722 value += array[iRow] * array[iRow]; 723 array[iRow] = 0.0; 730 724 } 731 725 temp.setNumElements(0); 732 726 double w = CoinMax(weights[i], value) * .1; 733 727 if (fabs(weights[i]  value) > w) { 734 735 736 728 printf("XX row %d (variable %d) old %g, true %g\n", i, pivotVariable[i], 729 weights[i], value); 730 weights[i] = value; // to reduce printout 737 731 } 738 732 } … … 740 734 #endif 741 735 norm_ = 0.0; 742 double * 736 double *COIN_RESTRICT work2 = input.denseVector(); 743 737 int numberNonZero = input.getNumElements(); 744 int * 745 738 int *COIN_RESTRICT which = input.getIndices(); 739 746 740 //compute norm 747 for (int i = 0; i < numberNonZero; i ++) {741 for (int i = 0; i < numberNonZero; i++) { 748 742 int iRow = which[i]; 749 743 double value = work2[iRow]; … … 752 746 // Do FT update 753 747 if (true) { 754 model_>factorization()>updateTwoColumnsFT(updateColumn, input);755 #if CLP_DEBUG >3748 model_>factorization()>updateTwoColumnsFT(updateColumn, input); 749 #if CLP_DEBUG > 3 756 750 printf("REGION after %d els\n", input.getNumElements()); 757 751 input.print(); … … 763 757 } 764 758 // pivot element 765 assert (model_>pivotRow()==model_>lastPivotRow());759 assert(model_>pivotRow() == model_>lastPivotRow()); 766 760 return updateColumn.denseVector()[model_>lastPivotRow()]; 767 761 } 768 void AbcDualRowSteepest::updateWeightsOnly(CoinIndexedVector & 769 { 770 if (mode_ ==2)762 void AbcDualRowSteepest::updateWeightsOnly(CoinIndexedVector &input) 763 { 764 if (mode_ == 2) 771 765 return; 772 766 norm_ = 0.0; 773 double * 767 double *COIN_RESTRICT work2 = input.denseVector(); 774 768 int numberNonZero = input.getNumElements(); 775 int * 776 769 int *COIN_RESTRICT which = input.getIndices(); 770 777 771 //compute norm 778 for (int i = 0; i < numberNonZero; i ++) {772 for (int i = 0; i < numberNonZero; i++) { 779 773 int iRow = which[i]; 780 774 double value = work2[iRow]; … … 784 778 } 785 779 // Actually updates weights 786 void AbcDualRowSteepest::updateWeights2(CoinIndexedVector & input,CoinIndexedVector &updateColumn)787 { 788 if (mode_ ==2)780 void AbcDualRowSteepest::updateWeights2(CoinIndexedVector &input, CoinIndexedVector &updateColumn) 781 { 782 if (mode_ == 2) 789 783 return; 790 double * 791 const double * 784 double *COIN_RESTRICT weights = weights_>denseVector(); 785 const double *COIN_RESTRICT work2 = input.denseVector(); 792 786 int pivotRow = model_>lastPivotRow(); 793 assert (model_>pivotRow()==model_>lastPivotRow());787 assert(model_>pivotRow() == model_>lastPivotRow()); 794 788 #ifdef CLP_DEBUG 795 if ( model_>logLevel ( ) > 4 && 796 fabs(norm_  weights[pivotRow]) > 1.0e2 * (1.0 + norm_)) 789 if (model_>logLevel() > 4 && fabs(norm_  weights[pivotRow]) > 1.0e2 * (1.0 + norm_)) 797 790 printf("on row %d, true weight %g, old %g\n", 798 791 pivotRow, sqrt(norm_), sqrt(weights[pivotRow])); 799 792 #endif 800 793 // pivot element 801 double alpha =model_>alpha();794 double alpha = model_>alpha(); 802 795 // could reinitialize here (could be expensive) 803 796 norm_ /= alpha * alpha; … … 806 799 double multiplier = 2.0 / alpha; 807 800 // look at updated column 808 const double * 801 const double *COIN_RESTRICT work = updateColumn.denseVector(); 809 802 int numberNonZero = updateColumn.getNumElements(); 810 const int * 811 double multiplier2 =1.0;803 const int *which = updateColumn.getIndices(); 804 double multiplier2 = 1.0; 812 805 #ifdef PRINT_WEIGHTS_FREQUENCY 813 806 xweights2++; 814 if ((xweights1 +xweights2+xweights3)%100==0) {807 if ((xweights1 + xweights2 + xweights3) % 100 == 0) { 815 808 printf("ZZ iteration %d sum weights %d  %d %d %d\n", 816 model_>numberIterations(),xweights1+xweights2+xweights3,xweights1,xweights2,xweights3);817 assert(abs(model_>numberIterations() (xweights1+xweights2+xweights3))<=1);818 } 819 #endif 820 if (model_>directionOut() <0) {809 model_>numberIterations(), xweights1 + xweights2 + xweights3, xweights1, xweights2, xweights3); 810 assert(abs(model_>numberIterations()  (xweights1 + xweights2 + xweights3)) <= 1); 811 } 812 #endif 813 if (model_>directionOut() < 0) { 821 814 #ifdef ABC_DEBUG 822 printf("Minus out in %d %d, out %d %d\n", model_>sequenceIn(),model_>directionIn(),823 model_>sequenceOut(),model_>directionOut());824 #endif 825 multiplier2 =1.0;815 printf("Minus out in %d %d, out %d %d\n", model_>sequenceIn(), model_>directionIn(), 816 model_>sequenceOut(), model_>directionOut()); 817 #endif 818 multiplier2 = 1.0; 826 819 } 827 820 for (int i = 0; i < numberNonZero; i++) { 828 821 int iRow = which[i]; 829 double theta = multiplier2 *work[iRow];822 double theta = multiplier2 * work[iRow]; 830 823 double devex = weights[iRow]; 831 824 double value = work2[iRow]; 832 devex += 825 devex += theta * (theta * norm_ + value * multiplier); 833 826 if (devex < DEVEX_TRY_NORM) 834 827 devex = DEVEX_TRY_NORM; … … 850 843 Computes change in objective function 851 844 */ 852 void 853 AbcDualRowSteepest::updatePrimalSolution( 854 CoinIndexedVector & primalUpdate, 855 double primalRatio) 856 { 857 double * COIN_RESTRICT work = primalUpdate.denseVector(); 845 void AbcDualRowSteepest::updatePrimalSolution( 846 CoinIndexedVector &primalUpdate, 847 double primalRatio) 848 { 849 double *COIN_RESTRICT work = primalUpdate.denseVector(); 858 850 int number = primalUpdate.getNumElements(); 859 int * 851 int *COIN_RESTRICT which = primalUpdate.getIndices(); 860 852 double tolerance = model_>currentPrimalTolerance(); 861 double * 862 double * 863 const double * 864 const double * 865 assert 853 double *COIN_RESTRICT infeas = infeasible_>denseVector(); 854 double *COIN_RESTRICT solutionBasic = model_>solutionBasic(); 855 const double *COIN_RESTRICT lowerBasic = model_>lowerBasic(); 856 const double *COIN_RESTRICT upperBasic = model_>upperBasic(); 857 assert(!primalUpdate.packedMode()); 866 858 for (int i = 0; i < number; i++) { 867 859 int iRow = which[i]; … … 879 871 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 880 872 if (lower == upper) 881 873 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 882 874 #endif 883 875 // store square in list 884 876 if (infeas[iRow]) 885 877 infeas[iRow] = value; // already there 886 878 else 887 879 infeasible_>quickAdd(iRow, value); 888 880 } else if (value > upper + tolerance) { 889 881 value = upper; … … 891 883 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 892 884 if (lower == upper) 893 885 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 894 886 #endif 895 887 // store square in list 896 888 if (infeas[iRow]) 897 889 infeas[iRow] = value; // already there 898 890 else 899 891 infeasible_>quickAdd(iRow, value); 900 892 } else { 901 893 // feasible  was it infeasible  if so set tiny 902 894 if (infeas[iRow]) 903 895 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 904 896 } 905 897 } 906 898 // Do pivot row 907 899 int iRow = model_>lastPivotRow(); 908 assert (model_>pivotRow()==model_>lastPivotRow());900 assert(model_>pivotRow() == model_>lastPivotRow()); 909 901 // feasible  was it infeasible  if so set tiny 910 902 //assert (infeas[iRow]); … … 913 905 primalUpdate.setNumElements(0); 914 906 } 915 #if ABC_PARALLEL ==2907 #if ABC_PARALLEL == 2 916 908 static void update(int first, int last, 917 const int * COIN_RESTRICT which, double * COIN_RESTRICT work, 918 const double * COIN_RESTRICT work2, double *COIN_RESTRICT weights,919 const double * COIN_RESTRICT lowerBasic,double *COIN_RESTRICT solutionBasic,920 const double *COIN_RESTRICT upperBasic,921 double multiplier,double multiplier2,922 double norm,double theta,double tolerance)923 { 924 if (last first>256) {925 int mid =(last+first)>>1;926 cilk_spawn update(first, mid,which,work,work2,weights,lowerBasic,solutionBasic,927 upperBasic,multiplier,multiplier2,norm,theta,tolerance);928 update(mid, last,which,work,work2,weights,lowerBasic,solutionBasic,929 upperBasic,multiplier,multiplier2,norm,theta,tolerance);909 const int *COIN_RESTRICT which, double *COIN_RESTRICT work, 910 const double *COIN_RESTRICT work2, double *COIN_RESTRICT weights, 911 const double *COIN_RESTRICT lowerBasic, double *COIN_RESTRICT solutionBasic, 912 const double *COIN_RESTRICT upperBasic, 913 double multiplier, double multiplier2, 914 double norm, double theta, double tolerance) 915 { 916 if (last  first > 256) { 917 int mid = (last + first) >> 1; 918 cilk_spawn update(first, mid, which, work, work2, weights, lowerBasic, solutionBasic, 919 upperBasic, multiplier, multiplier2, norm, theta, tolerance); 920 update(mid, last, which, work, work2, weights, lowerBasic, solutionBasic, 921 upperBasic, multiplier, multiplier2, norm, theta, tolerance); 930 922 cilk_sync; 931 923 } else { … … 933 925 int iRow = which[i]; 934 926 double updateValue = work[iRow]; 935 double thetaDevex = multiplier2 *updateValue;927 double thetaDevex = multiplier2 * updateValue; 936 928 double devex = weights[iRow]; 937 929 double valueDevex = work2[iRow]; 938 devex += 930 devex += thetaDevex * (thetaDevex * norm + valueDevex * multiplier); 939 931 if (devex < DEVEX_TRY_NORM) 940 932 devex = DEVEX_TRY_NORM; 941 933 weights[iRow] = devex; 942 934 double value = solutionBasic[iRow]; … … 947 939 solutionBasic[iRow] = value; 948 940 if (value < lower  tolerance) { 949 950 941 value = lower; 942 value *= value; 951 943 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 952 953 944 if (lower == upper) 945 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 954 946 #endif 955 947 } else if (value > upper + tolerance) { 956 957 948 value = upper; 949 value *= value; 958 950 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 959 960 951 if (lower == upper) 952 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 961 953 #endif 962 954 } else { 963 964 value=0.0;955 // feasible 956 value = 0.0; 965 957 } 966 958 // store square 967 work[iRow]=value; 968 } 969 } 970 } 971 #endif 972 void 973 AbcDualRowSteepest::updatePrimalSolutionAndWeights(CoinIndexedVector & weightsVector, 974 CoinIndexedVector & primalUpdate, 975 double theta) 959 work[iRow] = value; 960 } 961 } 962 } 963 #endif 964 void AbcDualRowSteepest::updatePrimalSolutionAndWeights(CoinIndexedVector &weightsVector, 965 CoinIndexedVector &primalUpdate, 966 double theta) 976 967 { 977 968 //printf("updatePrimal iteration %d  weightsVector has %d, primalUpdate has %d\n", 978 969 // model_>numberIterations(),weightsVector.getNumElements(),primalUpdate.getNumElements()); 979 double * 980 double * 970 double *COIN_RESTRICT weights = weights_>denseVector(); 971 double *COIN_RESTRICT work2 = weightsVector.denseVector(); 981 972 int pivotRow = model_>lastPivotRow(); 982 assert (model_>pivotRow()==model_>lastPivotRow());983 double * 973 assert(model_>pivotRow() == model_>lastPivotRow()); 974 double *COIN_RESTRICT work = primalUpdate.denseVector(); 984 975 int numberNonZero = primalUpdate.getNumElements(); 985 int * 976 int *COIN_RESTRICT which = primalUpdate.getIndices(); 986 977 double tolerance = model_>currentPrimalTolerance(); 987 double * 988 double * 989 const double * 990 const double * 991 assert 978 double *COIN_RESTRICT infeas = infeasible_>denseVector(); 979 double *COIN_RESTRICT solutionBasic = model_>solutionBasic(); 980 const double *COIN_RESTRICT lowerBasic = model_>lowerBasic(); 981 const double *COIN_RESTRICT upperBasic = model_>upperBasic(); 982 assert(!primalUpdate.packedMode()); 992 983 #ifdef CLP_DEBUG 993 if ( model_>logLevel ( ) > 4 && 994 fabs(norm_  weights[pivotRow]) > 1.0e2 * (1.0 + norm_)) 984 if (model_>logLevel() > 4 && fabs(norm_  weights[pivotRow]) > 1.0e2 * (1.0 + norm_)) 995 985 printf("on row %d, true weight %g, old %g\n", 996 986 pivotRow, sqrt(norm_), sqrt(weights[pivotRow])); 997 987 #endif 998 988 // pivot element 999 double alpha =model_>alpha();989 double alpha = model_>alpha(); 1000 990 // could reinitialize here (could be expensive) 1001 991 norm_ /= alpha * alpha; … … 1006 996 // only weights3 is used?  slightly faster to just update weights if going to invert 1007 997 xweights3++; 1008 if ((xweights1 +xweights2+xweights3)%1000==0) {998 if ((xweights1 + xweights2 + xweights3) % 1000 == 0) { 1009 999 printf("ZZ iteration %d sum weights %d  %d %d %d\n", 1010 model_>numberIterations(),xweights1+xweights2+xweights3,xweights1,xweights2,xweights3);1011 assert(abs(model_>numberIterations() (xweights1+xweights2+xweights3))<=1);1000 model_>numberIterations(), xweights1 + xweights2 + xweights3, xweights1, xweights2, xweights3); 1001 assert(abs(model_>numberIterations()  (xweights1 + xweights2 + xweights3)) <= 1); 1012 1002 } 1013 1003 #endif 1014 1004 // look at updated column 1015 double multiplier2 =1.0;1016 if (model_>directionOut() <0) {1005 double multiplier2 = 1.0; 1006 if (model_>directionOut() < 0) { 1017 1007 #ifdef ABC_DEBUG 1018 printf("Minus out in %d %d, out %d %d\n", model_>sequenceIn(),model_>directionIn(),1019 model_>sequenceOut(),model_>directionOut());1020 #endif 1021 multiplier2 =1.0;1022 } 1023 #if ABC_PARALLEL ==21024 update(0, numberNonZero,which,work,work2,weights,1025 lowerBasic,solutionBasic,upperBasic,1026 multiplier,multiplier2,norm_,theta,tolerance);1008 printf("Minus out in %d %d, out %d %d\n", model_>sequenceIn(), model_>directionIn(), 1009 model_>sequenceOut(), model_>directionOut()); 1010 #endif 1011 multiplier2 = 1.0; 1012 } 1013 #if ABC_PARALLEL == 2 1014 update(0, numberNonZero, which, work, work2, weights, 1015 lowerBasic, solutionBasic, upperBasic, 1016 multiplier, multiplier2, norm_, theta, tolerance); 1027 1017 for (int i = 0; i < numberNonZero; i++) { 1028 1018 int iRow = which[i]; 1029 double infeasibility =work[iRow];1030 work[iRow] =0.0;1019 double infeasibility = work[iRow]; 1020 work[iRow] = 0.0; 1031 1021 if (infeasibility) { 1032 1022 if (infeas[iRow]) 1033 1023 infeas[iRow] = infeasibility; // already there 1034 1024 else 1035 1025 infeasible_>quickAdd(iRow, infeasibility); 1036 1026 } else { 1037 1027 // feasible  was it infeasible  if so set tiny 1038 1028 if (infeas[iRow]) 1039 1029 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 1040 1030 } 1041 1031 } … … 1044 1034 int iRow = which[i]; 1045 1035 double updateValue = work[iRow]; 1046 work[iRow] =0.0;1047 double thetaDevex = multiplier2 *updateValue;1036 work[iRow] = 0.0; 1037 double thetaDevex = multiplier2 * updateValue; 1048 1038 double devex = weights[iRow]; 1049 1039 double valueDevex = work2[iRow]; 1050 devex += 1040 devex += thetaDevex * (thetaDevex * norm_ + valueDevex * multiplier); 1051 1041 if (devex < DEVEX_TRY_NORM) 1052 1042 devex = DEVEX_TRY_NORM; … … 1063 1053 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1064 1054 if (lower == upper) 1065 1055 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1066 1056 #endif 1067 1057 // store square in list 1068 1058 if (infeas[iRow]) 1069 1059 infeas[iRow] = value; // already there 1070 1060 else 1071 1061 infeasible_>quickAdd(iRow, value); 1072 1062 } else if (value > upper + tolerance) { 1073 1063 value = upper; … … 1075 1065 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1076 1066 if (lower == upper) 1077 1067 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1078 1068 #endif 1079 1069 // store square in list 1080 1070 if (infeas[iRow]) 1081 1071 infeas[iRow] = value; // already there 1082 1072 else 1083 1073 infeasible_>quickAdd(iRow, value); 1084 1074 } else { 1085 1075 // feasible  was it infeasible  if so set tiny 1086 1076 if (infeas[iRow]) 1087 1077 infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 1088 1078 } 1089 1079 } … … 1106 1096 // see 1/2 above primalUpdate.setNumElements(0); 1107 1097 } 1108 extern void parallelDual5(AbcSimplexFactorization * 1109 CoinIndexedVector **whichVector,1110 1111 1112 double *weights);1098 extern void parallelDual5(AbcSimplexFactorization *factorization, 1099 CoinIndexedVector **whichVector, 1100 int numberCpu, 1101 int whichCpu, 1102 double *weights); 1113 1103 /* Saves any weights round factorization as pivot rows may change 1114 1104 1) before factorization … … 1117 1107 4) restore weights 1118 1108 */ 1119 void 1120 AbcDualRowSteepest::saveWeights(AbcSimplex * model, int mode) 1109 void AbcDualRowSteepest::saveWeights(AbcSimplex *model, int mode) 1121 1110 { 1122 1111 model_ = model; 1123 1112 int numberRows = model_>numberRows(); 1124 1113 // faster to save indices in array numberTotal in length 1125 const int * 1126 double * 1114 const int *pivotVariable = model_>pivotVariable(); 1115 double *weights = weights_ ? weights_>denseVector() : NULL; 1127 1116 // see if we need to switch 1128 1117 //if (mode_==2&&mode==2) { 1129 1130 1118 // switch 1119 //mode=5; 1131 1120 //} 1132 1121 if (mode == 1) { 1133 if (weights_) {1122 if (weights_) { 1134 1123 // Check if size has changed 1135 1124 if (infeasible_>capacity() == numberRows) { 1136 1137 CoinAbcMemcpy(weights_>getIndices(),pivotVariable,numberRows);1138 1125 // change from row numbers to sequence numbers 1126 CoinAbcMemcpy(weights_>getIndices(), pivotVariable, numberRows); 1127 state_ = 1; 1139 1128 } else { 1140 1141 1142 1143 1144 1145 1146 1147 1148 } 1149 } 1150 } else if (mode != 3) {1129 // size has changed  clear everything 1130 delete weights_; 1131 weights_ = NULL; 1132 delete infeasible_; 1133 infeasible_ = NULL; 1134 delete savedWeights_; 1135 savedWeights_ = NULL; 1136 state_ = 1; 1137 } 1138 } 1139 } else if (mode != 3) { 1151 1140 // restore 1152 1141 if (!weights_  state_ == 1  mode == 5) { 1153 1142 // see if crash basis??? 1154 int numberBasicSlacks=0; 1155 bool initializeNorms=false; 1156 for (int i=0;i<numberRows;i++) { 1157 if (model_>getInternalStatus(i)==AbcSimplex::basic) 1158 numberBasicSlacks++; 1159 } 1160 if (numberBasicSlacks<numberRows&&!model_>numberIterations()) { 1161 char line[100]; 1162 if (numberBasicSlacks>numberRows(numberRows>>1)&&mode_==3&& 1163 (model_>moreSpecialOptions()&256)==0) { 1164 sprintf(line,"Steep initial basis has %d structurals out of %d  initializing norms\n",numberRowsnumberBasicSlacks,numberRows); 1165 initializeNorms=true; 1166 } else { 1167 sprintf(line,"Steep initial basis has %d structurals out of %d  too many\n",numberRowsnumberBasicSlacks,numberRows); 1168 } 1169 model_>messageHandler()>message(CLP_GENERAL2,*model_>messagesPointer()) 1170 << line << CoinMessageEol; 1143 int numberBasicSlacks = 0; 1144 bool initializeNorms = false; 1145 for (int i = 0; i < numberRows; i++) { 1146 if (model_>getInternalStatus(i) == AbcSimplex::basic) 1147 numberBasicSlacks++; 1148 } 1149 if (numberBasicSlacks < numberRows && !model_>numberIterations()) { 1150 char line[100]; 1151 if (numberBasicSlacks > numberRows  (numberRows >> 1) && mode_ == 3 && (model_>moreSpecialOptions() & 256) == 0) { 1152 sprintf(line, "Steep initial basis has %d structurals out of %d  initializing norms\n", numberRows  numberBasicSlacks, numberRows); 1153 initializeNorms = true; 1154 } else { 1155 sprintf(line, "Steep initial basis has %d structurals out of %d  too many\n", numberRows  numberBasicSlacks, numberRows); 1156 } 1157 model_>messageHandler()>message(CLP_GENERAL2, *model_>messagesPointer()) 1158 << line << CoinMessageEol; 1171 1159 } 1172 1160 // initialize weights … … 1178 1166 savedWeights_ = new CoinIndexedVector(); 1179 1167 // allow for dense multiple of 8! 1180 savedWeights_>reserve(numberRows +((initializeNorms) ? 8 : 0));1168 savedWeights_>reserve(numberRows + ((initializeNorms) ? 8 : 0)); 1181 1169 infeasible_ = new CoinIndexedVector(); 1182 1170 infeasible_>reserve(numberRows); 1183 1171 weights = weights_>denseVector(); 1184 if (mode == 5 (mode_!=1&&!initializeNorms)) {1185 1186 1187 1188 1172 if (mode == 5  (mode_ != 1 && !initializeNorms)) { 1173 // initialize to 1.0 (can we do better?) 1174 for (int i = 0; i < numberRows; i++) { 1175 weights[i] = 1.0; 1176 } 1189 1177 } else { 1190 double *COIN_RESTRICT array = savedWeights_>denseVector();1191 int *COIN_RESTRICT which = savedWeights_>getIndices();1192 1193 1178 double *COIN_RESTRICT array = savedWeights_>denseVector(); 1179 int *COIN_RESTRICT which = savedWeights_>getIndices(); 1180 // Later look at factorization to see which can be skipped 1181 // also use threads to parallelize 1194 1182 #if 0 1195 1183 // use infeasible to collect row counts … … 1213 1201 int nSlack=0; 1214 1202 #endif 1215 1216 1217 #if ABC_PARALLEL ==21218 if (model_>parallelMode()==0) {1219 #endif 1220 1221 1222 1223 1224 1225 1226 1227 1228 int k= which[j];1229 1230 1231 1203 //#undef ABC_PARALLEL 1204 // for now just with cilk 1205 #if ABC_PARALLEL == 2 1206 if (model_>parallelMode() == 0) { 1207 #endif 1208 for (int i = 0; i < numberRows; i++) { 1209 double value = 0.0; 1210 array[i] = 1.0; 1211 which[0] = i; 1212 savedWeights_>setNumElements(1); 1213 model_>factorization()>updateColumnTranspose(*savedWeights_); 1214 int number = savedWeights_>getNumElements(); 1215 for (int j = 0; j < number; j++) { 1216 int k = which[j]; 1217 value += array[k] * array[k]; 1218 array[k] = 0.0; 1219 } 1232 1220 #if 0 1233 1221 if (!counts[i]) { … … 1236 1224 } 1237 1225 #endif 1238 1239 1240 1241 #if ABC_PARALLEL ==21242 1243 1244 1245 int numberCpuMinusOne=CoinMin(model_>parallelMode(),3);1246 1247 CoinIndexedVector *whichVector[4];1248 for (int i=0;i<numberCpuMinusOne;i++) {1249 which[i]=model_>getAvailableArray();1250 whichVector[i]=model_>usefulArray(which[i]);1251 1252 1253 whichVector[numberCpuMinusOne]=savedWeights_;1254 parallelDual5(model_>factorization(),whichVector,numberCpuMinusOne+1,numberCpuMinusOne,weights);1255 for (int i=0;i<numberCpuMinusOne;i++) 1256 1257 1258 #endif 1226 savedWeights_>setNumElements(0); 1227 weights[i] = value; 1228 } 1229 #if ABC_PARALLEL == 2 1230 } else { 1231 // parallel 1232 // get arrays 1233 int numberCpuMinusOne = CoinMin(model_>parallelMode(), 3); 1234 int which[3]; 1235 CoinIndexedVector *whichVector[4]; 1236 for (int i = 0; i < numberCpuMinusOne; i++) { 1237 which[i] = model_>getAvailableArray(); 1238 whichVector[i] = model_>usefulArray(which[i]); 1239 assert(!whichVector[i]>packedMode()); 1240 } 1241 whichVector[numberCpuMinusOne] = savedWeights_; 1242 parallelDual5(model_>factorization(), whichVector, numberCpuMinusOne + 1, numberCpuMinusOne, weights); 1243 for (int i = 0; i < numberCpuMinusOne; i++) 1244 model_>clearArrays(which[i]); 1245 } 1246 #endif 1259 1247 #if 0 1260 1248 printf("Steep %d can be skipped\n",nSlack); … … 1262 1250 } 1263 1251 } else { 1264 int useArray =model_>getAvailableArrayPublic();1265 CoinIndexedVector * 1266 int * 1267 double * 1268 if (mode !=4) {1269 CoinAbcMemcpy(which,weights_>getIndices(),numberRows);1270 CoinAbcScatterTo(weights_>denseVector(),element,which,numberRows);1252 int useArray = model_>getAvailableArrayPublic(); 1253 CoinIndexedVector *usefulArray = model_>usefulArray(useArray); 1254 int *COIN_RESTRICT which = usefulArray>getIndices(); 1255 double *COIN_RESTRICT element = usefulArray>denseVector(); 1256 if (mode != 4) { 1257 CoinAbcMemcpy(which, weights_>getIndices(), numberRows); 1258 CoinAbcScatterTo(weights_>denseVector(), element, which, numberRows); 1271 1259 } else { 1272 1273 CoinAbcMemcpy(which,savedWeights_>getIndices(),numberRows);1274 CoinAbcScatterTo(savedWeights_>denseVector(),element,which,numberRows);1260 // restore weights 1261 CoinAbcMemcpy(which, savedWeights_>getIndices(), numberRows); 1262 CoinAbcScatterTo(savedWeights_>denseVector(), element, which, numberRows); 1275 1263 } 1276 1264 for (int i = 0; i < numberRows; i++) { 1277 1278 double thisWeight=element[iPivot];1279 1280 1281 1282 1283 1284 1285 1286 1265 int iPivot = pivotVariable[i]; 1266 double thisWeight = element[iPivot]; 1267 if (thisWeight) { 1268 if (thisWeight < DEVEX_TRY_NORM) 1269 thisWeight = DEVEX_TRY_NORM; // may need to check more 1270 weights[i] = thisWeight; 1271 } else { 1272 // odd 1273 weights[i] = 1.0; 1274 } 1287 1275 } 1288 1276 if (0) { 1289 1290 1291 assert (numberRows<=1000);1292 memcpy(w,weights,numberRows*sizeof(double));1293 memcpy(s,pivotVariable,numberRows*sizeof(int));1294 CoinSort_2(s,s+numberRows,w);1295 1296 for (int i=0;i<numberRows;i++) {1297 printf("%d seq %d weight %g\n",i,s[i],w[i]);1298 1299 1277 double w[1000]; 1278 int s[1000]; 1279 assert(numberRows <= 1000); 1280 memcpy(w, weights, numberRows * sizeof(double)); 1281 memcpy(s, pivotVariable, numberRows * sizeof(int)); 1282 CoinSort_2(s, s + numberRows, w); 1283 printf("Weights===========\n"); 1284 for (int i = 0; i < numberRows; i++) { 1285 printf("%d seq %d weight %g\n", i, s[i], w[i]); 1286 } 1287 printf("End weights===========\n"); 1300 1288 } 1301 1289 usefulArray>setNumElements(numberRows); … … 1305 1293 } 1306 1294 state_ = 0; 1307 CoinAbcMemcpy(savedWeights_>denseVector(), weights_>denseVector(),numberRows);1308 CoinAbcMemcpy(savedWeights_>getIndices(), pivotVariable,numberRows);1295 CoinAbcMemcpy(savedWeights_>denseVector(), weights_>denseVector(), numberRows); 1296 CoinAbcMemcpy(savedWeights_>getIndices(), pivotVariable, numberRows); 1309 1297 } 1310 1298 if (mode >= 2) { 1311 1299 infeasible_>clear(); 1312 1300 double tolerance = model_>currentPrimalTolerance(); 1313 const double * 1314 const double * 1315 const double * 1316 if ((model_>stateOfProblem() &VALUES_PASS2)==0) {1301 const double *COIN_RESTRICT solutionBasic = model_>solutionBasic(); 1302 const double *COIN_RESTRICT lowerBasic = model_>lowerBasic(); 1303 const double *COIN_RESTRICT upperBasic = model_>upperBasic(); 1304 if ((model_>stateOfProblem() & VALUES_PASS2) == 0) { 1317 1305 for (int iRow = 0; iRow < numberRows; iRow++) { 1318 1319 1320 1321 1322 1323 1306 double value = solutionBasic[iRow]; 1307 double lower = lowerBasic[iRow]; 1308 double upper = upperBasic[iRow]; 1309 if (value < lower  tolerance) { 1310 value = lower; 1311 value *= value; 1324 1312 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1325 1326 1327 #endif 1328 1329 1330 1331 1332 1313 if (lower == upper) 1314 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1315 #endif 1316 // store square in list 1317 infeasible_>quickAdd(iRow, value); 1318 } else if (value > upper + tolerance) { 1319 value = upper; 1320 value *= value; 1333 1321 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1334 1335 1336 #endif 1337 1338 1339 1322 if (lower == upper) 1323 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1324 #endif 1325 // store square in list 1326 infeasible_>quickAdd(iRow, value); 1327 } 1340 1328 } 1341 1329 } else { 1342 1330 // check values pass 1343 const double * 1344 int numberFake =0;1331 const double *fakeDjs = model_>fakeDjs(); 1332 int numberFake = 0; 1345 1333 for (int iRow = 0; iRow < numberRows; iRow++) { 1346 1347 1348 1349 1350 1351 1352 1353 1354 value=0.0;1355 1356 1357 1334 double value = solutionBasic[iRow]; 1335 double lower = lowerBasic[iRow]; 1336 double upper = upperBasic[iRow]; 1337 if (value < lower  tolerance) { 1338 value = lower; 1339 } else if (value > upper + tolerance) { 1340 value = upper; 1341 } else { 1342 value = 0.0; 1343 } 1344 if (value) { 1345 value *= value; 1358 1346 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1359 1360 1361 #endif 1362 1363 1364 int iSequence=pivotVariable[iRow];1365 double fakeDj=fakeDjs[iSequence];1366 if (fabs(fakeDj)>10.0*tolerance) 1367 1368 1347 if (lower == upper) 1348 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1349 #endif 1350 // store square in list 1351 infeasible_>quickAdd(iRow, value); 1352 int iSequence = pivotVariable[iRow]; 1353 double fakeDj = fakeDjs[iSequence]; 1354 if (fabs(fakeDj) > 10.0 * tolerance) 1355 numberFake++; 1356 } 1369 1357 } 1370 1358 if (numberFake) 1371 printf("%d fake basic djs\n",numberFake);1359 printf("%d fake basic djs\n", numberFake); 1372 1360 else 1373 model_>setStateOfProblem(model_>stateOfProblem()&~VALUES_PASS2);1361 model_>setStateOfProblem(model_>stateOfProblem() & ~VALUES_PASS2); 1374 1362 } 1375 1363 } 1376 1364 } 1377 1365 // Recompute infeasibilities 1378 void 1379 AbcDualRowSteepest::recomputeInfeasibilities() 1366 void AbcDualRowSteepest::recomputeInfeasibilities() 1380 1367 { 1381 1368 int numberRows = model_>numberRows(); 1382 1369 infeasible_>clear(); 1383 1370 double tolerance = model_>currentPrimalTolerance(); 1384 const double * 1385 const double * 1386 const double * 1371 const double *COIN_RESTRICT solutionBasic = model_>solutionBasic(); 1372 const double *COIN_RESTRICT lowerBasic = model_>lowerBasic(); 1373 const double *COIN_RESTRICT upperBasic = model_>upperBasic(); 1387 1374 for (int iRow = 0; iRow < numberRows; iRow++) { 1388 1375 double value = solutionBasic[iRow]; … … 1394 1381 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1395 1382 if (lower == upper) 1396 1383 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1397 1384 #endif 1398 1385 // store square in list … … 1403 1390 #ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER 1404 1391 if (lower == upper) 1405 1392 value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables 1406 1393 #endif 1407 1394 // store square in list … … 1413 1400 // Clone 1414 1401 // 1415 AbcDualRowPivot * 1402 AbcDualRowPivot *AbcDualRowSteepest::clone(bool CopyData) const 1416 1403 { 1417 1404 if (CopyData) { … … 1422 1409 } 1423 1410 // Gets rid of all arrays 1424 void 1425 AbcDualRowSteepest::clearArrays() 1411 void AbcDualRowSteepest::clearArrays() 1426 1412 { 1427 1413 if (persistence_ == normal) { … … 1436 1422 } 1437 1423 // Returns true if would not find any row 1438 bool 1439 AbcDualRowSteepest::looksOptimal() const 1424 bool AbcDualRowSteepest::looksOptimal() const 1440 1425 { 1441 1426 int iRow; … … 1445 1430 double error = CoinMin(1.0e2, model_>largestPrimalError()); 1446 1431 // allow tolerance at least slightly bigger than standard 1447 tolerance = tolerance + 1432 tolerance = tolerance + error; 1448 1433 // But cap 1449 1434 tolerance = CoinMin(1000.0, tolerance); 1450 1435 int numberRows = model_>numberRows(); 1451 1436 int numberInfeasible = 0; 1452 const double * 1453 const double * 1454 const double * 1437 const double *COIN_RESTRICT solutionBasic = model_>solutionBasic(); 1438 const double *COIN_RESTRICT lowerBasic = model_>lowerBasic(); 1439 const double *COIN_RESTRICT upperBasic = model_>upperBasic(); 1455 1440 for (iRow = 0; iRow < numberRows; iRow++) { 1456 1441 double value = solutionBasic[iRow]; … … 1465 1450 return (numberInfeasible == 0); 1466 1451 } 1452 1453 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 1454 */
Note: See TracChangeset
for help on using the changeset viewer.