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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/AbcMatrix.cpp
r2074 r2385 3 3 // Corporation and others, Copyright (C) 2012, FasterCoin. All Rights Reserved. 4 4 // This code is licensed under the terms of the Eclipse Public License (EPL). 5 6 7 5 8 6 #include <cstdio> … … 24 22 #include "mkl_spblas.h" 25 23 #endif 26 #if ABC_INSTRUMENT >124 #if ABC_INSTRUMENT > 1 27 25 extern int abcPricing[20]; 28 26 extern int abcPricingDense[20]; … … 37 35 // Default Constructor 38 36 // 39 AbcMatrix::AbcMatrix () 40 : matrix_(NULL), 41 model_(NULL), 42 rowStart_(NULL), 43 element_(NULL), 44 column_(NULL), 45 numberColumnBlocks_(0), 46 numberRowBlocks_(0), 37 AbcMatrix::AbcMatrix() 38 : matrix_(NULL) 39 , model_(NULL) 40 , rowStart_(NULL) 41 , element_(NULL) 42 , column_(NULL) 43 , numberColumnBlocks_(0) 44 , numberRowBlocks_(0) 45 , 47 46 #ifdef COUNT_COPY 48 countRealColumn_(NULL), 49 countStartLarge_(NULL), 50 countRow_(NULL), 51 countElement_(NULL), 52 smallestCount_(0), 53 largestCount_(0), 54 #endif 55 startFraction_(0.0), 56 endFraction_(1.0), 57 savedBestDj_(0.0), 58 originalWanted_(0), 59 currentWanted_(0), 60 savedBestSequence_(1), 61 minimumObjectsScan_(1), 62 minimumGoodReducedCosts_(1) 47 countRealColumn_(NULL) 48 , countStartLarge_(NULL) 49 , countRow_(NULL) 50 , countElement_(NULL) 51 , smallestCount_(0) 52 , largestCount_(0) 53 , 54 #endif 55 startFraction_(0.0) 56 , endFraction_(1.0) 57 , savedBestDj_(0.0) 58 , originalWanted_(0) 59 , currentWanted_(0) 60 , savedBestSequence_(1) 61 , minimumObjectsScan_(1) 62 , minimumGoodReducedCosts_(1) 63 63 { 64 64 } … … 67 67 // Copy constructor 68 68 // 69 AbcMatrix::AbcMatrix (const AbcMatrix &rhs)69 AbcMatrix::AbcMatrix(const AbcMatrix &rhs) 70 70 { 71 71 #ifndef COIN_SPARSE_MATRIX … … 74 74 matrix_ = new CoinPackedMatrix(*(rhs.matrix_), 0, 0); 75 75 #endif 76 model_ =rhs.model_;76 model_ = rhs.model_; 77 77 rowStart_ = NULL; 78 78 element_ = NULL; … … 85 85 #endif 86 86 numberColumnBlocks_ = rhs.numberColumnBlocks_; 87 CoinAbcMemcpy(startColumnBlock_, rhs.startColumnBlock_,numberColumnBlocks_+1);87 CoinAbcMemcpy(startColumnBlock_, rhs.startColumnBlock_, numberColumnBlocks_ + 1); 88 88 numberRowBlocks_ = rhs.numberRowBlocks_; 89 89 if (numberRowBlocks_) { 90 assert 91 int numberRows =model_>numberRows();90 assert(model_); 91 int numberRows = model_>numberRows(); 92 92 int numberElements = matrix_>getNumElements(); 93 memcpy(blockStart_, rhs.blockStart_,sizeof(blockStart_));94 rowStart_ =CoinCopyOfArray(rhs.rowStart_,numberRows*(numberRowBlocks_+2));95 element_ =CoinCopyOfArray(rhs.element_,numberElements);96 column_ =CoinCopyOfArray(rhs.column_,numberElements);93 memcpy(blockStart_, rhs.blockStart_, sizeof(blockStart_)); 94 rowStart_ = CoinCopyOfArray(rhs.rowStart_, numberRows * (numberRowBlocks_ + 2)); 95 element_ = CoinCopyOfArray(rhs.element_, numberElements); 96 column_ = CoinCopyOfArray(rhs.column_, numberElements); 97 97 #ifdef COUNT_COPY 98 98 smallestCount_ = rhs.smallestCount_; 99 99 largestCount_ = rhs.largestCount_; 100 int numberColumns=model_>numberColumns(); 101 countRealColumn_=CoinCopyOfArray(rhs.countRealColumn_,numberColumns); 102 memcpy(countStart_,rhs.countStart_,reinterpret_cast<char *>(&countRealColumn_) 103 reinterpret_cast<char *>(countStart_)); 104 int numberLarge = numberColumnscountStart_[MAX_COUNT]; 105 countStartLarge_=CoinCopyOfArray(rhs.countStartLarge_,numberLarge+1); 106 numberElements=countStartLarge_[numberLarge]; 107 countElement_=CoinCopyOfArray(rhs.countElement_,numberElements); 108 countRow_=CoinCopyOfArray(rhs.countRow_,numberElements); 100 int numberColumns = model_>numberColumns(); 101 countRealColumn_ = CoinCopyOfArray(rhs.countRealColumn_, numberColumns); 102 memcpy(countStart_, rhs.countStart_, reinterpret_cast< char * >(&countRealColumn_)  reinterpret_cast< char * >(countStart_)); 103 int numberLarge = numberColumns  countStart_[MAX_COUNT]; 104 countStartLarge_ = CoinCopyOfArray(rhs.countStartLarge_, numberLarge + 1); 105 numberElements = countStartLarge_[numberLarge]; 106 countElement_ = CoinCopyOfArray(rhs.countElement_, numberElements); 107 countRow_ = CoinCopyOfArray(rhs.countRow_, numberElements); 109 108 #endif 110 109 } 111 110 } 112 111 113 AbcMatrix::AbcMatrix (const CoinPackedMatrix &rhs)112 AbcMatrix::AbcMatrix(const CoinPackedMatrix &rhs) 114 113 { 115 114 #ifndef COIN_SPARSE_MATRIX … … 119 118 #endif 120 119 matrix_>cleanMatrix(); 121 model_ =NULL;120 model_ = NULL; 122 121 rowStart_ = NULL; 123 122 element_ = NULL; … … 131 130 largestCount_ = 0; 132 131 #endif 133 numberColumnBlocks_ =1;134 startColumnBlock_[0] =0;135 startColumnBlock_[1] =0;132 numberColumnBlocks_ = 1; 133 startColumnBlock_[0] = 0; 134 startColumnBlock_[1] = 0; 136 135 numberRowBlocks_ = 0; 137 136 startFraction_ = 0; … … 146 145 #ifdef ABC_SPRINT 147 146 /* Subset constructor (without gaps). */ 148 AbcMatrix::AbcMatrix (const AbcMatrix &wholeMatrix,149 int numberRows, const int *whichRows,150 int numberColumns, const int *whichColumns)147 AbcMatrix::AbcMatrix(const AbcMatrix &wholeMatrix, 148 int numberRows, const int *whichRows, 149 int numberColumns, const int *whichColumns) 151 150 { 152 151 #ifndef COIN_SPARSE_MATRIX 153 matrix_ = new CoinPackedMatrix(*wholeMatrix.matrix_, numberRows,whichRows,154 numberColumns,whichColumns);152 matrix_ = new CoinPackedMatrix(*wholeMatrix.matrix_, numberRows, whichRows, 153 numberColumns, whichColumns); 155 154 #else 156 155 matrix_ = new CoinPackedMatrix(rhs, 0, 0); … … 158 157 #endif 159 158 matrix_>cleanMatrix(); 160 model_ =NULL;159 model_ = NULL; 161 160 rowStart_ = NULL; 162 161 element_ = NULL; … … 170 169 largestCount_ = 0; 171 170 #endif 172 numberColumnBlocks_ =1;173 startColumnBlock_[0] =0;174 startColumnBlock_[1] =0;171 numberColumnBlocks_ = 1; 172 startColumnBlock_[0] = 0; 173 startColumnBlock_[1] = 0; 175 174 numberRowBlocks_ = 0; 176 175 startFraction_ = 0; … … 187 186 // Destructor 188 187 // 189 AbcMatrix::~AbcMatrix 188 AbcMatrix::~AbcMatrix() 190 189 { 191 190 delete matrix_; 192 delete 193 delete 194 delete 191 delete[] rowStart_; 192 delete[] element_; 193 delete[] column_; 195 194 #ifdef COUNT_COPY 196 delete 197 delete 198 delete 199 delete 195 delete[] countRealColumn_; 196 delete[] countStartLarge_; 197 delete[] countRow_; 198 delete[] countElement_; 200 199 #endif 201 200 } … … 205 204 // 206 205 AbcMatrix & 207 AbcMatrix::operator=(const AbcMatrix &rhs)206 AbcMatrix::operator=(const AbcMatrix &rhs) 208 207 { 209 208 if (this != &rhs) { … … 214 213 matrix_ = new CoinPackedMatrix(*(rhs.matrix_), 0, 0); 215 214 #endif 216 model_ =rhs.model_;217 delete 218 delete 219 delete 215 model_ = rhs.model_; 216 delete[] rowStart_; 217 delete[] element_; 218 delete[] column_; 220 219 #ifdef COUNT_COPY 221 delete 222 delete 223 delete 224 delete 220 delete[] countRealColumn_; 221 delete[] countStartLarge_; 222 delete[] countRow_; 223 delete[] countElement_; 225 224 #endif 226 225 rowStart_ = NULL; … … 234 233 #endif 235 234 numberColumnBlocks_ = rhs.numberColumnBlocks_; 236 CoinAbcMemcpy(startColumnBlock_, rhs.startColumnBlock_,numberColumnBlocks_+1);235 CoinAbcMemcpy(startColumnBlock_, rhs.startColumnBlock_, numberColumnBlocks_ + 1); 237 236 numberRowBlocks_ = rhs.numberRowBlocks_; 238 237 if (numberRowBlocks_) { 239 assert 240 int numberRows =model_>numberRows();238 assert(model_); 239 int numberRows = model_>numberRows(); 241 240 int numberElements = matrix_>getNumElements(); 242 memcpy(blockStart_, rhs.blockStart_,sizeof(blockStart_));243 rowStart_ =CoinCopyOfArray(rhs.rowStart_,numberRows*(numberRowBlocks_+2));244 element_ =CoinCopyOfArray(rhs.element_,numberElements);245 column_ =CoinCopyOfArray(rhs.column_,numberElements);241 memcpy(blockStart_, rhs.blockStart_, sizeof(blockStart_)); 242 rowStart_ = CoinCopyOfArray(rhs.rowStart_, numberRows * (numberRowBlocks_ + 2)); 243 element_ = CoinCopyOfArray(rhs.element_, numberElements); 244 column_ = CoinCopyOfArray(rhs.column_, numberElements); 246 245 #ifdef COUNT_COPY 247 246 smallestCount_ = rhs.smallestCount_; 248 247 largestCount_ = rhs.largestCount_; 249 int numberColumns=model_>numberColumns(); 250 countRealColumn_=CoinCopyOfArray(rhs.countRealColumn_,numberColumns); 251 memcpy(countStart_,rhs.countStart_,reinterpret_cast<char *>(&countRealColumn_) 252 reinterpret_cast<char *>(countStart_)); 253 int numberLarge = numberColumnscountStart_[MAX_COUNT]; 254 countStartLarge_=CoinCopyOfArray(rhs.countStartLarge_,numberLarge+1); 255 numberElements=countStartLarge_[numberLarge]; 256 countElement_=CoinCopyOfArray(rhs.countElement_,numberElements); 257 countRow_=CoinCopyOfArray(rhs.countRow_,numberElements); 248 int numberColumns = model_>numberColumns(); 249 countRealColumn_ = CoinCopyOfArray(rhs.countRealColumn_, numberColumns); 250 memcpy(countStart_, rhs.countStart_, reinterpret_cast< char * >(&countRealColumn_)  reinterpret_cast< char * >(countStart_)); 251 int numberLarge = numberColumns  countStart_[MAX_COUNT]; 252 countStartLarge_ = CoinCopyOfArray(rhs.countStartLarge_, numberLarge + 1); 253 numberElements = countStartLarge_[numberLarge]; 254 countElement_ = CoinCopyOfArray(rhs.countElement_, numberElements); 255 countRow_ = CoinCopyOfArray(rhs.countRow_, numberElements); 258 256 #endif 259 257 } … … 270 268 } 271 269 // Sets model 272 void 273 AbcMatrix::setModel(AbcSimplex * model) 274 { 275 model_=model; 276 int numberColumns=model_>numberColumns(); 277 bool needExtension=numberColumns>matrix_>getNumCols(); 270 void AbcMatrix::setModel(AbcSimplex *model) 271 { 272 model_ = model; 273 int numberColumns = model_>numberColumns(); 274 bool needExtension = numberColumns > matrix_>getNumCols(); 278 275 if (needExtension) { 279 276 CoinBigIndex lastElement = matrix_>getNumElements(); 280 matrix_>reserve(numberColumns, lastElement,true);281 CoinBigIndex * 282 for (int i =numberColumns;i>=0;i) {283 if (columnStart[i] ==0)284 columnStart[i]=lastElement;277 matrix_>reserve(numberColumns, lastElement, true); 278 CoinBigIndex *columnStart = matrix_>getMutableVectorStarts(); 279 for (int i = numberColumns; i >= 0; i) { 280 if (columnStart[i] == 0) 281 columnStart[i] = lastElement; 285 282 else 286 287 } 288 assert (lastElement==columnStart[numberColumns]);283 break; 284 } 285 assert(lastElement == columnStart[numberColumns]); 289 286 } 290 287 } … … 293 290 AbcMatrix::reverseOrderedCopy() const 294 291 { 295 CoinPackedMatrix * 292 CoinPackedMatrix *matrix = new CoinPackedMatrix(); 296 293 matrix>setExtraGap(0.0); 297 294 matrix>setExtraMajor(0.0); … … 301 298 /// returns number of elements in column part of basis, 302 299 CoinBigIndex 303 AbcMatrix::countBasis( const int *whichColumn,304 int &numberColumnBasic)305 { 306 const int * 300 AbcMatrix::countBasis(const int *whichColumn, 301 int &numberColumnBasic) 302 { 303 const int *COIN_RESTRICT columnLength = matrix_>getVectorLengths(); 307 304 CoinBigIndex numberElements = 0; 308 int numberRows =model_>numberRows();305 int numberRows = model_>numberRows(); 309 306 // just count  can be over so ignore zero problem 310 307 for (int i = 0; i < numberColumnBasic; i++) { 311 int iColumn = whichColumn[i] numberRows;308 int iColumn = whichColumn[i]  numberRows; 312 309 numberElements += columnLength[iColumn]; 313 310 } 314 311 return numberElements; 315 312 } 316 void 317 AbcMatrix::fillBasis(const int * COIN_RESTRICT whichColumn, 318 int & numberColumnBasic, 319 int * COIN_RESTRICT indexRowU, 320 int * COIN_RESTRICT start, 321 int * COIN_RESTRICT rowCount, 322 int * COIN_RESTRICT columnCount, 323 CoinFactorizationDouble * COIN_RESTRICT elementU) 324 { 325 const int * COIN_RESTRICT columnLength = matrix_>getVectorLengths(); 313 void AbcMatrix::fillBasis(const int *COIN_RESTRICT whichColumn, 314 int &numberColumnBasic, 315 int *COIN_RESTRICT indexRowU, 316 int *COIN_RESTRICT start, 317 int *COIN_RESTRICT rowCount, 318 int *COIN_RESTRICT columnCount, 319 CoinFactorizationDouble *COIN_RESTRICT elementU) 320 { 321 const int *COIN_RESTRICT columnLength = matrix_>getVectorLengths(); 326 322 CoinBigIndex numberElements = start[0]; 327 323 // fill 328 const CoinBigIndex * 329 const int * 330 const double * 331 int numberRows =model_>numberRows();324 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts(); 325 const int *COIN_RESTRICT row = matrix_>getIndices(); 326 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 327 int numberRows = model_>numberRows(); 332 328 for (int i = 0; i < numberColumnBasic; i++) { 333 int iColumn = whichColumn[i] numberRows;329 int iColumn = whichColumn[i]  numberRows; 334 330 int length = columnLength[iColumn]; 335 331 CoinBigIndex startThis = columnStart[iColumn]; … … 340 336 indexRowU[numberElements] = iRow; 341 337 rowCount[iRow]++; 342 assert 338 assert(elementByColumn[j]); 343 339 elementU[numberElements++] = elementByColumn[j]; 344 340 } 345 start[i +1] = numberElements;341 start[i + 1] = numberElements; 346 342 } 347 343 } 348 344 #ifdef ABC_LONG_FACTORIZATION 349 void 350 AbcMatrix::fillBasis(const int * COIN_RESTRICT whichColumn, 351 int & numberColumnBasic, 352 int * COIN_RESTRICT indexRowU, 353 int * COIN_RESTRICT start, 354 int * COIN_RESTRICT rowCount, 355 int * COIN_RESTRICT columnCount, 356 long double * COIN_RESTRICT elementU) 357 { 358 const int * COIN_RESTRICT columnLength = matrix_>getVectorLengths(); 345 void AbcMatrix::fillBasis(const int *COIN_RESTRICT whichColumn, 346 int &numberColumnBasic, 347 int *COIN_RESTRICT indexRowU, 348 int *COIN_RESTRICT start, 349 int *COIN_RESTRICT rowCount, 350 int *COIN_RESTRICT columnCount, 351 long double *COIN_RESTRICT elementU) 352 { 353 const int *COIN_RESTRICT columnLength = matrix_>getVectorLengths(); 359 354 CoinBigIndex numberElements = start[0]; 360 355 // fill 361 const CoinBigIndex * 362 const int * 363 const double * 364 int numberRows =model_>numberRows();356 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts(); 357 const int *COIN_RESTRICT row = matrix_>getIndices(); 358 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 359 int numberRows = model_>numberRows(); 365 360 for (int i = 0; i < numberColumnBasic; i++) { 366 int iColumn = whichColumn[i] numberRows;361 int iColumn = whichColumn[i]  numberRows; 367 362 int length = columnLength[iColumn]; 368 363 CoinBigIndex startThis = columnStart[iColumn]; … … 373 368 indexRowU[numberElements] = iRow; 374 369 rowCount[iRow]++; 375 assert 370 assert(elementByColumn[j]); 376 371 elementU[numberElements++] = elementByColumn[j]; 377 372 } 378 start[i +1] = numberElements;379 } 380 } 381 #endif 373 start[i + 1] = numberElements; 374 } 375 } 376 #endif 382 377 #if 0 383 378 /// Move largest in column to beginning … … 416 411 #endif 417 412 // Creates row copy 418 void 419 AbcMatrix::createRowCopy() 413 void AbcMatrix::createRowCopy() 420 414 { 421 415 #if ABC_PARALLEL 422 if (model_>parallelMode() ==0)423 #endif 424 numberRowBlocks_ =1;416 if (model_>parallelMode() == 0) 417 #endif 418 numberRowBlocks_ = 1; 425 419 #if ABC_PARALLEL 426 420 else 427 numberRowBlocks_ =CoinMin(NUMBER_ROW_BLOCKS,model_>numberCpus());428 #endif 429 int maximumRows =model_>maximumAbcNumberRows();430 int numberRows =model_>numberRows();431 int numberColumns =model_>numberColumns();421 numberRowBlocks_ = CoinMin(NUMBER_ROW_BLOCKS, model_>numberCpus()); 422 #endif 423 int maximumRows = model_>maximumAbcNumberRows(); 424 int numberRows = model_>numberRows(); 425 int numberColumns = model_>numberColumns(); 432 426 int numberElements = matrix_>getNumElements(); 433 assert 434 char * whichBlock_=new char[numberColumns];435 rowStart_ =new CoinBigIndex[numberRows*(numberRowBlocks_+2)];436 element_ = new double 437 column_ = new int 438 const CoinBigIndex * 439 memset(blockStart_, 0,sizeof(blockStart_));427 assert(!rowStart_); 428 char *whichBlock_ = new char[numberColumns]; 429 rowStart_ = new CoinBigIndex[numberRows * (numberRowBlocks_ + 2)]; 430 element_ = new double[numberElements]; 431 column_ = new int[numberElements]; 432 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts(); 433 memset(blockStart_, 0, sizeof(blockStart_)); 440 434 int ecount[10]; 441 assert (numberRowBlocks_<16);442 CoinAbcMemset0(ecount, 10);435 assert(numberRowBlocks_ < 16); 436 CoinAbcMemset0(ecount, 10); 443 437 // allocate to blocks (put a bit less in first as will be dealing with slacks) LATER 444 CoinBigIndex start =0;445 int block =0;446 CoinBigIndex work =(2*numberColumns+matrix_>getNumElements()+numberRowBlocks_1)/numberRowBlocks_;447 CoinBigIndex thisWork =work;448 for (int iColumn =0;iColumn<numberColumns;iColumn++) {449 CoinBigIndex end =columnStart[iColumn+1];450 assert (block>=0&&block<numberRowBlocks_);451 whichBlock_[iColumn] =static_cast<char>(block);452 thisWork = 2 +endstart;453 ecount[block] +=endstart;454 start =end;438 CoinBigIndex start = 0; 439 int block = 0; 440 CoinBigIndex work = (2 * numberColumns + matrix_>getNumElements() + numberRowBlocks_  1) / numberRowBlocks_; 441 CoinBigIndex thisWork = work; 442 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 443 CoinBigIndex end = columnStart[iColumn + 1]; 444 assert(block >= 0 && block < numberRowBlocks_); 445 whichBlock_[iColumn] = static_cast< char >(block); 446 thisWork = 2 + end  start; 447 ecount[block] += end  start; 448 start = end; 455 449 blockStart_[block]++; 456 if (thisWork <=0) {450 if (thisWork <= 0) { 457 451 block++; 458 thisWork =work;452 thisWork = work; 459 453 } 460 454 } … … 465 459 printf("\n"); 466 460 #endif 467 start =0;468 for (int i =0;i<numberRowBlocks_;i++) {469 int next =start+blockStart_[i];470 blockStart_[i] =start;471 start =next;472 } 473 blockStart_[numberRowBlocks_] =start;474 assert (start==numberColumns);475 const int * 476 const double * 461 start = 0; 462 for (int i = 0; i < numberRowBlocks_; i++) { 463 int next = start + blockStart_[i]; 464 blockStart_[i] = start; 465 start = next; 466 } 467 blockStart_[numberRowBlocks_] = start; 468 assert(start == numberColumns); 469 const int *COIN_RESTRICT row = matrix_>getIndices(); 470 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 477 471 // counts 478 CoinAbcMemset0(rowStart_, numberRows*(numberRowBlocks_+2));479 int * COIN_RESTRICT last = rowStart_+numberRows*(numberRowBlocks_+1);480 for (int iColumn =0;iColumn<numberColumns;iColumn++) {481 int block =whichBlock_[iColumn];472 CoinAbcMemset0(rowStart_, numberRows * (numberRowBlocks_ + 2)); 473 int *COIN_RESTRICT last = rowStart_ + numberRows * (numberRowBlocks_ + 1); 474 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 475 int block = whichBlock_[iColumn]; 482 476 blockStart_[block]++; 483 int base =(block+1)*numberRows;484 for (CoinBigIndex j =columnStart[iColumn];j<columnStart[iColumn+1];j++) {485 int iRow =row[j];486 rowStart_[base +iRow]++;477 int base = (block + 1) * numberRows; 478 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) { 479 int iRow = row[j]; 480 rowStart_[base + iRow]++; 487 481 last[iRow]++; 488 482 } 489 483 } 490 484 // back to real starts 491 for (int iBlock =numberRowBlocks_1;iBlock>=0;iBlock)492 blockStart_[iBlock +1]=blockStart_[iBlock];493 blockStart_[0] =0;494 CoinBigIndex put =0;495 for (int iRow =1;iRow<numberRows;iRow++) {496 int n =last[iRow1];497 last[iRow 1]=put;485 for (int iBlock = numberRowBlocks_  1; iBlock >= 0; iBlock) 486 blockStart_[iBlock + 1] = blockStart_[iBlock]; 487 blockStart_[0] = 0; 488 CoinBigIndex put = 0; 489 for (int iRow = 1; iRow < numberRows; iRow++) { 490 int n = last[iRow  1]; 491 last[iRow  1] = put; 498 492 put += n; 499 rowStart_[iRow] =put;500 } 501 int n =last[numberRows1];502 last[numberRows 1]=put;493 rowStart_[iRow] = put; 494 } 495 int n = last[numberRows  1]; 496 last[numberRows  1] = put; 503 497 put += n; 504 assert (put==numberElements);498 assert(put == numberElements); 505 499 //last[numberRows1]=put; 506 500 // starts 507 int base =0;508 for (int iBlock =0;iBlock<numberRowBlocks_;iBlock++) {509 for (int iRow =0;iRow<numberRows;iRow++) {510 rowStart_[base +numberRows+iRow]+=rowStart_[base+iRow];511 } 512 base += numberRows;513 } 514 for (int iColumn =0;iColumn<numberColumns;iColumn++) {515 int block =whichBlock_[iColumn];516 int where =blockStart_[block];501 int base = 0; 502 for (int iBlock = 0; iBlock < numberRowBlocks_; iBlock++) { 503 for (int iRow = 0; iRow < numberRows; iRow++) { 504 rowStart_[base + numberRows + iRow] += rowStart_[base + iRow]; 505 } 506 base += numberRows; 507 } 508 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 509 int block = whichBlock_[iColumn]; 510 int where = blockStart_[block]; 517 511 blockStart_[block]++; 518 int base =block*numberRows;519 for (CoinBigIndex j =columnStart[iColumn];j<columnStart[iColumn+1];j++) {520 int iRow =row[j];521 CoinBigIndex put =rowStart_[base+iRow];522 rowStart_[base +iRow]++;523 column_[put] =where+maximumRows;524 element_[put] =elementByColumn[j];512 int base = block * numberRows; 513 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) { 514 int iRow = row[j]; 515 CoinBigIndex put = rowStart_[base + iRow]; 516 rowStart_[base + iRow]++; 517 column_[put] = where + maximumRows; 518 element_[put] = elementByColumn[j]; 525 519 } 526 520 } 527 521 // back to real starts etc 528 base =numberRows*numberRowBlocks_;529 for (int iBlock =numberRowBlocks_1;iBlock>=0;iBlock){530 blockStart_[iBlock +1]=blockStart_[iBlock]+maximumRows;531 CoinAbcMemcpy(rowStart_ +base,rowStart_+basenumberRows,numberRows);522 base = numberRows * numberRowBlocks_; 523 for (int iBlock = numberRowBlocks_  1; iBlock >= 0; iBlock) { 524 blockStart_[iBlock + 1] = blockStart_[iBlock] + maximumRows; 525 CoinAbcMemcpy(rowStart_ + base, rowStart_ + base  numberRows, numberRows); 532 526 base = numberRows; 533 527 } 534 blockStart_[0] =0;//maximumRows;535 delete 536 // and move 537 CoinAbcMemcpy(rowStart_, last,numberRows);528 blockStart_[0] = 0; //maximumRows; 529 delete[] whichBlock_; 530 // and move 531 CoinAbcMemcpy(rowStart_, last, numberRows); 538 532 // All in useful 539 CoinAbcMemcpy(rowStart_ +(numberRowBlocks_+1)*numberRows,540 rowStart_+(numberRowBlocks_)*numberRows,numberRows);533 CoinAbcMemcpy(rowStart_ + (numberRowBlocks_ + 1) * numberRows, 534 rowStart_ + (numberRowBlocks_)*numberRows, numberRows); 541 535 #ifdef COUNT_COPY 542 536 // now blocked by element count 543 countRealColumn_ =new int[numberColumns];544 int counts [2*MAX_COUNT];545 memset(counts, 0,sizeof(counts));537 countRealColumn_ = new int[numberColumns]; 538 int counts[2 * MAX_COUNT]; 539 memset(counts, 0, sizeof(counts)); 546 540 //memset(countFirst_,0,sizeof(countFirst_)); 547 int numberLarge =0;548 for (int iColumn =0;iColumn<numberColumns;iColumn++) {549 int n =columnStart[iColumn+1]columnStart[iColumn];550 if (n <MAX_COUNT)541 int numberLarge = 0; 542 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 543 int n = columnStart[iColumn + 1]  columnStart[iColumn]; 544 if (n < MAX_COUNT) 551 545 counts[n]++; 552 546 else 553 547 numberLarge++; 554 548 } 555 CoinBigIndex numberExtra =3; // for alignment549 CoinBigIndex numberExtra = 3; // for alignment 556 550 #define SMALL_COUNT 1 557 for (int i =0;i<MAX_COUNT;i++) {558 int n =counts[i];559 if (n >=SMALL_COUNT) {551 for (int i = 0; i < MAX_COUNT; i++) { 552 int n = counts[i]; 553 if (n >= SMALL_COUNT) { 560 554 n &= 3; 561 int extra =(4n)&3;562 numberExtra += i*extra;555 int extra = (4  n) & 3; 556 numberExtra += i * extra; 563 557 } else { 564 558 // treat as large 565 numberLarge +=n;566 } 567 } 568 countElement_ = new double [numberElements+numberExtra];569 countRow_ =new int [numberElements+numberExtra];570 countStartLarge_ = new CoinBigIndex [numberLarge+1];571 countStartLarge_[numberLarge] =numberElements+numberExtra;559 numberLarge += n; 560 } 561 } 562 countElement_ = new double[numberElements + numberExtra]; 563 countRow_ = new int[numberElements + numberExtra]; 564 countStartLarge_ = new CoinBigIndex[numberLarge + 1]; 565 countStartLarge_[numberLarge] = numberElements + numberExtra; 572 566 //return; 573 CoinInt64 xx = reinterpret_cast< CoinInt64>(countElement_);574 int iBottom = static_cast< int>(xx & 31);575 int offset = iBottom >>3;576 CoinBigIndex firstElementLarge =0;567 CoinInt64 xx = reinterpret_cast< CoinInt64 >(countElement_); 568 int iBottom = static_cast< int >(xx & 31); 569 int offset = iBottom >> 3; 570 CoinBigIndex firstElementLarge = 0; 577 571 if (offset) 578 firstElementLarge =4offset;572 firstElementLarge = 4  offset; 579 573 //countStart_[0]=firstElementLarge; 580 int positionLarge =0;581 smallestCount_ =0;582 largestCount_ =0;583 for (int i =0;i<MAX_COUNT;i++) {584 int n =counts[i];585 countFirst_[i] =positionLarge;586 countStart_[i] =firstElementLarge;587 if (n >=SMALL_COUNT) {588 counts[i +MAX_COUNT]=1;589 if (smallestCount_ ==0)590 smallestCount_=i;591 largestCount_ =i;592 positionLarge +=n;593 firstElementLarge +=n*i;574 int positionLarge = 0; 575 smallestCount_ = 0; 576 largestCount_ = 0; 577 for (int i = 0; i < MAX_COUNT; i++) { 578 int n = counts[i]; 579 countFirst_[i] = positionLarge; 580 countStart_[i] = firstElementLarge; 581 if (n >= SMALL_COUNT) { 582 counts[i + MAX_COUNT] = 1; 583 if (smallestCount_ == 0) 584 smallestCount_ = i; 585 largestCount_ = i; 586 positionLarge += n; 587 firstElementLarge += n * i; 594 588 n &= 3; 595 int extra =(4n)&3;596 firstElementLarge += i*extra;597 } 598 counts[i] =0;589 int extra = (4  n) & 3; 590 firstElementLarge += i * extra; 591 } 592 counts[i] = 0; 599 593 } 600 594 largestCount_++; 601 countFirst_[MAX_COUNT] =positionLarge;602 countStart_[MAX_COUNT] =firstElementLarge;603 numberLarge =0;604 for (int iColumn =0;iColumn<numberColumns;iColumn++) {605 CoinBigIndex start =columnStart[iColumn];606 int n =columnStart[iColumn+1]start;595 countFirst_[MAX_COUNT] = positionLarge; 596 countStart_[MAX_COUNT] = firstElementLarge; 597 numberLarge = 0; 598 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 599 CoinBigIndex start = columnStart[iColumn]; 600 int n = columnStart[iColumn + 1]  start; 607 601 CoinBigIndex put; 608 602 int position; 609 if (n <MAX_COUNT&&counts[MAX_COUNT+n]!=0) {610 int iWhich =counts[n];611 position = countFirst_[n] +iWhich;612 int iBlock4 =iWhich&(~3);603 if (n < MAX_COUNT && counts[MAX_COUNT + n] != 0) { 604 int iWhich = counts[n]; 605 position = countFirst_[n] + iWhich; 606 int iBlock4 = iWhich & (~3); 613 607 iWhich &= 3; 614 put = countStart_[n] +iBlock4*n;608 put = countStart_[n] + iBlock4 * n; 615 609 put += iWhich; 616 610 counts[n]++; 617 for (int i =0;i<n;i++) {618 countRow_[put]=row[start+i];619 countElement_[put]=elementByColumn[start+i];620 put+=4;611 for (int i = 0; i < n; i++) { 612 countRow_[put] = row[start + i]; 613 countElement_[put] = elementByColumn[start + i]; 614 put += 4; 621 615 } 622 616 } else { 623 position = positionLarge +numberLarge;624 put =firstElementLarge;625 countStartLarge_[numberLarge] =put;626 firstElementLarge +=n;617 position = positionLarge + numberLarge; 618 put = firstElementLarge; 619 countStartLarge_[numberLarge] = put; 620 firstElementLarge += n; 627 621 numberLarge++; 628 CoinAbcMemcpy(countRow_ +put,row+start,n);629 CoinAbcMemcpy(countElement_ +put,elementByColumn+start,n);630 } 631 countRealColumn_[position] =iColumn+maximumRows;632 } 633 countStartLarge_[numberLarge] =firstElementLarge;634 assert (firstElementLarge<=numberElements+numberExtra);622 CoinAbcMemcpy(countRow_ + put, row + start, n); 623 CoinAbcMemcpy(countElement_ + put, elementByColumn + start, n); 624 } 625 countRealColumn_[position] = iColumn + maximumRows; 626 } 627 countStartLarge_[numberLarge] = firstElementLarge; 628 assert(firstElementLarge <= numberElements + numberExtra); 635 629 #endif 636 630 } 637 631 // Make all useful 638 void 639 AbcMatrix::makeAllUseful(CoinIndexedVector & /*spare*/) 640 { 641 int numberRows=model_>numberRows(); 642 CoinBigIndex * COIN_RESTRICT rowEnd = rowStart_+numberRows; 643 const CoinBigIndex * COIN_RESTRICT rowReallyEnd = rowStart_+2*numberRows; 644 for (int iRow=0;iRow<numberRows;iRow++) { 645 rowEnd[iRow]=rowReallyEnd[iRow]; 632 void AbcMatrix::makeAllUseful(CoinIndexedVector & /*spare*/) 633 { 634 int numberRows = model_>numberRows(); 635 CoinBigIndex *COIN_RESTRICT rowEnd = rowStart_ + numberRows; 636 const CoinBigIndex *COIN_RESTRICT rowReallyEnd = rowStart_ + 2 * numberRows; 637 for (int iRow = 0; iRow < numberRows; iRow++) { 638 rowEnd[iRow] = rowReallyEnd[iRow]; 646 639 } 647 640 } 648 641 #define SQRT_ARRAY 649 642 // Creates scales for column copy (rowCopy in model may be modified) 650 void 651 AbcMatrix::scale(int numberAlreadyScaled) 652 { 653 int numberRows=model_>numberRows(); 654 bool inBranchAndBound=(model_>specialOptions(),0x1000000)!=0; 655 bool doScaling=numberAlreadyScaled>=0; 643 void AbcMatrix::scale(int numberAlreadyScaled) 644 { 645 int numberRows = model_>numberRows(); 646 bool inBranchAndBound = (model_>specialOptions(), 0x1000000) != 0; 647 bool doScaling = numberAlreadyScaled >= 0; 656 648 if (!doScaling) 657 numberAlreadyScaled =0;658 if (numberAlreadyScaled ==numberRows)649 numberAlreadyScaled = 0; 650 if (numberAlreadyScaled == numberRows) 659 651 return; // no need to do anything 660 int numberColumns =model_>numberColumns();661 double * COIN_RESTRICT rowScale=model_>rowScale2();662 double * COIN_RESTRICT inverseRowScale=model_>inverseRowScale2();663 double * COIN_RESTRICT columnScale=model_>columnScale2();664 double * COIN_RESTRICT inverseColumnScale=model_>inverseColumnScale2();652 int numberColumns = model_>numberColumns(); 653 double *COIN_RESTRICT rowScale = model_>rowScale2(); 654 double *COIN_RESTRICT inverseRowScale = model_>inverseRowScale2(); 655 double *COIN_RESTRICT columnScale = model_>columnScale2(); 656 double *COIN_RESTRICT inverseColumnScale = model_>inverseColumnScale2(); 665 657 // we are going to mark bits we are interested in 666 int whichArray =model_>getAvailableArrayPublic();667 char * COIN_RESTRICT usefulColumn = reinterpret_cast<char *>(model_>usefulArray(whichArray)>getIndices());668 memset(usefulColumn, 1,numberColumns);669 const double * 670 const double * 671 const double * 672 const double * 658 int whichArray = model_>getAvailableArrayPublic(); 659 char *COIN_RESTRICT usefulColumn = reinterpret_cast< char * >(model_>usefulArray(whichArray)>getIndices()); 660 memset(usefulColumn, 1, numberColumns); 661 const double *COIN_RESTRICT rowLower = model_>rowLower(); 662 const double *COIN_RESTRICT rowUpper = model_>rowUpper(); 663 const double *COIN_RESTRICT columnLower = model_>columnLower(); 664 const double *COIN_RESTRICT columnUpper = model_>columnUpper(); 673 665 //#define LEAVE_FIXED 674 666 // mark empty and fixed columns 675 667 // get matrix data pointers 676 const int * 677 const CoinBigIndex * 678 double * 679 CoinPackedMatrix * 680 const int * 681 const CoinBigIndex * 682 const double * 683 assert (numberAlreadyScaled>=0&&numberAlreadyScaled<numberRows);668 const int *COIN_RESTRICT row = matrix_>getIndices(); 669 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts(); 670 double *COIN_RESTRICT elementByColumn = matrix_>getMutableElements(); 671 CoinPackedMatrix *COIN_RESTRICT rowCopy = reverseOrderedCopy(); 672 const int *column = rowCopy>getIndices(); 673 const CoinBigIndex *COIN_RESTRICT rowStart = rowCopy>getVectorStarts(); 674 const double *COIN_RESTRICT element = rowCopy>getElements(); 675 assert(numberAlreadyScaled >= 0 && numberAlreadyScaled < numberRows); 684 676 #ifndef LEAVE_FIXED 685 677 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 686 if (columnUpper[iColumn] == 687 usefulColumn[iColumn] =0;678 if (columnUpper[iColumn] == columnLower[iColumn]) 679 usefulColumn[iColumn] = 0; 688 680 } 689 681 #endif … … 692 684 if (!numberAlreadyScaled) { 693 685 // may be redundant 694 CoinFillN(rowScale, numberRows,1.0);695 CoinFillN(columnScale, numberColumns,1.0);686 CoinFillN(rowScale, numberRows, 1.0); 687 CoinFillN(columnScale, numberColumns, 1.0); 696 688 CoinBigIndex start = 0; 697 689 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 698 CoinBigIndex end = columnStart[iColumn +1];690 CoinBigIndex end = columnStart[iColumn + 1]; 699 691 if (usefulColumn[iColumn]) { 700 if (end>start) {701 702 703 704 705 706 707 usefulColumn[iColumn]=0;708 709 } 710 start =end;692 if (end > start) { 693 for (CoinBigIndex j = start; j < end; j++) { 694 double value = fabs(elementByColumn[j]); 695 overallLargest = CoinMax(overallLargest, value); 696 overallSmallest = CoinMin(overallSmallest, value); 697 } 698 } else { 699 usefulColumn[iColumn] = 0; 700 } 701 } 702 start = end; 711 703 } 712 704 } else { 713 705 CoinBigIndex start = rowStart[numberAlreadyScaled]; 714 706 for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) { 715 rowScale[iRow] =1.0;716 CoinBigIndex end = rowStart[iRow +1];707 rowScale[iRow] = 1.0; 708 CoinBigIndex end = rowStart[iRow + 1]; 717 709 for (CoinBigIndex j = start; j < end; j++) { 718 int iColumn=column[j];719 720 double value = fabs(elementByColumn[j])*columnScale[iColumn];721 722 723 724 } 725 } 726 } 727 if ((overallSmallest >= 0.5 && overallLargest <= 2.0) !doScaling) {710 int iColumn = column[j]; 711 if (usefulColumn[iColumn]) { 712 double value = fabs(elementByColumn[j]) * columnScale[iColumn]; 713 overallLargest = CoinMax(overallLargest, value); 714 overallSmallest = CoinMin(overallSmallest, value); 715 } 716 } 717 } 718 } 719 if ((overallSmallest >= 0.5 && overallLargest <= 2.0)  !doScaling) { 728 720 //printf("no scaling\n"); 729 721 delete rowCopy; 730 722 model_>clearArraysPublic(whichArray); 731 CoinFillN(inverseRowScale +numberAlreadyScaled,numberRowsnumberAlreadyScaled,1.0);732 if (!numberAlreadyScaled) 733 CoinFillN(inverseColumnScale, numberColumns,1.0);723 CoinFillN(inverseRowScale + numberAlreadyScaled, numberRows  numberAlreadyScaled, 1.0); 724 if (!numberAlreadyScaled) 725 CoinFillN(inverseColumnScale, numberColumns, 1.0); 734 726 //moveLargestToStart(); 735 return 727 return; 736 728 } 737 729 // need to scale 738 730 double largest; 739 731 double smallest; 740 int scalingMethod =model_>scalingFlag();732 int scalingMethod = model_>scalingFlag(); 741 733 if (scalingMethod == 4) { 742 734 // As auto … … 751 743 // if scalingMethod 3 then may change 752 744 bool extraDetails = (model_>logLevel() > 2); 753 bool secondTime =false;745 bool secondTime = false; 754 746 while (!finished) { 755 747 int numberPass = !numberAlreadyScaled ? 3 : 1; … … 757 749 overallSmallest = 1.0e20; 758 750 if (!secondTime) { 759 secondTime =true;751 secondTime = true; 760 752 } else { 761 CoinFillN (rowScale, numberRows, 1.0);762 CoinFillN (columnScale, numberColumns, 1.0);763 } 753 CoinFillN(rowScale, numberRows, 1.0); 754 CoinFillN(columnScale, numberColumns, 1.0); 755 } 764 756 if (scalingMethod == 1  scalingMethod == 3) { 765 757 // Maximum in each row 766 758 for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) { 767 768 for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {769 770 771 772 773 assert(largest < 1.0e40);774 775 776 759 largest = 1.0e10; 760 for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) { 761 int iColumn = column[j]; 762 if (usefulColumn[iColumn]) { 763 double value = fabs(element[j]); 764 largest = CoinMax(largest, value); 765 assert(largest < 1.0e40); 766 } 767 } 768 rowScale[iRow] = 1.0 / largest; 777 769 #ifdef COIN_DEVELOP 778 779 780 781 770 if (extraDetails) { 771 overallLargest = CoinMax(overallLargest, largest); 772 overallSmallest = CoinMin(overallSmallest, largest); 773 } 782 774 #endif 783 775 } 784 776 } else { 785 777 while (numberPass) { 786 787 788 789 790 791 792 793 for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {794 795 796 797 798 799 800 801 778 overallLargest = 0.0; 779 overallSmallest = 1.0e50; 780 numberPass; 781 // Geometric mean on row scales 782 for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) { 783 largest = 1.0e50; 784 smallest = 1.0e50; 785 for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) { 786 int iColumn = column[j]; 787 if (usefulColumn[iColumn]) { 788 double value = fabs(element[j]); 789 value *= columnScale[iColumn]; 790 largest = CoinMax(largest, value); 791 smallest = CoinMin(smallest, value); 792 } 793 } 802 794 #ifdef SQRT_ARRAY 803 795 rowScale[iRow] = smallest * largest; 804 796 #else 805 806 #endif 807 808 809 810 811 812 813 797 rowScale[iRow] = 1.0 / sqrt(smallest * largest); 798 #endif 799 if (extraDetails) { 800 overallLargest = CoinMax(largest * rowScale[iRow], overallLargest); 801 overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest); 802 } 803 } 804 if (model_>scalingFlag() == 5) 805 break; // just scale rows 814 806 #ifdef SQRT_ARRAY 815 816 #endif 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 j < columnStart[iColumn+1]; j++) {832 833 834 835 836 837 807 CoinAbcInverseSqrts(rowScale, numberRows); 808 #endif 809 if (!inBranchAndBound) 810 model_>messageHandler()>message(CLP_PACKEDSCALE_WHILE, *model_>messagesPointer()) 811 << overallSmallest 812 << overallLargest 813 << CoinMessageEol; 814 // skip last column round 815 if (numberPass == 1) 816 break; 817 // Geometric mean on column scales 818 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 819 if (usefulColumn[iColumn]) { 820 largest = 1.0e50; 821 smallest = 1.0e50; 822 for (CoinBigIndex j = columnStart[iColumn]; 823 j < columnStart[iColumn + 1]; j++) { 824 int iRow = row[j]; 825 double value = fabs(elementByColumn[j]); 826 value *= rowScale[iRow]; 827 largest = CoinMax(largest, value); 828 smallest = CoinMin(smallest, value); 829 } 838 830 #ifdef SQRT_ARRAY 839 831 columnScale[iColumn] = smallest * largest; 840 832 #else 841 842 #endif 843 844 833 columnScale[iColumn] = 1.0 / sqrt(smallest * largest); 834 #endif 835 } 836 } 845 837 #ifdef SQRT_ARRAY 846 838 CoinAbcInverseSqrts(columnScale, numberColumns); 847 839 #endif 848 840 } … … 853 845 double scaledDifference = difference * rowScale[iRow]; 854 846 if (scaledDifference > tolerance && scaledDifference < 1.0e4) { 855 856 857 858 859 847 // make gap larger 848 rowScale[iRow] *= 1.0e4 / scaledDifference; 849 rowScale[iRow] = CoinMax(1.0e10, CoinMin(1.0e10, rowScale[iRow])); 850 //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference, 851 // scaledDifference,difference*rowScale[iRow]); 860 852 } 861 853 } … … 865 857 overallSmallest = 1.0e50; 866 858 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 867 868 869 870 for (CoinBigIndex j = columnStart[iColumn];j < columnStart[iColumn+1]; j++) {871 872 873 874 875 876 877 878 859 if (usefulColumn[iColumn]) { 860 largest = 1.0e20; 861 smallest = 1.0e50; 862 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) { 863 int iRow = row[j]; 864 double value = fabs(elementByColumn[j] * rowScale[iRow]); 865 largest = CoinMax(largest, value); 866 smallest = CoinMin(smallest, value); 867 } 868 if (overallSmallest * largest > smallest) 869 overallSmallest = smallest / largest; 870 } 879 871 } 880 872 } … … 885 877 scalingMethod = 4; 886 878 } else { 887 assert 879 assert(scalingMethod == 4); 888 880 if (overallSmallest > 2.0 * savedOverallRatio) { 889 890 891 892 893 894 895 881 finished = true; // geometric was better 882 if (model_>scalingFlag() == 4) { 883 // If in Branch and bound change 884 if ((model_>specialOptions() & 1024) != 0) { 885 model_>scaling(2); 886 } 887 } 896 888 } else { 897 898 899 900 901 902 903 889 scalingMethod = 1; // redo equilibrium 890 if (model_>scalingFlag() == 4) { 891 // If in Branch and bound change 892 if ((model_>specialOptions() & 1024) != 0) { 893 model_>scaling(1); 894 } 895 } 904 896 } 905 897 #if 0 … … 928 920 overallLargest = CoinMin(100.0, overallLargest); 929 921 overallSmallest = 1.0e50; 930 char * COIN_RESTRICT usedRow = reinterpret_cast<char *>(inverseRowScale);922 char *COIN_RESTRICT usedRow = reinterpret_cast< char * >(inverseRowScale); 931 923 memset(usedRow, 0, numberRows); 932 924 //printf("scaling %d\n",model_>scalingFlag()); … … 934 926 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 935 927 if (usefulColumn[iColumn]) { 936 937 938 for (CoinBigIndex j = columnStart[iColumn];j < columnStart[iColumn+1]; j++) {939 940 941 942 943 944 945 946 928 largest = 1.0e20; 929 smallest = 1.0e50; 930 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) { 931 int iRow = row[j]; 932 usedRow[iRow] = 1; 933 double value = fabs(elementByColumn[j] * rowScale[iRow]); 934 largest = CoinMax(largest, value); 935 smallest = CoinMin(smallest, value); 936 } 937 columnScale[iColumn] = overallLargest / largest; 938 //columnScale[iColumn]=CoinMax(1.0e10,CoinMin(1.0e10,columnScale[iColumn])); 947 939 #ifdef RANDOMIZE 948 949 950 951 952 #endif 953 954 955 956 957 958 959 960 961 962 963 940 if (!numberAlreadyScaled) { 941 double value = 0.5  randomNumberGenerator_.randomDouble(); //between 0.5 to + 0.5 942 columnScale[iColumn] *= (1.0 + 0.1 * value); 943 } 944 #endif 945 double difference = columnUpper[iColumn]  columnLower[iColumn]; 946 if (difference < 1.0e5 * columnScale[iColumn]) { 947 // make gap larger 948 columnScale[iColumn] = difference / 1.0e5; 949 //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference, 950 // scaledDifference,difference*columnScale[iColumn]); 951 } 952 double value = smallest * columnScale[iColumn]; 953 if (overallSmallest > value) 954 overallSmallest = value; 955 //overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]); 964 956 } else { 965 966 957 assert(columnScale[iColumn] == 1.0); 958 //columnScale[iColumn]=1.0; 967 959 } 968 960 } 969 961 for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) { 970 962 if (!usedRow[iRow]) { 971 963 rowScale[iRow] = 1.0; 972 964 } 973 965 } 974 966 } 975 967 if (!inBranchAndBound) 976 model_>messageHandler()>message(CLP_PACKEDSCALE_FINAL, *model_>messagesPointer())977 << overallSmallest978 << overallLargest979 << CoinMessageEol;968 model_>messageHandler()>message(CLP_PACKEDSCALE_FINAL, *model_>messagesPointer()) 969 << overallSmallest 970 << overallLargest 971 << CoinMessageEol; 980 972 if (overallSmallest < 1.0e13) { 981 973 // Change factorization zero tolerance 982 974 double newTolerance = CoinMax(1.0e15 * (overallSmallest / 1.0e13), 983 975 1.0e18); 984 976 if (model_>factorization()>zeroTolerance() > newTolerance) 985 977 model_>factorization()>zeroTolerance(newTolerance); … … 987 979 model_>setZeroTolerance(newTolerance); 988 980 #ifndef NDEBUG 989 assert (newTolerance<0.0); // just so we can fix981 assert(newTolerance < 0.0); // just so we can fix 990 982 #endif 991 983 } 992 984 // make copy (could do faster by using previous values) 993 985 // could just do partial 994 CoinAbcReciprocal(inverseRowScale +numberAlreadyScaled,numberRowsnumberAlreadyScaled,995 rowScale+numberAlreadyScaled);996 if (!numberAlreadyScaled) 997 CoinAbcReciprocal(inverseColumnScale, numberColumns,columnScale);986 CoinAbcReciprocal(inverseRowScale + numberAlreadyScaled, numberRows  numberAlreadyScaled, 987 rowScale + numberAlreadyScaled); 988 if (!numberAlreadyScaled) 989 CoinAbcReciprocal(inverseColumnScale, numberColumns, columnScale); 998 990 // Do scaled copy //NO and move largest to start 999 991 CoinBigIndex start = 0; 1000 992 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 1001 double scale =columnScale[iColumn];1002 CoinBigIndex end = columnStart[iColumn +1];993 double scale = columnScale[iColumn]; 994 CoinBigIndex end = columnStart[iColumn + 1]; 1003 995 for (CoinBigIndex j = start; j < end; j++) { 1004 double value =elementByColumn[j];996 double value = elementByColumn[j]; 1005 997 int iRow = row[j]; 1006 if (iRow >=numberAlreadyScaled) {1007 value*= scale*rowScale[iRow];1008 elementByColumn[j]= value;1009 } 1010 } 1011 start =end;998 if (iRow >= numberAlreadyScaled) { 999 value *= scale * rowScale[iRow]; 1000 elementByColumn[j] = value; 1001 } 1002 } 1003 start = end; 1012 1004 } 1013 1005 delete rowCopy; … … 1041 1033 @pre <code>y</code> must be of size <code>numRows()</code> */ 1042 1034 //scaled versions 1043 void 1044 AbcMatrix::timesModifyExcludingSlacks(double scalar, 1045 const double * x, double * y) const 1035 void AbcMatrix::timesModifyExcludingSlacks(double scalar, 1036 const double *x, double *y) const 1046 1037 { 1047 1038 int numberTotal = model_>numberTotal(); 1048 1039 int maximumRows = model_>maximumAbcNumberRows(); 1049 1040 // get matrix data pointers 1050 const int * 1051 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;1052 const double * 1041 const int *COIN_RESTRICT row = matrix_>getIndices(); 1042 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 1043 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 1053 1044 for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) { 1054 1045 double value = x[iColumn]; 1055 1046 if (value) { 1056 1047 CoinBigIndex start = columnStart[iColumn]; 1057 CoinBigIndex end = columnStart[iColumn +1];1048 CoinBigIndex end = columnStart[iColumn + 1]; 1058 1049 value *= scalar; 1059 1050 for (CoinBigIndex j = start; j < end; j++) { 1060 1061 1051 int iRow = row[j]; 1052 y[iRow] += value * elementByColumn[j]; 1062 1053 } 1063 1054 } … … 1067 1058 @pre <code>x</code> must be of size <code>numColumns()+numRows()</code> 1068 1059 @pre <code>y</code> must be of size <code>numRows()</code> */ 1069 void 1070 AbcMatrix::timesModifyIncludingSlacks(double scalar, 1071 const double * x, double * y) const 1072 { 1073 int numberRows=model_>numberRows(); 1060 void AbcMatrix::timesModifyIncludingSlacks(double scalar, 1061 const double *x, double *y) const 1062 { 1063 int numberRows = model_>numberRows(); 1074 1064 int numberTotal = model_>numberTotal(); 1075 1065 int maximumRows = model_>maximumAbcNumberRows(); 1076 1066 // makes no sense for x==y?? 1077 assert (x!=y);1067 assert(x != y); 1078 1068 // For now just by column and assumes already scaled (use reallyScale) 1079 1069 // get matrix data pointers 1080 const int * 1081 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;1082 const double * 1083 if (scalar ==1.0) {1070 const int *COIN_RESTRICT row = matrix_>getIndices(); 1071 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 1072 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 1073 if (scalar == 1.0) { 1084 1074 // add 1085 if (x !=y) {1086 for (int i =0;i<numberRows;i++)1087 y[i]+=x[i];1075 if (x != y) { 1076 for (int i = 0; i < numberRows; i++) 1077 y[i] += x[i]; 1088 1078 } else { 1089 for (int i =0;i<numberRows;i++)1090 y[i]+=y[i];1079 for (int i = 0; i < numberRows; i++) 1080 y[i] += y[i]; 1091 1081 } 1092 1082 for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) { 1093 1083 double value = x[iColumn]; 1094 1084 if (value) { 1095 1096 CoinBigIndex next = columnStart[iColumn+1];1097 1098 1099 y[jRow]+=value*elementByColumn[j];1100 1085 CoinBigIndex start = columnStart[iColumn]; 1086 CoinBigIndex next = columnStart[iColumn + 1]; 1087 for (CoinBigIndex j = start; j < next; j++) { 1088 int jRow = row[j]; 1089 y[jRow] += value * elementByColumn[j]; 1090 } 1101 1091 } 1102 1092 } 1103 1093 } else { 1104 1094 // subtract 1105 if (x !=y) {1106 for (int i =0;i<numberRows;i++)1107 y[i]=x[i];1095 if (x != y) { 1096 for (int i = 0; i < numberRows; i++) 1097 y[i] = x[i]; 1108 1098 } else { 1109 for (int i =0;i<numberRows;i++)1110 y[i]=0.0;1099 for (int i = 0; i < numberRows; i++) 1100 y[i] = 0.0; 1111 1101 } 1112 1102 for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) { 1113 1103 double value = x[iColumn]; 1114 1104 if (value) { 1115 1116 CoinBigIndex next = columnStart[iColumn+1];1117 1118 1119 y[jRow]=value*elementByColumn[j];1120 1105 CoinBigIndex start = columnStart[iColumn]; 1106 CoinBigIndex next = columnStart[iColumn + 1]; 1107 for (CoinBigIndex j = start; j < next; j++) { 1108 int jRow = row[j]; 1109 y[jRow] = value * elementByColumn[j]; 1110 } 1121 1111 } 1122 1112 } … … 1126 1116 @pre <code>x</code> must be of size <code>numColumns()+numRows()</code> 1127 1117 @pre <code>y</code> must be of size <code>numRows()</code> */ 1128 void 1129 AbcMatrix::timesIncludingSlacks(double scalar, 1130 const double * x, double * y) const 1131 { 1132 int numberRows=model_>numberRows(); 1118 void AbcMatrix::timesIncludingSlacks(double scalar, 1119 const double *x, double *y) const 1120 { 1121 int numberRows = model_>numberRows(); 1133 1122 int numberTotal = model_>numberTotal(); 1134 1123 int maximumRows = model_>maximumAbcNumberRows(); 1135 1124 // For now just by column and assumes already scaled (use reallyScale) 1136 1125 // get matrix data pointers 1137 const int * 1138 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;1139 const double * 1140 if (scalar ==1.0) {1126 const int *COIN_RESTRICT row = matrix_>getIndices(); 1127 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 1128 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 1129 if (scalar == 1.0) { 1141 1130 // add 1142 if (x !=y) {1143 for (int i =0;i<numberRows;i++)1144 y[i]=x[i];1145 } 1131 if (x != y) { 1132 for (int i = 0; i < numberRows; i++) 1133 y[i] = x[i]; 1134 } 1146 1135 for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) { 1147 1136 double value = x[iColumn]; 1148 1137 if (value) { 1149 1150 CoinBigIndex next = columnStart[iColumn+1];1151 1152 1153 y[jRow]+=value*elementByColumn[j];1154 1138 CoinBigIndex start = columnStart[iColumn]; 1139 CoinBigIndex next = columnStart[iColumn + 1]; 1140 for (CoinBigIndex j = start; j < next; j++) { 1141 int jRow = row[j]; 1142 y[jRow] += value * elementByColumn[j]; 1143 } 1155 1144 } 1156 1145 } 1157 1146 } else { 1158 1147 // subtract 1159 if (x !=y) {1160 for (int i =0;i<numberRows;i++)1161 y[i]=x[i];1148 if (x != y) { 1149 for (int i = 0; i < numberRows; i++) 1150 y[i] = x[i]; 1162 1151 } else { 1163 for (int i =0;i<numberRows;i++)1164 y[i]=y[i];1152 for (int i = 0; i < numberRows; i++) 1153 y[i] = y[i]; 1165 1154 } 1166 1155 for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) { 1167 1156 double value = x[iColumn]; 1168 1157 if (value) { 1169 1170 CoinBigIndex next = columnStart[iColumn+1];1171 1172 1173 y[jRow]=value*elementByColumn[j];1174 1175 } 1176 } 1177 } 1178 } 1179 extern double parallelDual4(AbcSimplexDual * 1180 static double firstPass(AbcSimplex * model,int iBlock,1181 double upperTheta, int &freeSequence,1182 CoinPartitionedVector &tableauRow,1183 CoinPartitionedVector & candidateList) 1184 { 1185 int * 1186 double * 1187 double * COIN_RESTRICT arrayCandidate=candidateList.denseVector();1188 int * 1189 const double * 1190 const unsigned char * 1158 CoinBigIndex start = columnStart[iColumn]; 1159 CoinBigIndex next = columnStart[iColumn + 1]; 1160 for (CoinBigIndex j = start; j < next; j++) { 1161 int jRow = row[j]; 1162 y[jRow] = value * elementByColumn[j]; 1163 } 1164 } 1165 } 1166 } 1167 } 1168 extern double parallelDual4(AbcSimplexDual *); 1169 static double firstPass(AbcSimplex *model, int iBlock, 1170 double upperTheta, int &freeSequence, 1171 CoinPartitionedVector &tableauRow, 1172 CoinPartitionedVector &candidateList) 1173 { 1174 int *COIN_RESTRICT index = tableauRow.getIndices(); 1175 double *COIN_RESTRICT array = tableauRow.denseVector(); 1176 double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); 1177 int *COIN_RESTRICT indexCandidate = candidateList.getIndices(); 1178 const double *COIN_RESTRICT abcDj = model>djRegion(); 1179 const unsigned char *COIN_RESTRICT internalStatus = model>internalStatus(); 1191 1180 // do first pass to get possibles 1192 1181 double bestPossible = 0.0; … … 1194 1183 double tentativeTheta = 1.0e25; // try with this much smaller as guess 1195 1184 double acceptablePivot = model>currentAcceptablePivot(); 1196 double dualT =model>currentDualTolerance();1185 double dualT = model>currentDualTolerance(); 1197 1186 // fixed will have been taken out by now 1198 const double multiplier[] = { 1.0, 1.0 };1199 freeSequence =1;1200 int firstIn =model>abcMatrix()>blockStart(iBlock);1201 int numberNonZero =tableauRow.getNumElements(iBlock)+firstIn;1202 int numberRemaining =firstIn;1187 const double multiplier[] = { 1.0, 1.0 }; 1188 freeSequence = 1; 1189 int firstIn = model>abcMatrix()>blockStart(iBlock); 1190 int numberNonZero = tableauRow.getNumElements(iBlock) + firstIn; 1191 int numberRemaining = firstIn; 1203 1192 //first=tableauRow.getNumElements(); 1204 1193 // could pass in last and if numberNonZero==lastfirstIn scan as well … … 1206 1195 for (int i = firstIn; i < numberNonZero; i++) { 1207 1196 int iSequence = index[i]; 1208 double tableauValue =array[i];1209 unsigned char iStatus =internalStatus[iSequence]&7;1197 double tableauValue = array[i]; 1198 unsigned char iStatus = internalStatus[iSequence] & 7; 1210 1199 double mult = multiplier[iStatus]; 1211 1200 double alpha = tableauValue * mult; … … 1213 1202 double value = oldValue  tentativeTheta * alpha; 1214 1203 if (value < dualT) { 1215 1216 1217 1218 1219 1220 1221 1222 1204 bestPossible = CoinMax(bestPossible, alpha); 1205 value = oldValue  upperTheta * alpha; 1206 if (value < dualT && alpha >= acceptablePivot) { 1207 upperTheta = (oldValue  dualT) / alpha; 1208 } 1209 // add to list 1210 arrayCandidate[numberRemaining] = alpha; 1211 indexCandidate[numberRemaining++] = iSequence; 1223 1212 } 1224 1213 } … … 1226 1215 double badFree = 0.0; 1227 1216 double freeAlpha = model>currentAcceptablePivot(); 1228 int freeSequenceIn =model>freeSequenceIn();1217 int freeSequenceIn = model>freeSequenceIn(); 1229 1218 double currentDualTolerance = model>currentDualTolerance(); 1230 1219 for (int i = firstIn; i < numberNonZero; i++) { 1231 1220 int iSequence = index[i]; 1232 double tableauValue =array[i];1233 unsigned char iStatus =internalStatus[iSequence]&7;1234 if ((iStatus &4)==0) {1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1221 double tableauValue = array[i]; 1222 unsigned char iStatus = internalStatus[iSequence] & 7; 1223 if ((iStatus & 4) == 0) { 1224 double mult = multiplier[iStatus]; 1225 double alpha = tableauValue * mult; 1226 double oldValue = abcDj[iSequence] * mult; 1227 double value = oldValue  tentativeTheta * alpha; 1228 if (value < dualT) { 1229 bestPossible = CoinMax(bestPossible, alpha); 1230 value = oldValue  upperTheta * alpha; 1231 if (value < dualT && alpha >= acceptablePivot) { 1232 upperTheta = (oldValue  dualT) / alpha; 1233 } 1234 // add to list 1235 arrayCandidate[numberRemaining] = alpha; 1236 indexCandidate[numberRemaining++] = iSequence; 1237 } 1249 1238 } else { 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1239 bool keep; 1240 bestPossible = CoinMax(bestPossible, fabs(tableauValue)); 1241 double oldValue = abcDj[iSequence]; 1242 // If free has to be very large  should come in via dualRow 1243 //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e3) 1244 //break; 1245 if (oldValue > currentDualTolerance) { 1246 keep = true; 1247 } else if (oldValue < currentDualTolerance) { 1248 keep = true; 1249 } else { 1250 if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e5)) { 1251 keep = true; 1252 } else { 1253 keep = false; 1254 badFree = CoinMax(badFree, fabs(tableauValue)); 1255 } 1256 } 1257 if (keep) { 1269 1258 #ifdef PAN 1270 if (model>fakeSuperBasic(iSequence)>=0) {1271 #endif 1272 if (iSequence==freeSequenceIn)1273 tableauValue=COIN_DBL_MAX;1274 1275 1276 1277 1278 1259 if (model>fakeSuperBasic(iSequence) >= 0) { 1260 #endif 1261 if (iSequence == freeSequenceIn) 1262 tableauValue = COIN_DBL_MAX; 1263 // free  choose largest 1264 if (fabs(tableauValue) > fabs(freeAlpha)) { 1265 freeAlpha = tableauValue; 1266 freeSequence = iSequence; 1267 } 1279 1268 #ifdef PAN 1280 1281 #endif 1282 1269 } 1270 #endif 1271 } 1283 1272 } 1284 1273 } … … 1286 1275 //firstInX=numberNonZerofirstIn; 1287 1276 //lastInX=1;//numberRemaininglastInX; 1288 tableauRow.setNumElementsPartition(iBlock, numberNonZerofirstIn);1289 candidateList.setNumElementsPartition(iBlock, numberRemainingfirstIn);1277 tableauRow.setNumElementsPartition(iBlock, numberNonZero  firstIn); 1278 candidateList.setNumElementsPartition(iBlock, numberRemaining  firstIn); 1290 1279 return upperTheta; 1291 1280 } 1292 1281 // gets sorted tableau row and a possible value of theta 1293 double 1294 AbcMatrix::dualColumn1Row(int iBlock, double upperTheta, int & 1295 const CoinIndexedVector &update,1296 CoinPartitionedVector &tableauRow,1297 CoinPartitionedVector &candidateList) const1282 double 1283 AbcMatrix::dualColumn1Row(int iBlock, double upperTheta, int &freeSequence, 1284 const CoinIndexedVector &update, 1285 CoinPartitionedVector &tableauRow, 1286 CoinPartitionedVector &candidateList) const 1298 1287 { 1299 1288 int maximumRows = model_>maximumAbcNumberRows(); 1300 int number =update.getNumElements();1301 const double * COIN_RESTRICT pi=update.denseVector();1302 const int * 1303 const CoinBigIndex * 1304 int numberRows =model_>numberRows();1305 const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;1289 int number = update.getNumElements(); 1290 const double *COIN_RESTRICT pi = update.denseVector(); 1291 const int *COIN_RESTRICT piIndex = update.getIndices(); 1292 const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_; 1293 int numberRows = model_>numberRows(); 1294 const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_; 1306 1295 // count down 1307 1296 int nColumns; 1308 int firstIn =blockStart_[iBlock];1309 int first =firstIn;1297 int firstIn = blockStart_[iBlock]; 1298 int first = firstIn; 1310 1299 if (!first) 1311 first =maximumRows;1312 int last =blockStart_[iBlock+1];1313 nColumns =lastfirst;1314 int target =nColumns;1315 rowStart += iBlock *numberRows;1316 rowEnd = rowStart +numberRows;1317 for (int i =0;i<number;i++) {1318 int iRow =piIndex[i];1300 first = maximumRows; 1301 int last = blockStart_[iBlock + 1]; 1302 nColumns = last  first; 1303 int target = nColumns; 1304 rowStart += iBlock * numberRows; 1305 rowEnd = rowStart + numberRows; 1306 for (int i = 0; i < number; i++) { 1307 int iRow = piIndex[i]; 1319 1308 CoinBigIndex end = rowEnd[iRow]; 1320 target = end rowStart[iRow];1321 if (target <0)1309 target = end  rowStart[iRow]; 1310 if (target < 0) 1322 1311 break; 1323 1312 } 1324 if (target >0) {1313 if (target > 0) { 1325 1314 //printf("going to few %d ops %d\n",number,nColumnstarget); 1326 1315 return dualColumn1RowFew(iBlock, upperTheta, freeSequence, 1327 1328 } 1329 const unsigned char * 1330 int * 1331 double * 1332 const double * 1333 const int * 1334 for (int i =0;i<number;i++) {1335 int iRow =piIndex[i];1316 update, tableauRow, candidateList); 1317 } 1318 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 1319 int *COIN_RESTRICT index = tableauRow.getIndices(); 1320 double *COIN_RESTRICT array = tableauRow.denseVector(); 1321 const double *COIN_RESTRICT element = element_; 1322 const int *COIN_RESTRICT column = column_; 1323 for (int i = 0; i < number; i++) { 1324 int iRow = piIndex[i]; 1336 1325 double piValue = pi[iRow]; 1337 1326 CoinBigIndex end = rowEnd[iRow]; 1338 for (CoinBigIndex j =rowStart[iRow];j<end;j++) {1327 for (CoinBigIndex j = rowStart[iRow]; j < end; j++) { 1339 1328 int iColumn = column[j]; 1340 assert (iColumn>=first&&iColumn<last);1329 assert(iColumn >= first && iColumn < last); 1341 1330 double value = element[j]; 1342 array[iColumn] += piValue*value;1331 array[iColumn] += piValue * value; 1343 1332 } 1344 1333 } 1345 1334 int numberNonZero; 1346 1335 int numberRemaining; 1347 if (iBlock ==0) {1336 if (iBlock == 0) { 1348 1337 #if 1 1349 upperTheta =static_cast<AbcSimplexDual *>(model_)>dualColumn1A();1350 numberNonZero =tableauRow.getNumElements(0);1351 numberRemaining =candidateList.getNumElements(0);1338 upperTheta = static_cast< AbcSimplexDual * >(model_)>dualColumn1A(); 1339 numberNonZero = tableauRow.getNumElements(0); 1340 numberRemaining = candidateList.getNumElements(0); 1352 1341 #else 1353 numberNonZero =0;1354 for (int i =0;i<number;i++) {1355 int iRow =piIndex[i];1356 unsigned char type =internalStatus[iRow];1357 if ((type &7)<6) {1358 index[numberNonZero]=iRow;1359 double piValue=pi[iRow];1360 array[numberNonZero++]=piValue;1361 } 1362 } 1363 numberRemaining =0;1342 numberNonZero = 0; 1343 for (int i = 0; i < number; i++) { 1344 int iRow = piIndex[i]; 1345 unsigned char type = internalStatus[iRow]; 1346 if ((type & 7) < 6) { 1347 index[numberNonZero] = iRow; 1348 double piValue = pi[iRow]; 1349 array[numberNonZero++] = piValue; 1350 } 1351 } 1352 numberRemaining = 0; 1364 1353 #endif 1365 1354 } else { 1366 numberNonZero =firstIn;1367 numberRemaining =firstIn;1368 } 1369 double * COIN_RESTRICT arrayCandidate=candidateList.denseVector();1370 int * 1355 numberNonZero = firstIn; 1356 numberRemaining = firstIn; 1357 } 1358 double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); 1359 int *COIN_RESTRICT indexCandidate = candidateList.getIndices(); 1371 1360 //printf("first %d last %d firstIn %d lastIn %d\n", 1372 1361 // first,last,firstIn,lastIn); 1373 const double * 1362 const double *COIN_RESTRICT abcDj = model_>djRegion(); 1374 1363 // do first pass to get possibles 1375 1364 double bestPossible = 0.0; … … 1377 1366 double tentativeTheta = 1.0e25; // try with this much smaller as guess 1378 1367 double acceptablePivot = model_>currentAcceptablePivot(); 1379 double dualT =model_>currentDualTolerance();1380 const double multiplier[] = { 1.0, 1.0 };1368 double dualT = model_>currentDualTolerance(); 1369 const double multiplier[] = { 1.0, 1.0 }; 1381 1370 double zeroTolerance = model_>zeroTolerance(); 1382 freeSequence =1;1371 freeSequence = 1; 1383 1372 if (model_>ordinaryVariables()) { 1384 1373 for (int iSequence = first; iSequence < last; iSequence++) { 1385 double tableauValue =array[iSequence];1374 double tableauValue = array[iSequence]; 1386 1375 if (tableauValue) { 1387 array[iSequence]=0.0;1388 if (fabs(tableauValue)>zeroTolerance) {1389 unsigned char iStatus=internalStatus[iSequence]&7;1390 if (iStatus<4) {1391 index[numberNonZero]=iSequence;1392 array[numberNonZero++]=tableauValue;1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1376 array[iSequence] = 0.0; 1377 if (fabs(tableauValue) > zeroTolerance) { 1378 unsigned char iStatus = internalStatus[iSequence] & 7; 1379 if (iStatus < 4) { 1380 index[numberNonZero] = iSequence; 1381 array[numberNonZero++] = tableauValue; 1382 double mult = multiplier[iStatus]; 1383 double alpha = tableauValue * mult; 1384 double oldValue = abcDj[iSequence] * mult; 1385 double value = oldValue  tentativeTheta * alpha; 1386 if (value < dualT) { 1387 bestPossible = CoinMax(bestPossible, alpha); 1388 value = oldValue  upperTheta * alpha; 1389 if (value < dualT && alpha >= acceptablePivot) { 1390 upperTheta = (oldValue  dualT) / alpha; 1391 } 1392 // add to list 1393 arrayCandidate[numberRemaining] = alpha; 1394 indexCandidate[numberRemaining++] = iSequence; 1395 } 1396 } 1397 } 1409 1398 } 1410 1399 } … … 1412 1401 double badFree = 0.0; 1413 1402 double freeAlpha = model_>currentAcceptablePivot(); 1414 int freeSequenceIn =model_>freeSequenceIn();1403 int freeSequenceIn = model_>freeSequenceIn(); 1415 1404 //printf("block %d freeSequence %d acceptable %g\n",iBlock,freeSequenceIn,freeAlpha); 1416 1405 double currentDualTolerance = model_>currentDualTolerance(); 1417 1406 for (int iSequence = first; iSequence < last; iSequence++) { 1418 double tableauValue =array[iSequence];1407 double tableauValue = array[iSequence]; 1419 1408 if (tableauValue) { 1420 array[iSequence]=0.0;1421 if (fabs(tableauValue)>zeroTolerance) {1422 unsigned char iStatus=internalStatus[iSequence]&7;1423 if (iStatus<6) {1424 if ((iStatus&4)==0) {1425 index[numberNonZero]=iSequence;1426 array[numberNonZero++]=tableauValue;1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 index[numberNonZero]=iSequence;1444 array[numberNonZero++]=tableauValue;1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1409 array[iSequence] = 0.0; 1410 if (fabs(tableauValue) > zeroTolerance) { 1411 unsigned char iStatus = internalStatus[iSequence] & 7; 1412 if (iStatus < 6) { 1413 if ((iStatus & 4) == 0) { 1414 index[numberNonZero] = iSequence; 1415 array[numberNonZero++] = tableauValue; 1416 double mult = multiplier[iStatus]; 1417 double alpha = tableauValue * mult; 1418 double oldValue = abcDj[iSequence] * mult; 1419 double value = oldValue  tentativeTheta * alpha; 1420 if (value < dualT) { 1421 bestPossible = CoinMax(bestPossible, alpha); 1422 value = oldValue  upperTheta * alpha; 1423 if (value < dualT && alpha >= acceptablePivot) { 1424 upperTheta = (oldValue  dualT) / alpha; 1425 } 1426 // add to list 1427 arrayCandidate[numberRemaining] = alpha; 1428 indexCandidate[numberRemaining++] = iSequence; 1429 } 1430 } else { 1431 bool keep; 1432 index[numberNonZero] = iSequence; 1433 array[numberNonZero++] = tableauValue; 1434 bestPossible = CoinMax(bestPossible, fabs(tableauValue)); 1435 double oldValue = abcDj[iSequence]; 1436 // If free has to be very large  should come in via dualRow 1437 //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e3) 1438 //break; 1439 // may be fake super basic 1440 if (oldValue > currentDualTolerance) { 1441 keep = true; 1442 } else if (oldValue < currentDualTolerance) { 1443 keep = true; 1444 } else { 1445 if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e5)) { 1446 keep = true; 1447 } else { 1448 keep = false; 1449 badFree = CoinMax(badFree, fabs(tableauValue)); 1450 } 1451 } 1463 1452 #if 0 1464 1453 if (iSequence==freeSequenceIn) 1465 1454 assert (keep); 1466 1455 #endif 1467 1456 if (keep) { 1468 1457 #ifdef PAN 1469 if (model_>fakeSuperBasic(iSequence)>=0) {1470 #endif 1471 if (iSequence==freeSequenceIn)1472 tableauValue=COIN_DBL_MAX;1473 1474 1475 1476 1477 1458 if (model_>fakeSuperBasic(iSequence) >= 0) { 1459 #endif 1460 if (iSequence == freeSequenceIn) 1461 tableauValue = COIN_DBL_MAX; 1462 // free  choose largest 1463 if (fabs(tableauValue) > fabs(freeAlpha)) { 1464 freeAlpha = tableauValue; 1465 freeSequence = iSequence; 1466 } 1478 1467 #ifdef PAN 1479 1480 #endif 1481 1482 1483 1484 1468 } 1469 #endif 1470 } 1471 } 1472 } 1473 } 1485 1474 } 1486 1475 } … … 1496 1485 xxInfo[4][iBlock]=numberRemainingfirstIn; 1497 1486 #endif 1498 tableauRow.setNumElementsPartition(iBlock, numberNonZerofirstIn);1499 candidateList.setNumElementsPartition(iBlock, numberRemainingfirstIn);1487 tableauRow.setNumElementsPartition(iBlock, numberNonZero  firstIn); 1488 candidateList.setNumElementsPartition(iBlock, numberRemaining  firstIn); 1500 1489 return upperTheta; 1501 1490 } 1502 1491 // gets sorted tableau row and a possible value of theta 1503 double 1504 AbcMatrix::dualColumn1Row2(double upperTheta, int & 1505 const CoinIndexedVector &update,1506 CoinPartitionedVector &tableauRow,1507 CoinPartitionedVector &candidateList) const1492 double 1493 AbcMatrix::dualColumn1Row2(double upperTheta, int &freeSequence, 1494 const CoinIndexedVector &update, 1495 CoinPartitionedVector &tableauRow, 1496 CoinPartitionedVector &candidateList) const 1508 1497 { 1509 1498 //int first=model_>maximumAbcNumberRows(); 1510 assert(update.getNumElements() ==2);1511 const double * COIN_RESTRICT pi=update.denseVector();1512 const int * 1513 int * 1514 double * 1515 const CoinBigIndex * 1516 int numberRows =model_>numberRows();1517 const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;1518 const double * 1519 const int * 1499 assert(update.getNumElements() == 2); 1500 const double *COIN_RESTRICT pi = update.denseVector(); 1501 const int *COIN_RESTRICT piIndex = update.getIndices(); 1502 int *COIN_RESTRICT index = tableauRow.getIndices(); 1503 double *COIN_RESTRICT array = tableauRow.denseVector(); 1504 const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_; 1505 int numberRows = model_>numberRows(); 1506 const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_; 1507 const double *COIN_RESTRICT element = element_; 1508 const int *COIN_RESTRICT column = column_; 1520 1509 int iRow0 = piIndex[0]; 1521 1510 int iRow1 = piIndex[1]; 1522 1511 CoinBigIndex end0 = rowEnd[iRow0]; 1523 1512 CoinBigIndex end1 = rowEnd[iRow1]; 1524 if (end0 rowStart[iRow0]>end1rowStart[iRow1]) {1525 int temp =iRow0;1526 iRow0 =iRow1;1527 iRow1 =temp;1513 if (end0  rowStart[iRow0] > end1  rowStart[iRow1]) { 1514 int temp = iRow0; 1515 iRow0 = iRow1; 1516 iRow1 = temp; 1528 1517 } 1529 1518 CoinBigIndex start = rowStart[iRow0]; 1530 1519 CoinBigIndex end = rowEnd[iRow0]; 1531 1520 double piValue = pi[iRow0]; 1532 double * COIN_RESTRICT arrayCandidate=candidateList.denseVector();1521 double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); 1533 1522 int numberNonZero; 1534 numberNonZero =tableauRow.getNumElements(0);1535 int n =numberNonZero;1536 for (CoinBigIndex j =start;j<end;j++) {1523 numberNonZero = tableauRow.getNumElements(0); 1524 int n = numberNonZero; 1525 for (CoinBigIndex j = start; j < end; j++) { 1537 1526 int iSequence = column[j]; 1538 1527 double value = element[j]; 1539 arrayCandidate[iSequence] = piValue*value;1540 index[n++] =iSequence;1528 arrayCandidate[iSequence] = piValue * value; 1529 index[n++] = iSequence; 1541 1530 } 1542 1531 start = rowStart[iRow1]; 1543 1532 end = rowEnd[iRow1]; 1544 1533 piValue = pi[iRow1]; 1545 for (CoinBigIndex j =start;j<end;j++) {1534 for (CoinBigIndex j = start; j < end; j++) { 1546 1535 int iSequence = column[j]; 1547 1536 double value = element[j]; 1548 1537 value *= piValue; 1549 1538 if (!arrayCandidate[iSequence]) { 1550 arrayCandidate[iSequence] = value;1551 index[n++] =iSequence;1552 } else { 1553 arrayCandidate[iSequence] += value; 1539 arrayCandidate[iSequence] = value; 1540 index[n++] = iSequence; 1541 } else { 1542 arrayCandidate[iSequence] += value; 1554 1543 } 1555 1544 } 1556 1545 // pack down 1557 const unsigned char * 1546 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 1558 1547 double zeroTolerance = model_>zeroTolerance(); 1559 while (numberNonZero <n) {1560 int iSequence =index[numberNonZero];1561 double value =arrayCandidate[iSequence];1562 arrayCandidate[iSequence] =0.0;1563 unsigned char iStatus =internalStatus[iSequence]&7;1564 if (fabs(value) >zeroTolerance&&iStatus<6) {1565 index[numberNonZero] =iSequence;1566 array[numberNonZero++] =value;1548 while (numberNonZero < n) { 1549 int iSequence = index[numberNonZero]; 1550 double value = arrayCandidate[iSequence]; 1551 arrayCandidate[iSequence] = 0.0; 1552 unsigned char iStatus = internalStatus[iSequence] & 7; 1553 if (fabs(value) > zeroTolerance && iStatus < 6) { 1554 index[numberNonZero] = iSequence; 1555 array[numberNonZero++] = value; 1567 1556 } else { 1568 1557 // kill 1569 1558 n; 1570 index[numberNonZero] =index[n];1571 } 1572 } 1573 tableauRow.setNumElementsPartition(0, numberNonZero);1574 return firstPass(model_, 0, upperTheta, freeSequence,1575 1576 candidateList);1559 index[numberNonZero] = index[n]; 1560 } 1561 } 1562 tableauRow.setNumElementsPartition(0, numberNonZero); 1563 return firstPass(model_, 0, upperTheta, freeSequence, 1564 tableauRow, 1565 candidateList); 1577 1566 } 1578 1567 //static int ixxxxxx=1; 1579 1568 // gets sorted tableau row and a possible value of theta 1580 double 1581 AbcMatrix::dualColumn1RowFew(int iBlock, double upperTheta, int & 1582 const CoinIndexedVector &update,1583 CoinPartitionedVector &tableauRow,1584 CoinPartitionedVector &candidateList) const1569 double 1570 AbcMatrix::dualColumn1RowFew(int iBlock, double upperTheta, int &freeSequence, 1571 const CoinIndexedVector &update, 1572 CoinPartitionedVector &tableauRow, 1573 CoinPartitionedVector &candidateList) const 1585 1574 { 1586 1575 //int first=model_>maximumAbcNumberRows(); 1587 1576 int number = update.getNumElements(); 1588 const double * COIN_RESTRICT pi=update.denseVector();1589 const int * 1590 int * 1591 double * 1592 const CoinBigIndex * 1593 int numberRows =model_>numberRows();1594 const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;1595 const double * 1596 const int * 1597 double * COIN_RESTRICT arrayCandidate=candidateList.denseVector();1577 const double *COIN_RESTRICT pi = update.denseVector(); 1578 const int *COIN_RESTRICT piIndex = update.getIndices(); 1579 int *COIN_RESTRICT index = tableauRow.getIndices(); 1580 double *COIN_RESTRICT array = tableauRow.denseVector(); 1581 const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_; 1582 int numberRows = model_>numberRows(); 1583 const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_; 1584 const double *COIN_RESTRICT element = element_; 1585 const int *COIN_RESTRICT column = column_; 1586 double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); 1598 1587 int numberNonZero; 1599 assert (iBlock>=0);1600 const unsigned char * 1601 int firstIn =blockStart_[iBlock];1602 if (iBlock ==0) {1603 numberNonZero =0;1604 for (int i =0;i<number;i++) {1605 int iRow =piIndex[i];1606 unsigned char type =internalStatus[iRow];1607 if ((type &7)<6) {1608 index[numberNonZero]=iRow;1609 double piValue=pi[iRow];1610 array[numberNonZero++]=piValue;1588 assert(iBlock >= 0); 1589 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 1590 int firstIn = blockStart_[iBlock]; 1591 if (iBlock == 0) { 1592 numberNonZero = 0; 1593 for (int i = 0; i < number; i++) { 1594 int iRow = piIndex[i]; 1595 unsigned char type = internalStatus[iRow]; 1596 if ((type & 7) < 6) { 1597 index[numberNonZero] = iRow; 1598 double piValue = pi[iRow]; 1599 array[numberNonZero++] = piValue; 1611 1600 } 1612 1601 } 1613 1602 } else { 1614 numberNonZero =firstIn;1615 } 1616 int n =numberNonZero;1617 rowStart += iBlock *numberRows;1618 rowEnd = rowStart +numberRows;1619 for (int i =0;i<number;i++) {1620 int iRow =piIndex[i];1603 numberNonZero = firstIn; 1604 } 1605 int n = numberNonZero; 1606 rowStart += iBlock * numberRows; 1607 rowEnd = rowStart + numberRows; 1608 for (int i = 0; i < number; i++) { 1609 int iRow = piIndex[i]; 1621 1610 double piValue = pi[iRow]; 1622 1611 CoinBigIndex end = rowEnd[iRow]; 1623 for (CoinBigIndex j =rowStart[iRow];j<end;j++) {1612 for (CoinBigIndex j = rowStart[iRow]; j < end; j++) { 1624 1613 int iColumn = column[j]; 1625 double value = element[j] *piValue;1626 double oldValue =arrayCandidate[iColumn];1614 double value = element[j] * piValue; 1615 double oldValue = arrayCandidate[iColumn]; 1627 1616 value += oldValue; 1628 1617 if (!oldValue) { 1629 arrayCandidate[iColumn]= value; 1630 index[n++]=iColumn;1631 } else if (value) { 1632 arrayCandidate[iColumn] = value; 1618 arrayCandidate[iColumn] = value; 1619 index[n++] = iColumn; 1620 } else if (value) { 1621 arrayCandidate[iColumn] = value; 1633 1622 } else { 1634 1623 arrayCandidate[iColumn] = COIN_INDEXED_REALLY_TINY_ELEMENT; 1635 1624 } 1636 1625 } … … 1638 1627 // pack down 1639 1628 double zeroTolerance = model_>zeroTolerance(); 1640 while (numberNonZero <n) {1641 int iSequence =index[numberNonZero];1642 double value =arrayCandidate[iSequence];1643 arrayCandidate[iSequence] =0.0;1644 unsigned char iStatus =internalStatus[iSequence]&7;1645 if (fabs(value) >zeroTolerance&&iStatus<6) {1646 index[numberNonZero] =iSequence;1647 array[numberNonZero++] =value;1629 while (numberNonZero < n) { 1630 int iSequence = index[numberNonZero]; 1631 double value = arrayCandidate[iSequence]; 1632 arrayCandidate[iSequence] = 0.0; 1633 unsigned char iStatus = internalStatus[iSequence] & 7; 1634 if (fabs(value) > zeroTolerance && iStatus < 6) { 1635 index[numberNonZero] = iSequence; 1636 array[numberNonZero++] = value; 1648 1637 } else { 1649 1638 // kill 1650 1639 n; 1651 index[numberNonZero] =index[n];1652 } 1653 } 1654 tableauRow.setNumElementsPartition(iBlock, numberNonZerofirstIn);1655 upperTheta = firstPass(model_, iBlock, upperTheta, freeSequence,1656 1657 candidateList);1640 index[numberNonZero] = index[n]; 1641 } 1642 } 1643 tableauRow.setNumElementsPartition(iBlock, numberNonZero  firstIn); 1644 upperTheta = firstPass(model_, iBlock, upperTheta, freeSequence, 1645 tableauRow, 1646 candidateList); 1658 1647 return upperTheta; 1659 1648 } 1660 1649 // gets sorted tableau row and a possible value of theta 1661 double 1662 AbcMatrix::dualColumn1Row1(double upperTheta, int & 1663 const CoinIndexedVector &update,1664 CoinPartitionedVector &tableauRow,1665 CoinPartitionedVector &candidateList) const1666 { 1667 assert(update.getNumElements() ==1);1650 double 1651 AbcMatrix::dualColumn1Row1(double upperTheta, int &freeSequence, 1652 const CoinIndexedVector &update, 1653 CoinPartitionedVector &tableauRow, 1654 CoinPartitionedVector &candidateList) const 1655 { 1656 assert(update.getNumElements() == 1); 1668 1657 int iRow = update.getIndices()[0]; 1669 1658 double piValue = update.denseVector()[iRow]; 1670 int * 1671 double * 1672 const CoinBigIndex * 1673 int numberRows =model_>numberRows();1674 const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;1659 int *COIN_RESTRICT index = tableauRow.getIndices(); 1660 double *COIN_RESTRICT array = tableauRow.denseVector(); 1661 const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_; 1662 int numberRows = model_>numberRows(); 1663 const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_; 1675 1664 CoinBigIndex start = rowStart[iRow]; 1676 1665 CoinBigIndex end = rowEnd[iRow]; 1677 const double * 1678 const int * 1666 const double *COIN_RESTRICT element = element_; 1667 const int *COIN_RESTRICT column = column_; 1679 1668 int numberNonZero; 1680 1669 int numberRemaining; 1681 numberNonZero =tableauRow.getNumElements(0);1670 numberNonZero = tableauRow.getNumElements(0); 1682 1671 numberRemaining = candidateList.getNumElements(0); 1683 double * COIN_RESTRICT arrayCandidate=candidateList.denseVector();1684 int * 1685 const double * 1686 const unsigned char * 1672 double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); 1673 int *COIN_RESTRICT indexCandidate = candidateList.getIndices(); 1674 const double *COIN_RESTRICT abcDj = model_>djRegion(); 1675 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 1687 1676 // do first pass to get possibles 1688 1677 double bestPossible = 0.0; … … 1690 1679 double tentativeTheta = 1.0e25; // try with this much smaller as guess 1691 1680 double acceptablePivot = model_>currentAcceptablePivot(); 1692 double dualT =model_>currentDualTolerance();1693 const double multiplier[] = { 1.0, 1.0 };1681 double dualT = model_>currentDualTolerance(); 1682 const double multiplier[] = { 1.0, 1.0 }; 1694 1683 double zeroTolerance = model_>zeroTolerance(); 1695 freeSequence =1;1684 freeSequence = 1; 1696 1685 if (model_>ordinaryVariables()) { 1697 for (CoinBigIndex j =start;j<end;j++) {1686 for (CoinBigIndex j = start; j < end; j++) { 1698 1687 int iSequence = column[j]; 1699 1688 double value = element[j]; 1700 double tableauValue = piValue *value;1701 if (fabs(tableauValue) >zeroTolerance) {1702 unsigned char iStatus=internalStatus[iSequence]&7;1703 if (iStatus<4) {1704 index[numberNonZero]=iSequence;1705 array[numberNonZero++]=tableauValue;1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1689 double tableauValue = piValue * value; 1690 if (fabs(tableauValue) > zeroTolerance) { 1691 unsigned char iStatus = internalStatus[iSequence] & 7; 1692 if (iStatus < 4) { 1693 index[numberNonZero] = iSequence; 1694 array[numberNonZero++] = tableauValue; 1695 double mult = multiplier[iStatus]; 1696 double alpha = tableauValue * mult; 1697 double oldValue = abcDj[iSequence] * mult; 1698 double value = oldValue  tentativeTheta * alpha; 1699 if (value < dualT) { 1700 bestPossible = CoinMax(bestPossible, alpha); 1701 value = oldValue  upperTheta * alpha; 1702 if (value < dualT && alpha >= acceptablePivot) { 1703 upperTheta = (oldValue  dualT) / alpha; 1704 } 1705 // add to list 1706 arrayCandidate[numberRemaining] = alpha; 1707 indexCandidate[numberRemaining++] = iSequence; 1708 } 1709 } 1721 1710 } 1722 1711 } … … 1724 1713 double badFree = 0.0; 1725 1714 double freeAlpha = model_>currentAcceptablePivot(); 1726 int freeSequenceIn =model_>freeSequenceIn();1715 int freeSequenceIn = model_>freeSequenceIn(); 1727 1716 double currentDualTolerance = model_>currentDualTolerance(); 1728 for (CoinBigIndex j =start;j<end;j++) {1717 for (CoinBigIndex j = start; j < end; j++) { 1729 1718 int iSequence = column[j]; 1730 1719 double value = element[j]; 1731 double tableauValue = piValue *value;1732 if (fabs(tableauValue) >zeroTolerance) {1733 unsigned char iStatus=internalStatus[iSequence]&7;1734 if (iStatus<6) {1735 if ((iStatus&4)==0) {1736 index[numberNonZero]=iSequence;1737 array[numberNonZero++]=tableauValue;1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 index[numberNonZero]=iSequence;1755 array[numberNonZero++]=tableauValue;1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1720 double tableauValue = piValue * value; 1721 if (fabs(tableauValue) > zeroTolerance) { 1722 unsigned char iStatus = internalStatus[iSequence] & 7; 1723 if (iStatus < 6) { 1724 if ((iStatus & 4) == 0) { 1725 index[numberNonZero] = iSequence; 1726 array[numberNonZero++] = tableauValue; 1727 double mult = multiplier[iStatus]; 1728 double alpha = tableauValue * mult; 1729 double oldValue = abcDj[iSequence] * mult; 1730 double value = oldValue  tentativeTheta * alpha; 1731 if (value < dualT) { 1732 bestPossible = CoinMax(bestPossible, alpha); 1733 value = oldValue  upperTheta * alpha; 1734 if (value < dualT && alpha >= acceptablePivot) { 1735 upperTheta = (oldValue  dualT) / alpha; 1736 } 1737 // add to list 1738 arrayCandidate[numberRemaining] = alpha; 1739 indexCandidate[numberRemaining++] = iSequence; 1740 } 1741 } else { 1742 bool keep; 1743 index[numberNonZero] = iSequence; 1744 array[numberNonZero++] = tableauValue; 1745 bestPossible = CoinMax(bestPossible, fabs(tableauValue)); 1746 double oldValue = abcDj[iSequence]; 1747 // If free has to be very large  should come in via dualRow 1748 //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e3) 1749 //break; 1750 if (oldValue > currentDualTolerance) { 1751 keep = true; 1752 } else if (oldValue < currentDualTolerance) { 1753 keep = true; 1754 } else { 1755 if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e5)) { 1756 keep = true; 1757 } else { 1758 keep = false; 1759 badFree = CoinMax(badFree, fabs(tableauValue)); 1760 } 1761 } 1762 if (keep) { 1774 1763 #ifdef PAN 1775 if (model_>fakeSuperBasic(iSequence)>=0) {1776 #endif 1777 if (iSequence==freeSequenceIn)1778 tableauValue=COIN_DBL_MAX;1779 1780 1781 1782 1783 1764 if (model_>fakeSuperBasic(iSequence) >= 0) { 1765 #endif 1766 if (iSequence == freeSequenceIn) 1767 tableauValue = COIN_DBL_MAX; 1768 // free  choose largest 1769 if (fabs(tableauValue) > fabs(freeAlpha)) { 1770 freeAlpha = tableauValue; 1771 freeSequence = iSequence; 1772 } 1784 1773 #ifdef PAN 1785 1786 #endif 1787 1788 1789 1790 } 1791 } 1792 } 1793 tableauRow.setNumElementsPartition(0, numberNonZero);1794 candidateList.setNumElementsPartition(0, numberRemaining);1774 } 1775 #endif 1776 } 1777 } 1778 } 1779 } 1780 } 1781 } 1782 tableauRow.setNumElementsPartition(0, numberNonZero); 1783 candidateList.setNumElementsPartition(0, numberRemaining); 1795 1784 return upperTheta; 1796 1785 } 1797 1786 //#define PARALLEL2 1798 1787 #ifdef PARALLEL2 1799 #undef cilk_for 1788 #undef cilk_for 1800 1789 #undef cilk_spawn 1801 1790 #undef cilk_sync … … 1823 1812 } 1824 1813 #endif 1825 void 1826 AbcMatrix::rebalance() const 1814 void AbcMatrix::rebalance() const 1827 1815 { 1828 1816 int maximumRows = model_>maximumAbcNumberRows(); … … 1836 1824 */ 1837 1825 #if ABC_PARALLEL 1838 int howOften =CoinMax(model_>factorization()>maximumPivots(),200);1839 if ((model_>numberIterations() %howOften)==0!startColumnBlock_[1]) {1840 int numberCpus =model_>numberCpus();1841 if (numberCpus >1) {1842 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;1843 const unsigned char * 1844 int numberRows =model_>numberRows();1845 int total =0;1846 for (int iSequence =0;iSequence<numberRows;iSequence++) {1847 unsigned char iStatus=internalStatus[iSequence]&7;1848 if (iStatus<4) 1849 total+=3;1850 } 1851 int totalSlacks =total;1826 int howOften = CoinMax(model_>factorization()>maximumPivots(), 200); 1827 if ((model_>numberIterations() % howOften) == 0  !startColumnBlock_[1]) { 1828 int numberCpus = model_>numberCpus(); 1829 if (numberCpus > 1) { 1830 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 1831 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 1832 int numberRows = model_>numberRows(); 1833 int total = 0; 1834 for (int iSequence = 0; iSequence < numberRows; iSequence++) { 1835 unsigned char iStatus = internalStatus[iSequence] & 7; 1836 if (iStatus < 4) 1837 total += 3; 1838 } 1839 int totalSlacks = total; 1852 1840 CoinBigIndex end = columnStart[maximumRows]; 1853 for (int iSequence =maximumRows;iSequence<numberTotal;iSequence++) {1854 CoinBigIndex start=end;1855 end = columnStart[iSequence+1];1856 unsigned char iStatus=internalStatus[iSequence]&7;1857 if (iStatus<4) 1858 total+=5+2*(endstart);1859 1860 total+=1;1861 } 1862 int chunk = (total +numberCpus1)/numberCpus;1863 total =totalSlacks;1864 int iCpu =0;1865 startColumnBlock_[0] =0;1841 for (int iSequence = maximumRows; iSequence < numberTotal; iSequence++) { 1842 CoinBigIndex start = end; 1843 end = columnStart[iSequence + 1]; 1844 unsigned char iStatus = internalStatus[iSequence] & 7; 1845 if (iStatus < 4) 1846 total += 5 + 2 * (end  start); 1847 else 1848 total += 1; 1849 } 1850 int chunk = (total + numberCpus  1) / numberCpus; 1851 total = totalSlacks; 1852 int iCpu = 0; 1853 startColumnBlock_[0] = 0; 1866 1854 end = columnStart[maximumRows]; 1867 for (int iSequence =maximumRows;iSequence<numberTotal;iSequence++) {1868 CoinBigIndex start=end;1869 end = columnStart[iSequence+1];1870 unsigned char iStatus=internalStatus[iSequence]&7;1871 if (iStatus<4) 1872 total+=5+2*(endstart);1873 1874 total+=1;1875 if (total>chunk) {1876 1877 total=0;1878 startColumnBlock_[iCpu]=iSequence;1879 1880 } 1881 assert(iCpu <numberCpus);1855 for (int iSequence = maximumRows; iSequence < numberTotal; iSequence++) { 1856 CoinBigIndex start = end; 1857 end = columnStart[iSequence + 1]; 1858 unsigned char iStatus = internalStatus[iSequence] & 7; 1859 if (iStatus < 4) 1860 total += 5 + 2 * (end  start); 1861 else 1862 total += 1; 1863 if (total > chunk) { 1864 iCpu++; 1865 total = 0; 1866 startColumnBlock_[iCpu] = iSequence; 1867 } 1868 } 1869 assert(iCpu < numberCpus); 1882 1870 iCpu++; 1883 for (; iCpu<=numberCpus;iCpu++)1884 startColumnBlock_[iCpu]=numberTotal;1885 numberColumnBlocks_ =numberCpus;1871 for (; iCpu <= numberCpus; iCpu++) 1872 startColumnBlock_[iCpu] = numberTotal; 1873 numberColumnBlocks_ = numberCpus; 1886 1874 } else { 1887 numberColumnBlocks_ =1;1888 startColumnBlock_[0] =0;1889 startColumnBlock_[1] =numberTotal;1875 numberColumnBlocks_ = 1; 1876 startColumnBlock_[0] = 0; 1877 startColumnBlock_[1] = numberTotal; 1890 1878 } 1891 1879 } 1892 1880 #else 1893 startColumnBlock_[0] =0;1894 startColumnBlock_[1] =numberTotal;1895 #endif 1896 } 1897 double 1898 AbcMatrix::dualColumn1(const CoinIndexedVector & 1899 CoinPartitionedVector & tableauRow,CoinPartitionedVector &candidateList) const1881 startColumnBlock_[0] = 0; 1882 startColumnBlock_[1] = numberTotal; 1883 #endif 1884 } 1885 double 1886 AbcMatrix::dualColumn1(const CoinIndexedVector &update, 1887 CoinPartitionedVector &tableauRow, CoinPartitionedVector &candidateList) const 1900 1888 { 1901 1889 int numberTotal = model_>numberTotal(); … … 1904 1892 tableauRow.setPackedMode(true); 1905 1893 candidateList.setPackedMode(true); 1906 int number =update.getNumElements();1894 int number = update.getNumElements(); 1907 1895 double upperTheta; 1908 if (rowStart_ &&number<3) {1909 #if ABC_INSTRUMENT >11896 if (rowStart_ && number < 3) { 1897 #if ABC_INSTRUMENT > 1 1910 1898 { 1911 int n =update.getNumElements();1912 abcPricing[n <19 ? n : 19]++;1899 int n = update.getNumElements(); 1900 abcPricing[n < 19 ? n : 19]++; 1913 1901 } 1914 1902 #endif … … 1916 1904 // do slacks first 1917 1905 int starts[2]; 1918 starts[0] =0;1919 starts[1] =numberTotal;1920 tableauRow.setPartitions(1, starts);1921 candidateList.setPartitions(1, starts);1922 upperTheta = static_cast< AbcSimplexDual *>(model_)>dualColumn1A();1906 starts[0] = 0; 1907 starts[1] = numberTotal; 1908 tableauRow.setPartitions(1, starts); 1909 candidateList.setPartitions(1, starts); 1910 upperTheta = static_cast< AbcSimplexDual * >(model_)>dualColumn1A(); 1923 1911 //assert (upperTheta>100*model_>dualTolerance()model_>sequenceIn()>=0model_>lastFirstFree()>=0); 1924 int freeSequence =1;1912 int freeSequence = 1; 1925 1913 // worth using row copy 1926 assert 1927 if (number ==2) {1928 upperTheta =dualColumn1Row2(upperTheta,freeSequence,update,tableauRow,candidateList);1914 assert(number); 1915 if (number == 2) { 1916 upperTheta = dualColumn1Row2(upperTheta, freeSequence, update, tableauRow, candidateList); 1929 1917 } else { 1930 upperTheta =dualColumn1Row1(upperTheta,freeSequence,update,tableauRow,candidateList);1931 } 1932 if (freeSequence >=0) {1933 int numberNonZero =tableauRow.getNumElements(0);1934 const int * 1935 const double * 1918 upperTheta = dualColumn1Row1(upperTheta, freeSequence, update, tableauRow, candidateList); 1919 } 1920 if (freeSequence >= 0) { 1921 int numberNonZero = tableauRow.getNumElements(0); 1922 const int *COIN_RESTRICT index = tableauRow.getIndices(); 1923 const double *COIN_RESTRICT array = tableauRow.denseVector(); 1936 1924 // search for free coming in 1937 double freeAlpha =0.0;1938 int bestSequence =model_>sequenceIn();1939 if (bestSequence >=0)1940 freeAlpha=model_>alpha();1925 double freeAlpha = 0.0; 1926 int bestSequence = model_>sequenceIn(); 1927 if (bestSequence >= 0) 1928 freeAlpha = model_>alpha(); 1941 1929 index = tableauRow.getIndices(); 1942 1930 array = tableauRow.denseVector(); 1943 // free variable  search 1944 for (int k =0;k<numberNonZero;k++) {1945 if (freeSequence==index[k]) {1946 if (fabs(freeAlpha)<fabs(array[k])) {1947 freeAlpha=array[k];1948 1949 1950 1951 } 1952 if (model_>sequenceIn() <0fabs(freeAlpha)>fabs(model_>alpha())) {1953 1954 1955 1956 1931 // free variable  search 1932 for (int k = 0; k < numberNonZero; k++) { 1933 if (freeSequence == index[k]) { 1934 if (fabs(freeAlpha) < fabs(array[k])) { 1935 freeAlpha = array[k]; 1936 } 1937 break; 1938 } 1939 } 1940 if (model_>sequenceIn() < 0  fabs(freeAlpha) > fabs(model_>alpha())) { 1941 double oldValue = model_>djRegion()[freeSequence]; 1942 model_>setSequenceIn(freeSequence); 1943 model_>setAlpha(freeAlpha); 1944 model_>setTheta(oldValue / freeAlpha); 1957 1945 } 1958 1946 } … … 1960 1948 // three or more 1961 1949 // need to do better job on dividing up (but wait until vector or by row) 1962 upperTheta = parallelDual4(static_cast< AbcSimplexDual *>(model_));1950 upperTheta = parallelDual4(static_cast< AbcSimplexDual * >(model_)); 1963 1951 } 1964 1952 //tableauRow.compact(); … … 1968 1956 #endif 1969 1957 candidateList.computeNumberElements(); 1970 int numberRemaining =candidateList.getNumElements();1971 if (!numberRemaining &&model_>sequenceIn()<0) {1958 int numberRemaining = candidateList.getNumElements(); 1959 if (!numberRemaining && model_>sequenceIn() < 0) { 1972 1960 return COIN_DBL_MAX; // Looks infeasible 1973 1961 } else { … … 1975 1963 } 1976 1964 } 1977 #define _mm256_broadcast_sd(x) static_cast< __m256d> (__builtin_ia32_vbroadcastsd256(x))1965 #define _mm256_broadcast_sd(x) static_cast< __m256d >(__builtin_ia32_vbroadcastsd256(x)) 1978 1966 #define _mm256_load_pd(x) *(__m256d *)(x) 1979 #define _mm256_store_pd (s, x) *((__m256d *)s)=x 1980 void 1981 AbcMatrix::dualColumn1Part(int iBlock, int & sequenceIn, double & upperThetaResult, 1982 const CoinIndexedVector & update, 1983 CoinPartitionedVector & tableauRow,CoinPartitionedVector & candidateList) const 1984 { 1985 double upperTheta=upperThetaResult; 1967 #define _mm256_store_pd (s, x) * ((__m256d *)s) = x 1968 void AbcMatrix::dualColumn1Part(int iBlock, int &sequenceIn, double &upperThetaResult, 1969 const CoinIndexedVector &update, 1970 CoinPartitionedVector &tableauRow, CoinPartitionedVector &candidateList) const 1971 { 1972 double upperTheta = upperThetaResult; 1986 1973 #if 0 1987 1974 double time0=CoinCpuTime(); 1988 1975 #endif 1989 1976 int maximumRows = model_>maximumAbcNumberRows(); 1990 int firstIn =startColumnBlock_[iBlock];1991 int last = startColumnBlock_[iBlock +1];1977 int firstIn = startColumnBlock_[iBlock]; 1978 int last = startColumnBlock_[iBlock + 1]; 1992 1979 int numberNonZero; 1993 1980 int numberRemaining; 1994 1981 int first; 1995 if (firstIn ==0) {1996 upperTheta =static_cast<AbcSimplexDual *>(model_)>dualColumn1A();1997 numberNonZero =tableauRow.getNumElements(0);1982 if (firstIn == 0) { 1983 upperTheta = static_cast< AbcSimplexDual * >(model_)>dualColumn1A(); 1984 numberNonZero = tableauRow.getNumElements(0); 1998 1985 numberRemaining = candidateList.getNumElements(0); 1999 first =maximumRows;1986 first = maximumRows; 2000 1987 } else { 2001 first =firstIn;2002 numberNonZero =first;2003 numberRemaining =first;2004 } 2005 sequenceIn =1;1988 first = firstIn; 1989 numberNonZero = first; 1990 numberRemaining = first; 1991 } 1992 sequenceIn = 1; 2006 1993 // get matrix data pointers 2007 const int * 2008 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;2009 const double * 2010 const double * COIN_RESTRICT pi=update.denseVector();1994 const int *COIN_RESTRICT row = matrix_>getIndices(); 1995 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 1996 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 1997 const double *COIN_RESTRICT pi = update.denseVector(); 2011 1998 //const int * COIN_RESTRICT piIndex = update.getIndices(); 2012 int * 2013 double * 1999 int *COIN_RESTRICT index = tableauRow.getIndices(); 2000 double *COIN_RESTRICT array = tableauRow.denseVector(); 2014 2001 // pivot elements 2015 double * COIN_RESTRICT arrayCandidate=candidateList.denseVector();2002 double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); 2016 2003 // indices 2017 int * 2018 const double * 2019 const unsigned char * 2004 int *COIN_RESTRICT indexCandidate = candidateList.getIndices(); 2005 const double *COIN_RESTRICT abcDj = model_>djRegion(); 2006 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 2020 2007 // do first pass to get possibles 2021 2008 double bestPossible = 0.0; … … 2023 2010 double tentativeTheta = 1.0e25; // try with this much smaller as guess 2024 2011 double acceptablePivot = model_>currentAcceptablePivot(); 2025 double dualT =model_>currentDualTolerance();2026 const double multiplier[] = { 1.0, 1.0 };2012 double dualT = model_>currentDualTolerance(); 2013 const double multiplier[] = { 1.0, 1.0 }; 2027 2014 double zeroTolerance = model_>zeroTolerance(); 2028 2015 if (model_>ordinaryVariables()) { 2029 2016 #ifdef COUNT_COPY 2030 if (first >maximumRowslast<model_>numberTotal()false) {2017 if (first > maximumRows  last < model_>numberTotal()  false) { 2031 2018 #endif 2032 2019 #if 1 2033 for (int iSequence = first; iSequence < last; iSequence++) {2034 unsigned char iStatus=internalStatus[iSequence]&7;2035 if (iStatus<4) {2036 2037 CoinBigIndex next = columnStart[iSequence+1];2038 2039 2040 2041 2042 2043 if (fabs(tableauValue)>zeroTolerance) {2020 for (int iSequence = first; iSequence < last; iSequence++) { 2021 unsigned char iStatus = internalStatus[iSequence] & 7; 2022 if (iStatus < 4) { 2023 CoinBigIndex start = columnStart[iSequence]; 2024 CoinBigIndex next = columnStart[iSequence + 1]; 2025 double tableauValue = 0.0; 2026 for (CoinBigIndex j = start; j < next; j++) { 2027 int jRow = row[j]; 2028 tableauValue += pi[jRow] * elementByColumn[j]; 2029 } 2030 if (fabs(tableauValue) > zeroTolerance) { 2044 2031 #else 2045 cilk_for (int iSequence = first; iSequence < last; iSequence++) { 2046 unsigned char iStatus=internalStatus[iSequence]&7; 2047 if (iStatus<4) { 2048 CoinBigIndex start = columnStart[iSequence]; 2049 CoinBigIndex next = columnStart[iSequence+1]; 2050 double tableauValue = 0.0; 2051 for (CoinBigIndex j = start; j < next; j++) { 2052 int jRow = row[j]; 2053 tableauValue += pi[jRow] * elementByColumn[j]; 2054 } 2055 array[iSequence]=tableauValue; 2032 cilk_for(int iSequence = first; iSequence < last; iSequence++) 2033 { 2034 unsigned char iStatus = internalStatus[iSequence] & 7; 2035 if (iStatus < 4) { 2036 CoinBigIndex start = columnStart[iSequence]; 2037 CoinBigIndex next = columnStart[iSequence + 1]; 2038 double tableauValue = 0.0; 2039 for (CoinBigIndex j = start; j < next; j++) { 2040 int jRow = row[j]; 2041 tableauValue += pi[jRow] * elementByColumn[j]; 2042 } 2043 array[iSequence] = tableauValue; 2056 2044 } 2057 2045 } 2058 2046 //printf("time %.6g\n",CoinCpuTime()time0); 2059 2047 for (int iSequence = first; iSequence < last; iSequence++) { 2060 double tableauValue =array[iSequence];2048 double tableauValue = array[iSequence]; 2061 2049 if (tableauValue) { 2062 array[iSequence]=0.0;2063 if (fabs(tableauValue)>zeroTolerance) {2064 unsigned char iStatus=internalStatus[iSequence]&7;2065 #endif 2066 index[numberNonZero]=iSequence;2067 array[numberNonZero++]=tableauValue;2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082 2083 }2084 }2050 array[iSequence] = 0.0; 2051 if (fabs(tableauValue) > zeroTolerance) { 2052 unsigned char iStatus = internalStatus[iSequence] & 7; 2053 #endif 2054 index[numberNonZero] = iSequence; 2055 array[numberNonZero++] = tableauValue; 2056 double mult = multiplier[iStatus]; 2057 double alpha = tableauValue * mult; 2058 double oldValue = abcDj[iSequence] * mult; 2059 double value = oldValue  tentativeTheta * alpha; 2060 if (value < dualT) { 2061 bestPossible = CoinMax(bestPossible, alpha); 2062 value = oldValue  upperTheta * alpha; 2063 if (value < dualT && alpha >= acceptablePivot) { 2064 upperTheta = (oldValue  dualT) / alpha; 2065 } 2066 // add to list 2067 arrayCandidate[numberRemaining] = alpha; 2068 indexCandidate[numberRemaining++] = iSequence; 2069 } 2070 } 2071 } 2072 } 2085 2073 #ifdef COUNT_COPY 2086 2074 } else { 2087 const double * 2088 const int * 2089 double temp[4] __attribute__ ((aligned(32)));2090 memset(temp, 0,sizeof(temp));2091 for (int iCount =smallestCount_;iCount<largestCount_;iCount++) {2092 int iStart=countFirst_[iCount];2093 int number=countFirst_[iCount+1]iStart;2094 2095 2096 const int * COIN_RESTRICT blockRow = row+countStart_[iCount];2097 const double * COIN_RESTRICT blockElement = element+countStart_[iCount];2098 const int * COIN_RESTRICT sequence = countRealColumn_+iStart;2099 2100 int numberBlocks=number>>2;2101 2102 for (int iBlock=0;iBlock<numberBlocks;iBlock++) {2103 double tableauValue0=0.0;2104 double tableauValue1=0.0;2105 double tableauValue2=0.0;2106 double tableauValue3=0.0;2107 for (int j=0;j<iCount;j++) {2108 tableauValue0 += pi[blockRow[0]]*blockElement[0];2109 tableauValue1 += pi[blockRow[1]]*blockElement[1];2110 tableauValue2 += pi[blockRow[2]]*blockElement[2];2111 tableauValue3 += pi[blockRow[3]]*blockElement[3];2112 blockRow+=4;2113 blockElement+=4;2114 2115 if (fabs(tableauValue0)>zeroTolerance) {2116 int iSequence=sequence[0];2117 unsigned char iStatus=internalStatus[iSequence]&7;2118 if (iStatus<4) {2119 index[numberNonZero]=iSequence;2120 array[numberNonZero++]=tableauValue0;2121 2122 2123 2124 2125 2126 2127 2128 2129 2130 2131 2132 2133 2134 2135 2136 2137 if (fabs(tableauValue1)>zeroTolerance) {2138 int iSequence=sequence[1];2139 unsigned char iStatus=internalStatus[iSequence]&7;2140 if (iStatus<4) {2141 index[numberNonZero]=iSequence;2142 array[numberNonZero++]=tableauValue1;2143 2144 2145 2146 2147 2148 2149 2150 2151 2152 2153 2154 2155 2156 2157 2158 2159 if (fabs(tableauValue2)>zeroTolerance) {2160 int iSequence=sequence[2];2161 unsigned char iStatus=internalStatus[iSequence]&7;2162 if (iStatus<4) {2163 index[numberNonZero]=iSequence;2164 array[numberNonZero++]=tableauValue2;2165 2166 2167 2168 2169 2170 2171 2172 2173 2174 2175 2176 2177 2178 2179 2180 2181 if (fabs(tableauValue3)>zeroTolerance) {2182 int iSequence=sequence[3];2183 unsigned char iStatus=internalStatus[iSequence]&7;2184 if (iStatus<4) {2185 index[numberNonZero]=iSequence;2186 array[numberNonZero++]=tableauValue3;2187 2188 2189 2190 2191 2192 2193 2194 2195 2196 2197 2198 2199 2200 2201 2202 2203 sequence+=4;2204 2205 for (int i=0;i<number;i++) {2206 int iSequence=sequence[i];2207 unsigned char iStatus=internalStatus[iSequence]&7;2208 if (iStatus<4) {2209 double tableauValue=0.0;2210 for (int j=0;j<iCount;j++) {2211 int iRow=blockRow[4*j];2212 tableauValue += pi[iRow]*blockElement[4*j];2213 2214 if (fabs(tableauValue)>zeroTolerance) {2215 index[numberNonZero]=iSequence;2216 array[numberNonZero++]=tableauValue;2217 2218 2219 2220 2221 2222 2223 2224 2225 2226 2227 2228 2229 2230 2231 2232 2233 2234 2235 2236 } 2237 int numberColumns =model_>numberColumns();2075 const double *COIN_RESTRICT element = countElement_; 2076 const int *COIN_RESTRICT row = countRow_; 2077 double temp[4] __attribute__((aligned(32))); 2078 memset(temp, 0, sizeof(temp)); 2079 for (int iCount = smallestCount_; iCount < largestCount_; iCount++) { 2080 int iStart = countFirst_[iCount]; 2081 int number = countFirst_[iCount + 1]  iStart; 2082 if (!number) 2083 continue; 2084 const int *COIN_RESTRICT blockRow = row + countStart_[iCount]; 2085 const double *COIN_RESTRICT blockElement = element + countStart_[iCount]; 2086 const int *COIN_RESTRICT sequence = countRealColumn_ + iStart; 2087 // really need to sort and swap to avoid tests 2088 int numberBlocks = number >> 2; 2089 number &= 3; 2090 for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { 2091 double tableauValue0 = 0.0; 2092 double tableauValue1 = 0.0; 2093 double tableauValue2 = 0.0; 2094 double tableauValue3 = 0.0; 2095 for (int j = 0; j < iCount; j++) { 2096 tableauValue0 += pi[blockRow[0]] * blockElement[0]; 2097 tableauValue1 += pi[blockRow[1]] * blockElement[1]; 2098 tableauValue2 += pi[blockRow[2]] * blockElement[2]; 2099 tableauValue3 += pi[blockRow[3]] * blockElement[3]; 2100 blockRow += 4; 2101 blockElement += 4; 2102 } 2103 if (fabs(tableauValue0) > zeroTolerance) { 2104 int iSequence = sequence[0]; 2105 unsigned char iStatus = internalStatus[iSequence] & 7; 2106 if (iStatus < 4) { 2107 index[numberNonZero] = iSequence; 2108 array[numberNonZero++] = tableauValue0; 2109 double mult = multiplier[iStatus]; 2110 double alpha = tableauValue0 * mult; 2111 double oldValue = abcDj[iSequence] * mult; 2112 double value = oldValue  tentativeTheta * alpha; 2113 if (value < dualT) { 2114 bestPossible = CoinMax(bestPossible, alpha); 2115 value = oldValue  upperTheta * alpha; 2116 if (value < dualT && alpha >= acceptablePivot) { 2117 upperTheta = (oldValue  dualT) / alpha; 2118 } 2119 // add to list 2120 arrayCandidate[numberRemaining] = alpha; 2121 indexCandidate[numberRemaining++] = iSequence; 2122 } 2123 } 2124 } 2125 if (fabs(tableauValue1) > zeroTolerance) { 2126 int iSequence = sequence[1]; 2127 unsigned char iStatus = internalStatus[iSequence] & 7; 2128 if (iStatus < 4) { 2129 index[numberNonZero] = iSequence; 2130 array[numberNonZero++] = tableauValue1; 2131 double mult = multiplier[iStatus]; 2132 double alpha = tableauValue1 * mult; 2133 double oldValue = abcDj[iSequence] * mult; 2134 double value = oldValue  tentativeTheta * alpha; 2135 if (value < dualT) { 2136 bestPossible = CoinMax(bestPossible, alpha); 2137 value = oldValue  upperTheta * alpha; 2138 if (value < dualT && alpha >= acceptablePivot) { 2139 upperTheta = (oldValue  dualT) / alpha; 2140 } 2141 // add to list 2142 arrayCandidate[numberRemaining] = alpha; 2143 indexCandidate[numberRemaining++] = iSequence; 2144 } 2145 } 2146 } 2147 if (fabs(tableauValue2) > zeroTolerance) { 2148 int iSequence = sequence[2]; 2149 unsigned char iStatus = internalStatus[iSequence] & 7; 2150 if (iStatus < 4) { 2151 index[numberNonZero] = iSequence; 2152 array[numberNonZero++] = tableauValue2; 2153 double mult = multiplier[iStatus]; 2154 double alpha = tableauValue2 * mult; 2155 double oldValue = abcDj[iSequence] * mult; 2156 double value = oldValue  tentativeTheta * alpha; 2157 if (value < dualT) { 2158 bestPossible = CoinMax(bestPossible, alpha); 2159 value = oldValue  upperTheta * alpha; 2160 if (value < dualT && alpha >= acceptablePivot) { 2161 upperTheta = (oldValue  dualT) / alpha; 2162 } 2163 // add to list 2164 arrayCandidate[numberRemaining] = alpha; 2165 indexCandidate[numberRemaining++] = iSequence; 2166 } 2167 } 2168 } 2169 if (fabs(tableauValue3) > zeroTolerance) { 2170 int iSequence = sequence[3]; 2171 unsigned char iStatus = internalStatus[iSequence] & 7; 2172 if (iStatus < 4) { 2173 index[numberNonZero] = iSequence; 2174 array[numberNonZero++] = tableauValue3; 2175 double mult = multiplier[iStatus]; 2176 double alpha = tableauValue3 * mult; 2177 double oldValue = abcDj[iSequence] * mult; 2178 double value = oldValue  tentativeTheta * alpha; 2179 if (value < dualT) { 2180 bestPossible = CoinMax(bestPossible, alpha); 2181 value = oldValue  upperTheta * alpha; 2182 if (value < dualT && alpha >= acceptablePivot) { 2183 upperTheta = (oldValue  dualT) / alpha; 2184 } 2185 // add to list 2186 arrayCandidate[numberRemaining] = alpha; 2187 indexCandidate[numberRemaining++] = iSequence; 2188 } 2189 } 2190 } 2191 sequence += 4; 2192 } 2193 for (int i = 0; i < number; i++) { 2194 int iSequence = sequence[i]; 2195 unsigned char iStatus = internalStatus[iSequence] & 7; 2196 if (iStatus < 4) { 2197 double tableauValue = 0.0; 2198 for (int j = 0; j < iCount; j++) { 2199 int iRow = blockRow[4 * j]; 2200 tableauValue += pi[iRow] * blockElement[4 * j]; 2201 } 2202 if (fabs(tableauValue) > zeroTolerance) { 2203 index[numberNonZero] = iSequence; 2204 array[numberNonZero++] = tableauValue; 2205 double mult = multiplier[iStatus]; 2206 double alpha = tableauValue * mult; 2207 double oldValue = abcDj[iSequence] * mult; 2208 double value = oldValue  tentativeTheta * alpha; 2209 if (value < dualT) { 2210 bestPossible = CoinMax(bestPossible, alpha); 2211 value = oldValue  upperTheta * alpha; 2212 if (value < dualT && alpha >= acceptablePivot) { 2213 upperTheta = (oldValue  dualT) / alpha; 2214 } 2215 // add to list 2216 arrayCandidate[numberRemaining] = alpha; 2217 indexCandidate[numberRemaining++] = iSequence; 2218 } 2219 } 2220 } 2221 blockRow++; 2222 blockElement++; 2223 } 2224 } 2225 int numberColumns = model_>numberColumns(); 2238 2226 // large ones 2239 const CoinBigIndex * COIN_RESTRICT largeStart = countStartLarge_countFirst_[MAX_COUNT];2240 for (int i =countFirst_[MAX_COUNT];i<numberColumns;i++) {2241 int iSequence=countRealColumn_[i];2242 unsigned char iStatus=internalStatus[iSequence]&7;2243 if (iStatus<4) {2244 2245 CoinBigIndex next = largeStart[i+1];2246 2247 2248 2249 2250 2251 if (fabs(tableauValue)>zeroTolerance) {2252 index[numberNonZero]=iSequence;2253 array[numberNonZero++]=tableauValue;2254 2255 2256 2257 2258 2259 2260 2261 2262 2263 2264 2265 2266 2267 2268 2269 2227 const CoinBigIndex *COIN_RESTRICT largeStart = countStartLarge_  countFirst_[MAX_COUNT]; 2228 for (int i = countFirst_[MAX_COUNT]; i < numberColumns; i++) { 2229 int iSequence = countRealColumn_[i]; 2230 unsigned char iStatus = internalStatus[iSequence] & 7; 2231 if (iStatus < 4) { 2232 CoinBigIndex start = largeStart[i]; 2233 CoinBigIndex next = largeStart[i + 1]; 2234 double tableauValue = 0.0; 2235 for (CoinBigIndex j = start; j < next; j++) { 2236 int jRow = row[j]; 2237 tableauValue += pi[jRow] * element[j]; 2238 } 2239 if (fabs(tableauValue) > zeroTolerance) { 2240 index[numberNonZero] = iSequence; 2241 array[numberNonZero++] = tableauValue; 2242 double mult = multiplier[iStatus]; 2243 double alpha = tableauValue * mult; 2244 double oldValue = abcDj[iSequence] * mult; 2245 double value = oldValue  tentativeTheta * alpha; 2246 if (value < dualT) { 2247 bestPossible = CoinMax(bestPossible, alpha); 2248 value = oldValue  upperTheta * alpha; 2249 if (value < dualT && alpha >= acceptablePivot) { 2250 upperTheta = (oldValue  dualT) / alpha; 2251 } 2252 // add to list 2253 arrayCandidate[numberRemaining] = alpha; 2254 indexCandidate[numberRemaining++] = iSequence; 2255 } 2256 } 2257 } 2270 2258 } 2271 2259 } … … 2274 2262 double badFree = 0.0; 2275 2263 double freePivot = model_>currentAcceptablePivot(); 2276 int freeSequenceIn =model_>freeSequenceIn();2264 int freeSequenceIn = model_>freeSequenceIn(); 2277 2265 double currentDualTolerance = model_>currentDualTolerance(); 2278 2266 for (int iSequence = first; iSequence < last; iSequence++) { 2279 unsigned char iStatus =internalStatus[iSequence]&7;2280 if (iStatus <6) {2281 2282 CoinBigIndex next = columnStart[iSequence+1];2283 2284 2285 2286 2287 2288 if (fabs(tableauValue)>zeroTolerance) {2289 if ((iStatus&4)==0) {2290 index[numberNonZero]=iSequence;2291 array[numberNonZero++]=tableauValue;2292 2293 2294 2295 2296 2297 2298 2299 2300 2301 2302 2303 2304 2305 2306 2307 2308 index[numberNonZero]=iSequence;2309 array[numberNonZero++]=tableauValue;2310 2311 2312 2313 2314 2315 2316 2317 2318 2319 2320 2321 2322 2323 2324 2325 2326 2327 2267 unsigned char iStatus = internalStatus[iSequence] & 7; 2268 if (iStatus < 6) { 2269 CoinBigIndex start = columnStart[iSequence]; 2270 CoinBigIndex next = columnStart[iSequence + 1]; 2271 double tableauValue = 0.0; 2272 for (CoinBigIndex j = start; j < next; j++) { 2273 int jRow = row[j]; 2274 tableauValue += pi[jRow] * elementByColumn[j]; 2275 } 2276 if (fabs(tableauValue) > zeroTolerance) { 2277 if ((iStatus & 4) == 0) { 2278 index[numberNonZero] = iSequence; 2279 array[numberNonZero++] = tableauValue; 2280 double mult = multiplier[iStatus]; 2281 double alpha = tableauValue * mult; 2282 double oldValue = abcDj[iSequence] * mult; 2283 double value = oldValue  tentativeTheta * alpha; 2284 if (value < dualT) { 2285 bestPossible = CoinMax(bestPossible, alpha); 2286 value = oldValue  upperTheta * alpha; 2287 if (value < dualT && alpha >= acceptablePivot) { 2288 upperTheta = (oldValue  dualT) / alpha; 2289 } 2290 // add to list 2291 arrayCandidate[numberRemaining] = alpha; 2292 indexCandidate[numberRemaining++] = iSequence; 2293 } 2294 } else { 2295 bool keep; 2296 index[numberNonZero] = iSequence; 2297 array[numberNonZero++] = tableauValue; 2298 bestPossible = CoinMax(bestPossible, fabs(tableauValue)); 2299 double oldValue = abcDj[iSequence]; 2300 // If free has to be very large  should come in via dualRow 2301 //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e3) 2302 //break; 2303 if (oldValue > currentDualTolerance) { 2304 keep = true; 2305 } else if (oldValue < currentDualTolerance) { 2306 keep = true; 2307 } else { 2308 if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e5)) { 2309 keep = true; 2310 } else { 2311 keep = false; 2312 badFree = CoinMax(badFree, fabs(tableauValue)); 2313 } 2314 } 2315 if (keep) { 2328 2316 #ifdef PAN 2329 if (model_>fakeSuperBasic(iSequence)>=0) {2330 #endif 2331 if (iSequence==freeSequenceIn)2332 tableauValue=COIN_DBL_MAX;2333 2334 2335 2336 2337 2317 if (model_>fakeSuperBasic(iSequence) >= 0) { 2318 #endif 2319 if (iSequence == freeSequenceIn) 2320 tableauValue = COIN_DBL_MAX; 2321 // free  choose largest 2322 if (fabs(tableauValue) > fabs(freePivot)) { 2323 freePivot = tableauValue; 2324 sequenceIn = iSequence; 2325 } 2338 2326 #ifdef PAN 2339 2340 #endif 2341 2342 2343 2344 } 2345 } 2346 } 2327 } 2328 #endif 2329 } 2330 } 2331 } 2332 } 2333 } 2334 } 2347 2335 // adjust 2348 2336 numberNonZero = firstIn; 2349 2337 numberRemaining = firstIn; 2350 tableauRow.setNumElementsPartition(iBlock, numberNonZero);2351 candidateList.setNumElementsPartition(iBlock, numberRemaining);2352 if (!numberRemaining &&model_>sequenceIn()<0) {2353 2354 upperThetaResult =COIN_DBL_MAX; // Looks infeasible2338 tableauRow.setNumElementsPartition(iBlock, numberNonZero); 2339 candidateList.setNumElementsPartition(iBlock, numberRemaining); 2340 if (!numberRemaining && model_>sequenceIn() < 0) { 2341 2342 upperThetaResult = COIN_DBL_MAX; // Looks infeasible 2355 2343 } else { 2356 upperThetaResult=upperTheta; 2357 } 2358 } 2359 // gets tableau row  returns number of slacks in block 2360 int 2361 AbcMatrix::primalColumnRow(int iBlock,bool doByRow,const CoinIndexedVector & updates, 2362 CoinPartitionedVector & tableauRow) const 2344 upperThetaResult = upperTheta; 2345 } 2346 } 2347 // gets tableau row  returns number of slacks in block 2348 int AbcMatrix::primalColumnRow(int iBlock, bool doByRow, const CoinIndexedVector &updates, 2349 CoinPartitionedVector &tableauRow) const 2363 2350 { 2364 2351 int maximumRows = model_>maximumAbcNumberRows(); 2365 int first =tableauRow.startPartition(iBlock);2366 int last =tableauRow.startPartition(iBlock+1);2367 if (doByRow) { 2368 assert(first ==blockStart_[iBlock]);2369 assert(last ==blockStart_[iBlock+1]);2352 int first = tableauRow.startPartition(iBlock); 2353 int last = tableauRow.startPartition(iBlock + 1); 2354 if (doByRow) { 2355 assert(first == blockStart_[iBlock]); 2356 assert(last == blockStart_[iBlock + 1]); 2370 2357 } else { 2371 assert(first ==startColumnBlock_[iBlock]);2372 assert(last ==startColumnBlock_[iBlock+1]);2373 } 2374 int numberSlacks =updates.getNumElements();2375 int numberNonZero =0;2376 const double * COIN_RESTRICT pi=updates.denseVector();2377 const int * 2378 int * COIN_RESTRICT index = tableauRow.getIndices()+first;2379 double * COIN_RESTRICT array = tableauRow.denseVector()+first;2380 const unsigned char * 2358 assert(first == startColumnBlock_[iBlock]); 2359 assert(last == startColumnBlock_[iBlock + 1]); 2360 } 2361 int numberSlacks = updates.getNumElements(); 2362 int numberNonZero = 0; 2363 const double *COIN_RESTRICT pi = updates.denseVector(); 2364 const int *COIN_RESTRICT piIndex = updates.getIndices(); 2365 int *COIN_RESTRICT index = tableauRow.getIndices() + first; 2366 double *COIN_RESTRICT array = tableauRow.denseVector() + first; 2367 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 2381 2368 if (!iBlock) { 2382 numberNonZero =numberSlacks;2383 for (int i =0;i<numberSlacks;i++) {2384 int iRow =piIndex[i];2385 index[i] =iRow;2386 array[i] =pi[iRow]; // ? what if small or basic2387 } 2388 first =maximumRows;2369 numberNonZero = numberSlacks; 2370 for (int i = 0; i < numberSlacks; i++) { 2371 int iRow = piIndex[i]; 2372 index[i] = iRow; 2373 array[i] = pi[iRow]; // ? what if small or basic 2374 } 2375 first = maximumRows; 2389 2376 } 2390 2377 double zeroTolerance = model_>zeroTolerance(); 2391 2378 if (doByRow) { 2392 int numberRows =model_>numberRows();2393 const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_+iBlock*numberRows;2394 const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows;2395 const double * 2396 const int * 2397 if (numberSlacks >1) {2398 double * 2399 for (int i =0;i<numberSlacks;i++) {2400 int iRow=piIndex[i];2401 2402 2403 for (CoinBigIndex j=rowStart[iRow];j<end;j++) {2404 2405 arrayTemp[iColumn] += element[j]*piValue;2406 2407 } 2408 for (int iSequence =first;iSequence<last;iSequence++) {2409 2410 2411 arrayTemp[iSequence]=0.0;2412 if (fabs(tableauValue)>zeroTolerance) {2413 unsigned char iStatus=internalStatus[iSequence]&7;2414 if (iStatus<6) {2415 index[numberNonZero]=iSequence;2416 array[numberNonZero++]=tableauValue;2417 2418 2419 2379 int numberRows = model_>numberRows(); 2380 const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_ + iBlock * numberRows; 2381 const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows; 2382 const double *COIN_RESTRICT element = element_; 2383 const int *COIN_RESTRICT column = column_; 2384 if (numberSlacks > 1) { 2385 double *COIN_RESTRICT arrayTemp = tableauRow.denseVector(); 2386 for (int i = 0; i < numberSlacks; i++) { 2387 int iRow = piIndex[i]; 2388 double piValue = pi[iRow]; 2389 CoinBigIndex end = rowEnd[iRow]; 2390 for (CoinBigIndex j = rowStart[iRow]; j < end; j++) { 2391 int iColumn = column[j]; 2392 arrayTemp[iColumn] += element[j] * piValue; 2393 } 2394 } 2395 for (int iSequence = first; iSequence < last; iSequence++) { 2396 double tableauValue = arrayTemp[iSequence]; 2397 if (tableauValue) { 2398 arrayTemp[iSequence] = 0.0; 2399 if (fabs(tableauValue) > zeroTolerance) { 2400 unsigned char iStatus = internalStatus[iSequence] & 7; 2401 if (iStatus < 6) { 2402 index[numberNonZero] = iSequence; 2403 array[numberNonZero++] = tableauValue; 2404 } 2405 } 2406 } 2420 2407 } 2421 2408 } else { 2422 int iRow =piIndex[0];2409 int iRow = piIndex[0]; 2423 2410 double piValue = pi[iRow]; 2424 2411 CoinBigIndex end = rowEnd[iRow]; 2425 for (CoinBigIndex j =rowStart[iRow];j<end;j++) {2426 2427 double tableauValue = element[j]*piValue;2428 if (fabs(tableauValue)>zeroTolerance) {2429 unsigned char iStatus=internalStatus[iSequence]&7;2430 if (iStatus<6) {2431 index[numberNonZero]=iSequence;2432 array[numberNonZero++]=tableauValue;2433 2434 2412 for (CoinBigIndex j = rowStart[iRow]; j < end; j++) { 2413 int iSequence = column[j]; 2414 double tableauValue = element[j] * piValue; 2415 if (fabs(tableauValue) > zeroTolerance) { 2416 unsigned char iStatus = internalStatus[iSequence] & 7; 2417 if (iStatus < 6) { 2418 index[numberNonZero] = iSequence; 2419 array[numberNonZero++] = tableauValue; 2420 } 2421 } 2435 2422 } 2436 2423 } 2437 2424 } else { 2438 2425 // get matrix data pointers 2439 const int * 2440 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;2441 const double * 2442 const double * COIN_RESTRICT pi=updates.denseVector();2426 const int *COIN_RESTRICT row = matrix_>getIndices(); 2427 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 2428 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 2429 const double *COIN_RESTRICT pi = updates.denseVector(); 2443 2430 for (int iSequence = first; iSequence < last; iSequence++) { 2444 unsigned char iStatus =internalStatus[iSequence]&7;2445 if (iStatus <6) {2446 2447 CoinBigIndex next = columnStart[iSequence+1];2448 2449 2450 2451 2452 2453 if (fabs(tableauValue)>zeroTolerance) {2454 index[numberNonZero]=iSequence;2455 array[numberNonZero++]=tableauValue;2456 2457 } 2458 } 2459 } 2460 tableauRow.setNumElementsPartition(iBlock, numberNonZero);2431 unsigned char iStatus = internalStatus[iSequence] & 7; 2432 if (iStatus < 6) { 2433 CoinBigIndex start = columnStart[iSequence]; 2434 CoinBigIndex next = columnStart[iSequence + 1]; 2435 double tableauValue = 0.0; 2436 for (CoinBigIndex j = start; j < next; j++) { 2437 int jRow = row[j]; 2438 tableauValue += pi[jRow] * elementByColumn[j]; 2439 } 2440 if (fabs(tableauValue) > zeroTolerance) { 2441 index[numberNonZero] = iSequence; 2442 array[numberNonZero++] = tableauValue; 2443 } 2444 } 2445 } 2446 } 2447 tableauRow.setNumElementsPartition(iBlock, numberNonZero); 2461 2448 return numberSlacks; 2462 2449 } 2463 2450 // Get sequenceIn when Dantzig 2464 int 2465 AbcMatrix::pivotColumnDantzig(const CoinIndexedVector & updates, 2466 CoinPartitionedVector & spare) const 2467 { 2468 int maximumRows = model_>maximumAbcNumberRows(); 2469 // rebalance 2470 rebalance(); 2471 spare.setPackedMode(true); 2472 bool useRowCopy=false; 2473 if (rowStart_) { 2474 // see if worth using row copy 2475 int number=updates.getNumElements(); 2476 assert (number); 2477 useRowCopy = (number<<2)<maximumRows; 2478 } 2479 const int * starts; 2480 if (useRowCopy) 2481 starts=blockStart(); 2482 else 2483 starts=startColumnBlock(); 2451 int AbcMatrix::pivotColumnDantzig(const CoinIndexedVector &updates, 2452 CoinPartitionedVector &spare) const 2453 { 2454 int maximumRows = model_>maximumAbcNumberRows(); 2455 // rebalance 2456 rebalance(); 2457 spare.setPackedMode(true); 2458 bool useRowCopy = false; 2459 if (rowStart_) { 2460 // see if worth using row copy 2461 int number = updates.getNumElements(); 2462 assert(number); 2463 useRowCopy = (number << 2) < maximumRows; 2464 } 2465 const int *starts; 2466 if (useRowCopy) 2467 starts = blockStart(); 2468 else 2469 starts = startColumnBlock(); 2484 2470 #if ABC_PARALLEL 2485 2471 #define NUMBER_BLOCKS NUMBER_ROW_BLOCKS 2486 int numberBlocks =CoinMin(NUMBER_BLOCKS,model_>numberCpus());2472 int numberBlocks = CoinMin(NUMBER_BLOCKS, model_>numberCpus()); 2487 2473 #else 2488 2474 #define NUMBER_BLOCKS 1 2489 int numberBlocks =NUMBER_BLOCKS;2475 int numberBlocks = NUMBER_BLOCKS; 2490 2476 #endif 2491 2477 #if ABC_PARALLEL 2492 if (model_>parallelMode()==0) 2493 numberBlocks=1; 2494 #endif 2495 spare.setPartitions(numberBlocks,starts); 2496 int which[NUMBER_BLOCKS]; 2497 double best[NUMBER_BLOCKS]; 2498 for (int i=0;i<numberBlocks1;i++) 2499 which[i]=cilk_spawn pivotColumnDantzig(i,useRowCopy,updates,spare,best[i]); 2500 which[numberBlocks1]=pivotColumnDantzig(numberBlocks1,useRowCopy,updates, 2501 spare,best[numberBlocks1]); 2502 cilk_sync; 2503 int bestSequence=1; 2504 double bestValue=model_>dualTolerance(); 2505 for (int i=0;i<numberBlocks;i++) { 2506 if (best[i]>bestValue) { 2507 bestValue=best[i]; 2508 bestSequence=which[i]; 2509 } 2510 } 2511 return bestSequence; 2512 } 2513 // Get sequenceIn when Dantzig (One block) 2514 int 2515 AbcMatrix::pivotColumnDantzig(int iBlock, bool doByRow,const CoinIndexedVector & updates, 2516 CoinPartitionedVector & spare, 2517 double & bestValue) const 2518 { 2519 // could rewrite to do more directly 2520 int numberSlacks=primalColumnRow(iBlock,doByRow,updates,spare); 2521 double * COIN_RESTRICT reducedCost = model_>djRegion(); 2522 int first=spare.startPartition(iBlock); 2523 int last=spare.startPartition(iBlock+1); 2524 int * COIN_RESTRICT index = spare.getIndices()+first; 2525 double * COIN_RESTRICT array = spare.denseVector()+first; 2526 int numberNonZero=spare.getNumElements(iBlock); 2527 int bestSequence=1; 2528 double bestDj=model_>dualTolerance(); 2529 double bestFreeDj = model_>dualTolerance(); 2530 int bestFreeSequence = 1; 2531 // redo LOTS so signs for slacks and columns same 2532 if (!first) { 2533 first = model_>maximumAbcNumberRows(); 2534 for (int i=0;i<numberSlacks;i++) { 2535 int iSequence = index[i]; 2536 double value = reducedCost[iSequence]; 2537 value += array[i]; 2538 array[i] = 0.0; 2539 reducedCost[iSequence] = value; 2540 } 2541 #ifndef CLP_PRIMAL_SLACK_MULTIPLIER 2478 if (model_>parallelMode() == 0) 2479 numberBlocks = 1; 2480 #endif 2481 spare.setPartitions(numberBlocks, starts); 2482 int which[NUMBER_BLOCKS]; 2483 double best[NUMBER_BLOCKS]; 2484 for (int i = 0; i < numberBlocks  1; i++) 2485 which[i] = cilk_spawn pivotColumnDantzig(i, useRowCopy, updates, spare, best[i]); 2486 which[numberBlocks  1] = pivotColumnDantzig(numberBlocks  1, useRowCopy, updates, 2487 spare, best[numberBlocks  1]); 2488 cilk_sync; 2489 int bestSequence = 1; 2490 double bestValue = model_>dualTolerance(); 2491 for (int i = 0; i < numberBlocks; i++) { 2492 if (best[i] > bestValue) { 2493 bestValue = best[i]; 2494 bestSequence = which[i]; 2495 } 2496 } 2497 return bestSequence; 2498 } 2499 // Get sequenceIn when Dantzig (One block) 2500 int AbcMatrix::pivotColumnDantzig(int iBlock, bool doByRow, const CoinIndexedVector &updates, 2501 CoinPartitionedVector &spare, 2502 double &bestValue) const 2503 { 2504 // could rewrite to do more directly 2505 int numberSlacks = primalColumnRow(iBlock, doByRow, updates, spare); 2506 double *COIN_RESTRICT reducedCost = model_>djRegion(); 2507 int first = spare.startPartition(iBlock); 2508 int last = spare.startPartition(iBlock + 1); 2509 int *COIN_RESTRICT index = spare.getIndices() + first; 2510 double *COIN_RESTRICT array = spare.denseVector() + first; 2511 int numberNonZero = spare.getNumElements(iBlock); 2512 int bestSequence = 1; 2513 double bestDj = model_>dualTolerance(); 2514 double bestFreeDj = model_>dualTolerance(); 2515 int bestFreeSequence = 1; 2516 // redo LOTS so signs for slacks and columns same 2517 if (!first) { 2518 first = model_>maximumAbcNumberRows(); 2519 for (int i = 0; i < numberSlacks; i++) { 2520 int iSequence = index[i]; 2521 double value = reducedCost[iSequence]; 2522 value += array[i]; 2523 array[i] = 0.0; 2524 reducedCost[iSequence] = value; 2525 } 2526 #ifndef CLP_PRIMAL_SLACK_MULTIPLIER 2542 2527 #define CLP_PRIMAL_SLACK_MULTIPLIER 1.0 2543 2528 #endif 2544 int numberRows=model_>numberRows(); 2545 for (int iSequence=0 ; iSequence < numberRows; iSequence++) { 2546 // check flagged variable 2547 if (!model_>flagged(iSequence)) { 2548 double value = reducedCost[iSequence] * CLP_PRIMAL_SLACK_MULTIPLIER; 2549 AbcSimplex::Status status = model_>getInternalStatus(iSequence); 2550 2551 switch(status) { 2552 2553 case AbcSimplex::basic: 2554 case AbcSimplex::isFixed: 2555 break; 2556 case AbcSimplex::isFree: 2557 case AbcSimplex::superBasic: 2558 if (fabs(value) > bestFreeDj) { 2559 bestFreeDj = fabs(value); 2560 bestFreeSequence = iSequence; 2561 } 2562 break; 2563 case AbcSimplex::atUpperBound: 2564 if (value > bestDj) { 2565 bestDj = value; 2566 bestSequence = iSequence; 2567 } 2568 break; 2569 case AbcSimplex::atLowerBound: 2570 if (value < bestDj) { 2571 bestDj = value; 2572 bestSequence = iSequence; 2573 } 2574 } 2575 } 2576 } 2577 } 2578 for (int i=numberSlacks;i<numberNonZero;i++) { 2579 int iSequence = index[i]; 2580 double value = reducedCost[iSequence]; 2581 value += array[i]; 2582 array[i] = 0.0; 2583 reducedCost[iSequence] = value; 2584 } 2585 for (int iSequence=first ; iSequence < last; iSequence++) { 2586 // check flagged variable 2587 if (!model_>flagged(iSequence)) { 2588 double value = reducedCost[iSequence]; 2589 AbcSimplex::Status status = model_>getInternalStatus(iSequence); 2590 2591 switch(status) { 2592 2593 case AbcSimplex::basic: 2594 case AbcSimplex::isFixed: 2595 break; 2596 case AbcSimplex::isFree: 2597 case AbcSimplex::superBasic: 2598 if (fabs(value) > bestFreeDj) { 2599 bestFreeDj = fabs(value); 2600 bestFreeSequence = iSequence; 2601 } 2602 break; 2603 case AbcSimplex::atUpperBound: 2604 if (value > bestDj) { 2605 bestDj = value; 2606 bestSequence = iSequence; 2607 } 2608 break; 2609 case AbcSimplex::atLowerBound: 2610 if (value < bestDj) { 2611 bestDj = value; 2612 bestSequence = iSequence; 2613 } 2614 } 2615 } 2616 } 2617 spare.setNumElementsPartition(iBlock,0); 2618 // bias towards free 2619 if (bestFreeSequence >= 0 && bestFreeDj > 0.1 * bestDj) { 2620 bestSequence = bestFreeSequence; 2621 bestDj=10.0*bestFreeDj; 2622 } 2623 bestValue=bestDj; 2624 return bestSequence; 2625 } 2626 // gets tableau row and dj row  returns number of slacks in block 2627 int 2628 AbcMatrix::primalColumnRowAndDjs(int iBlock,const CoinIndexedVector & updateTableau, 2629 const CoinIndexedVector & updateDjs, 2630 CoinPartitionedVector & tableauRow) const 2529 int numberRows = model_>numberRows(); 2530 for (int iSequence = 0; iSequence < numberRows; iSequence++) { 2531 // check flagged variable 2532 if (!model_>flagged(iSequence)) { 2533 double value = reducedCost[iSequence] * CLP_PRIMAL_SLACK_MULTIPLIER; 2534 AbcSimplex::Status status = model_>getInternalStatus(iSequence); 2535 2536 switch (status) { 2537 2538 case AbcSimplex::basic: 2539 case AbcSimplex::isFixed: 2540 break; 2541 case AbcSimplex::isFree: 2542 case AbcSimplex::superBasic: 2543 if (fabs(value) > bestFreeDj) { 2544 bestFreeDj = fabs(value); 2545 bestFreeSequence = iSequence; 2546 } 2547 break; 2548 case AbcSimplex::atUpperBound: 2549 if (value > bestDj) { 2550 bestDj = value; 2551 bestSequence = iSequence; 2552 } 2553 break; 2554 case AbcSimplex::atLowerBound: 2555 if (value < bestDj) { 2556 bestDj = value; 2557 bestSequence = iSequence; 2558 } 2559 } 2560 } 2561 } 2562 } 2563 for (int i = numberSlacks; i < numberNonZero; i++) { 2564 int iSequence = index[i]; 2565 double value = reducedCost[iSequence]; 2566 value += array[i]; 2567 array[i] = 0.0; 2568 reducedCost[iSequence] = value; 2569 } 2570 for (int iSequence = first; iSequence < last; iSequence++) { 2571 // check flagged variable 2572 if (!model_>flagged(iSequence)) { 2573 double value = reducedCost[iSequence]; 2574 AbcSimplex::Status status = model_>getInternalStatus(iSequence); 2575 2576 switch (status) { 2577 2578 case AbcSimplex::basic: 2579 case AbcSimplex::isFixed: 2580 break; 2581 case AbcSimplex::isFree: 2582 case AbcSimplex::superBasic: 2583 if (fabs(value) > bestFreeDj) { 2584 bestFreeDj = fabs(value); 2585 bestFreeSequence = iSequence; 2586 } 2587 break; 2588 case AbcSimplex::atUpperBound: 2589 if (value > bestDj) { 2590 bestDj = value; 2591 bestSequence = iSequence; 2592 } 2593 break; 2594 case AbcSimplex::atLowerBound: 2595 if (value < bestDj) { 2596 bestDj = value; 2597 bestSequence = iSequence; 2598 } 2599 } 2600 } 2601 } 2602 spare.setNumElementsPartition(iBlock, 0); 2603 // bias towards free 2604 if (bestFreeSequence >= 0 && bestFreeDj > 0.1 * bestDj) { 2605 bestSequence = bestFreeSequence; 2606 bestDj = 10.0 * bestFreeDj; 2607 } 2608 bestValue = bestDj; 2609 return bestSequence; 2610 } 2611 // gets tableau row and dj row  returns number of slacks in block 2612 int AbcMatrix::primalColumnRowAndDjs(int iBlock, const CoinIndexedVector &updateTableau, 2613 const CoinIndexedVector &updateDjs, 2614 CoinPartitionedVector &tableauRow) const 2631 2615 { 2632 2616 int maximumRows = model_>maximumAbcNumberRows(); 2633 int first =tableauRow.startPartition(iBlock);2634 int last =tableauRow.startPartition(iBlock+1);2635 assert(first ==startColumnBlock_[iBlock]);2636 assert(last ==startColumnBlock_[iBlock+1]);2637 int numberSlacks =updateTableau.getNumElements();2638 int numberNonZero =0;2639 const double * COIN_RESTRICT piTableau=updateTableau.denseVector();2640 const double * COIN_RESTRICT pi=updateDjs.denseVector();2641 int * COIN_RESTRICT index = tableauRow.getIndices()+first;2642 double * COIN_RESTRICT array = tableauRow.denseVector()+first;2643 const unsigned char * 2644 double * 2617 int first = tableauRow.startPartition(iBlock); 2618 int last = tableauRow.startPartition(iBlock + 1); 2619 assert(first == startColumnBlock_[iBlock]); 2620 assert(last == startColumnBlock_[iBlock + 1]); 2621 int numberSlacks = updateTableau.getNumElements(); 2622 int numberNonZero = 0; 2623 const double *COIN_RESTRICT piTableau = updateTableau.denseVector(); 2624 const double *COIN_RESTRICT pi = updateDjs.denseVector(); 2625 int *COIN_RESTRICT index = tableauRow.getIndices() + first; 2626 double *COIN_RESTRICT array = tableauRow.denseVector() + first; 2627 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 2628 double *COIN_RESTRICT reducedCost = model_>djRegion(); 2645 2629 if (!iBlock) { 2646 const int * 2647 numberNonZero =numberSlacks;2648 for (int i =0;i<numberSlacks;i++) {2649 int iRow =piIndexTableau[i];2650 index[i] =iRow;2651 array[i] =piTableau[iRow]; // ? what if small or basic2652 } 2653 const int * 2654 int number =updateDjs.getNumElements();2655 for (int i =0;i<number;i++) {2656 int iRow =piIndex[i];2630 const int *COIN_RESTRICT piIndexTableau = updateTableau.getIndices(); 2631 numberNonZero = numberSlacks; 2632 for (int i = 0; i < numberSlacks; i++) { 2633 int iRow = piIndexTableau[i]; 2634 index[i] = iRow; 2635 array[i] = piTableau[iRow]; // ? what if small or basic 2636 } 2637 const int *COIN_RESTRICT piIndex = updateDjs.getIndices(); 2638 int number = updateDjs.getNumElements(); 2639 for (int i = 0; i < number; i++) { 2640 int iRow = piIndex[i]; 2657 2641 reducedCost[iRow] = pi[iRow]; // sign? 2658 2642 } 2659 first =maximumRows;2643 first = maximumRows; 2660 2644 } 2661 2645 double zeroTolerance = model_>zeroTolerance(); 2662 2646 // get matrix data pointers 2663 const int * 2664 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;2665 const double * 2647 const int *COIN_RESTRICT row = matrix_>getIndices(); 2648 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 2649 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 2666 2650 for (int iSequence = first; iSequence < last; iSequence++) { 2667 unsigned char iStatus =internalStatus[iSequence]&7;2668 if (iStatus <6) {2651 unsigned char iStatus = internalStatus[iSequence] & 7; 2652 if (iStatus < 6) { 2669 2653 CoinBigIndex start = columnStart[iSequence]; 2670 CoinBigIndex next = columnStart[iSequence +1];2654 CoinBigIndex next = columnStart[iSequence + 1]; 2671 2655 double tableauValue = 0.0; 2672 double djValue =reducedCost[iSequence];2656 double djValue = reducedCost[iSequence]; 2673 2657 for (CoinBigIndex j = start; j < next; j++) { 2674 2675 2676 2677 } 2678 reducedCost[iSequence] =djValue;2679 if (fabs(tableauValue) >zeroTolerance) {2680 index[numberNonZero]=iSequence;2681 array[numberNonZero++]=tableauValue;2682 } 2683 } 2684 } 2685 tableauRow.setNumElementsPartition(iBlock, numberNonZero);2658 int jRow = row[j]; 2659 tableauValue += piTableau[jRow] * elementByColumn[j]; 2660 djValue = pi[jRow] * elementByColumn[j]; // sign? 2661 } 2662 reducedCost[iSequence] = djValue; 2663 if (fabs(tableauValue) > zeroTolerance) { 2664 index[numberNonZero] = iSequence; 2665 array[numberNonZero++] = tableauValue; 2666 } 2667 } 2668 } 2669 tableauRow.setNumElementsPartition(iBlock, numberNonZero); 2686 2670 return numberSlacks; 2687 2671 } … … 2690 2674 otherwise use updateForDjs 2691 2675 */ 2692 int 2693 AbcMatrix::primalColumnDouble(int iBlock,CoinPartitionedVector & updateForTableauRow, 2694 CoinPartitionedVector & updateForDjs, 2695 const CoinIndexedVector & updateForWeights, 2696 CoinPartitionedVector & spareColumn1, 2697 double * infeasibilities, 2698 double referenceIn, double devex, 2699 // Array for exact devex to say what is in reference framework 2700 unsigned int * reference, 2701 double * weights, double scaleFactor) const 2676 int AbcMatrix::primalColumnDouble(int iBlock, CoinPartitionedVector &updateForTableauRow, 2677 CoinPartitionedVector &updateForDjs, 2678 const CoinIndexedVector &updateForWeights, 2679 CoinPartitionedVector &spareColumn1, 2680 double *infeasibilities, 2681 double referenceIn, double devex, 2682 // Array for exact devex to say what is in reference framework 2683 unsigned int *reference, 2684 double *weights, double scaleFactor) const 2702 2685 { 2703 2686 int maximumRows = model_>maximumAbcNumberRows(); 2704 int first =startColumnBlock_[iBlock];2705 int last =startColumnBlock_[iBlock+1];2706 const double * COIN_RESTRICT piTableau=updateForTableauRow.denseVector();2707 const double * COIN_RESTRICT pi=updateForDjs.denseVector();2708 const double * COIN_RESTRICT piWeights=updateForWeights.denseVector();2709 const unsigned char * 2710 double * 2687 int first = startColumnBlock_[iBlock]; 2688 int last = startColumnBlock_[iBlock + 1]; 2689 const double *COIN_RESTRICT piTableau = updateForTableauRow.denseVector(); 2690 const double *COIN_RESTRICT pi = updateForDjs.denseVector(); 2691 const double *COIN_RESTRICT piWeights = updateForWeights.denseVector(); 2692 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 2693 double *COIN_RESTRICT reducedCost = model_>djRegion(); 2711 2694 double tolerance = model_>currentDualTolerance(); 2712 2695 // use spareColumn to track new infeasibilities 2713 int * COIN_RESTRICT newInfeas = spareColumn1.getIndices()+first;2714 int numberNewInfeas =0;2696 int *COIN_RESTRICT newInfeas = spareColumn1.getIndices() + first; 2697 int numberNewInfeas = 0; 2715 2698 // we can't really trust infeasibilities if there is dual error 2716 2699 // this coding has to mimic coding in checkDualSolution 2717 2700 double error = CoinMin(1.0e2, model_>largestDualError()); 2718 2701 // allow tolerance at least slightly bigger than standard 2719 tolerance = tolerance + 2720 double mult[2] ={1.0,1.0};2721 bool updateWeights =devex!=0.0;2722 #define isReference(i) (((reference[i >>5] >> (i & 31)) & 1) != 0)2702 tolerance = tolerance + error; 2703 double mult[2] = { 1.0, 1.0 }; 2704 bool updateWeights = devex != 0.0; 2705 #define isReference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0) 2723 2706 // Scale factor is other way round so tableau row is 1.0* while dj update is scaleFactor* 2724 2707 if (!iBlock) { 2725 int numberSlacks =updateForTableauRow.getNumElements();2726 const int * 2708 int numberSlacks = updateForTableauRow.getNumElements(); 2709 const int *COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices(); 2727 2710 if (!scaleFactor) { 2728 const int * 2729 int number =updateForDjs.getNumElements();2730 for (int i =0;i<number;i++) {2731 int iRow=piIndex[i];2732 int iStatus=internalStatus[iRow]&7;2733 double value=reducedCost[iRow];2734 2735 if (iStatus<4) {2736 reducedCost[iRow]=value;2737 2738 if (value<0.0) {2739 2740 newInfeas[numberNewInfeas++]=iRow;2741 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 2742 infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;2711 const int *COIN_RESTRICT piIndex = updateForDjs.getIndices(); 2712 int number = updateForDjs.getNumElements(); 2713 for (int i = 0; i < number; i++) { 2714 int iRow = piIndex[i]; 2715 int iStatus = internalStatus[iRow] & 7; 2716 double value = reducedCost[iRow]; 2717 value += pi[iRow]; 2718 if (iStatus < 4) { 2719 reducedCost[iRow] = value; 2720 value *= mult[iStatus]; 2721 if (value < 0.0) { 2722 if (!infeasibilities[iRow]) 2723 newInfeas[numberNewInfeas++] = iRow; 2724 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 2725 infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER; 2743 2726 #else 2744 infeasibilities[iRow]=value*value;2745 #endif 2746 2747 2748 infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;2749 2750 2751 } 2752 } else if (scaleFactor !=COIN_DBL_MAX) {2753 for (int i =0;i<numberSlacks;i++) {2754 int iRow=piIndexTableau[i];2755 int iStatus=internalStatus[iRow]&7;2756 if (iStatus<4) {2757 double value=reducedCost[iRow];2758 value += scaleFactor*piTableau[iRow]; // sign?2759 reducedCost[iRow]=value;2760 2761 if (value<0.0) {2762 2763 newInfeas[numberNewInfeas++]=iRow;2764 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 2765 infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;2727 infeasibilities[iRow] = value * value; 2728 #endif 2729 } else { 2730 if (infeasibilities[iRow]) 2731 infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 2732 } 2733 } 2734 } 2735 } else if (scaleFactor != COIN_DBL_MAX) { 2736 for (int i = 0; i < numberSlacks; i++) { 2737 int iRow = piIndexTableau[i]; 2738 int iStatus = internalStatus[iRow] & 7; 2739 if (iStatus < 4) { 2740 double value = reducedCost[iRow]; 2741 value += scaleFactor * piTableau[iRow]; // sign? 2742 reducedCost[iRow] = value; 2743 value *= mult[iStatus]; 2744 if (value < 0.0) { 2745 if (!infeasibilities[iRow]) 2746 newInfeas[numberNewInfeas++] = iRow; 2747 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 2748 infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER; 2766 2749 #else 2767 infeasibilities[iRow]=value*value;2768 #endif 2769 2770 2771 infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;2772 2773 2750 infeasibilities[iRow] = value * value; 2751 #endif 2752 } else { 2753 if (infeasibilities[iRow]) 2754 infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 2755 } 2756 } 2774 2757 } 2775 2758 } 2776 2759 if (updateWeights) { 2777 for (int i =0;i<numberSlacks;i++) {2778 int iRow=piIndexTableau[i];2779 double modification=piWeights[iRow];2780 2781 2782 2783 2784 2785 2786 2787 2788 2789 2790 2791 2792 2793 2794 2795 2796 2797 } 2798 } 2799 first =maximumRows;2760 for (int i = 0; i < numberSlacks; i++) { 2761 int iRow = piIndexTableau[i]; 2762 double modification = piWeights[iRow]; 2763 double thisWeight = weights[iRow]; 2764 double pivot = piTableau[iRow]; 2765 double pivotSquared = pivot * pivot; 2766 thisWeight += pivotSquared * devex  pivot * modification; 2767 if (thisWeight < DEVEX_TRY_NORM) { 2768 if (referenceIn < 0.0) { 2769 // steepest 2770 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); 2771 } else { 2772 // exact 2773 thisWeight = referenceIn * pivotSquared; 2774 if (isReference(iRow)) 2775 thisWeight += 1.0; 2776 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); 2777 } 2778 } 2779 weights[iRow] = thisWeight; 2780 } 2781 } 2782 first = maximumRows; 2800 2783 } 2801 2784 double zeroTolerance = model_>zeroTolerance(); 2802 double freeTolerance = FREE_ACCEPT * tolerance;; 2785 double freeTolerance = FREE_ACCEPT * tolerance; 2786 ; 2803 2787 // get matrix data pointers 2804 const int * 2805 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;2806 const double * 2807 bool updateDjs =scaleFactor!=COIN_DBL_MAX;2788 const int *COIN_RESTRICT row = matrix_>getIndices(); 2789 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 2790 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 2791 bool updateDjs = scaleFactor != COIN_DBL_MAX; 2808 2792 for (int iSequence = first; iSequence < last; iSequence++) { 2809 unsigned char iStatus =internalStatus[iSequence]&7;2810 if (iStatus <6) {2793 unsigned char iStatus = internalStatus[iSequence] & 7; 2794 if (iStatus < 6) { 2811 2795 CoinBigIndex start = columnStart[iSequence]; 2812 CoinBigIndex next = columnStart[iSequence +1];2796 CoinBigIndex next = columnStart[iSequence + 1]; 2813 2797 double tableauValue = 0.0; 2814 double djValue =reducedCost[iSequence];2798 double djValue = reducedCost[iSequence]; 2815 2799 if (!scaleFactor) { 2816 2817 2818 2819 2820 2800 for (CoinBigIndex j = start; j < next; j++) { 2801 int jRow = row[j]; 2802 tableauValue += piTableau[jRow] * elementByColumn[j]; 2803 djValue += pi[jRow] * elementByColumn[j]; 2804 } 2821 2805 } else { 2822 2823 2824 2825 2826 djValue += tableauValue*scaleFactor; // sign?2806 for (CoinBigIndex j = start; j < next; j++) { 2807 int jRow = row[j]; 2808 tableauValue += piTableau[jRow] * elementByColumn[j]; 2809 } 2810 djValue += tableauValue * scaleFactor; // sign? 2827 2811 } 2828 2812 if (updateDjs) { 2829 reducedCost[iSequence]=djValue;2830 if (iStatus<4) {2831 2832 if (djValue<0.0) {2833 2834 newInfeas[numberNewInfeas++]=iSequence;2835 infeasibilities[iSequence]=djValue*djValue;2836 2837 2838 infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;2839 2840 2841 2842 2843 2844 2845 newInfeas[numberNewInfeas++]=iSequence;2846 2847 2848 2849 2850 infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;2851 2852 2853 } 2854 if (fabs(tableauValue) >zeroTolerance&&updateWeights) {2855 double modification=0.0;2856 2857 2858 2859 2860 2861 2862 2863 2864 2865 2866 2867 2868 2869 2870 2871 2872 2873 2874 2875 2876 2877 } 2878 } 2879 } 2880 spareColumn1.setTempNumElementsPartition(iBlock, numberNewInfeas);2881 int bestSequence =1;2813 reducedCost[iSequence] = djValue; 2814 if (iStatus < 4) { 2815 djValue *= mult[iStatus]; 2816 if (djValue < 0.0) { 2817 if (!infeasibilities[iSequence]) 2818 newInfeas[numberNewInfeas++] = iSequence; 2819 infeasibilities[iSequence] = djValue * djValue; 2820 } else { 2821 if (infeasibilities[iSequence]) 2822 infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT; 2823 } 2824 } else { 2825 if (fabs(djValue) > freeTolerance) { 2826 // we are going to bias towards free (but only if reasonable) 2827 djValue *= FREE_BIAS; 2828 if (!infeasibilities[iSequence]) 2829 newInfeas[numberNewInfeas++] = iSequence; 2830 // store square in list 2831 infeasibilities[iSequence] = djValue * djValue; 2832 } else { 2833 if (infeasibilities[iSequence]) 2834 infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT; 2835 } 2836 } 2837 } 2838 if (fabs(tableauValue) > zeroTolerance && updateWeights) { 2839 double modification = 0.0; 2840 for (CoinBigIndex j = start; j < next; j++) { 2841 int jRow = row[j]; 2842 modification += piWeights[jRow] * elementByColumn[j]; 2843 } 2844 double thisWeight = weights[iSequence]; 2845 double pivot = tableauValue; 2846 double pivotSquared = pivot * pivot; 2847 thisWeight += pivotSquared * devex  pivot * modification; 2848 if (thisWeight < DEVEX_TRY_NORM) { 2849 if (referenceIn < 0.0) { 2850 // steepest 2851 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); 2852 } else { 2853 // exact 2854 thisWeight = referenceIn * pivotSquared; 2855 if (isReference(iSequence)) 2856 thisWeight += 1.0; 2857 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); 2858 } 2859 } 2860 weights[iSequence] = thisWeight; 2861 } 2862 } 2863 } 2864 spareColumn1.setTempNumElementsPartition(iBlock, numberNewInfeas); 2865 int bestSequence = 1; 2882 2866 #if 0 2883 2867 if (!iBlock) … … 2900 2884 otherwise use updateForDjs 2901 2885 */ 2902 int 2903 AbcMatrix::primalColumnSparseDouble(int iBlock,CoinPartitionedVector & updateForTableauRow, 2904 CoinPartitionedVector & updateForDjs, 2905 const CoinIndexedVector & updateForWeights, 2906 CoinPartitionedVector & spareColumn1, 2907 double * infeasibilities, 2908 double referenceIn, double devex, 2909 // Array for exact devex to say what is in reference framework 2910 unsigned int * reference, 2911 double * weights, double scaleFactor) const 2886 int AbcMatrix::primalColumnSparseDouble(int iBlock, CoinPartitionedVector &updateForTableauRow, 2887 CoinPartitionedVector &updateForDjs, 2888 const CoinIndexedVector &updateForWeights, 2889 CoinPartitionedVector &spareColumn1, 2890 double *infeasibilities, 2891 double referenceIn, double devex, 2892 // Array for exact devex to say what is in reference framework 2893 unsigned int *reference, 2894 double *weights, double scaleFactor) const 2912 2895 { 2913 2896 int maximumRows = model_>maximumAbcNumberRows(); 2914 int first =blockStart_[iBlock];2897 int first = blockStart_[iBlock]; 2915 2898 //int last=blockStart_[iBlock+1]; 2916 const double * COIN_RESTRICT piTableau=updateForTableauRow.denseVector();2899 const double *COIN_RESTRICT piTableau = updateForTableauRow.denseVector(); 2917 2900 //const double * COIN_RESTRICT pi=updateForDjs.denseVector(); 2918 const double * COIN_RESTRICT piWeights=updateForWeights.denseVector();2919 const unsigned char * 2920 double * 2901 const double *COIN_RESTRICT piWeights = updateForWeights.denseVector(); 2902 const unsigned char *COIN_RESTRICT internalStatus = model_>internalStatus(); 2903 double *COIN_RESTRICT reducedCost = model_>djRegion(); 2921 2904 double tolerance = model_>currentDualTolerance(); 2922 2905 // use spareColumn to track new infeasibilities 2923 int * COIN_RESTRICT newInfeas = spareColumn1.getIndices()+first;2924 int numberNewInfeas =0;2906 int *COIN_RESTRICT newInfeas = spareColumn1.getIndices() + first; 2907 int numberNewInfeas = 0; 2925 2908 // we can't really trust infeasibilities if there is dual error 2926 2909 // this coding has to mimic coding in checkDualSolution 2927 2910 double error = CoinMin(1.0e2, model_>largestDualError()); 2928 2911 // allow tolerance at least slightly bigger than standard 2929 tolerance = tolerance + 2930 double mult[2] ={1.0,1.0};2931 bool updateWeights =devex!=0.0;2932 int numberSlacks =updateForTableauRow.getNumElements();2933 const int * 2934 #define isReference(i) (((reference[i >>5] >> (i & 31)) & 1) != 0)2912 tolerance = tolerance + error; 2913 double mult[2] = { 1.0, 1.0 }; 2914 bool updateWeights = devex != 0.0; 2915 int numberSlacks = updateForTableauRow.getNumElements(); 2916 const int *COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices(); 2917 #define isReference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0) 2935 2918 // Scale factor is other way round so tableau row is 1.0* while dj update is scaleFactor* 2936 assert 2919 assert(scaleFactor); 2937 2920 if (!iBlock) { 2938 if (scaleFactor !=COIN_DBL_MAX) {2939 for (int i =0;i<numberSlacks;i++) {2940 int iRow=piIndexTableau[i];2941 int iStatus=internalStatus[iRow]&7;2942 if (iStatus<4) {2943 double value=reducedCost[iRow];2944 value += scaleFactor*piTableau[iRow]; // sign?2945 reducedCost[iRow]=value;2946 2947 if (value<0.0) {2948 2949 newInfeas[numberNewInfeas++]=iRow;2950 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 2951 infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;2921 if (scaleFactor != COIN_DBL_MAX) { 2922 for (int i = 0; i < numberSlacks; i++) { 2923 int iRow = piIndexTableau[i]; 2924 int iStatus = internalStatus[iRow] & 7; 2925 if (iStatus < 4) { 2926 double value = reducedCost[iRow]; 2927 value += scaleFactor * piTableau[iRow]; // sign? 2928 reducedCost[iRow] = value; 2929 value *= mult[iStatus]; 2930 if (value < 0.0) { 2931 if (!infeasibilities[iRow]) 2932 newInfeas[numberNewInfeas++] = iRow; 2933 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER 2934 infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER; 2952 2935 #else 2953 infeasibilities[iRow]=value*value;2954 #endif 2955 2956 2957 infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;2958 2959 2936 infeasibilities[iRow] = value * value; 2937 #endif 2938 } else { 2939 if (infeasibilities[iRow]) 2940 infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT; 2941 } 2942 } 2960 2943 } 2961 2944 } 2962 2945 if (updateWeights) { 2963 for (int i =0;i<numberSlacks;i++) {2964 int iRow=piIndexTableau[i];2965 double modification=piWeights[iRow];2966 2967 2968 2969 2970 2971 2972 2973 2974 2975 2976 2977 2978 2979 2980 2981 2982 2983 } 2984 } 2985 first =maximumRows;2946 for (int i = 0; i < numberSlacks; i++) { 2947 int iRow = piIndexTableau[i]; 2948 double modification = piWeights[iRow]; 2949 double thisWeight = weights[iRow]; 2950 double pivot = piTableau[iRow]; 2951 double pivotSquared = pivot * pivot; 2952 thisWeight += pivotSquared * devex  pivot * modification; 2953 if (thisWeight < DEVEX_TRY_NORM) { 2954 if (referenceIn < 0.0) { 2955 // steepest 2956 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); 2957 } else { 2958 // exact 2959 thisWeight = referenceIn * pivotSquared; 2960 if (isReference(iRow)) 2961 thisWeight += 1.0; 2962 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); 2963 } 2964 } 2965 weights[iRow] = thisWeight; 2966 } 2967 } 2968 first = maximumRows; 2986 2969 } 2987 2970 double zeroTolerance = model_>zeroTolerance(); 2988 double freeTolerance = FREE_ACCEPT * tolerance;; 2971 double freeTolerance = FREE_ACCEPT * tolerance; 2972 ; 2989 2973 // get matrix data pointers 2990 const int * 2991 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;2992 const double * 2974 const int *COIN_RESTRICT row = matrix_>getIndices(); 2975 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 2976 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 2993 2977 // get tableau row 2994 int * COIN_RESTRICT index = updateForTableauRow.getIndices()+first;2995 double * 2996 int numberRows =model_>numberRows();2997 const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_+iBlock*numberRows;2998 const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows;2999 const double * 3000 const int * 3001 int numberInRow =0;3002 if (numberSlacks >2) {3003 for (int i =0;i<numberSlacks;i++) {3004 int iRow =piIndexTableau[i];2978 int *COIN_RESTRICT index = updateForTableauRow.getIndices() + first; 2979 double *COIN_RESTRICT array = updateForTableauRow.denseVector(); 2980 int numberRows = model_>numberRows(); 2981 const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_ + iBlock * numberRows; 2982 const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows; 2983 const double *COIN_RESTRICT element = element_; 2984 const int *COIN_RESTRICT column = column_; 2985 int numberInRow = 0; 2986 if (numberSlacks > 2) { 2987 for (int i = 0; i < numberSlacks; i++) { 2988 int iRow = piIndexTableau[i]; 3005 2989 double piValue = piTableau[iRow]; 3006 2990 CoinBigIndex end = rowEnd[iRow]; 3007 for (CoinBigIndex j =rowStart[iRow];j<end;j++) {3008 3009 double value = element[j]*piValue;3010 double oldValue=array[iSequence];3011 3012 3013 array[iSequence]= value; 3014 index[numberInRow++]=iSequence;3015 } else if (value) { 3016 array[iSequence] = value; 3017 3018 3019 3020 } 3021 } 3022 } else if (numberSlacks ==2) {2991 for (CoinBigIndex j = rowStart[iRow]; j < end; j++) { 2992 int iSequence = column[j]; 2993 double value = element[j] * piValue; 2994 double oldValue = array[iSequence]; 2995 value += oldValue; 2996 if (!oldValue) { 2997 array[iSequence] = value; 2998 index[numberInRow++] = iSequence; 2999 } else if (value) { 3000 array[iSequence] = value; 3001 } else { 3002 array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT; 3003 } 3004 } 3005 } 3006 } else if (numberSlacks == 2) { 3023 3007 int iRow0 = piIndexTableau[0]; 3024 3008 int iRow1 = piIndexTableau[1]; 3025 3009 CoinBigIndex end0 = rowEnd[iRow0]; 3026 3010 CoinBigIndex end1 = rowEnd[iRow1]; 3027 if (end0 rowStart[iRow0]>end1rowStart[iRow1]) {3028 int temp =iRow0;3029 iRow0 =iRow1;3030 iRow1 =temp;3011 if (end0  rowStart[iRow0] > end1  rowStart[iRow1]) { 3012 int temp = iRow0; 3013 iRow0 = iRow1; 3014 iRow1 = temp; 3031 3015 } 3032 3016 CoinBigIndex start = rowStart[iRow0]; 3033 3017 CoinBigIndex end = rowEnd[iRow0]; 3034 3018 double piValue = piTableau[iRow0]; 3035 for (CoinBigIndex j =start;j<end;j++) {3019 for (CoinBigIndex j = start; j < end; j++) { 3036 3020 int iSequence = column[j]; 3037 3021 double value = element[j]; 3038 array[iSequence] = piValue*value;3039 index[numberInRow++] =iSequence;3022 array[iSequence] = piValue * value; 3023 index[numberInRow++] = iSequence; 3040 3024 } 3041 3025 start = rowStart[iRow1]; 3042 3026 end = rowEnd[iRow1]; 3043 3027 piValue = piTableau[iRow1]; 3044 for (CoinBigIndex j =start;j<end;j++) {3028 for (CoinBigIndex j = start; j < end; j++) { 3045 3029 int iSequence = column[j]; 3046 3030 double value = element[j]; 3047 3031 value *= piValue; 3048 3032 if (!array[iSequence]) { 3049 array[iSequence]= value; 3050 index[numberInRow++]=iSequence;3051 } else { 3052 3053 3054 array[iSequence]= value; 3055 3056 3033 array[iSequence] = value; 3034 index[numberInRow++] = iSequence; 3035 } else { 3036 value += array[iSequence]; 3037 if (value) 3038 array[iSequence] = value; 3039 else 3040 array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT; 3057 3041 } 3058 3042 } … … 3062 3046 CoinBigIndex end = rowEnd[iRow0]; 3063 3047 double piValue = piTableau[iRow0]; 3064 for (CoinBigIndex j =start;j<end;j++) {3048 for (CoinBigIndex j = start; j < end; j++) { 3065 3049 int iSequence = column[j]; 3066 3050 double value = element[j]; 3067 array[iSequence] = piValue*value;3068 index[numberInRow++] =iSequence;3069 } 3070 } 3071 bool updateDjs =scaleFactor!=COIN_DBL_MAX;3072 for (int iLook = 0 ;iLook<numberInRow;iLook++) {3073 int iSequence =index[iLook];3074 unsigned char iStatus =internalStatus[iSequence]&7;3075 if (iStatus <6) {3051 array[iSequence] = piValue * value; 3052 index[numberInRow++] = iSequence; 3053 } 3054 } 3055 bool updateDjs = scaleFactor != COIN_DBL_MAX; 3056 for (int iLook = 0; iLook < numberInRow; iLook++) { 3057 int iSequence = index[iLook]; 3058 unsigned char iStatus = internalStatus[iSequence] & 7; 3059 if (iStatus < 6) { 3076 3060 double tableauValue = array[iSequence]; 3077 array[iSequence] =0.0;3078 double djValue =reducedCost[iSequence];3079 djValue += tableauValue *scaleFactor; // sign?3061 array[iSequence] = 0.0; 3062 double djValue = reducedCost[iSequence]; 3063 djValue += tableauValue * scaleFactor; // sign? 3080 3064 if (updateDjs) { 3081 reducedCost[iSequence]=djValue;3082 if (iStatus<4) {3083 3084 if (djValue<0.0) {3085 3086 newInfeas[numberNewInfeas++]=iSequence;3087 infeasibilities[iSequence]=djValue*djValue;3088 3089 3090 infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;3091 3092 3093 3094 3095 3096 3097 newInfeas[numberNewInfeas++]=iSequence;3098 3099 3100 3101 3102 infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;3103 3104 3105 } 3106 if (fabs(tableauValue) >zeroTolerance&&updateWeights) {3107 3108 CoinBigIndex next = columnStart[iSequence+1];3109 double modification=0.0;3110 3111 3112 3113 3114 3115 3116 3117 3118 3119 3120 3121 3122 3123 3124 3125 3126 3127 3128 3129 3130 3065 reducedCost[iSequence] = djValue; 3066 if (iStatus < 4) { 3067 djValue *= mult[iStatus]; 3068 if (djValue < 0.0) { 3069 if (!infeasibilities[iSequence]) 3070 newInfeas[numberNewInfeas++] = iSequence; 3071 infeasibilities[iSequence] = djValue * djValue; 3072 } else { 3073 if (infeasibilities[iSequence]) 3074 infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT; 3075 } 3076 } else { 3077 if (fabs(djValue) > freeTolerance) { 3078 // we are going to bias towards free (but only if reasonable) 3079 djValue *= FREE_BIAS; 3080 if (!infeasibilities[iSequence]) 3081 newInfeas[numberNewInfeas++] = iSequence; 3082 // store square in list 3083 infeasibilities[iSequence] = djValue * djValue; 3084 } else { 3085 if (infeasibilities[iSequence]) 3086 infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT; 3087 } 3088 } 3089 } 3090 if (fabs(tableauValue) > zeroTolerance && updateWeights) { 3091 CoinBigIndex start = columnStart[iSequence]; 3092 CoinBigIndex next = columnStart[iSequence + 1]; 3093 double modification = 0.0; 3094 for (CoinBigIndex j = start; j < next; j++) { 3095 int jRow = row[j]; 3096 modification += piWeights[jRow] * elementByColumn[j]; 3097 } 3098 double thisWeight = weights[iSequence]; 3099 double pivot = tableauValue; 3100 double pivotSquared = pivot * pivot; 3101 thisWeight += pivotSquared * devex  pivot * modification; 3102 if (thisWeight < DEVEX_TRY_NORM) { 3103 if (referenceIn < 0.0) { 3104 // steepest 3105 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); 3106 } else { 3107 // exact 3108 thisWeight = referenceIn * pivotSquared; 3109 if (isReference(iSequence)) 3110 thisWeight += 1.0; 3111 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); 3112 } 3113 } 3114 weights[iSequence] = thisWeight; 3131 3115 } 3132 3116 } else { 3133 array[iSequence] =0.0;3134 } 3135 } 3136 spareColumn1.setTempNumElementsPartition(iBlock, numberNewInfeas);3137 int bestSequence =1;3117 array[iSequence] = 0.0; 3118 } 3119 } 3120 spareColumn1.setTempNumElementsPartition(iBlock, numberNewInfeas); 3121 int bestSequence = 1; 3138 3122 #if 0 3139 3123 if (!iBlock) … … 3167 3151 Returns best 3168 3152 */ 3169 int 3170 AbcMatrix::primalColumnDouble(CoinPartitionedVector & updateForTableauRow, 3171 CoinPartitionedVector & updateForDjs, 3172 const CoinIndexedVector & updateForWeights, 3173 CoinPartitionedVector & spareColumn1, 3174 CoinIndexedVector & infeasible, 3175 double referenceIn, double devex, 3176 // Array for exact devex to say what is in reference framework 3177 unsigned int * reference, 3178 double * weights, double scaleFactor) const 3153 int AbcMatrix::primalColumnDouble(CoinPartitionedVector &updateForTableauRow, 3154 CoinPartitionedVector &updateForDjs, 3155 const CoinIndexedVector &updateForWeights, 3156 CoinPartitionedVector &spareColumn1, 3157 CoinIndexedVector &infeasible, 3158 double referenceIn, double devex, 3159 // Array for exact devex to say what is in reference framework 3160 unsigned int *reference, 3161 double *weights, double scaleFactor) const 3179 3162 { 3180 3163 //int maximumRows = model_>maximumAbcNumberRows(); 3181 3182 3164 // rebalance 3165 rebalance(); 3183 3166 #ifdef PRICE_IN_ABC_MATRIX 3184 3185 #endif 3186 double * infeasibilities=infeasible.denseVector();3187 int bestSequence=1;3188 3189 int maximumRows=model_>maximumAbcNumberRows();3190 int number=updateForTableauRow.getNumElements();3167 int which[NUMBER_BLOCKS]; 3168 #endif 3169 double *infeasibilities = infeasible.denseVector(); 3170 int bestSequence = 1; 3171 // see if worth using row copy 3172 int maximumRows = model_>maximumAbcNumberRows(); 3173 int number = updateForTableauRow.getNumElements(); 3191 3174 #ifdef GCC_4_9 3192 assert(number);3175 assert(number); 3193 3176 #else 3194 3195 3196 3197 #endif 3198 bool useRowCopy = (gotRowCopy()&&(number<<2)<maximumRows);3199 3200 3201 useRowCopy=false; // look again later3202 const int *starts;3203 3204 3205 starts=blockStart_;3206 numberBlocks=numberRowBlocks_;3207 3208 starts=startColumnBlock_;3209 numberBlocks=numberColumnBlocks_;3210 3211 3212 for (int i=0;i<numberBlocks;i++)3177 if (!number) { 3178 printf("Null tableau row!\n"); 3179 } 3180 #endif 3181 bool useRowCopy = (gotRowCopy() && (number << 2) < maximumRows); 3182 //useRowCopy=false; 3183 if (!scaleFactor) 3184 useRowCopy = false; // look again later 3185 const int *starts; 3186 int numberBlocks; 3187 if (useRowCopy) { 3188 starts = blockStart_; 3189 numberBlocks = numberRowBlocks_; 3190 } else { 3191 starts = startColumnBlock_; 3192 numberBlocks = numberColumnBlocks_; 3193 } 3194 if (useRowCopy) { 3195 for (int i = 0; i < numberBlocks; i++) 3213 3196 #ifdef PRICE_IN_ABC_MATRIX 3214 which[i]=3215 #endif 3216 cilk_spawn primalColumnSparseDouble(i,updateForTableauRow,updateForDjs,updateForWeights,3217 3218 infeasibilities,referenceIn,devex,reference,weights,scaleFactor);3219 3220 3221 for (int i=0;i<numberBlocks;i++)3197 which[i] = 3198 #endif 3199 cilk_spawn primalColumnSparseDouble(i, updateForTableauRow, updateForDjs, updateForWeights, 3200 spareColumn1, 3201 infeasibilities, referenceIn, devex, reference, weights, scaleFactor); 3202 cilk_sync; 3203 } else { 3204 for (int i = 0; i < numberBlocks; i++) 3222 3205 #ifdef PRICE_IN_ABC_MATRIX 3223 which[i]=3224 #endif 3225 cilk_spawn primalColumnDouble(i,updateForTableauRow,updateForDjs,updateForWeights,3226 3227 infeasibilities,referenceIn,devex,reference,weights,scaleFactor);3228 3229 3206 which[i] = 3207 #endif 3208 cilk_spawn primalColumnDouble(i, updateForTableauRow, updateForDjs, updateForWeights, 3209 spareColumn1, 3210 infeasibilities, referenceIn, devex, reference, weights, scaleFactor); 3211 cilk_sync; 3212 } 3230 3213 #ifdef PRICE_IN_ABC_MATRIX 3231 double bestValue=model_>dualTolerance();3232 int sequenceIn[8]={1,1,1,1,1,1,1,1};3233 #endif 3234 assert (numberColumnBlocks_<=8);3214 double bestValue = model_>dualTolerance(); 3215 int sequenceIn[8] = { 1, 1, 1, 1, 1, 1, 1, 1 }; 3216 #endif 3217 assert(numberColumnBlocks_ <= 8); 3235 3218 #ifdef PRICE_IN_ABC_MATRIX 3236 3237 int nPossible=0;3238 #endif 3239 3240 3241 int numberInfeas=infeasible.getNumElements();3242 int *COIN_RESTRICT infeas = infeasible.getIndices();3243 const int *COIN_RESTRICT newInfeasAll = spareColumn1.getIndices();3244 for (int i=0;i<numberColumnBlocks_;i++) {3245 const int * COIN_RESTRICT newInfeas = newInfeasAll+starts[i];3246 int numberNewInfeas=spareColumn1.getNumElements(i);3247 spareColumn1.setTempNumElementsPartition(i,0);3248 for (int j=0;j<numberNewInfeas;j++)3249 infeas[numberInfeas++]=newInfeas[j];3219 double weightedDj[8]; 3220 int nPossible = 0; 3221 #endif 3222 //const double * reducedCost = model_>djRegion(); 3223 // use spareColumn to track new infeasibilities 3224 int numberInfeas = infeasible.getNumElements(); 3225 int *COIN_RESTRICT infeas = infeasible.getIndices(); 3226 const int *COIN_RESTRICT newInfeasAll = spareColumn1.getIndices(); 3227 for (int i = 0; i < numberColumnBlocks_; i++) { 3228 const int *COIN_RESTRICT newInfeas = newInfeasAll + starts[i]; 3229 int numberNewInfeas = spareColumn1.getNumElements(i); 3230 spareColumn1.setTempNumElementsPartition(i, 0); 3231 for (int j = 0; j < numberNewInfeas; j++) 3232 infeas[numberInfeas++] = newInfeas[j]; 3250 3233 #ifdef PRICE_IN_ABC_MATRIX 3251 if (which[i]>=0) {3252 int iSequence=which[i];3253 double value=fabs(infeasibilities[iSequence]/weights[iSequence]);3254 if (value>bestValue) {3255 bestValue=value;3256 bestSequence=which[i];3257 3258 weightedDj[nPossible]=value;3259 sequenceIn[nPossible++]=iSequence;3260 3261 #endif 3262 3263 3234 if (which[i] >= 0) { 3235 int iSequence = which[i]; 3236 double value = fabs(infeasibilities[iSequence] / weights[iSequence]); 3237 if (value > bestValue) { 3238 bestValue = value; 3239 bestSequence = which[i]; 3240 } 3241 weightedDj[nPossible] = value; 3242 sequenceIn[nPossible++] = iSequence; 3243 } 3244 #endif 3245 } 3246 infeasible.setNumElements(numberInfeas); 3264 3247 #ifdef PRICE_IN_ABC_MATRIX 3265 CoinSort_2(weightedDj,weightedDj+nPossible,sequenceIn); 3266 //model_>setMultipleSequenceIn(sequenceIn); 3267 #endif 3268 return bestSequence; 3269 3248 CoinSort_2(weightedDj, weightedDj + nPossible, sequenceIn); 3249 //model_>setMultipleSequenceIn(sequenceIn); 3250 #endif 3251 return bestSequence; 3270 3252 } 3271 3253 // Partial pricing 3272 void 3273 AbcMatrix::partialPricing(double startFraction, double endFraction, 3274 int & bestSequence, int & numberWanted) 3254 void AbcMatrix::partialPricing(double startFraction, double endFraction, 3255 int &bestSequence, int &numberWanted) 3275 3256 { 3276 3257 numberWanted = currentWanted_; 3277 3258 int maximumRows = model_>maximumAbcNumberRows(); 3278 3259 // get matrix data pointers 3279 const int * 3280 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;3281 const double * 3282 int numberColumns =model_>numberColumns();3283 int start = static_cast< int>(startFraction * numberColumns);3284 int end = CoinMin(static_cast< int>(endFraction * numberColumns + 1), numberColumns);3260 const int *COIN_RESTRICT row = matrix_>getIndices(); 3261 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 3262 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 3263 int numberColumns = model_>numberColumns(); 3264 int start = static_cast< int >(startFraction * numberColumns); 3265 int end = CoinMin(static_cast< int >(endFraction * numberColumns + 1), numberColumns); 3285 3266 // adjust 3286 3267 start += maximumRows; 3287 3268 end += maximumRows; 3288 3269 double tolerance = model_>currentDualTolerance(); 3289 double * 3290 const double * 3291 const double * 3270 double *reducedCost = model_>djRegion(); 3271 const double *duals = model_>dualRowSolution(); 3272 const double *cost = model_>costRegion(); 3292 3273 double bestDj; 3293 3274 if (bestSequence >= 0) … … 3303 3284 double value; 3304 3285 AbcSimplex::Status status = model_>getInternalStatus(iSequence); 3305 switch (status) {3286 switch (status) { 3306 3287 case AbcSimplex::basic: 3307 3288 case AbcSimplex::isFixed: 3308 3289 break; 3309 3290 case AbcSimplex::isFree: 3310 3291 case AbcSimplex::superBasic: 3311 3312 3313 j < columnStart[iSequence+1]; j++) {3314 3315 3316 3317 3318 3319 3320 3321 3322 3323 3324 3325 3326 3327 3328 3329 3330 3331 3332 3333 3292 value = cost[iSequence]; 3293 for (CoinBigIndex j = columnStart[iSequence]; 3294 j < columnStart[iSequence + 1]; j++) { 3295 int jRow = row[j]; 3296 value = duals[jRow] * elementByColumn[j]; 3297 } 3298 value = fabs(value); 3299 if (value > FREE_ACCEPT * tolerance) { 3300 numberWanted; 3301 // we are going to bias towards free (but only if reasonable) 3302 value *= FREE_BIAS; 3303 if (value > bestDj) { 3304 // check flagged variable and correct dj 3305 if (!model_>flagged(iSequence)) { 3306 bestDj = value; 3307 bestSequence = iSequence; 3308 } else { 3309 // just to make sure we don't exit before got something 3310 numberWanted++; 3311 } 3312 } 3313 } 3314 break; 3334 3315 case AbcSimplex::atUpperBound: 3335 3336 3337 3338 j < columnStart[iSequence+1]; j++) {3339 3340 3341 3342 3343 3344 3345 3346 3347 3348 3349 3350 3351 3352 3353 3354 3355 3316 value = cost[iSequence]; 3317 // scaled 3318 for (CoinBigIndex j = columnStart[iSequence]; 3319 j < columnStart[iSequence + 1]; j++) { 3320 int jRow = row[j]; 3321 value = duals[jRow] * elementByColumn[j]; 3322 } 3323 if (value > tolerance) { 3324 numberWanted; 3325 if (value > bestDj) { 3326 // check flagged variable and correct dj 3327 if (!model_>flagged(iSequence)) { 3328 bestDj = value; 3329 bestSequence = iSequence; 3330 } else { 3331 // just to make sure we don't exit before got something 3332 numberWanted++; 3333 } 3334 } 3335 } 3336 break; 3356 3337 case AbcSimplex::atLowerBound: 3357 3358 3359 j < columnStart[iSequence+1]; j++) {3360 3361 3362 3363 3364 3365 3366 3367 3368 3369 3370 3371 3372 3373 3374 3375 3376 3377 3338 value = cost[iSequence]; 3339 for (CoinBigIndex j = columnStart[iSequence]; 3340 j < columnStart[iSequence + 1]; j++) { 3341 int jRow = row[j]; 3342 value = duals[jRow] * elementByColumn[j]; 3343 } 3344 value = value; 3345 if (value > tolerance) { 3346 numberWanted; 3347 if (value > bestDj) { 3348 // check flagged variable and correct dj 3349 if (!model_>flagged(iSequence)) { 3350 bestDj = value; 3351 bestSequence = iSequence; 3352 } else { 3353 // just to make sure we don't exit before got something 3354 numberWanted++; 3355 } 3356 } 3357 } 3358 break; 3378 3359 } 3379 3360 } … … 3389 3370 double value = cost[bestSequence]; 3390 3371 for (CoinBigIndex j = columnStart[bestSequence]; 3391 j < columnStart[bestSequence+1]; j++) {3372 j < columnStart[bestSequence + 1]; j++) { 3392 3373 int jRow = row[j]; 3393 3374 value = duals[jRow] * elementByColumn[j]; … … 3402 3383 just for indices Already in z. 3403 3384 Note  z always packed mode */ 3404 void 3405 AbcMatrix::subsetTransposeTimes(const CoinIndexedVector & x, 3406 CoinIndexedVector & z) const 3385 void AbcMatrix::subsetTransposeTimes(const CoinIndexedVector &x, 3386 CoinIndexedVector &z) const 3407 3387 { 3408 3388 int maximumRows = model_>maximumAbcNumberRows(); 3409 3389 // get matrix data pointers 3410 const int * 3411 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;3412 const double * 3413 const double * COIN_RESTRICT other=x.denseVector();3390 const int *COIN_RESTRICT row = matrix_>getIndices(); 3391 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 3392 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 3393 const double *COIN_RESTRICT other = x.denseVector(); 3414 3394 int numberNonZero = z.getNumElements(); 3415 int * 3416 double * 3417 int numberRows =model_>numberRows();3418 for (int i =0;i<numberNonZero;i++) {3419 int iSequence =index[i];3420 if (iSequence >=numberRows) {3395 int *COIN_RESTRICT index = z.getIndices(); 3396 double *COIN_RESTRICT array = z.denseVector(); 3397 int numberRows = model_>numberRows(); 3398 for (int i = 0; i < numberNonZero; i++) { 3399 int iSequence = index[i]; 3400 if (iSequence >= numberRows) { 3421 3401 CoinBigIndex start = columnStart[iSequence]; 3422 CoinBigIndex next = columnStart[iSequence +1];3402 CoinBigIndex next = columnStart[iSequence + 1]; 3423 3403 double value = 0.0; 3424 3404 for (CoinBigIndex j = start; j < next; j++) { 3425 int jRow = row[j]; 3426 value = other[jRow] * elementByColumn[j]; 3427 } 3428 array[i]=value; 3429 } else { 3430 array[i]=other[iSequence]; 3431 } 3432 } 3433 } 3434 // Return <code>x *A</code> in <code>z</code> 3435 void 3436 AbcMatrix::transposeTimes(const CoinIndexedVector & x, 3437 CoinIndexedVector & z) const 3405 int jRow = row[j]; 3406 value = other[jRow] * elementByColumn[j]; 3407 } 3408 array[i] = value; 3409 } else { 3410 array[i] = other[iSequence]; 3411 } 3412 } 3413 } 3414 // Return <code>x *A</code> in <code>z</code> 3415 void AbcMatrix::transposeTimes(const CoinIndexedVector &x, 3416 CoinIndexedVector &z) const 3438 3417 { 3439 3418 int maximumRows = model_>maximumAbcNumberRows(); 3440 3419 // get matrix data pointers 3441 const int * 3442 const CoinBigIndex * COIN_RESTRICT columnStart = matrix_>getVectorStarts()maximumRows;3443 const double * 3444 const double * COIN_RESTRICT other=x.denseVector();3420 const int *COIN_RESTRICT row = matrix_>getIndices(); 3421 const CoinBigIndex *COIN_RESTRICT columnStart = matrix_>getVectorStarts()  maximumRows; 3422 const double *COIN_RESTRICT elementByColumn = matrix_>getElements(); 3423 const double *COIN_RESTRICT other = x.denseVector(); 3445 3424 int numberNonZero = 0; 3446 int * 3447 double * 3448 int numberColumns =model_>numberColumns();3449 double zeroTolerance =model_>zeroTolerance();3450 for (int iSequence =maximumRows;iSequence<maximumRows+numberColumns;iSequence++) {3425 int *COIN_RESTRICT index = z.getIndices(); 3426 double *COIN_RESTRICT array = z.denseVector(); 3427 int numberColumns = model_>numberColumns(); 3428 double zeroTolerance = model_>zeroTolerance(); 3429 for (int iSequence = maximumRows; iSequence < maximumRows + numberColumns; iSequence++) { 3451 3430 CoinBigIndex start = columnStart[iSequence]; 3452 CoinBigIndex next = columnStart[iSequence +1];3431 CoinBigIndex next = columnStart[iSequence + 1]; 3453 3432 double value = 0.0; 3454 3433 for (CoinBigIndex j = start; j < next; j++) { … … 3456 3435 value = other[jRow] * elementByColumn[j]; 3457 3436 } 3458 if (fabs(value) >zeroTolerance) {3437 if (fabs(value) > zeroTolerance) { 3459 3438 // TEMP 3460 array[iSequence maximumRows]=value;3461 index[numberNonZero++] =iSequencemaximumRows;3439