00001
00002
00003
00004 #include "CoinPragma.hpp"
00005 #include "ClpNetworkBasis.hpp"
00006 #include "CoinHelperFunctions.hpp"
00007 #include "ClpSimplex.hpp"
00008 #include "ClpMatrixBase.hpp"
00009 #include "CoinIndexedVector.hpp"
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 ClpNetworkBasis::ClpNetworkBasis ()
00020 {
00021 slackValue_=-1.0;
00022 numberRows_=0;
00023 numberColumns_=0;
00024 parent_ = NULL;
00025 descendant_ = NULL;
00026 pivot_ = NULL;
00027 rightSibling_ = NULL;
00028 leftSibling_ = NULL;
00029 sign_ = NULL;
00030 stack_ = NULL;
00031 permute_ = NULL;
00032 permuteBack_ = NULL;
00033 stack2_=NULL;
00034 depth_ = NULL;
00035 mark_ = NULL;
00036 model_=NULL;
00037 }
00038
00039 ClpNetworkBasis::ClpNetworkBasis(const ClpSimplex * model,
00040 int numberRows, const double * pivotRegion,
00041 const int * permuteBack,
00042 const int * startColumn,
00043 const int * numberInColumn,
00044 const int * indexRow, const double * element)
00045 {
00046 slackValue_=-1.0;
00047 numberRows_=numberRows;
00048 numberColumns_=numberRows;
00049 parent_ = new int [ numberRows_+1];
00050 descendant_ = new int [ numberRows_+1];
00051 pivot_ = new int [ numberRows_+1];
00052 rightSibling_ = new int [ numberRows_+1];
00053 leftSibling_ = new int [ numberRows_+1];
00054 sign_ = new double [ numberRows_+1];
00055 stack_ = new int [ numberRows_+1];
00056 stack2_ = new int[numberRows_+1];
00057 depth_ = new int[numberRows_+1];
00058 mark_ = new char[numberRows_+1];
00059 permute_ = new int [numberRows_ + 1];
00060 permuteBack_ = new int [numberRows_ + 1];
00061 int i;
00062 for (i=0;i<numberRows_+1;i++) {
00063 parent_[i]=-1;
00064 descendant_[i]=-1;
00065 pivot_[i]=-1;
00066 rightSibling_[i]=-1;
00067 leftSibling_[i]=-1;
00068 sign_[i]=-1.0;
00069 stack_[i]=-1;
00070 permute_[i]=i;
00071 permuteBack_[i]=i;
00072 stack2_[i]=-1;
00073 depth_[i]=-1;
00074 mark_[i]=0;
00075 }
00076 mark_[numberRows_]=1;
00077
00078
00079
00080
00081 int lastPivot=numberRows_;
00082 for (i=0;i<numberRows_;i++) {
00083 int iPivot = permuteBack[i];
00084 lastPivot=iPivot;
00085 double sign;
00086 if (pivotRegion[i]>0.0)
00087 sign = 1.0;
00088 else
00089 sign =-1.0;
00090 int other;
00091 if (numberInColumn[i]>0) {
00092 int iRow = indexRow[startColumn[i]];
00093 other = permuteBack[iRow];
00094
00095 } else {
00096 other = numberRows_;
00097 }
00098 sign_[iPivot] = sign;
00099 int iParent = other;
00100 parent_[iPivot] = other;
00101 if (descendant_[iParent]>=0) {
00102
00103 int iRight = descendant_[iParent];
00104 rightSibling_[iPivot]=iRight;
00105 leftSibling_[iRight]=iPivot;
00106 } else {
00107 rightSibling_[iPivot]=-1;
00108 }
00109 descendant_[iParent] = iPivot;
00110 leftSibling_[iPivot]=-1;
00111 }
00112
00113 int iPivot = numberRows_;
00114 int nStack = 1;
00115 stack_[0]=descendant_[numberRows_];
00116 depth_[numberRows_]=-1;
00117 while (nStack) {
00118
00119 int iNext = stack_[--nStack];
00120 if (iNext>=0) {
00121 depth_[iNext] = nStack;
00122 iPivot = iNext;
00123 int iRight = rightSibling_[iNext];
00124 stack_[nStack++] = iRight;
00125 if (descendant_[iNext]>=0)
00126 stack_[nStack++] = descendant_[iNext];
00127 }
00128 }
00129 model_=model;
00130 check();
00131 }
00132
00133
00134
00135
00136 ClpNetworkBasis::ClpNetworkBasis (const ClpNetworkBasis & rhs)
00137 {
00138 slackValue_=rhs.slackValue_;
00139 numberRows_=rhs.numberRows_;
00140 numberColumns_=rhs.numberColumns_;
00141 if (rhs.parent_) {
00142 parent_ = new int [numberRows_+1];
00143 memcpy(parent_,rhs.parent_,(numberRows_+1)*sizeof(int));
00144 } else {
00145 parent_ = NULL;
00146 }
00147 if (rhs.descendant_) {
00148 descendant_ = new int [numberRows_+1];
00149 memcpy(descendant_,rhs.descendant_,(numberRows_+1)*sizeof(int));
00150 } else {
00151 descendant_ = NULL;
00152 }
00153 if (rhs.pivot_) {
00154 pivot_ = new int [numberRows_+1];
00155 memcpy(pivot_,rhs.pivot_,(numberRows_+1)*sizeof(int));
00156 } else {
00157 pivot_ = NULL;
00158 }
00159 if (rhs.rightSibling_) {
00160 rightSibling_ = new int [numberRows_+1];
00161 memcpy(rightSibling_,rhs.rightSibling_,(numberRows_+1)*sizeof(int));
00162 } else {
00163 rightSibling_ = NULL;
00164 }
00165 if (rhs.leftSibling_) {
00166 leftSibling_ = new int [numberRows_+1];
00167 memcpy(leftSibling_,rhs.leftSibling_,(numberRows_+1)*sizeof(int));
00168 } else {
00169 leftSibling_ = NULL;
00170 }
00171 if (rhs.sign_) {
00172 sign_ = new double [numberRows_+1];
00173 memcpy(sign_,rhs.sign_,(numberRows_+1)*sizeof(double));
00174 } else {
00175 sign_ = NULL;
00176 }
00177 if (rhs.stack_) {
00178 stack_ = new int [numberRows_+1];
00179 memcpy(stack_,rhs.stack_,(numberRows_+1)*sizeof(int));
00180 } else {
00181 stack_ = NULL;
00182 }
00183 if (rhs.permute_) {
00184 permute_ = new int [numberRows_+1];
00185 memcpy(permute_,rhs.permute_,(numberRows_+1)*sizeof(int));
00186 } else {
00187 permute_ = NULL;
00188 }
00189 if (rhs.permuteBack_) {
00190 permuteBack_ = new int [numberRows_+1];
00191 memcpy(permuteBack_,rhs.permuteBack_,(numberRows_+1)*sizeof(int));
00192 } else {
00193 permuteBack_ = NULL;
00194 }
00195 if (rhs.stack2_) {
00196 stack2_ = new int [numberRows_+1];
00197 memcpy(stack2_,rhs.stack2_,(numberRows_+1)*sizeof(int));
00198 } else {
00199 stack2_ = NULL;
00200 }
00201 if (rhs.depth_) {
00202 depth_ = new int [numberRows_+1];
00203 memcpy(depth_,rhs.depth_,(numberRows_+1)*sizeof(int));
00204 } else {
00205 depth_ = NULL;
00206 }
00207 if (rhs.mark_) {
00208 mark_ = new char [numberRows_+1];
00209 memcpy(mark_,rhs.mark_,(numberRows_+1)*sizeof(char));
00210 } else {
00211 mark_ = NULL;
00212 }
00213 model_=rhs.model_;
00214 }
00215
00216
00217
00218
00219 ClpNetworkBasis::~ClpNetworkBasis ()
00220 {
00221 delete [] parent_;
00222 delete [] descendant_;
00223 delete [] pivot_;
00224 delete [] rightSibling_;
00225 delete [] leftSibling_;
00226 delete [] sign_;
00227 delete [] stack_;
00228 delete [] permute_;
00229 delete [] permuteBack_;
00230 delete [] stack2_;
00231 delete [] depth_;
00232 delete [] mark_;
00233 }
00234
00235
00236
00237
00238 ClpNetworkBasis &
00239 ClpNetworkBasis::operator=(const ClpNetworkBasis& rhs)
00240 {
00241 if (this != &rhs) {
00242 delete [] parent_;
00243 delete [] descendant_;
00244 delete [] pivot_;
00245 delete [] rightSibling_;
00246 delete [] leftSibling_;
00247 delete [] sign_;
00248 delete [] stack_;
00249 delete [] permute_;
00250 delete [] permuteBack_;
00251 delete [] stack2_;
00252 delete [] depth_;
00253 delete [] mark_;
00254 slackValue_=rhs.slackValue_;
00255 numberRows_=rhs.numberRows_;
00256 numberColumns_=rhs.numberColumns_;
00257 if (rhs.parent_) {
00258 parent_ = new int [numberRows_+1];
00259 memcpy(parent_,rhs.parent_,(numberRows_+1)*sizeof(int));
00260 } else {
00261 parent_ = NULL;
00262 }
00263 if (rhs.descendant_) {
00264 descendant_ = new int [numberRows_+1];
00265 memcpy(descendant_,rhs.descendant_,(numberRows_+1)*sizeof(int));
00266 } else {
00267 descendant_ = NULL;
00268 }
00269 if (rhs.pivot_) {
00270 pivot_ = new int [numberRows_+1];
00271 memcpy(pivot_,rhs.pivot_,(numberRows_+1)*sizeof(int));
00272 } else {
00273 pivot_ = NULL;
00274 }
00275 if (rhs.rightSibling_) {
00276 rightSibling_ = new int [numberRows_+1];
00277 memcpy(rightSibling_,rhs.rightSibling_,(numberRows_+1)*sizeof(int));
00278 } else {
00279 rightSibling_ = NULL;
00280 }
00281 if (rhs.leftSibling_) {
00282 leftSibling_ = new int [numberRows_+1];
00283 memcpy(leftSibling_,rhs.leftSibling_,(numberRows_+1)*sizeof(int));
00284 } else {
00285 leftSibling_ = NULL;
00286 }
00287 if (rhs.sign_) {
00288 sign_ = new double [numberRows_+1];
00289 memcpy(sign_,rhs.sign_,(numberRows_+1)*sizeof(double));
00290 } else {
00291 sign_ = NULL;
00292 }
00293 if (rhs.stack_) {
00294 stack_ = new int [numberRows_+1];
00295 memcpy(stack_,rhs.stack_,(numberRows_+1)*sizeof(int));
00296 } else {
00297 stack_ = NULL;
00298 }
00299 if (rhs.permute_) {
00300 permute_ = new int [numberRows_+1];
00301 memcpy(permute_,rhs.permute_,(numberRows_+1)*sizeof(int));
00302 } else {
00303 permute_ = NULL;
00304 }
00305 if (rhs.permuteBack_) {
00306 permuteBack_ = new int [numberRows_+1];
00307 memcpy(permuteBack_,rhs.permuteBack_,(numberRows_+1)*sizeof(int));
00308 } else {
00309 permuteBack_ = NULL;
00310 }
00311 if (rhs.stack2_) {
00312 stack2_ = new int [numberRows_+1];
00313 memcpy(stack2_,rhs.stack2_,(numberRows_+1)*sizeof(int));
00314 } else {
00315 stack2_ = NULL;
00316 }
00317 if (rhs.depth_) {
00318 depth_ = new int [numberRows_+1];
00319 memcpy(depth_,rhs.depth_,(numberRows_+1)*sizeof(int));
00320 } else {
00321 depth_ = NULL;
00322 }
00323 if (rhs.mark_) {
00324 mark_ = new char [numberRows_+1];
00325 memcpy(mark_,rhs.mark_,(numberRows_+1)*sizeof(char));
00326 } else {
00327 mark_ = NULL;
00328 }
00329 }
00330 return *this;
00331 }
00332
00333 void ClpNetworkBasis::check()
00334 {
00335
00336 int iPivot = numberRows_;
00337 int nStack = 1;
00338 stack_[0]=descendant_[numberRows_];
00339 depth_[numberRows_]=-1;
00340 while (nStack) {
00341
00342 int iNext = stack_[--nStack];
00343 if (iNext>=0) {
00344
00345 depth_[iNext] = nStack;
00346 iPivot = iNext;
00347 int iRight = rightSibling_[iNext];
00348 stack_[nStack++] = iRight;
00349 if (descendant_[iNext]>=0)
00350 stack_[nStack++] = descendant_[iNext];
00351 }
00352 }
00353 }
00354
00355 void ClpNetworkBasis::print()
00356 {
00357 int i;
00358 printf(" parent descendant left right sign depth\n");
00359 for (i=0;i<numberRows_+1;i++)
00360 printf("%4d %7d %8d %7d %7d %5g %7d\n",
00361 i,parent_[i],descendant_[i],leftSibling_[i],rightSibling_[i],
00362 sign_[i],depth_[i]);
00363 }
00364
00365
00366
00367 int
00368 ClpNetworkBasis::replaceColumn ( CoinIndexedVector * regionSparse,
00369 int pivotRow)
00370 {
00371
00372
00373
00374 assert (!regionSparse->getNumElements());
00375 model_->unpack(regionSparse, model_->sequenceIn());
00376
00377
00378
00379 int * indices = regionSparse->getIndices();
00380 int iRow0 = indices[0];
00381 int iRow1;
00382 if (regionSparse->getNumElements()==2)
00383 iRow1 = indices[1];
00384 else
00385 iRow1 = numberRows_;
00386 double sign = -regionSparse->denseVector()[iRow0];
00387 regionSparse->clear();
00388
00389 model_->unpack(regionSparse,model_->pivotVariable()[pivotRow]);
00390 int jRow0 = indices[0];
00391 int jRow1;
00392 if (regionSparse->getNumElements()==2)
00393 jRow1 = indices[1];
00394 else
00395 jRow1 = numberRows_;
00396 regionSparse->clear();
00397
00398
00399 #ifdef FULL_DEBUG
00400 printf ("irow %d %d, jrow %d %d\n",
00401 iRow0,iRow1,jRow0,jRow1);
00402 #endif
00403 if (parent_[jRow0]==jRow1) {
00404 int newPivot = jRow0;
00405 if (newPivot!=pivotRow) {
00406 #ifdef FULL_DEBUG
00407 printf("pivot row of %d permuted to %d\n",pivotRow,newPivot);
00408 #endif
00409 pivotRow = newPivot;
00410 }
00411 } else {
00412
00413 int newPivot = jRow1;
00414 if (newPivot!=pivotRow) {
00415 #ifdef FULL_DEBUG
00416 printf("pivot row of %d permuted to %d\n",pivotRow,newPivot);
00417 #endif
00418 pivotRow = newPivot;
00419 }
00420 }
00421 bool extraPrint = (model_->numberIterations()>-3)&&
00422 (model_->logLevel()>10);
00423 if (extraPrint)
00424 print();
00425 #ifdef FULL_DEBUG
00426 printf("In %d (region= %g, stored %g) %d (%g) pivoting on %d (%g)\n",
00427 iRow1, sign, sign_[iRow1], iRow0, sign_[iRow0] ,pivotRow,sign_[pivotRow]);
00428 #endif
00429
00430 int kRow = -1;
00431 int jRow = iRow1;
00432 while (jRow!=numberRows_) {
00433 if (jRow==pivotRow) {
00434 kRow = iRow1;
00435 break;
00436 } else {
00437 jRow = parent_[jRow];
00438 }
00439 }
00440 if (kRow<0) {
00441 jRow = iRow0;
00442 while (jRow!=numberRows_) {
00443 if (jRow==pivotRow) {
00444 kRow = iRow0;
00445 break;
00446 } else {
00447 jRow = parent_[jRow];
00448 }
00449 }
00450 }
00451
00452 if (iRow0==kRow) {
00453 iRow0 = iRow1;
00454 iRow1 = kRow;
00455 sign = -sign;
00456 }
00457
00458
00459
00460 int nStack=1;
00461 stack_[0]=iRow0;
00462 while (kRow!=pivotRow) {
00463 stack_[nStack++] = kRow;
00464 if (sign*sign_[kRow]<0.0) {
00465 sign_[kRow]= -sign_[kRow];
00466 } else {
00467 sign = -sign;
00468 }
00469 kRow = parent_[kRow];
00470
00471 }
00472 stack_[nStack++]=pivotRow;
00473 if (sign*sign_[pivotRow]<0.0) {
00474 sign_[pivotRow]= -sign_[pivotRow];
00475 } else {
00476 sign = -sign;
00477 }
00478 int iParent = parent_[pivotRow];
00479 while (nStack>1) {
00480 int iLeft;
00481 int iRight;
00482 kRow = stack_[--nStack];
00483 int newParent = stack_[nStack-1];
00484 #ifdef FULL_DEBUG
00485 printf("row %d, old parent %d, new parent %d, pivotrow %d\n",kRow,
00486 iParent,newParent,pivotRow);
00487 #endif
00488 int i1 = permuteBack_[pivotRow];
00489 int i2 = permuteBack_[kRow];
00490 permuteBack_[pivotRow]=i2;
00491 permuteBack_[kRow]=i1;
00492
00493 permute_[i1]=kRow;
00494 permute_[i2]=pivotRow;
00495 pivotRow = kRow;
00496
00497 iLeft = leftSibling_[kRow];
00498 iRight = rightSibling_[kRow];
00499 if (iLeft<0) {
00500
00501 if (iRight>=0) {
00502 leftSibling_[iRight]=iLeft;
00503 descendant_[iParent]=iRight;
00504 } else {
00505 #ifdef FULL_DEBUG
00506 printf("Saying %d (old parent of %d) has no descendants\n",
00507 iParent, kRow);
00508 #endif
00509 descendant_[iParent]=-1;
00510 }
00511 } else {
00512
00513 rightSibling_[iLeft] = iRight;
00514 if (iRight>=0)
00515 leftSibling_[iRight]=iLeft;
00516 }
00517 leftSibling_[kRow]=-1;
00518 rightSibling_[kRow]=-1;
00519
00520
00521
00522 if (descendant_[newParent]>=0) {
00523
00524 int jRight = descendant_[newParent];
00525 rightSibling_[kRow]=jRight;
00526 leftSibling_[jRight]=kRow;
00527 } else {
00528 rightSibling_[kRow]=-1;
00529 }
00530 descendant_[newParent] = kRow;
00531 leftSibling_[kRow]=-1;
00532 parent_[kRow]=newParent;
00533
00534 iParent = kRow;
00535 }
00536
00537
00538 {
00539 int iPivot = stack_[1];
00540 int iDepth=depth_[parent_[iPivot]];
00541 iDepth ++;
00542 int nStack = 1;
00543 stack_[0]=iPivot;
00544 while (nStack) {
00545
00546 int iNext = stack_[--nStack];
00547 if (iNext>=0) {
00548
00549 depth_[iNext]=nStack+iDepth;
00550 stack_[nStack++] = rightSibling_[iNext];
00551 if (descendant_[iNext]>=0)
00552 stack_[nStack++] = descendant_[iNext];
00553 }
00554 }
00555 }
00556 if (extraPrint)
00557 print();
00558
00559 return 0;
00560 }
00561
00562
00563 double
00564 ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse,
00565 CoinIndexedVector * regionSparse2,
00566 int pivotRow)
00567 {
00568 regionSparse->clear ( );
00569 double *region = regionSparse->denseVector ( );
00570 double *region2 = regionSparse2->denseVector ( );
00571 int *regionIndex2 = regionSparse2->getIndices ( );
00572 int numberNonZero = regionSparse2->getNumElements ( );
00573 int *regionIndex = regionSparse->getIndices ( );
00574 int i;
00575 bool doTwo = (numberNonZero==2);
00576 int i0=-1;
00577 int i1=-1;
00578 if (doTwo) {
00579 i0 = regionIndex2[0];
00580 i1 = regionIndex2[1];
00581 }
00582 double returnValue=0.0;
00583 bool packed = regionSparse2->packedMode();
00584 if (packed) {
00585 if (doTwo&®ion2[0]*region2[1]<0.0) {
00586 region[i0] = region2[0];
00587 region2[0]=0.0;
00588 region[i1] = region2[1];
00589 region2[1]=0.0;
00590 int iDepth0 = depth_[i0];
00591 int iDepth1 = depth_[i1];
00592 if (iDepth1>iDepth0) {
00593 int temp = i0;
00594 i0 = i1;
00595 i1 = temp;
00596 temp = iDepth0;
00597 iDepth0 = iDepth1;
00598 iDepth1 = temp;
00599 }
00600 numberNonZero=0;
00601 if (pivotRow<0) {
00602 while (iDepth0>iDepth1) {
00603 double pivotValue = region[i0];
00604
00605 int iBack = permuteBack_[i0];
00606 region2[numberNonZero] = pivotValue*sign_[i0];
00607 regionIndex2[numberNonZero++] = iBack;
00608 int otherRow = parent_[i0];
00609 region[i0] =0.0;
00610 region[otherRow] += pivotValue;
00611 iDepth0--;
00612 i0 = otherRow;
00613 }
00614 while (i0!=i1) {
00615 double pivotValue = region[i0];
00616
00617 int iBack = permuteBack_[i0];
00618 region2[numberNonZero] = pivotValue*sign_[i0];
00619 regionIndex2[numberNonZero++] = iBack;
00620 int otherRow = parent_[i0];
00621 region[i0] =0.0;
00622 region[otherRow] += pivotValue;
00623 i0 = otherRow;
00624 double pivotValue1 = region[i1];
00625
00626 int iBack1 = permuteBack_[i1];
00627 region2[numberNonZero] = pivotValue1*sign_[i1];
00628 regionIndex2[numberNonZero++] = iBack1;
00629 int otherRow1 = parent_[i1];
00630 region[i1] =0.0;
00631 region[otherRow1] += pivotValue1;
00632 i1 = otherRow1;
00633 }
00634 } else {
00635 while (iDepth0>iDepth1) {
00636 double pivotValue = region[i0];
00637
00638 int iBack = permuteBack_[i0];
00639 double value = pivotValue*sign_[i0];
00640 region2[numberNonZero] = value;
00641 regionIndex2[numberNonZero++] = iBack;
00642 if (iBack==pivotRow)
00643 returnValue = value;
00644 int otherRow = parent_[i0];
00645 region[i0] =0.0;
00646 region[otherRow] += pivotValue;
00647 iDepth0--;
00648 i0 = otherRow;
00649 }
00650 while (i0!=i1) {
00651 double pivotValue = region[i0];
00652
00653 int iBack = permuteBack_[i0];
00654 double value = pivotValue*sign_[i0];
00655 region2[numberNonZero] = value;
00656 regionIndex2[numberNonZero++] = iBack;
00657 if (iBack==pivotRow)
00658 returnValue = value;
00659 int otherRow = parent_[i0];
00660 region[i0] =0.0;
00661 region[otherRow] += pivotValue;
00662 i0 = otherRow;
00663 double pivotValue1 = region[i1];
00664
00665 int iBack1 = permuteBack_[i1];
00666 value = pivotValue1*sign_[i1];
00667 region2[numberNonZero] = value;
00668 regionIndex2[numberNonZero++] = iBack1;
00669 if (iBack1==pivotRow)
00670 returnValue = value;
00671 int otherRow1 = parent_[i1];
00672 region[i1] =0.0;
00673 region[otherRow1] += pivotValue1;
00674 i1 = otherRow1;
00675 }
00676 }
00677 } else {
00678
00679
00680 int greatestDepth=-1;
00681
00682 for (i=0;i<numberNonZero;i++) {
00683 int j = regionIndex2[i];
00684 double value = region2[i];
00685 region2[i]=0.0;
00686 region[j]=value;
00687 regionIndex[i]=j;
00688 int iDepth = depth_[j];
00689 if (iDepth>greatestDepth)
00690 greatestDepth = iDepth;
00691
00692 while (!mark_[j]) {
00693 int iNext = stack2_[iDepth];
00694 stack2_[iDepth]=j;
00695 stack_[j]=iNext;
00696 mark_[j]=1;
00697 iDepth--;
00698 j=parent_[j];
00699 }
00700 }
00701 numberNonZero=0;
00702 if (pivotRow<0) {
00703 for (;greatestDepth>=0; greatestDepth--) {
00704 int iPivot = stack2_[greatestDepth];
00705 stack2_[greatestDepth]=-1;
00706 while (iPivot>=0) {
00707 mark_[iPivot]=0;
00708 double pivotValue = region[iPivot];
00709 if (pivotValue) {
00710
00711 int iBack = permuteBack_[iPivot];
00712 region2[numberNonZero] = pivotValue*sign_[iPivot];
00713 regionIndex2[numberNonZero++] = iBack;
00714 int otherRow = parent_[iPivot];
00715 region[iPivot] =0.0;
00716 region[otherRow] += pivotValue;
00717 }
00718 iPivot = stack_[iPivot];
00719 }
00720 }
00721 } else {
00722 for (;greatestDepth>=0; greatestDepth--) {
00723 int iPivot = stack2_[greatestDepth];
00724 stack2_[greatestDepth]=-1;
00725 while (iPivot>=0) {
00726 mark_[iPivot]=0;
00727 double pivotValue = region[iPivot];
00728 if (pivotValue) {
00729
00730 int iBack = permuteBack_[iPivot];
00731 double value = pivotValue*sign_[iPivot];
00732 region2[numberNonZero] = value;
00733 regionIndex2[numberNonZero++] = iBack;
00734 if (iBack==pivotRow)
00735 returnValue = value;
00736 int otherRow = parent_[iPivot];
00737 region[iPivot] =0.0;
00738 region[otherRow] += pivotValue;
00739 }
00740 iPivot = stack_[iPivot];
00741 }
00742 }
00743 }
00744 }
00745 } else {
00746 if (doTwo&®ion2[i0]*region2[i1]<0.0) {
00747
00748 region[i0] = region2[i0];
00749 region2[i0]=0.0;
00750 region[i1] = region2[i1];
00751 region2[i1]=0.0;
00752 int iDepth0 = depth_[i0];
00753 int iDepth1 = depth_[i1];
00754 if (iDepth1>iDepth0) {
00755 int temp = i0;
00756 i0 = i1;
00757 i1 = temp;
00758 temp = iDepth0;
00759 iDepth0 = iDepth1;
00760 iDepth1 = temp;
00761 }
00762 numberNonZero=0;
00763 while (iDepth0>iDepth1) {
00764 double pivotValue = region[i0];
00765
00766 int iBack = permuteBack_[i0];
00767 regionIndex2[numberNonZero++] = iBack;
00768 int otherRow = parent_[i0];
00769 region2[iBack] = pivotValue*sign_[i0];
00770 region[i0] =0.0;
00771 region[otherRow] += pivotValue;
00772 iDepth0--;
00773 i0 = otherRow;
00774 }
00775 while (i0!=i1) {
00776 double pivotValue = region[i0];
00777
00778 int iBack = permuteBack_[i0];
00779 regionIndex2[numberNonZero++] = iBack;
00780 int otherRow = parent_[i0];
00781 region2[iBack] = pivotValue*sign_[i0];
00782 region[i0] =0.0;
00783 region[otherRow] += pivotValue;
00784 i0 = otherRow;
00785 double pivotValue1 = region[i1];
00786
00787 int iBack1 = permuteBack_[i1];
00788 regionIndex2[numberNonZero++] = iBack1;
00789 int otherRow1 = parent_[i1];
00790 region2[iBack1] = pivotValue1*sign_[i1];
00791 region[i1] =0.0;
00792 region[otherRow1] += pivotValue1;
00793 i1 = otherRow1;
00794 }
00795 } else {
00796
00797
00798 int greatestDepth=-1;
00799
00800 for (i=0;i<numberNonZero;i++) {
00801 int j = regionIndex2[i];
00802 double value = region2[j];
00803 region2[j]=0.0;
00804 region[j]=value;
00805 regionIndex[i]=j;
00806 int iDepth = depth_[j];
00807 if (iDepth>greatestDepth)
00808 greatestDepth = iDepth;
00809
00810 while (!mark_[j]) {
00811 int iNext = stack2_[iDepth];
00812 stack2_[iDepth]=j;
00813 stack_[j]=iNext;
00814 mark_[j]=1;
00815 iDepth--;
00816 j=parent_[j];
00817 }
00818 }
00819 numberNonZero=0;
00820 for (;greatestDepth>=0; greatestDepth--) {
00821 int iPivot = stack2_[greatestDepth];
00822 stack2_[greatestDepth]=-1;
00823 while (iPivot>=0) {
00824 mark_[iPivot]=0;
00825 double pivotValue = region[iPivot];
00826 if (pivotValue) {
00827
00828 int iBack = permuteBack_[iPivot];
00829 regionIndex2[numberNonZero++] = iBack;
00830 int otherRow = parent_[iPivot];
00831 region2[iBack] = pivotValue*sign_[iPivot];
00832 region[iPivot] =0.0;
00833 region[otherRow] += pivotValue;
00834 }
00835 iPivot = stack_[iPivot];
00836 }
00837 }
00838 }
00839 if (pivotRow>=0)
00840 returnValue = region2[pivotRow];
00841 }
00842 region[numberRows_]=0.0;
00843 regionSparse2->setNumElements(numberNonZero);
00844 #ifdef FULL_DEBUG
00845 {
00846 int i;
00847 for (i=0;i<numberRows_;i++) {
00848 assert(!mark_[i]);
00849 assert (stack2_[i]==-1);
00850 }
00851 }
00852 #endif
00853 return returnValue;
00854 }
00855
00856
00857
00858
00859
00860 int
00861 ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse,
00862 double region2[] ) const
00863 {
00864 regionSparse->clear ( );
00865 double *region = regionSparse->denseVector ( );
00866 int numberNonZero = 0;
00867 int *regionIndex = regionSparse->getIndices ( );
00868 int i;
00869
00870
00871 int greatestDepth=-1;
00872 for (i=0;i<numberRows_;i++) {
00873 double value = region2[i];
00874 if (value) {
00875 region2[i]=0.0;
00876 region[i]=value;
00877 regionIndex[numberNonZero++]=i;
00878 int j=i;
00879 int iDepth = depth_[j];
00880 if (iDepth>greatestDepth)
00881 greatestDepth = iDepth;
00882
00883 while (!mark_[j]) {
00884 int iNext = stack2_[iDepth];
00885 stack2_[iDepth]=j;
00886 stack_[j]=iNext;
00887 mark_[j]=1;
00888 iDepth--;
00889 j=parent_[j];
00890 }
00891 }
00892 }
00893 numberNonZero=0;
00894 for (;greatestDepth>=0; greatestDepth--) {
00895 int iPivot = stack2_[greatestDepth];
00896 stack2_[greatestDepth]=-1;
00897 while (iPivot>=0) {
00898 mark_[iPivot]=0;
00899 double pivotValue = region[iPivot];
00900 if (pivotValue) {
00901
00902 int iBack = permuteBack_[iPivot];
00903 numberNonZero++;
00904 int otherRow = parent_[iPivot];
00905 region2[iBack] = pivotValue*sign_[iPivot];
00906 region[iPivot] =0.0;
00907 region[otherRow] += pivotValue;
00908 }
00909 iPivot = stack_[iPivot];
00910 }
00911 }
00912 region[numberRows_]=0.0;
00913 return numberNonZero;
00914 }
00915
00916
00917
00918
00919
00920
00921 int
00922 ClpNetworkBasis::updateColumnTranspose ( CoinIndexedVector * regionSparse,
00923 double region2[] ) const
00924 {
00925
00926
00927 double *region = regionSparse->denseVector ( );
00928 int *regionIndex = regionSparse->getIndices ( );
00929 int i;
00930 int numberNonZero=0;
00931 memcpy(region,region2,numberRows_*sizeof(double));
00932 for (i=0;i<numberRows_;i++) {
00933 double value = region[i];
00934 if (value) {
00935 int k = permute_[i];
00936 region[i]=0.0;
00937 region2[k]=value;
00938 regionIndex[numberNonZero++]=k;
00939 mark_[k]=1;
00940 }
00941 }
00942
00943
00944
00945 int greatestDepth=-1;
00946 int smallestDepth=numberRows_;
00947 for (i=0;i<numberNonZero;i++) {
00948 int j = regionIndex[i];
00949
00950 int iDepth = depth_[j];
00951 smallestDepth = min(iDepth,smallestDepth) ;
00952 greatestDepth = max(iDepth,greatestDepth) ;
00953 int jNext = stack2_[iDepth];
00954 stack2_[iDepth]=j;
00955 stack_[j]=jNext;
00956
00957 int iChild = descendant_[j];
00958 while (iChild>=0) {
00959 if (!mark_[iChild]) {
00960 regionIndex[numberNonZero++] = iChild;
00961 mark_[iChild]=1;
00962 }
00963 iChild = rightSibling_[iChild];
00964 }
00965 }
00966 numberNonZero=0;
00967 region2[numberRows_]=0.0;
00968 int iDepth;
00969 for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
00970 int iPivot = stack2_[iDepth];
00971 stack2_[iDepth]=-1;
00972 while (iPivot>=0) {
00973 mark_[iPivot]=0;
00974 double pivotValue = region2[iPivot];
00975 int otherRow = parent_[iPivot];
00976 double otherValue = region2[otherRow];
00977 pivotValue = sign_[iPivot]*pivotValue+otherValue;
00978 region2[iPivot]=pivotValue;
00979 if (pivotValue)
00980 numberNonZero++;
00981 iPivot = stack_[iPivot];
00982 }
00983 }
00984 return numberNonZero;
00985 }
00986
00987 int
00988 ClpNetworkBasis::updateColumnTranspose ( CoinIndexedVector * regionSparse,
00989 CoinIndexedVector * regionSparse2) const
00990 {
00991
00992
00993 regionSparse->clear ( );
00994 double *region = regionSparse->denseVector ( );
00995 double *region2 = regionSparse2->denseVector ( );
00996 int *regionIndex2 = regionSparse2->getIndices ( );
00997 int numberNonZero2 = regionSparse2->getNumElements ( );
00998 int *regionIndex = regionSparse->getIndices ( );
00999 int i;
01000 int numberNonZero=0;
01001 bool packed = regionSparse2->packedMode();
01002 if (packed) {
01003 for (i=0;i<numberNonZero2;i++) {
01004 int k = regionIndex2[i];
01005 int j = permute_[k];
01006 double value = region2[i];
01007 region2[i]=0.0;
01008 region[j]=value;
01009 mark_[j]=1;
01010 regionIndex[numberNonZero++]=j;
01011 }
01012
01013
01014 int greatestDepth=-1;
01015 int smallestDepth=numberRows_;
01016
01017 for (i=0;i<numberNonZero2;i++) {
01018 int j = regionIndex[i];
01019 regionIndex2[i]=j;
01020
01021 int iDepth = depth_[j];
01022 smallestDepth = min(iDepth,smallestDepth) ;
01023 greatestDepth = max(iDepth,greatestDepth) ;
01024 int jNext = stack2_[iDepth];
01025 stack2_[iDepth]=j;
01026 stack_[j]=jNext;
01027
01028 int iChild = descendant_[j];
01029 while (iChild>=0) {
01030 if (!mark_[iChild]) {
01031 regionIndex2[numberNonZero++] = iChild;
01032 mark_[iChild]=1;
01033 }
01034 iChild = rightSibling_[iChild];
01035 }
01036 }
01037 for (;i<numberNonZero;i++) {
01038 int j = regionIndex2[i];
01039
01040 int iDepth = depth_[j];
01041 smallestDepth = min(iDepth,smallestDepth) ;
01042 greatestDepth = max(iDepth,greatestDepth) ;
01043 int jNext = stack2_[iDepth];
01044 stack2_[iDepth]=j;
01045 stack_[j]=jNext;
01046
01047 int iChild = descendant_[j];
01048 while (iChild>=0) {
01049 if (!mark_[iChild]) {
01050 regionIndex2[numberNonZero++] = iChild;
01051 mark_[iChild]=1;
01052 }
01053 iChild = rightSibling_[iChild];
01054 }
01055 }
01056 numberNonZero2=0;
01057 region[numberRows_]=0.0;
01058 int iDepth;
01059 for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
01060 int iPivot = stack2_[iDepth];
01061 stack2_[iDepth]=-1;
01062 while (iPivot>=0) {
01063 mark_[iPivot]=0;
01064 double pivotValue = region[iPivot];
01065 int otherRow = parent_[iPivot];
01066 double otherValue = region[otherRow];
01067 pivotValue = sign_[iPivot]*pivotValue+otherValue;
01068 region[iPivot]=pivotValue;
01069 if (pivotValue) {
01070 region2[numberNonZero2]=pivotValue;
01071 regionIndex2[numberNonZero2++]=iPivot;
01072 }
01073 iPivot = stack_[iPivot];
01074 }
01075 }
01076
01077 for (i=0;i<numberNonZero2;i++) {
01078 int k = regionIndex2[i];
01079 region[k]=0.0;
01080 }
01081 } else {
01082 for (i=0;i<numberNonZero2;i++) {
01083 int k = regionIndex2[i];
01084 int j = permute_[k];
01085 double value = region2[k];
01086 region2[k]=0.0;
01087 region[j]=value;
01088 mark_[j]=1;
01089 regionIndex[numberNonZero++]=j;
01090 }
01091
01092
01093
01094 int greatestDepth=-1;
01095 int smallestDepth=numberRows_;
01096
01097 for (i=0;i<numberNonZero2;i++) {
01098 int j = regionIndex[i];
01099 double value = region[j];
01100 region[j]=0.0;
01101 region2[j]=value;
01102 regionIndex2[i]=j;
01103
01104 int iDepth = depth_[j];
01105 smallestDepth = min(iDepth,smallestDepth) ;
01106 greatestDepth = max(iDepth,greatestDepth) ;
01107 int jNext = stack2_[iDepth];
01108 stack2_[iDepth]=j;
01109 stack_[j]=jNext;
01110
01111 int iChild = descendant_[j];
01112 while (iChild>=0) {
01113 if (!mark_[iChild]) {
01114 regionIndex2[numberNonZero++] = iChild;
01115 mark_[iChild]=1;
01116 }
01117 iChild = rightSibling_[iChild];
01118 }
01119 }
01120 for (;i<numberNonZero;i++) {
01121 int j = regionIndex2[i];
01122
01123 int iDepth = depth_[j];
01124 smallestDepth = min(iDepth,smallestDepth) ;
01125 greatestDepth = max(iDepth,greatestDepth) ;
01126 int jNext = stack2_[iDepth];
01127 stack2_[iDepth]=j;
01128 stack_[j]=jNext;
01129
01130 int iChild = descendant_[j];
01131 while (iChild>=0) {
01132 if (!mark_[iChild]) {
01133 regionIndex2[numberNonZero++] = iChild;
01134 mark_[iChild]=1;
01135 }
01136 iChild = rightSibling_[iChild];
01137 }
01138 }
01139 numberNonZero2=0;
01140 region2[numberRows_]=0.0;
01141 int iDepth;
01142 for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
01143 int iPivot = stack2_[iDepth];
01144 stack2_[iDepth]=-1;
01145 while (iPivot>=0) {
01146 mark_[iPivot]=0;
01147 double pivotValue = region2[iPivot];
01148 int otherRow = parent_[iPivot];
01149 double otherValue = region2[otherRow];
01150 pivotValue = sign_[iPivot]*pivotValue+otherValue;
01151 region2[iPivot]=pivotValue;
01152 if (pivotValue)
01153 regionIndex2[numberNonZero2++]=iPivot;
01154 iPivot = stack_[iPivot];
01155 }
01156 }
01157 }
01158 regionSparse2->setNumElements(numberNonZero2);
01159 #ifdef FULL_DEBUG
01160 {
01161 int i;
01162 for (i=0;i<numberRows_;i++) {
01163 assert(!mark_[i]);
01164 assert (stack2_[i]==-1);
01165 }
01166 }
01167 #endif
01168 return numberNonZero2;
01169 }