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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpNetworkBasis.cpp
r2271 r2385 11 11 #include "CoinIndexedVector.hpp" 12 12 13 14 13 //############################################################################# 15 14 // Constructors / Destructor / Assignment … … 19 18 // Default Constructor 20 19 // 21 ClpNetworkBasis::ClpNetworkBasis 20 ClpNetworkBasis::ClpNetworkBasis() 22 21 { 23 22 #ifndef COIN_FAST_CODE 24 23 slackValue_ = 1.0; 25 24 #endif 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 25 numberRows_ = 0; 26 numberColumns_ = 0; 27 parent_ = NULL; 28 descendant_ = NULL; 29 pivot_ = NULL; 30 rightSibling_ = NULL; 31 leftSibling_ = NULL; 32 sign_ = NULL; 33 stack_ = NULL; 34 permute_ = NULL; 35 permuteBack_ = NULL; 36 stack2_ = NULL; 37 depth_ = NULL; 38 mark_ = NULL; 39 model_ = NULL; 41 40 } 42 41 // Constructor from CoinFactorization 43 ClpNetworkBasis::ClpNetworkBasis(const ClpSimplex * 44 int numberRows, const CoinFactorizationDouble *pivotRegion,45 const int *permuteBack,46 const int *startColumn,47 const int *numberInColumn,48 const int *indexRow, const CoinFactorizationDouble * /*element*/)42 ClpNetworkBasis::ClpNetworkBasis(const ClpSimplex *model, 43 int numberRows, const CoinFactorizationDouble *pivotRegion, 44 const int *permuteBack, 45 const int *startColumn, 46 const int *numberInColumn, 47 const int *indexRow, const CoinFactorizationDouble * /*element*/) 49 48 { 50 49 #ifndef COIN_FAST_CODE 51 50 slackValue_ = 1.0; 52 51 #endif 53 54 55 parent_ = new int [ numberRows_+1];56 descendant_ = new int [ numberRows_+1];57 pivot_ = new int [ numberRows_+1];58 rightSibling_ = new int [ numberRows_+1];59 leftSibling_ = new int [ numberRows_+1];60 sign_ = new double [ numberRows_+1];61 stack_ = new int [ numberRows_+1];62 stack2_ = new int[numberRows_+1];63 depth_ = new int[numberRows_+1];64 mark_ = new char[numberRows_+1];65 permute_ = new int[numberRows_ + 1];66 permuteBack_ = new int[numberRows_ + 1];67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 52 numberRows_ = numberRows; 53 numberColumns_ = numberRows; 54 parent_ = new int[numberRows_ + 1]; 55 descendant_ = new int[numberRows_ + 1]; 56 pivot_ = new int[numberRows_ + 1]; 57 rightSibling_ = new int[numberRows_ + 1]; 58 leftSibling_ = new int[numberRows_ + 1]; 59 sign_ = new double[numberRows_ + 1]; 60 stack_ = new int[numberRows_ + 1]; 61 stack2_ = new int[numberRows_ + 1]; 62 depth_ = new int[numberRows_ + 1]; 63 mark_ = new char[numberRows_ + 1]; 64 permute_ = new int[numberRows_ + 1]; 65 permuteBack_ = new int[numberRows_ + 1]; 66 int i; 67 for (i = 0; i < numberRows_ + 1; i++) { 68 parent_[i] = 1; 69 descendant_[i] = 1; 70 pivot_[i] = 1; 71 rightSibling_[i] = 1; 72 leftSibling_[i] = 1; 73 sign_[i] = 1.0; 74 stack_[i] = 1; 75 permute_[i] = i; 76 permuteBack_[i] = i; 77 stack2_[i] = 1; 78 depth_[i] = 1; 79 mark_[i] = 0; 80 } 81 mark_[numberRows_] = 1; 82 // pivotColumnBack gives order of pivoting into basis 83 // so pivotColumnback[0] is first slack in basis and 84 // it pivots on row permuteBack[0] 85 // a known root is given by permuteBack[numberRows_1] 86 for (i = 0; i < numberRows_; i++) { 87 int iPivot = permuteBack[i]; 88 double sign; 89 if (pivotRegion[i] > 0.0) 90 sign = 1.0; 91 else 92 sign = 1.0; 93 int other; 94 if (numberInColumn[i] > 0) { 95 int iRow = indexRow[startColumn[i]]; 96 other = permuteBack[iRow]; 97 //assert (parent_[other]!=1); 98 } else { 99 other = numberRows_; 100 } 101 sign_[iPivot] = sign; 102 int iParent = other; 103 parent_[iPivot] = other; 104 if (descendant_[iParent] >= 0) { 105 // we have a sibling 106 int iRight = descendant_[iParent]; 107 rightSibling_[iPivot] = iRight; 108 leftSibling_[iRight] = iPivot; 109 } else { 110 rightSibling_[iPivot] = 1; 111 } 112 descendant_[iParent] = iPivot; 113 leftSibling_[iPivot] = 1; 114 } 115 // do depth 116 int nStack = 1; 117 stack_[0] = descendant_[numberRows_]; 118 depth_[numberRows_] = 1; // root 119 while (nStack) { 120 // take off 121 int iNext = stack_[nStack]; 122 if (iNext >= 0) { 123 depth_[iNext] = nStack; 124 int iRight = rightSibling_[iNext]; 125 stack_[nStack++] = iRight; 126 if (descendant_[iNext] >= 0) 127 stack_[nStack++] = descendant_[iNext]; 128 } 129 } 130 model_ = model; 131 check(); 133 132 } 134 133 … … 136 135 // Copy constructor 137 136 // 138 ClpNetworkBasis::ClpNetworkBasis (const ClpNetworkBasis &rhs)137 ClpNetworkBasis::ClpNetworkBasis(const ClpNetworkBasis &rhs) 139 138 { 140 139 #ifndef COIN_FAST_CODE 141 140 slackValue_ = rhs.slackValue_; 142 141 #endif 143 144 145 146 parent_ = new int [numberRows_+1];147 148 149 150 151 152 descendant_ = new int [numberRows_+1];153 154 155 156 157 158 pivot_ = new int [numberRows_+1];159 160 161 162 163 164 rightSibling_ = new int [numberRows_+1];165 166 167 168 169 170 leftSibling_ = new int [numberRows_+1];171 172 173 174 175 176 sign_ = new double [numberRows_+1];177 178 179 180 181 182 stack_ = new int [numberRows_+1];183 184 185 186 187 188 permute_ = new int [numberRows_+1];189 190 191 192 193 194 permuteBack_ = new int [numberRows_+1];195 196 197 198 199 200 stack2_ = new int [numberRows_+1];201 202 203 204 205 206 depth_ = new int [numberRows_+1];207 208 209 210 211 212 mark_ = new char [numberRows_+1];213 214 215 216 217 142 numberRows_ = rhs.numberRows_; 143 numberColumns_ = rhs.numberColumns_; 144 if (rhs.parent_) { 145 parent_ = new int[numberRows_ + 1]; 146 CoinMemcpyN(rhs.parent_, (numberRows_ + 1), parent_); 147 } else { 148 parent_ = NULL; 149 } 150 if (rhs.descendant_) { 151 descendant_ = new int[numberRows_ + 1]; 152 CoinMemcpyN(rhs.descendant_, (numberRows_ + 1), descendant_); 153 } else { 154 descendant_ = NULL; 155 } 156 if (rhs.pivot_) { 157 pivot_ = new int[numberRows_ + 1]; 158 CoinMemcpyN(rhs.pivot_, (numberRows_ + 1), pivot_); 159 } else { 160 pivot_ = NULL; 161 } 162 if (rhs.rightSibling_) { 163 rightSibling_ = new int[numberRows_ + 1]; 164 CoinMemcpyN(rhs.rightSibling_, (numberRows_ + 1), rightSibling_); 165 } else { 166 rightSibling_ = NULL; 167 } 168 if (rhs.leftSibling_) { 169 leftSibling_ = new int[numberRows_ + 1]; 170 CoinMemcpyN(rhs.leftSibling_, (numberRows_ + 1), leftSibling_); 171 } else { 172 leftSibling_ = NULL; 173 } 174 if (rhs.sign_) { 175 sign_ = new double[numberRows_ + 1]; 176 CoinMemcpyN(rhs.sign_, (numberRows_ + 1), sign_); 177 } else { 178 sign_ = NULL; 179 } 180 if (rhs.stack_) { 181 stack_ = new int[numberRows_ + 1]; 182 CoinMemcpyN(rhs.stack_, (numberRows_ + 1), stack_); 183 } else { 184 stack_ = NULL; 185 } 186 if (rhs.permute_) { 187 permute_ = new int[numberRows_ + 1]; 188 CoinMemcpyN(rhs.permute_, (numberRows_ + 1), permute_); 189 } else { 190 permute_ = NULL; 191 } 192 if (rhs.permuteBack_) { 193 permuteBack_ = new int[numberRows_ + 1]; 194 CoinMemcpyN(rhs.permuteBack_, (numberRows_ + 1), permuteBack_); 195 } else { 196 permuteBack_ = NULL; 197 } 198 if (rhs.stack2_) { 199 stack2_ = new int[numberRows_ + 1]; 200 CoinMemcpyN(rhs.stack2_, (numberRows_ + 1), stack2_); 201 } else { 202 stack2_ = NULL; 203 } 204 if (rhs.depth_) { 205 depth_ = new int[numberRows_ + 1]; 206 CoinMemcpyN(rhs.depth_, (numberRows_ + 1), depth_); 207 } else { 208 depth_ = NULL; 209 } 210 if (rhs.mark_) { 211 mark_ = new char[numberRows_ + 1]; 212 CoinMemcpyN(rhs.mark_, (numberRows_ + 1), mark_); 213 } else { 214 mark_ = NULL; 215 } 216 model_ = rhs.model_; 218 217 } 219 218 … … 221 220 // Destructor 222 221 // 223 ClpNetworkBasis::~ClpNetworkBasis 222 ClpNetworkBasis::~ClpNetworkBasis() 224 223 { 225 delete[] parent_;226 delete[] descendant_;227 delete[] pivot_;228 delete[] rightSibling_;229 delete[] leftSibling_;230 delete[] sign_;231 delete[] stack_;232 delete[] permute_;233 delete[] permuteBack_;234 delete[] stack2_;235 delete[] depth_;236 delete[] mark_;224 delete[] parent_; 225 delete[] descendant_; 226 delete[] pivot_; 227 delete[] rightSibling_; 228 delete[] leftSibling_; 229 delete[] sign_; 230 delete[] stack_; 231 delete[] permute_; 232 delete[] permuteBack_; 233 delete[] stack2_; 234 delete[] depth_; 235 delete[] mark_; 237 236 } 238 237 … … 241 240 // 242 241 ClpNetworkBasis & 243 ClpNetworkBasis::operator=(const ClpNetworkBasis &rhs)242 ClpNetworkBasis::operator=(const ClpNetworkBasis &rhs) 244 243 { 245 246 delete[] parent_;247 delete[] descendant_;248 delete[] pivot_;249 delete[] rightSibling_;250 delete[] leftSibling_;251 delete[] sign_;252 delete[] stack_;253 delete[] permute_;254 delete[] permuteBack_;255 delete[] stack2_;256 delete[] depth_;257 delete[] mark_;244 if (this != &rhs) { 245 delete[] parent_; 246 delete[] descendant_; 247 delete[] pivot_; 248 delete[] rightSibling_; 249 delete[] leftSibling_; 250 delete[] sign_; 251 delete[] stack_; 252 delete[] permute_; 253 delete[] permuteBack_; 254 delete[] stack2_; 255 delete[] depth_; 256 delete[] mark_; 258 257 #ifndef COIN_FAST_CODE 259 258 slackValue_ = rhs.slackValue_; 260 259 #endif 261 262 263 264 parent_ = new int [numberRows_+1];265 266 267 268 269 270 descendant_ = new int [numberRows_+1];271 272 273 274 275 276 pivot_ = new int [numberRows_+1];277 278 279 280 281 282 rightSibling_ = new int [numberRows_+1];283 284 285 286 287 288 leftSibling_ = new int [numberRows_+1];289 290 291 292 293 294 sign_ = new double [numberRows_+1];295 296 297 298 299 300 stack_ = new int [numberRows_+1];301 302 303 304 305 306 permute_ = new int [numberRows_+1];307 308 309 310 311 312 permuteBack_ = new int [numberRows_+1];313 314 315 316 317 318 stack2_ = new int [numberRows_+1];319 320 321 322 323 324 depth_ = new int [numberRows_+1];325 326 327 328 329 330 mark_ = new char [numberRows_+1];331 332 333 334 335 336 260 numberRows_ = rhs.numberRows_; 261 numberColumns_ = rhs.numberColumns_; 262 if (rhs.parent_) { 263 parent_ = new int[numberRows_ + 1]; 264 CoinMemcpyN(rhs.parent_, (numberRows_ + 1), parent_); 265 } else { 266 parent_ = NULL; 267 } 268 if (rhs.descendant_) { 269 descendant_ = new int[numberRows_ + 1]; 270 CoinMemcpyN(rhs.descendant_, (numberRows_ + 1), descendant_); 271 } else { 272 descendant_ = NULL; 273 } 274 if (rhs.pivot_) { 275 pivot_ = new int[numberRows_ + 1]; 276 CoinMemcpyN(rhs.pivot_, (numberRows_ + 1), pivot_); 277 } else { 278 pivot_ = NULL; 279 } 280 if (rhs.rightSibling_) { 281 rightSibling_ = new int[numberRows_ + 1]; 282 CoinMemcpyN(rhs.rightSibling_, (numberRows_ + 1), rightSibling_); 283 } else { 284 rightSibling_ = NULL; 285 } 286 if (rhs.leftSibling_) { 287 leftSibling_ = new int[numberRows_ + 1]; 288 CoinMemcpyN(rhs.leftSibling_, (numberRows_ + 1), leftSibling_); 289 } else { 290 leftSibling_ = NULL; 291 } 292 if (rhs.sign_) { 293 sign_ = new double[numberRows_ + 1]; 294 CoinMemcpyN(rhs.sign_, (numberRows_ + 1), sign_); 295 } else { 296 sign_ = NULL; 297 } 298 if (rhs.stack_) { 299 stack_ = new int[numberRows_ + 1]; 300 CoinMemcpyN(rhs.stack_, (numberRows_ + 1), stack_); 301 } else { 302 stack_ = NULL; 303 } 304 if (rhs.permute_) { 305 permute_ = new int[numberRows_ + 1]; 306 CoinMemcpyN(rhs.permute_, (numberRows_ + 1), permute_); 307 } else { 308 permute_ = NULL; 309 } 310 if (rhs.permuteBack_) { 311 permuteBack_ = new int[numberRows_ + 1]; 312 CoinMemcpyN(rhs.permuteBack_, (numberRows_ + 1), permuteBack_); 313 } else { 314 permuteBack_ = NULL; 315 } 316 if (rhs.stack2_) { 317 stack2_ = new int[numberRows_ + 1]; 318 CoinMemcpyN(rhs.stack2_, (numberRows_ + 1), stack2_); 319 } else { 320 stack2_ = NULL; 321 } 322 if (rhs.depth_) { 323 depth_ = new int[numberRows_ + 1]; 324 CoinMemcpyN(rhs.depth_, (numberRows_ + 1), depth_); 325 } else { 326 depth_ = NULL; 327 } 328 if (rhs.mark_) { 329 mark_ = new char[numberRows_ + 1]; 330 CoinMemcpyN(rhs.mark_, (numberRows_ + 1), mark_); 331 } else { 332 mark_ = NULL; 333 } 334 } 335 return *this; 337 336 } 338 337 // checks looks okay 339 338 void ClpNetworkBasis::check() 340 339 { 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 340 // check depth 341 int nStack = 1; 342 stack_[0] = descendant_[numberRows_]; 343 depth_[numberRows_] = 1; // root 344 while (nStack) { 345 // take off 346 int iNext = stack_[nStack]; 347 if (iNext >= 0) { 348 //assert (depth_[iNext]==nStack); 349 depth_[iNext] = nStack; 350 int iRight = rightSibling_[iNext]; 351 stack_[nStack++] = iRight; 352 if (descendant_[iNext] >= 0) 353 stack_[nStack++] = descendant_[iNext]; 354 } 355 } 357 356 } 358 357 // prints 359 358 void ClpNetworkBasis::print() 360 359 { 361 362 363 364 365 366 360 int i; 361 printf(" parent descendant left right sign depth\n"); 362 for (i = 0; i < numberRows_ + 1; i++) 363 printf("%4d %7d %8d %7d %7d %5g %7d\n", 364 i, parent_[i], descendant_[i], leftSibling_[i], rightSibling_[i], 365 sign_[i], depth_[i]); 367 366 } 368 367 /* Replaces one Column to basis, 369 368 returns 0=OK 370 369 */ 371 int 372 ClpNetworkBasis::replaceColumn ( CoinIndexedVector * regionSparse, 373 int pivotRow) 370 int ClpNetworkBasis::replaceColumn(CoinIndexedVector *regionSparse, 371 int pivotRow) 374 372 { 375 376 377 378 assert(!regionSparse>getNumElements());379 380 381 382 383 int *indices = regionSparse>getIndices();384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 373 // When things have settled down then redo this to make more elegant 374 // I am sure lots of loops can be combined 375 // regionSparse is empty 376 assert(!regionSparse>getNumElements()); 377 model_>unpack(regionSparse, model_>sequenceIn()); 378 // arc given by pivotRow is leaving basis 379 //int kParent = parent_[pivotRow]; 380 // arc coming in has these two nodes 381 int *indices = regionSparse>getIndices(); 382 int iRow0 = indices[0]; 383 int iRow1; 384 if (regionSparse>getNumElements() == 2) 385 iRow1 = indices[1]; 386 else 387 iRow1 = numberRows_; 388 double sign = regionSparse>denseVector()[iRow0]; 389 regionSparse>clear(); 390 // and outgoing 391 model_>unpack(regionSparse, model_>pivotVariable()[pivotRow]); 392 int jRow0 = indices[0]; 393 int jRow1; 394 if (regionSparse>getNumElements() == 2) 395 jRow1 = indices[1]; 396 else 397 jRow1 = numberRows_; 398 regionSparse>clear(); 399 // get correct pivotRow 400 //#define FULL_DEBUG 403 401 #ifdef FULL_DEBUG 404 printf("irow %d %d, jrow %d %d\n",405 402 printf("irow %d %d, jrow %d %d\n", 403 iRow0, iRow1, jRow0, jRow1); 406 404 #endif 407 408 409 405 if (parent_[jRow0] == jRow1) { 406 int newPivot = jRow0; 407 if (newPivot != pivotRow) { 410 408 #ifdef FULL_DEBUG 411 409 printf("pivot row of %d permuted to %d\n", pivotRow, newPivot); 412 410 #endif 413 414 415 416 417 418 411 pivotRow = newPivot; 412 } 413 } else { 414 //assert (parent_[jRow1]==jRow0); 415 int newPivot = jRow1; 416 if (newPivot != pivotRow) { 419 417 #ifdef FULL_DEBUG 420 418 printf("pivot row of %d permuted to %d\n", pivotRow, newPivot); 421 419 #endif 422 pivotRow = newPivot; 423 } 424 } 425 bool extraPrint = (model_>numberIterations() > 3) && 426 (model_>logLevel() > 10); 427 if (extraPrint) 428 print(); 420 pivotRow = newPivot; 421 } 422 } 423 bool extraPrint = (model_>numberIterations() > 3) && (model_>logLevel() > 10); 424 if (extraPrint) 425 print(); 429 426 #ifdef FULL_DEBUG 430 431 iRow1, sign, sign_[iRow1], iRow0, sign_[iRow0], pivotRow, sign_[pivotRow]);427 printf("In %d (region= %g, stored %g) %d (%g) pivoting on %d (%g)\n", 428 iRow1, sign, sign_[iRow1], iRow0, sign_[iRow0], pivotRow, sign_[pivotRow]); 432 429 #endif 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 int newParent = stack_[nStack1];430 // see which path outgoing pivot is on 431 int kRow = 1; 432 int jRow = iRow1; 433 while (jRow != numberRows_) { 434 if (jRow == pivotRow) { 435 kRow = iRow1; 436 break; 437 } else { 438 jRow = parent_[jRow]; 439 } 440 } 441 if (kRow < 0) { 442 jRow = iRow0; 443 while (jRow != numberRows_) { 444 if (jRow == pivotRow) { 445 kRow = iRow0; 446 break; 447 } else { 448 jRow = parent_[jRow]; 449 } 450 } 451 } 452 //assert (kRow>=0); 453 if (iRow0 == kRow) { 454 iRow0 = iRow1; 455 iRow1 = kRow; 456 sign = sign; 457 } 458 // pivot row is on path from iRow1 back to root 459 // get stack of nodes to change 460 // Also get precursors for cleaning order 461 int nStack = 1; 462 stack_[0] = iRow0; 463 while (kRow != pivotRow) { 464 stack_[nStack++] = kRow; 465 if (sign * sign_[kRow] < 0.0) { 466 sign_[kRow] = sign_[kRow]; 467 } else { 468 sign = sign; 469 } 470 kRow = parent_[kRow]; 471 //sign *= sign_[kRow]; 472 } 473 stack_[nStack++] = pivotRow; 474 if (sign * sign_[pivotRow] < 0.0) { 475 sign_[pivotRow] = sign_[pivotRow]; 476 } else { 477 sign = sign; 478 } 479 int iParent = parent_[pivotRow]; 480 while (nStack > 1) { 481 int iLeft; 482 int iRight; 483 kRow = stack_[nStack]; 484 int newParent = stack_[nStack  1]; 488 485 #ifdef FULL_DEBUG 489 490 486 printf("row %d, old parent %d, new parent %d, pivotrow %d\n", kRow, 487 iParent, newParent, pivotRow); 491 488 #endif 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 489 int i1 = permuteBack_[pivotRow]; 490 int i2 = permuteBack_[kRow]; 491 permuteBack_[pivotRow] = i2; 492 permuteBack_[kRow] = i1; 493 // do Btran permutation 494 permute_[i1] = kRow; 495 permute_[i2] = pivotRow; 496 pivotRow = kRow; 497 // Take out of old parent 498 iLeft = leftSibling_[kRow]; 499 iRight = rightSibling_[kRow]; 500 if (iLeft < 0) { 501 // take out of tree 502 if (iRight >= 0) { 503 leftSibling_[iRight] = iLeft; 504 descendant_[iParent] = iRight; 505 } else { 509 506 #ifdef FULL_DEBUG 510 511 507 printf("Saying %d (old parent of %d) has no descendants\n", 508 iParent, kRow); 512 509 #endif 513 514 515 516 517 518 519 520 521 522 510 descendant_[iParent] = 1; 511 } 512 } else { 513 // take out of tree 514 rightSibling_[iLeft] = iRight; 515 if (iRight >= 0) 516 leftSibling_[iRight] = iLeft; 517 } 518 leftSibling_[kRow] = 1; 519 rightSibling_[kRow] = 1; 523 520 524 525 526 527 528 529 530 531 532 533 534 535 536 521 // now insert new one 522 // make this descendant of that 523 if (descendant_[newParent] >= 0) { 524 // we will have a sibling 525 int jRight = descendant_[newParent]; 526 rightSibling_[kRow] = jRight; 527 leftSibling_[jRight] = kRow; 528 } else { 529 rightSibling_[kRow] = 1; 530 } 531 descendant_[newParent] = kRow; 532 leftSibling_[kRow] = 1; 533 parent_[kRow] = newParent; 537 534 538 539 540 541 542 543 int iPivot= stack_[1];544 545 iDepth++; //correct for this one546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 535 iParent = kRow; 536 } 537 // now redo all depths from stack_[1] 538 // This must be possible to combine  but later 539 { 540 int iPivot = stack_[1]; 541 int iDepth = depth_[parent_[iPivot]]; //depth of parent 542 iDepth++; //correct for this one 543 int nStack = 1; 544 stack_[0] = iPivot; 545 while (nStack) { 546 // take off 547 int iNext = stack_[nStack]; 548 if (iNext >= 0) { 549 // add stack level 550 depth_[iNext] = nStack + iDepth; 551 stack_[nStack++] = rightSibling_[iNext]; 552 if (descendant_[iNext] >= 0) 553 stack_[nStack++] = descendant_[iNext]; 554 } 555 } 556 } 557 if (extraPrint) 558 print(); 559 //check(); 560 return 0; 564 561 } 565 562 566 563 /* Updates one column (FTRAN) from region2 */ 567 564 double 568 ClpNetworkBasis::updateColumn ( CoinIndexedVector *regionSparse,569 CoinIndexedVector *regionSparse2,570 565 ClpNetworkBasis::updateColumn(CoinIndexedVector *regionSparse, 566 CoinIndexedVector *regionSparse2, 567 int pivotRow) 571 568 { 572 regionSparse>clear ( ); 573 double *region = regionSparse>denseVector ( ); 574 double *region2 = regionSparse2>denseVector ( ); 575 int *regionIndex2 = regionSparse2>getIndices ( ); 576 int numberNonZero = regionSparse2>getNumElements ( ); 577 int *regionIndex = regionSparse>getIndices ( ); 578 int i; 579 bool doTwo = (numberNonZero == 2); 580 int i0 = 1; 581 int i1 = 1; 582 if (doTwo) { 583 i0 = regionIndex2[0]; 584 i1 = regionIndex2[1]; 585 } 586 double returnValue = 0.0; 587 bool packed = regionSparse2>packedMode(); 588 if (packed) { 589 if (doTwo && region2[0]*region2[1] < 0.0) { 590 region[i0] = region2[0]; 591 region2[0] = 0.0; 592 region[i1] = region2[1]; 593 region2[1] = 0.0; 594 int iDepth0 = depth_[i0]; 595 int iDepth1 = depth_[i1]; 596 if (iDepth1 > iDepth0) { 597 int temp = i0; 598 i0 = i1; 599 i1 = temp; 600 temp = iDepth0; 601 iDepth0 = iDepth1; 602 iDepth1 = temp; 603 } 604 numberNonZero = 0; 605 if (pivotRow < 0) { 606 while (iDepth0 > iDepth1) { 607 double pivotValue = region[i0]; 608 // put back now ? 609 int iBack = permuteBack_[i0]; 610 region2[numberNonZero] = pivotValue * sign_[i0]; 611 regionIndex2[numberNonZero++] = iBack; 612 int otherRow = parent_[i0]; 613 region[i0] = 0.0; 614 region[otherRow] += pivotValue; 615 iDepth0; 616 i0 = otherRow; 617 } 618 while (i0 != i1) { 619 double pivotValue = region[i0]; 620 // put back now ? 621 int iBack = permuteBack_[i0]; 622 region2[numberNonZero] = pivotValue * sign_[i0]; 623 regionIndex2[numberNonZero++] = iBack; 624 int otherRow = parent_[i0]; 625 region[i0] = 0.0; 626 region[otherRow] += pivotValue; 627 i0 = otherRow; 628 double pivotValue1 = region[i1]; 629 // put back now ? 630 int iBack1 = permuteBack_[i1]; 631 region2[numberNonZero] = pivotValue1 * sign_[i1]; 632 regionIndex2[numberNonZero++] = iBack1; 633 int otherRow1 = parent_[i1]; 634 region[i1] = 0.0; 635 region[otherRow1] += pivotValue1; 636 i1 = otherRow1; 637 } 638 } else { 639 while (iDepth0 > iDepth1) { 640 double pivotValue = region[i0]; 641 // put back now ? 642 int iBack = permuteBack_[i0]; 643 double value = pivotValue * sign_[i0]; 644 region2[numberNonZero] = value; 645 regionIndex2[numberNonZero++] = iBack; 646 if (iBack == pivotRow) 647 returnValue = value; 648 int otherRow = parent_[i0]; 649 region[i0] = 0.0; 650 region[otherRow] += pivotValue; 651 iDepth0; 652 i0 = otherRow; 653 } 654 while (i0 != i1) { 655 double pivotValue = region[i0]; 656 // put back now ? 657 int iBack = permuteBack_[i0]; 658 double value = pivotValue * sign_[i0]; 659 region2[numberNonZero] = value; 660 regionIndex2[numberNonZero++] = iBack; 661 if (iBack == pivotRow) 662 returnValue = value; 663 int otherRow = parent_[i0]; 664 region[i0] = 0.0; 665 region[otherRow] += pivotValue; 666 i0 = otherRow; 667 double pivotValue1 = region[i1]; 668 // put back now ? 669 int iBack1 = permuteBack_[i1]; 670 value = pivotValue1 * sign_[i1]; 671 region2[numberNonZero] = value; 672 regionIndex2[numberNonZero++] = iBack1; 673 if (iBack1 == pivotRow) 674 returnValue = value; 675 int otherRow1 = parent_[i1]; 676 region[i1] = 0.0; 677 region[otherRow1] += pivotValue1; 678 i1 = otherRow1; 679 } 680 } 681 } else { 682 // set up linked lists at each depth 683 // stack2 is start, stack is next 684 int greatestDepth = 1; 685 //mark_[numberRows_]=1; 686 for (i = 0; i < numberNonZero; i++) { 687 int j = regionIndex2[i]; 688 double value = region2[i]; 689 region2[i] = 0.0; 690 region[j] = value; 691 regionIndex[i] = j; 692 int iDepth = depth_[j]; 693 if (iDepth > greatestDepth) 694 greatestDepth = iDepth; 695 // and back until marked 696 while (!mark_[j]) { 697 int iNext = stack2_[iDepth]; 698 stack2_[iDepth] = j; 699 stack_[j] = iNext; 700 mark_[j] = 1; 701 iDepth; 702 j = parent_[j]; 703 } 704 } 705 numberNonZero = 0; 706 if (pivotRow < 0) { 707 for (; greatestDepth >= 0; greatestDepth) { 708 int iPivot = stack2_[greatestDepth]; 709 stack2_[greatestDepth] = 1; 710 while (iPivot >= 0) { 711 mark_[iPivot] = 0; 712 double pivotValue = region[iPivot]; 713 if (pivotValue) { 714 // put back now ? 715 int iBack = permuteBack_[iPivot]; 716 region2[numberNonZero] = pivotValue * sign_[iPivot]; 717 regionIndex2[numberNonZero++] = iBack; 718 int otherRow = parent_[iPivot]; 719 region[iPivot] = 0.0; 720 region[otherRow] += pivotValue; 721 } 722 iPivot = stack_[iPivot]; 723 } 724 } 725 } else { 726 for (; greatestDepth >= 0; greatestDepth) { 727 int iPivot = stack2_[greatestDepth]; 728 stack2_[greatestDepth] = 1; 729 while (iPivot >= 0) { 730 mark_[iPivot] = 0; 731 double pivotValue = region[iPivot]; 732 if (pivotValue) { 733 // put back now ? 734 int iBack = permuteBack_[iPivot]; 735 double value = pivotValue * sign_[iPivot]; 736 region2[numberNonZero] = value; 737 regionIndex2[numberNonZero++] = iBack; 738 if (iBack == pivotRow) 739 returnValue = value; 740 int otherRow = parent_[iPivot]; 741 region[iPivot] = 0.0; 742 region[otherRow] += pivotValue; 743 } 744 iPivot = stack_[iPivot]; 745 } 746 } 747 } 569 regionSparse>clear(); 570 double *region = regionSparse>denseVector(); 571 double *region2 = regionSparse2>denseVector(); 572 int *regionIndex2 = regionSparse2>getIndices(); 573 int numberNonZero = regionSparse2>getNumElements(); 574 int *regionIndex = regionSparse>getIndices(); 575 int i; 576 bool doTwo = (numberNonZero == 2); 577 int i0 = 1; 578 int i1 = 1; 579 if (doTwo) { 580 i0 = regionIndex2[0]; 581 i1 = regionIndex2[1]; 582 } 583 double returnValue = 0.0; 584 bool packed = regionSparse2>packedMode(); 585 if (packed) { 586 if (doTwo && region2[0] * region2[1] < 0.0) { 587 region[i0] = region2[0]; 588 region2[0] = 0.0; 589 region[i1] = region2[1]; 590 region2[1] = 0.0; 591 int iDepth0 = depth_[i0]; 592 int iDepth1 = depth_[i1]; 593 if (iDepth1 > iDepth0) { 594 int temp = i0; 595 i0 = i1; 596 i1 = temp; 597 temp = iDepth0; 598 iDepth0 = iDepth1; 599 iDepth1 = temp; 600 } 601 numberNonZero = 0; 602 if (pivotRow < 0) { 603 while (iDepth0 > iDepth1) { 604 double pivotValue = region[i0]; 605 // put back now ? 606 int iBack = permuteBack_[i0]; 607 region2[numberNonZero] = pivotValue * sign_[i0]; 608 regionIndex2[numberNonZero++] = iBack; 609 int otherRow = parent_[i0]; 610 region[i0] = 0.0; 611 region[otherRow] += pivotValue; 612 iDepth0; 613 i0 = otherRow; 614 } 615 while (i0 != i1) { 616 double pivotValue = region[i0]; 617 // put back now ? 618 int iBack = permuteBack_[i0]; 619 region2[numberNonZero] = pivotValue * sign_[i0]; 620 regionIndex2[numberNonZero++] = iBack; 621 int otherRow = parent_[i0]; 622 region[i0] = 0.0; 623 region[otherRow] += pivotValue; 624 i0 = otherRow; 625 double pivotValue1 = region[i1]; 626 // put back now ? 627 int iBack1 = permuteBack_[i1]; 628 region2[numberNonZero] = pivotValue1 * sign_[i1]; 629 regionIndex2[numberNonZero++] = iBack1; 630 int otherRow1 = parent_[i1]; 631 region[i1] = 0.0; 632 region[otherRow1] += pivotValue1; 633 i1 = otherRow1; 634 } 635 } else { 636 while (iDepth0 > iDepth1) { 637 double pivotValue = region[i0]; 638 // put back now ? 639 int iBack = permuteBack_[i0]; 640 double value = pivotValue * sign_[i0]; 641 region2[numberNonZero] = value; 642 regionIndex2[numberNonZero++] = iBack; 643 if (iBack == pivotRow) 644 returnValue = value; 645 int otherRow = parent_[i0]; 646 region[i0] = 0.0; 647 region[otherRow] += pivotValue; 648 iDepth0; 649 i0 = otherRow; 650 } 651 while (i0 != i1) { 652 double pivotValue = region[i0]; 653 // put back now ? 654 int iBack = permuteBack_[i0]; 655 double value = pivotValue * sign_[i0]; 656 region2[numberNonZero] = value; 657 regionIndex2[numberNonZero++] = iBack; 658 if (iBack == pivotRow) 659 returnValue = value; 660 int otherRow = parent_[i0]; 661 region[i0] = 0.0; 662 region[otherRow] += pivotValue; 663 i0 = otherRow; 664 double pivotValue1 = region[i1]; 665 // put back now ? 666 int iBack1 = permuteBack_[i1]; 667 value = pivotValue1 * sign_[i1]; 668 region2[numberNonZero] = value; 669 regionIndex2[numberNonZero++] = iBack1; 670 if (iBack1 == pivotRow) 671 returnValue = value; 672 int otherRow1 = parent_[i1]; 673 region[i1] = 0.0; 674 region[otherRow1] += pivotValue1; 675 i1 = otherRow1; 676 } 677 } 678 } else { 679 // set up linked lists at each depth 680 // stack2 is start, stack is next 681 int greatestDepth = 1; 682 //mark_[numberRows_]=1; 683 for (i = 0; i < numberNonZero; i++) { 684 int j = regionIndex2[i]; 685 double value = region2[i]; 686 region2[i] = 0.0; 687 region[j] = value; 688 regionIndex[i] = j; 689 int iDepth = depth_[j]; 690 if (iDepth > greatestDepth) 691 greatestDepth = iDepth; 692 // and back until marked 693 while (!mark_[j]) { 694 int iNext = stack2_[iDepth]; 695 stack2_[iDepth] = j; 696 stack_[j] = iNext; 697 mark_[j] = 1; 698 iDepth; 699 j = parent_[j]; 700 } 701 } 702 numberNonZero = 0; 703 if (pivotRow < 0) { 704 for (; greatestDepth >= 0; greatestDepth) { 705 int iPivot = stack2_[greatestDepth]; 706 stack2_[greatestDepth] = 1; 707 while (iPivot >= 0) { 708 mark_[iPivot] = 0; 709 double pivotValue = region[iPivot]; 710 if (pivotValue) { 711 // put back now ? 712 int iBack = permuteBack_[iPivot]; 713 region2[numberNonZero] = pivotValue * sign_[iPivot]; 714 regionIndex2[numberNonZero++] = iBack; 715 int otherRow = parent_[iPivot]; 716 region[iPivot] = 0.0; 717 region[otherRow] += pivotValue; 718 } 719 iPivot = stack_[iPivot]; 748 720 } 749 } else { 750 if (doTwo && region2[i0]*region2[i1] < 0.0) { 751 // If just + 1 then could go backwards on depth until join 752 region[i0] = region2[i0]; 753 region2[i0] = 0.0; 754 region[i1] = region2[i1]; 755 region2[i1] = 0.0; 756 int iDepth0 = depth_[i0]; 757 int iDepth1 = depth_[i1]; 758 if (iDepth1 > iDepth0) { 759 int temp = i0; 760 i0 = i1; 761 i1 = temp; 762 temp = iDepth0; 763 iDepth0 = iDepth1; 764 iDepth1 = temp; 765 } 766 numberNonZero = 0; 767 while (iDepth0 > iDepth1) { 768 double pivotValue = region[i0]; 769 // put back now ? 770 int iBack = permuteBack_[i0]; 771 regionIndex2[numberNonZero++] = iBack; 772 int otherRow = parent_[i0]; 773 region2[iBack] = pivotValue * sign_[i0]; 774 region[i0] = 0.0; 775 region[otherRow] += pivotValue; 776 iDepth0; 777 i0 = otherRow; 778 } 779 while (i0 != i1) { 780 double pivotValue = region[i0]; 781 // put back now ? 782 int iBack = permuteBack_[i0]; 783 regionIndex2[numberNonZero++] = iBack; 784 int otherRow = parent_[i0]; 785 region2[iBack] = pivotValue * sign_[i0]; 786 region[i0] = 0.0; 787 region[otherRow] += pivotValue; 788 i0 = otherRow; 789 double pivotValue1 = region[i1]; 790 // put back now ? 791 int iBack1 = permuteBack_[i1]; 792 regionIndex2[numberNonZero++] = iBack1; 793 int otherRow1 = parent_[i1]; 794 region2[iBack1] = pivotValue1 * sign_[i1]; 795 region[i1] = 0.0; 796 region[otherRow1] += pivotValue1; 797 i1 = otherRow1; 798 } 799 } else { 800 // set up linked lists at each depth 801 // stack2 is start, stack is next 802 int greatestDepth = 1; 803 //mark_[numberRows_]=1; 804 for (i = 0; i < numberNonZero; i++) { 805 int j = regionIndex2[i]; 806 double value = region2[j]; 807 region2[j] = 0.0; 808 region[j] = value; 809 regionIndex[i] = j; 810 int iDepth = depth_[j]; 811 if (iDepth > greatestDepth) 812 greatestDepth = iDepth; 813 // and back until marked 814 while (!mark_[j]) { 815 int iNext = stack2_[iDepth]; 816 stack2_[iDepth] = j; 817 stack_[j] = iNext; 818 mark_[j] = 1; 819 iDepth; 820 j = parent_[j]; 821 } 822 } 823 numberNonZero = 0; 824 for (; greatestDepth >= 0; greatestDepth) { 825 int iPivot = stack2_[greatestDepth]; 826 stack2_[greatestDepth] = 1; 827 while (iPivot >= 0) { 828 mark_[iPivot] = 0; 829 double pivotValue = region[iPivot]; 830 if (pivotValue) { 831 // put back now ? 832 int iBack = permuteBack_[iPivot]; 833 regionIndex2[numberNonZero++] = iBack; 834 int otherRow = parent_[iPivot]; 835 region2[iBack] = pivotValue * sign_[iPivot]; 836 region[iPivot] = 0.0; 837 region[otherRow] += pivotValue; 838 } 839 iPivot = stack_[iPivot]; 840 } 841 } 721 } 722 } else { 723 for (; greatestDepth >= 0; greatestDepth) { 724 int iPivot = stack2_[greatestDepth]; 725 stack2_[greatestDepth] = 1; 726 while (iPivot >= 0) { 727 mark_[iPivot] = 0; 728 double pivotValue = region[iPivot]; 729 if (pivotValue) { 730 // put back now ? 731 int iBack = permuteBack_[iPivot]; 732 double value = pivotValue * sign_[iPivot]; 733 region2[numberNonZero] = value; 734 regionIndex2[numberNonZero++] = iBack; 735 if (iBack == pivotRow) 736 returnValue = value; 737 int otherRow = parent_[iPivot]; 738 region[iPivot] = 0.0; 739 region[otherRow] += pivotValue; 740 } 741 iPivot = stack_[iPivot]; 842 742 } 843 if (pivotRow >= 0) 844 returnValue = region2[pivotRow]; 845 } 846 region[numberRows_] = 0.0; 847 regionSparse2>setNumElements(numberNonZero); 743 } 744 } 745 } 746 } else { 747 if (doTwo && region2[i0] * region2[i1] < 0.0) { 748 // If just + 1 then could go backwards on depth until join 749 region[i0] = region2[i0]; 750 region2[i0] = 0.0; 751 region[i1] = region2[i1]; 752 region2[i1] = 0.0; 753 int iDepth0 = depth_[i0]; 754 int iDepth1 = depth_[i1]; 755 if (iDepth1 > iDepth0) { 756 int temp = i0; 757 i0 = i1; 758 i1 = temp; 759 temp = iDepth0; 760 iDepth0 = iDepth1; 761 iDepth1 = temp; 762 } 763 numberNonZero = 0; 764 while (iDepth0 > iDepth1) { 765 double pivotValue = region[i0]; 766 // put back now ? 767 int iBack = permuteBack_[i0]; 768 regionIndex2[numberNonZero++] = iBack; 769 int otherRow = parent_[i0]; 770 region2[iBack] = pivotValue * sign_[i0]; 771 region[i0] = 0.0; 772 region[otherRow] += pivotValue; 773 iDepth0; 774 i0 = otherRow; 775 } 776 while (i0 != i1) { 777 double pivotValue = region[i0]; 778 // put back now ? 779 int iBack = permuteBack_[i0]; 780 regionIndex2[numberNonZero++] = iBack; 781 int otherRow = parent_[i0]; 782 region2[iBack] = pivotValue * sign_[i0]; 783 region[i0] = 0.0; 784 region[otherRow] += pivotValue; 785 i0 = otherRow; 786 double pivotValue1 = region[i1]; 787 // put back now ? 788 int iBack1 = permuteBack_[i1]; 789 regionIndex2[numberNonZero++] = iBack1; 790 int otherRow1 = parent_[i1]; 791 region2[iBack1] = pivotValue1 * sign_[i1]; 792 region[i1] = 0.0; 793 region[otherRow1] += pivotValue1; 794 i1 = otherRow1; 795 } 796 } else { 797 // set up linked lists at each depth 798 // stack2 is start, stack is next 799 int greatestDepth = 1; 800 //mark_[numberRows_]=1; 801 for (i = 0; i < numberNonZero; i++) { 802 int j = regionIndex2[i]; 803 double value = region2[j]; 804 region2[j] = 0.0; 805 region[j] = value; 806 regionIndex[i] = j; 807 int iDepth = depth_[j]; 808 if (iDepth > greatestDepth) 809 greatestDepth = iDepth; 810 // and back until marked 811 while (!mark_[j]) { 812 int iNext = stack2_[iDepth]; 813 stack2_[iDepth] = j; 814 stack_[j] = iNext; 815 mark_[j] = 1; 816 iDepth; 817 j = parent_[j]; 818 } 819 } 820 numberNonZero = 0; 821 for (; greatestDepth >= 0; greatestDepth) { 822 int iPivot = stack2_[greatestDepth]; 823 stack2_[greatestDepth] = 1; 824 while (iPivot >= 0) { 825 mark_[iPivot] = 0; 826 double pivotValue = region[iPivot]; 827 if (pivotValue) { 828 // put back now ? 829 int iBack = permuteBack_[iPivot]; 830 regionIndex2[numberNonZero++] = iBack; 831 int otherRow = parent_[iPivot]; 832 region2[iBack] = pivotValue * sign_[iPivot]; 833 region[iPivot] = 0.0; 834 region[otherRow] += pivotValue; 835 } 836 iPivot = stack_[iPivot]; 837 } 838 } 839 } 840 if (pivotRow >= 0) 841 returnValue = region2[pivotRow]; 842 } 843 region[numberRows_] = 0.0; 844 regionSparse2>setNumElements(numberNonZero); 848 845 #ifdef FULL_DEBUG 849 850 851 852 853 assert(stack2_[i] == 1);854 855 846 { 847 int i; 848 for (i = 0; i < numberRows_; i++) { 849 assert(!mark_[i]); 850 assert(stack2_[i] == 1); 851 } 852 } 856 853 #endif 857 854 return returnValue; 858 855 } 859 856 /* Updates one column (FTRAN) to/from array … … 862 859 have got code working using this simple method  thank you! 863 860 (the only exception is if you know input is dense e.g. rhs) */ 864 int 865 ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse, 866 double region2[] ) const 861 int ClpNetworkBasis::updateColumn(CoinIndexedVector *regionSparse, 862 double region2[]) const 867 863 { 868 regionSparse>clear ();869 double *region = regionSparse>denseVector ();870 871 int *regionIndex = regionSparse>getIndices ();872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 864 regionSparse>clear(); 865 double *region = regionSparse>denseVector(); 866 int numberNonZero = 0; 867 int *regionIndex = regionSparse>getIndices(); 868 int i; 869 // set up linked lists at each depth 870 // stack2 is start, stack is next 871 int greatestDepth = 1; 872 for (i = 0; i < numberRows_; i++) { 873 double value = region2[i]; 874 if (value) { 875 region2[i] = 0.0; 876 region[i] = value; 877 regionIndex[numberNonZero++] = i; 878 int j = i; 879 int iDepth = depth_[j]; 880 if (iDepth > greatestDepth) 881 greatestDepth = iDepth; 882 // and back until marked 883 while (!mark_[j]) { 884 int iNext = stack2_[iDepth]; 885 stack2_[iDepth] = j; 886 stack_[j] = iNext; 887 mark_[j] = 1; 888 iDepth; 889 j = parent_[j]; 890 } 891 } 892 } 893 numberNonZero = 0; 894 for (; greatestDepth >= 0; greatestDepth) { 895 int iPivot = stack2_[greatestDepth]; 896 stack2_[greatestDepth] = 1; 897 while (iPivot >= 0) { 898 mark_[iPivot] = 0; 899 double pivotValue = region[iPivot]; 900 if (pivotValue) { 901 // put back now ? 902 int iBack = permuteBack_[iPivot]; 903 numberNonZero++; 904 int otherRow = parent_[iPivot]; 905 region2[iBack] = pivotValue * sign_[iPivot]; 906 region[iPivot] = 0.0; 907 region[otherRow] += pivotValue; 908 } 909 iPivot = stack_[iPivot]; 910 } 911 } 912 region[numberRows_] = 0.0; 913 return numberNonZero; 918 914 } 919 915 /* Updates one column transpose (BTRAN) … … 923 919 (the only exception is if you know input is dense e.g. dense objective) 924 920 returns number of nonzeros */ 925 int 926 ClpNetworkBasis::updateColumnTranspose ( CoinIndexedVector * regionSparse, 927 double region2[] ) const 921 int ClpNetworkBasis::updateColumnTranspose(CoinIndexedVector *regionSparse, 922 double region2[]) const 928 923 { 929 930 931 double *region = regionSparse>denseVector ();932 int *regionIndex = regionSparse>getIndices ();933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 smallestDepth = CoinMin(iDepth, smallestDepth);956 greatestDepth = CoinMax(iDepth, greatestDepth);957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 924 // permute in after copying 925 // so will end up in right place 926 double *region = regionSparse>denseVector(); 927 int *regionIndex = regionSparse>getIndices(); 928 int i; 929 int numberNonZero = 0; 930 CoinMemcpyN(region2, numberRows_, region); 931 for (i = 0; i < numberRows_; i++) { 932 double value = region[i]; 933 if (value) { 934 int k = permute_[i]; 935 region[i] = 0.0; 936 region2[k] = value; 937 regionIndex[numberNonZero++] = k; 938 mark_[k] = 1; 939 } 940 } 941 // copy back 942 // set up linked lists at each depth 943 // stack2 is start, stack is next 944 int greatestDepth = 1; 945 int smallestDepth = numberRows_; 946 for (i = 0; i < numberNonZero; i++) { 947 int j = regionIndex[i]; 948 // add in 949 int iDepth = depth_[j]; 950 smallestDepth = CoinMin(iDepth, smallestDepth); 951 greatestDepth = CoinMax(iDepth, greatestDepth); 952 int jNext = stack2_[iDepth]; 953 stack2_[iDepth] = j; 954 stack_[j] = jNext; 955 // and put all descendants on list 956 int iChild = descendant_[j]; 957 while (iChild >= 0) { 958 if (!mark_[iChild]) { 959 regionIndex[numberNonZero++] = iChild; 960 mark_[iChild] = 1; 961 } 962 iChild = rightSibling_[iChild]; 963 } 964 } 965 numberNonZero = 0; 966 region2[numberRows_] = 0.0; 967 int iDepth; 968 for (iDepth = smallestDepth; iDepth <= greatestDepth; iDepth++) { 969 int iPivot = stack2_[iDepth]; 970 stack2_[iDepth] = 1; 971 while (iPivot >= 0) { 972 mark_[iPivot] = 0; 973 double pivotValue = region2[iPivot]; 974 int otherRow = parent_[iPivot]; 975 double otherValue = region2[otherRow]; 976 pivotValue = sign_[iPivot] * pivotValue + otherValue; 977 region2[iPivot] = pivotValue; 978 if (pivotValue) 979 numberNonZero++; 980 iPivot = stack_[iPivot]; 981 } 982 } 983 return numberNonZero; 989 984 } 990 985 /* Updates one column (BTRAN) from region2 */ 991 int 992 ClpNetworkBasis::updateColumnTranspose ( CoinIndexedVector * regionSparse, 993 CoinIndexedVector * regionSparse2) const 986 int ClpNetworkBasis::updateColumnTranspose(CoinIndexedVector *regionSparse, 987 CoinIndexedVector *regionSparse2) const 994 988 { 995 996 997 regionSparse>clear ();998 double *region = regionSparse>denseVector ();999 double *region2 = regionSparse2>denseVector ();1000 int *regionIndex2 = regionSparse2>getIndices ();1001 int numberNonZero2 = regionSparse2>getNumElements ();1002 int *regionIndex = regionSparse>getIndices ();1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 smallestDepth = CoinMin(iDepth, smallestDepth);1027 greatestDepth = CoinMax(iDepth, greatestDepth);1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 smallestDepth = CoinMin(iDepth, smallestDepth);1046 greatestDepth = CoinMax(iDepth, greatestDepth);1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 smallestDepth = CoinMin(iDepth, smallestDepth);1110 greatestDepth = CoinMax(iDepth, greatestDepth);1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 smallestDepth = CoinMin(iDepth, smallestDepth);1129 greatestDepth = CoinMax(iDepth, greatestDepth);1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 989 // permute in  presume small number so copy back 990 // so will end up in right place 991 regionSparse>clear(); 992 double *region = regionSparse>denseVector(); 993 double *region2 = regionSparse2>denseVector(); 994 int *regionIndex2 = regionSparse2>getIndices(); 995 int numberNonZero2 = regionSparse2>getNumElements(); 996 int *regionIndex = regionSparse>getIndices(); 997 int i; 998 int numberNonZero = 0; 999 bool packed = regionSparse2>packedMode(); 1000 if (packed) { 1001 for (i = 0; i < numberNonZero2; i++) { 1002 int k = regionIndex2[i]; 1003 int j = permute_[k]; 1004 double value = region2[i]; 1005 region2[i] = 0.0; 1006 region[j] = value; 1007 mark_[j] = 1; 1008 regionIndex[numberNonZero++] = j; 1009 } 1010 // set up linked lists at each depth 1011 // stack2 is start, stack is next 1012 int greatestDepth = 1; 1013 int smallestDepth = numberRows_; 1014 //mark_[numberRows_]=1; 1015 for (i = 0; i < numberNonZero2; i++) { 1016 int j = regionIndex[i]; 1017 regionIndex2[i] = j; 1018 // add in 1019 int iDepth = depth_[j]; 1020 smallestDepth = CoinMin(iDepth, smallestDepth); 1021 greatestDepth = CoinMax(iDepth, greatestDepth); 1022 int jNext = stack2_[iDepth]; 1023 stack2_[iDepth] = j; 1024 stack_[j] = jNext; 1025 // and put all descendants on list 1026 int iChild = descendant_[j]; 1027 while (iChild >= 0) { 1028 if (!mark_[iChild]) { 1029 regionIndex2[numberNonZero++] = iChild; 1030 mark_[iChild] = 1; 1031 } 1032 iChild = rightSibling_[iChild]; 1033 } 1034 } 1035 for (; i < numberNonZero; i++) { 1036 int j = regionIndex2[i]; 1037 // add in 1038 int iDepth = depth_[j]; 1039 smallestDepth = CoinMin(iDepth, smallestDepth); 1040 greatestDepth = CoinMax(iDepth, greatestDepth); 1041 int jNext = stack2_[iDepth]; 1042 stack2_[iDepth] = j; 1043 stack_[j] = jNext; 1044 // and put all descendants on list 1045 int iChild = descendant_[j]; 1046 while (iChild >= 0) { 1047 if (!mark_[iChild]) { 1048 regionIndex2[numberNonZero++] = iChild; 1049 mark_[iChild] = 1; 1050 } 1051 iChild = rightSibling_[iChild]; 1052 } 1053 } 1054 numberNonZero2 = 0; 1055 region[numberRows_] = 0.0; 1056 int iDepth; 1057 for (iDepth = smallestDepth; iDepth <= greatestDepth; iDepth++) { 1058 int iPivot = stack2_[iDepth]; 1059 stack2_[iDepth] = 1; 1060 while (iPivot >= 0) { 1061 mark_[iPivot] = 0; 1062 double pivotValue = region[iPivot]; 1063 int otherRow = parent_[iPivot]; 1064 double otherValue = region[otherRow]; 1065 pivotValue = sign_[iPivot] * pivotValue + otherValue; 1066 region[iPivot] = pivotValue; 1067 if (pivotValue) { 1068 region2[numberNonZero2] = pivotValue; 1069 regionIndex2[numberNonZero2++] = iPivot; 1070 } 1071 iPivot = stack_[iPivot]; 1072 } 1073 } 1074 // zero out 1075 for (i = 0; i < numberNonZero2; i++) { 1076 int k = regionIndex2[i]; 1077 region[k] = 0.0; 1078 } 1079 } else { 1080 for (i = 0; i < numberNonZero2; i++) { 1081 int k = regionIndex2[i]; 1082 int j = permute_[k]; 1083 double value = region2[k]; 1084 region2[k] = 0.0; 1085 region[j] = value; 1086 mark_[j] = 1; 1087 regionIndex[numberNonZero++] = j; 1088 } 1089 // copy back 1090 // set up linked lists at each depth 1091 // stack2 is start, stack is next 1092 int greatestDepth = 1; 1093 int smallestDepth = numberRows_; 1094 //mark_[numberRows_]=1; 1095 for (i = 0; i < numberNonZero2; i++) { 1096 int j = regionIndex[i]; 1097 double value = region[j]; 1098 region[j] = 0.0; 1099 region2[j] = value; 1100 regionIndex2[i] = j; 1101 // add in 1102 int iDepth = depth_[j]; 1103 smallestDepth = CoinMin(iDepth, smallestDepth); 1104 greatestDepth = CoinMax(iDepth, greatestDepth); 1105 int jNext = stack2_[iDepth]; 1106 stack2_[iDepth] = j; 1107 stack_[j] = jNext; 1108 // and put all descendants on list 1109 int iChild = descendant_[j]; 1110 while (iChild >= 0) { 1111 if (!mark_[iChild]) { 1112 regionIndex2[numberNonZero++] = iChild; 1113 mark_[iChild] = 1; 1114 } 1115 iChild = rightSibling_[iChild]; 1116 } 1117 } 1118 for (; i < numberNonZero; i++) { 1119 int j = regionIndex2[i]; 1120 // add in 1121 int iDepth = depth_[j]; 1122 smallestDepth = CoinMin(iDepth, smallestDepth); 1123 greatestDepth = CoinMax(iDepth, greatestDepth); 1124 int jNext = stack2_[iDepth]; 1125 stack2_[iDepth] = j; 1126 stack_[j] = jNext; 1127 // and put all descendants on list 1128 int iChild = descendant_[j]; 1129 while (iChild >= 0) { 1130 if (!mark_[iChild]) { 1131 regionIndex2[numberNonZero++] = iChild; 1132 mark_[iChild] = 1; 1133 } 1134 iChild = rightSibling_[iChild]; 1135 } 1136 } 1137 numberNonZero2 = 0; 1138 region2[numberRows_] = 0.0; 1139 int iDepth; 1140 for (iDepth = smallestDepth; iDepth <= greatestDepth; iDepth++) { 1141 int iPivot = stack2_[iDepth]; 1142 stack2_[iDepth] = 1; 1143 while (iPivot >= 0) { 1144 mark_[iPivot] = 0; 1145 double pivotValue = region2[iPivot]; 1146 int otherRow = parent_[iPivot]; 1147 double otherValue = region2[otherRow]; 1148 pivotValue = sign_[iPivot] * pivotValue + otherValue; 1149 region2[iPivot] = pivotValue; 1150 if (pivotValue) 1151 regionIndex2[numberNonZero2++] = iPivot; 1152 iPivot = stack_[iPivot]; 1153 } 1154 } 1155 } 1156 regionSparse2>setNumElements(numberNonZero2); 1163 1157 #ifdef FULL_DEBUG 1164 1165 1166 1167 1168 assert(stack2_[i] == 1);1169 1170 1158 { 1159 int i; 1160 for (i = 0; i < numberRows_; i++) { 1161 assert(!mark_[i]); 1162 assert(stack2_[i] == 1); 1163 } 1164 } 1171 1165 #endif 1172 1166 return numberNonZero2; 1173 1167 } 1168 1169 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 1170 */
Note: See TracChangeset
for help on using the changeset viewer.