00001
00002
00003
00004 #if defined(_MSC_VER)
00005
00006 # pragma warning(disable:4786)
00007 #endif
00008
00009 #include <cassert>
00010 #include <cstdio>
00011
00012 #include "CoinFactorization.hpp"
00013 #include "CoinIndexedVector.hpp"
00014 #include "CoinHelperFunctions.hpp"
00015 #include <stdio.h>
00016 #include <iostream>
00017 #if DENSE_CODE==1
00018
00019 extern "C" int dgetrs_(const char *trans, const int *n, const int *nrhs,
00020 const double *a, const int *lda, const int *ipiv, double *b,
00021 const int * ldb, int *info);
00022 #endif
00023
00024
00025 #define BITS_PER_CHECK 8
00026 #define CHECK_SHIFT 3
00027 typedef unsigned char CoinCheckZero;
00028
00029
00030
00031
00032 void
00033 CoinFactorization::updateColumnLDensish ( CoinIndexedVector * regionSparse ,
00034 int * regionIndex)
00035 const
00036 {
00037 double *region = regionSparse->denseVector ( );
00038 int number = regionSparse->getNumElements ( );
00039 int numberNonZero;
00040 double tolerance = zeroTolerance_;
00041
00042 numberNonZero = 0;
00043 int k;
00044 int i , iPivot;
00045
00046 const CoinBigIndex *startColumn = startColumnL_;
00047 const int *indexRow = indexRowL_;
00048 const double *element = elementL_;
00049 int last = baseL_ + numberL_;
00050 int smallestIndex = numberRowsExtra_;
00051
00052 for (k=0;k<number;k++) {
00053 iPivot=regionIndex[k];
00054 if (iPivot<baseL_)
00055 regionIndex[numberNonZero++]=iPivot;
00056 else
00057 smallestIndex = min(iPivot,smallestIndex);
00058 }
00059
00060 for ( i = smallestIndex; i < last; i++ ) {
00061 double pivotValue = region[i];
00062 CoinBigIndex start = startColumn[i];
00063 CoinBigIndex end = startColumn[i + 1];
00064
00065 if ( fabs(pivotValue) > tolerance ) {
00066 CoinBigIndex j;
00067 for ( j = start; j < end; j ++ ) {
00068 int iRow0 = indexRow[j];
00069 double result0 = region[iRow0];
00070 double value0 = element[j];
00071
00072 region[iRow0] = result0 - value0 * pivotValue;
00073 }
00074 regionIndex[numberNonZero++] = i;
00075 } else {
00076 region[i] = 0.0;
00077 }
00078 }
00079 regionSparse->setNumElements ( numberNonZero );
00080 }
00081
00082 void
00083 CoinFactorization::updateColumnLSparsish ( CoinIndexedVector * regionSparse,
00084 int * regionIndex)
00085 const
00086 {
00087 double *region = regionSparse->denseVector ( );
00088 int number = regionSparse->getNumElements ( );
00089 int numberNonZero;
00090 double tolerance = zeroTolerance_;
00091
00092 numberNonZero = 0;
00093 int k;
00094 int i , iPivot;
00095
00096 const CoinBigIndex *startColumn = startColumnL_;
00097 const int *indexRow = indexRowL_;
00098 const double *element = elementL_;
00099 int last = baseL_ + numberL_;
00100 #ifdef DENSE_CODE2
00101 if (numberDense_) {
00102
00103 last = min(last,numberRows_-numberDense_);
00104 }
00105 #endif
00106
00107
00108
00109 int * stack = sparse_;
00110 int * list = stack + maximumRowsExtra_;
00111 int * next = list + maximumRowsExtra_;
00112 CoinCheckZero * mark = (CoinCheckZero *) (next + maximumRowsExtra_);
00113 int nMarked=0;
00114 int smallestIndex = numberRowsExtra_;
00115
00116 for (k=0;k<number;k++) {
00117 iPivot=regionIndex[k];
00118 if (iPivot<baseL_) {
00119 regionIndex[numberNonZero++]=iPivot;
00120 } else {
00121 smallestIndex = min(iPivot,smallestIndex);
00122 int iWord = iPivot>>CHECK_SHIFT;
00123 int iBit = iPivot-(iWord<<CHECK_SHIFT);
00124 if (mark[iWord]) {
00125 mark[iWord] |= 1<<iBit;
00126 } else {
00127 mark[iWord] = 1<<iBit;
00128 stack[nMarked++]=iWord;
00129 }
00130 }
00131 }
00132
00133
00134 int jLast = (smallestIndex+BITS_PER_CHECK-1)>>CHECK_SHIFT;
00135 jLast = min((jLast<<CHECK_SHIFT),last);
00136 for ( i = smallestIndex; i < jLast; i++ ) {
00137 double pivotValue = region[i];
00138 CoinBigIndex start = startColumn[i];
00139 CoinBigIndex end = startColumn[i + 1];
00140
00141 if ( fabs(pivotValue) > tolerance ) {
00142 CoinBigIndex j;
00143 for ( j = start; j < end; j ++ ) {
00144 int iRow0 = indexRow[j];
00145 double result0 = region[iRow0];
00146 double value0 = element[j];
00147 region[iRow0] = result0 - value0 * pivotValue;
00148 int iWord = iRow0>>CHECK_SHIFT;
00149 int iBit = iRow0-(iWord<<CHECK_SHIFT);
00150 if (mark[iWord]) {
00151 mark[iWord] |= 1<<iBit;
00152 } else {
00153 mark[iWord] = 1<<iBit;
00154 stack[nMarked++]=iWord;
00155 }
00156 }
00157 regionIndex[numberNonZero++] = i;
00158 } else {
00159 region[i] = 0.0;
00160 }
00161 }
00162
00163 int kLast = last>>CHECK_SHIFT;
00164 if (jLast<last) {
00165
00166 for (k=(jLast>>CHECK_SHIFT);k<kLast;k++) {
00167 unsigned int iMark = mark[k];
00168 if (iMark) {
00169
00170 i = k<<CHECK_SHIFT;
00171 int iLast = i+BITS_PER_CHECK;
00172 for ( ; i < iLast; i++ ) {
00173 double pivotValue = region[i];
00174 CoinBigIndex start = startColumn[i];
00175 CoinBigIndex end = startColumn[i + 1];
00176
00177 if ( fabs(pivotValue) > tolerance ) {
00178 CoinBigIndex j;
00179 for ( j = start; j < end; j ++ ) {
00180 int iRow0 = indexRow[j];
00181 double result0 = region[iRow0];
00182 double value0 = element[j];
00183 region[iRow0] = result0 - value0 * pivotValue;
00184 int iWord = iRow0>>CHECK_SHIFT;
00185 int iBit = iRow0-(iWord<<CHECK_SHIFT);
00186 if (mark[iWord]) {
00187 mark[iWord] |= 1<<iBit;
00188 } else {
00189 mark[iWord] = 1<<iBit;
00190 stack[nMarked++]=iWord;
00191 }
00192 }
00193 regionIndex[numberNonZero++] = i;
00194 } else {
00195 region[i] = 0.0;
00196 }
00197 }
00198 mark[k]=0;
00199 }
00200 }
00201 i = kLast<<CHECK_SHIFT;
00202 }
00203 for ( ; i < last; i++ ) {
00204 double pivotValue = region[i];
00205 CoinBigIndex start = startColumn[i];
00206 CoinBigIndex end = startColumn[i + 1];
00207
00208 if ( fabs(pivotValue) > tolerance ) {
00209 CoinBigIndex j;
00210 for ( j = start; j < end; j ++ ) {
00211 int iRow0 = indexRow[j];
00212 double result0 = region[iRow0];
00213 double value0 = element[j];
00214 region[iRow0] = result0 - value0 * pivotValue;
00215 int iWord = iRow0>>CHECK_SHIFT;
00216 int iBit = iRow0-(iWord<<CHECK_SHIFT);
00217 if (mark[iWord]) {
00218 mark[iWord] |= 1<<iBit;
00219 } else {
00220 mark[iWord] = 1<<iBit;
00221 stack[nMarked++]=iWord;
00222 }
00223 }
00224 regionIndex[numberNonZero++] = i;
00225 } else {
00226 region[i] = 0.0;
00227 }
00228 }
00229
00230 mark[smallestIndex>>CHECK_SHIFT]=0;
00231 mark[kLast]=0;
00232 regionSparse->setNumElements ( numberNonZero );
00233 }
00234
00235 void
00236 CoinFactorization::updateColumnLSparse ( CoinIndexedVector * regionSparse ,
00237 int * regionIndex)
00238 const
00239 {
00240 double *region = regionSparse->denseVector ( );
00241 int number = regionSparse->getNumElements ( );
00242 int numberNonZero;
00243 double tolerance = zeroTolerance_;
00244
00245 numberNonZero = 0;
00246 int j, k;
00247 int i , iPivot;
00248
00249 const CoinBigIndex *startColumn = startColumnL_;
00250 const int *indexRow = indexRowL_;
00251 const double *element = elementL_;
00252
00253
00254 int * stack = sparse_;
00255 int * list = stack + maximumRowsExtra_;
00256 int * next = list + maximumRowsExtra_;
00257 char * mark = (char *) (next + maximumRowsExtra_);
00258 int nList;
00259 #ifdef COIN_DEBUG
00260 for (i=0;i<maximumRowsExtra_;i++) {
00261 assert (!mark[i]);
00262 }
00263 #endif
00264 int nStack;
00265 nList=0;
00266 for (k=0;k<number;k++) {
00267 iPivot=regionIndex[k];
00268 if (iPivot>=baseL_) {
00269 if(!mark[iPivot]) {
00270 stack[0]=iPivot;
00271 int j=startColumn[iPivot+1]-1;
00272 if (j>=startColumn[iPivot]) {
00273 int kPivot=indexRow[j];
00274
00275 next[0] =j-1;
00276
00277 if (!mark[kPivot]) {
00278 stack[1]=kPivot;
00279 mark[kPivot]=2;
00280 next[1]=startColumn[kPivot+1]-1;
00281 nStack=2;
00282 } else {
00283 nStack=1;
00284 }
00285 while (nStack) {
00286 int kPivot,j;
00287
00288 kPivot=stack[--nStack];
00289 j=next[nStack];
00290 if (j<startColumn[kPivot]) {
00291
00292 list[nList++]=kPivot;
00293 mark[kPivot]=1;
00294 } else {
00295 kPivot=indexRow[j];
00296
00297 next[nStack++] --;
00298 if (!mark[kPivot]) {
00299
00300 stack[nStack]=kPivot;
00301 mark[kPivot]=2;
00302 next[nStack++]=startColumn[kPivot+1]-1;
00303 }
00304 }
00305 }
00306 } else {
00307
00308 list[nList++]=iPivot;
00309 mark[iPivot]=1;
00310 }
00311 }
00312 } else {
00313
00314 regionIndex[numberNonZero++]=iPivot;
00315 }
00316 }
00317 for (i=nList-1;i>=0;i--) {
00318 iPivot = list[i];
00319 mark[iPivot]=0;
00320 double pivotValue = region[iPivot];
00321 if ( fabs ( pivotValue ) > tolerance ) {
00322 regionIndex[numberNonZero++]=iPivot;
00323 for ( j = startColumn[iPivot]; j < startColumn[iPivot+1]; j ++ ) {
00324 int iRow = indexRow[j];
00325 double value = element[j];
00326 region[iRow] -= value * pivotValue;
00327 }
00328 } else {
00329 region[iPivot]=0.0;
00330 }
00331 }
00332 regionSparse->setNumElements ( numberNonZero );
00333 }
00334
00335 void
00336 CoinFactorization::updateColumnL ( CoinIndexedVector * regionSparse,
00337 int * regionIndex) const
00338 {
00339 if (numberL_) {
00340 int number = regionSparse->getNumElements ( );
00341 int goSparse;
00342
00343 if (sparseThreshold_>0) {
00344 if (ftranAverageAfterL_) {
00345 int newNumber = (int) (number*ftranAverageAfterL_);
00346 if (newNumber< sparseThreshold_&&(numberL_<<2)>newNumber)
00347 goSparse = 2;
00348 else if (newNumber< sparseThreshold2_&&(numberL_<<1)>newNumber)
00349 goSparse = 1;
00350 else
00351 goSparse = 0;
00352 } else {
00353 if (number<sparseThreshold_&&(numberL_<<2)>number)
00354 goSparse = 2;
00355 else
00356 goSparse = 0;
00357 }
00358 } else {
00359 goSparse=0;
00360 }
00361 switch (goSparse) {
00362 case 0:
00363 updateColumnLDensish(regionSparse,regionIndex);
00364 break;
00365 case 1:
00366 updateColumnLSparsish(regionSparse,regionIndex);
00367 break;
00368 case 2:
00369 updateColumnLSparse(regionSparse,regionIndex);
00370 break;
00371 }
00372 }
00373 #ifdef DENSE_CODE
00374 if (numberDense_) {
00375
00376 int lastSparse = numberRows_-numberDense_;
00377 int number = regionSparse->getNumElements();
00378 double *region = regionSparse->denseVector ( );
00379 int i=0;
00380 bool doDense=false;
00381 while (i<number) {
00382 int iRow = regionIndex[i];
00383 if (iRow>=lastSparse) {
00384 doDense=true;
00385 regionIndex[i] = regionIndex[--number];
00386 } else {
00387 i++;
00388 }
00389 }
00390 if (doDense) {
00391
00392
00393
00394 char trans = 'N';
00395 int ione=1;
00396 int info;
00397 dgetrs_(&trans,&numberDense_,&ione,denseArea_,&numberDense_,
00398 densePermute_,region+lastSparse,&numberDense_,&info);
00399 for (i=lastSparse;i<numberRows_;i++) {
00400 double value = region[i];
00401 if (value) {
00402 if (fabs(value)>=1.0e-15)
00403 regionIndex[number++] = i;
00404 else
00405 region[i]=0.0;
00406 }
00407 }
00408 regionSparse->setNumElements(number);
00409 }
00410 }
00411 #endif
00412 }
00413
00414 int CoinFactorization::checkPivot(double saveFromU,
00415 double oldPivot) const
00416 {
00417 int status;
00418 if ( fabs ( saveFromU ) > 1.0e-8 ) {
00419 double checkTolerance;
00420
00421 if ( numberRowsExtra_ < numberRows_ + 2 ) {
00422 checkTolerance = 1.0e-5;
00423 } else if ( numberRowsExtra_ < numberRows_ + 10 ) {
00424 checkTolerance = 1.0e-6;
00425 } else if ( numberRowsExtra_ < numberRows_ + 50 ) {
00426 checkTolerance = 1.0e-8;
00427 } else {
00428 checkTolerance = 1.0e-10;
00429 }
00430 checkTolerance *= relaxCheck_;
00431 if ( fabs ( 1.0 - fabs ( saveFromU / oldPivot ) ) < checkTolerance ) {
00432 status = 0;
00433 } else {
00434 #if COIN_DEBUG
00435 std::cout <<"inaccurate pivot "<< oldPivot << " "
00436 << saveFromU << std::endl;
00437 #endif
00438 if ( fabs ( fabs ( oldPivot ) - fabs ( saveFromU ) ) < 1.0e-12 ||
00439 fabs ( 1.0 - fabs ( saveFromU / oldPivot ) ) < 1.0e-8 ) {
00440 status = 1;
00441 } else {
00442 status = 2;
00443 }
00444 }
00445 } else {
00446
00447 status = 2;
00448 #if COIN_DEBUG
00449 std::cout <<"inaccurate pivot "<< saveFromU / oldPivot
00450 << " " << saveFromU << std::endl;
00451 #endif
00452 }
00453 return status;
00454 }
00455
00456
00457
00458 int
00459 CoinFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
00460 int pivotRow,
00461 double pivotCheck ,
00462 bool checkBeforeModifying)
00463 {
00464 CoinBigIndex *startColumn;
00465 int *indexRow;
00466 double *element;
00467
00468
00469 if ( numberColumnsExtra_ >= maximumColumnsExtra_ ) {
00470 return 5;
00471 }
00472 if ( lengthAreaU_ < startColumnU_[maximumColumnsExtra_] ) {
00473 return 3;
00474 }
00475
00476 int realPivotRow;
00477 realPivotRow = pivotColumn_[pivotRow];
00478
00479 double *region = regionSparse->denseVector ( );
00480
00481 element = elementU_;
00482
00483
00484
00485 if (!numberPivots_)
00486 checkBeforeModifying=true;
00487
00488 totalElements_ -= numberInColumn_[realPivotRow];
00489 double oldPivot = pivotRegion_[realPivotRow];
00490
00491 pivotCheck = pivotCheck / oldPivot;
00492 #if COIN_DEBUG>1
00493 int checkNumber=1000000;
00494
00495 if (numberR_>=checkNumber) {
00496 printf("pivot row %d, check %g - alpha region:\n",
00497 realPivotRow,pivotCheck);
00498
00499
00500
00501
00502
00503 }
00504 #endif
00505 pivotRegion_[realPivotRow] = 0.0;
00506 CoinBigIndex i;
00507
00508 CoinBigIndex saveEnd = startColumnU_[realPivotRow]
00509 + numberInColumn_[realPivotRow];
00510
00511 numberInColumn_[realPivotRow] = 0;
00512
00513 CoinBigIndex *startRow = startRowU_;
00514 CoinBigIndex start;
00515 CoinBigIndex end;
00516
00517 start = startRow[realPivotRow];
00518 end = start + numberInRow_[realPivotRow];
00519 int numberNonZero = 0;
00520 int *indexColumn = indexColumnU_;
00521 CoinBigIndex *convertRowToColumn = convertRowToColumnU_;
00522 int *regionIndex = regionSparse->getIndices ( );
00523
00524 #if COIN_DEBUG>1
00525 if (numberR_>=checkNumber)
00526 printf("Before btranu\n");
00527 #endif
00528 int smallestIndex=numberRowsExtra_;
00529 if (!checkBeforeModifying) {
00530 for ( i = start; i < end ; i ++ ) {
00531 int iColumn = indexColumn[i];
00532 smallestIndex = min(smallestIndex,iColumn);
00533 CoinBigIndex j = convertRowToColumn[i];
00534
00535 region[iColumn] = element[j];
00536 #if COIN_DEBUG>1
00537 if (numberR_>=checkNumber)
00538 printf("%d %g\n",iColumn,region[iColumn]);
00539 #endif
00540 element[j] = 0.0;
00541 regionIndex[numberNonZero++] = iColumn;
00542 }
00543 } else {
00544 for ( i = start; i < end ; i ++ ) {
00545 int iColumn = indexColumn[i];
00546 smallestIndex = min(smallestIndex,iColumn);
00547 CoinBigIndex j = convertRowToColumn[i];
00548
00549 region[iColumn] = element[j];
00550 #if COIN_DEBUG>1
00551 if (numberR_>=checkNumber)
00552 printf("%d %g\n",iColumn,region[iColumn]);
00553 #endif
00554 regionIndex[numberNonZero++] = iColumn;
00555 }
00556 }
00557
00558 regionSparse->setNumElements ( numberNonZero );
00559 updateColumnTransposeU ( regionSparse, smallestIndex );
00560 numberNonZero = regionSparse->getNumElements ( );
00561
00562 double saveFromU = 0.0;
00563
00564 CoinBigIndex startU = startColumnU_[numberColumnsExtra_];
00565 int *indexU = &indexRowU_[startU];
00566 double *elementU = &elementU_[startU];
00567
00568
00569
00570 if (checkBeforeModifying) {
00571 double tolerance = zeroTolerance_;
00572 int number = numberInColumn_[numberColumnsExtra_];
00573
00574 for ( i = 0; i < number; i++ ) {
00575 int iRow = indexU[i];
00576
00577
00578
00579 if ( fabs ( elementU[i] ) > tolerance ) {
00580 if ( iRow != realPivotRow ) {
00581 saveFromU -= elementU[i] * region[iRow];
00582 } else {
00583 saveFromU += elementU[i];
00584 }
00585 }
00586 }
00587
00588 int status = checkPivot(saveFromU,pivotCheck);
00589 if (status) {
00590
00591 pivotRegion_[realPivotRow] = oldPivot;
00592 number = saveEnd-startColumnU_[realPivotRow];
00593 totalElements_ += number;
00594 numberInColumn_[realPivotRow]=number;
00595 regionSparse->clear();
00596 return status;
00597 } else {
00598
00599 for ( i = start; i < end ; i ++ ) {
00600 CoinBigIndex j = convertRowToColumn[i];
00601 element[j] = 0.0;
00602 }
00603 }
00604 }
00605
00606
00607 for ( i = startColumnU_[realPivotRow]; i < saveEnd ; i ++ ) {
00608 element[i] = 0.0;
00609 }
00610
00611
00612 startColumn = startColumnR_;
00613 indexRow = indexRowR_;
00614 element = elementR_;
00615 CoinBigIndex l = lengthR_;
00616 int number = numberR_;
00617
00618 startColumn[number] = l;
00619 number++;
00620 startColumn[number] = l + numberNonZero;
00621 numberR_ = number;
00622 lengthR_ = l + numberNonZero;
00623 totalElements_ += numberNonZero;
00624 if ( lengthR_ >= lengthAreaR_ ) {
00625
00626 regionSparse->clear();
00627 return 3;
00628 }
00629 #if COIN_DEBUG>1
00630 if (numberR_>=checkNumber)
00631 printf("After btranu\n");
00632 #endif
00633 for ( i = 0; i < numberNonZero; i++ ) {
00634 int iRow = regionIndex[i];
00635 #if COIN_DEBUG>1
00636 if (numberR_>=checkNumber)
00637 printf("%d %g\n",iRow,region[iRow]);
00638 #endif
00639
00640 indexRow[l] = iRow;
00641 element[l] = region[iRow];
00642 l++;
00643 }
00644
00645 int next = nextRow_[realPivotRow];
00646 int last = lastRow_[realPivotRow];
00647
00648 nextRow_[last] = next;
00649 lastRow_[next] = last;
00650 numberInRow_[realPivotRow]=0;
00651 #if COIN_DEBUG
00652 nextRow_[realPivotRow] = 777777;
00653 lastRow_[realPivotRow] = 777777;
00654 #endif
00655
00656 permute_[numberRowsExtra_] = realPivotRow;
00657 permuteBack_[numberRowsExtra_] = -1;
00658
00659 permute_[numberRowsExtra_ + 1] = 0;
00660
00661 pivotColumn_[pivotRow] = numberRowsExtra_;
00662 pivotColumnBack_[numberRowsExtra_] = pivotRow;
00663 startColumn = startColumnU_;
00664 indexRow = indexRowU_;
00665 element = elementU_;
00666
00667 numberU_++;
00668 number = numberInColumn_[numberColumnsExtra_];
00669
00670 totalElements_ += number;
00671 lengthU_ += number;
00672 if ( lengthU_ >= lengthAreaU_ ) {
00673
00674 regionSparse->clear();
00675 return 3;
00676 }
00677
00678 saveFromU = 0.0;
00679
00680
00681
00682
00683 double tolerance = zeroTolerance_;
00684
00685 #if COIN_DEBUG>1
00686 if (numberR_>=checkNumber)
00687 printf("On U\n");
00688 #endif
00689 for ( i = 0; i < number; i++ ) {
00690 int iRow = indexU[i];
00691 #if COIN_DEBUG>1
00692 if (numberR_>=checkNumber)
00693 printf("%d %g\n",iRow,elementU[i]);
00694 #endif
00695
00696 if ( fabs ( elementU[i] ) > tolerance ) {
00697 if ( iRow != realPivotRow ) {
00698 int next = nextRow_[iRow];
00699 int numberInRow = numberInRow_[iRow];
00700 CoinBigIndex space;
00701 CoinBigIndex put = startRow[iRow] + numberInRow;
00702
00703 space = startRow[next] - put;
00704 if ( space <= 0 ) {
00705 getRowSpaceIterate ( iRow, numberInRow + 4 );
00706 put = startRow[iRow] + numberInRow;
00707 }
00708 indexColumn[put] = numberColumnsExtra_;
00709 convertRowToColumn[put] = i + startU;
00710 numberInRow_[iRow] = numberInRow + 1;
00711 saveFromU = saveFromU - elementU[i] * region[iRow];
00712 } else {
00713
00714 saveFromU += elementU[i];
00715 elementU[i] = 0.0;
00716 }
00717 } else {
00718 elementU[i] = 0.0;
00719 }
00720 }
00721
00722 last = lastRow_[maximumRowsExtra_];
00723 nextRow_[last] = numberRowsExtra_;
00724 lastRow_[maximumRowsExtra_] = numberRowsExtra_;
00725 lastRow_[numberRowsExtra_] = last;
00726 nextRow_[numberRowsExtra_] = maximumRowsExtra_;
00727 startRow[numberRowsExtra_] = startRow[maximumRowsExtra_];
00728 numberInRow_[numberRowsExtra_] = 0;
00729
00730 next = nextColumn_[maximumColumnsExtra_];
00731 lastColumn_[next] = numberColumnsExtra_;
00732 nextColumn_[maximumColumnsExtra_] = numberColumnsExtra_;
00733 nextColumn_[numberColumnsExtra_] = next;
00734 lastColumn_[numberColumnsExtra_] = maximumColumnsExtra_;
00735
00736 int status = checkPivot(saveFromU,pivotCheck);
00737
00738 if (status!=2) {
00739
00740 double pivotValue = 1.0 / saveFromU;
00741
00742 pivotRegion_[numberRowsExtra_] = pivotValue;
00743
00744 for ( i = 0; i < number; i++ ) {
00745 elementU[i] *= pivotValue;
00746 }
00747 numberRowsExtra_++;
00748 numberColumnsExtra_++;
00749 numberGoodU_++;
00750 numberPivots_++;
00751 }
00752 if ( numberRowsExtra_ > numberRows_ + 50 ) {
00753 CoinBigIndex extra = factorElements_ >> 1;
00754
00755 if ( numberRowsExtra_ > numberRows_ + 100 + numberRows_ / 500 ) {
00756 if ( extra < 2 * numberRows_ ) {
00757 extra = 2 * numberRows_;
00758 }
00759 } else {
00760 if ( extra < 5 * numberRows_ ) {
00761 extra = 5 * numberRows_;
00762 }
00763 }
00764 CoinBigIndex added = totalElements_ - factorElements_;
00765
00766 if ( added > extra && added > ( factorElements_ ) << 1 && !status
00767 && 3*totalElements_ > 2*(lengthAreaU_+lengthAreaL_)) {
00768 status = 3;
00769 if ( messageLevel_ & 4 ) {
00770 std::cout << "Factorization has "<< totalElements_
00771 << ", basis had " << factorElements_ <<std::endl;
00772 }
00773 }
00774 }
00775 if (numberInColumnPlus_&&status<2) {
00776
00777 double * elementR = elementR_ + lengthAreaR_;
00778 int * indexRowR = indexRowR_ + lengthAreaR_;
00779 int * startR = startColumnR_+maximumPivots_+1;
00780 int pivotRow = numberRowsExtra_-1;
00781 for ( i = 0; i < numberNonZero; i++ ) {
00782 int iRow = regionIndex[i];
00783 assert (pivotRow>iRow);
00784 next = nextColumn_[iRow];
00785 int space;
00786 if (next!=maximumColumnsExtra_)
00787 space = startR[next]-startR[iRow];
00788 else
00789 space = lengthAreaR_-startR[iRow];
00790 int numberInR = numberInColumnPlus_[iRow];
00791 if (space>numberInR) {
00792
00793 int put=startR[iRow]+numberInR;
00794 numberInColumnPlus_[iRow]=numberInR+1;
00795 indexRowR[put]=pivotRow;
00796 elementR[put]=region[iRow];
00797
00798 if (next==maximumColumnsExtra_)
00799 startR[maximumColumnsExtra_] = min(put + 4,lengthAreaR_);
00800 } else {
00801
00802 if (!getColumnSpaceIterateR(iRow,region[iRow],pivotRow)) {
00803
00804 delete [] numberInColumnPlus_;
00805 numberInColumnPlus_=NULL;
00806 regionSparse->clear();
00807 break;
00808 }
00809 }
00810 region[iRow]=0.0;
00811 }
00812 regionSparse->setNumElements(0);
00813 } else {
00814 regionSparse->clear();
00815 }
00816 return status;
00817 }
00818
00819
00820 int
00821 CoinFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
00822 CoinIndexedVector * regionSparse2 )
00823 const
00824 {
00825
00826 regionSparse->clear ( );
00827 double *region = regionSparse->denseVector ( );
00828 double * vector = regionSparse2->denseVector();
00829 int * index = regionSparse2->getIndices();
00830 int number = regionSparse2->getNumElements();
00831 int i;
00832 const int * pivotColumn = pivotColumn_;
00833
00834
00835 int *regionIndex = regionSparse->getIndices ( );
00836 int numberNonZero = number;
00837 int iRow;
00838 bool packed = regionSparse2->packedMode();
00839 if (packed) {
00840 for ( i = 0; i < number; i ++ ) {
00841 iRow = index[i];
00842 double value = vector[i];
00843 vector[i]=0.0;
00844 iRow=pivotColumn[iRow];
00845 region[iRow] = value;
00846 regionIndex[i] = iRow;
00847 }
00848 } else {
00849 for ( i = 0; i < number; i ++ ) {
00850 iRow = index[i];
00851 double value = vector[iRow];
00852 vector[iRow]=0.0;
00853 iRow=pivotColumn[iRow];
00854 region[iRow] = value;
00855 regionIndex[i] = iRow;
00856 }
00857 }
00858 regionSparse->setNumElements ( numberNonZero );
00859
00860
00861 int j;
00862 double *pivotRegion = pivotRegion_;
00863
00864 if (collectStatistics_) {
00865 numberBtranCounts_++;
00866 btranCountInput_ += (double) numberNonZero;
00867 }
00868
00869 int smallestIndex=numberRowsExtra_;
00870 for ( j = 0; j < numberNonZero; j++ ) {
00871 int iRow = regionIndex[j];
00872 smallestIndex = min (smallestIndex,iRow);
00873 region[iRow] *= pivotRegion[iRow];
00874 }
00875 updateColumnTransposeU ( regionSparse,smallestIndex );
00876 if (collectStatistics_)
00877 btranCountAfterU_ += (double) regionSparse->getNumElements();
00878
00879
00880
00881 double countBefore = btranCountAfterR_;
00882 updateColumnTransposeR ( regionSparse );
00883
00884
00885 numberNonZero = regionSparse->getNumElements ( );
00886
00887 updateColumnTransposeL ( regionSparse );
00888
00889 if (collectStatistics_)
00890 if (numberNonZero<=numberRows_)
00891 btranCountAfterL_ += (double) regionSparse->getNumElements();
00892 else
00893 btranCountAfterL_ += btranCountAfterR_-countBefore;
00894 number = regionSparse->getNumElements ( );
00895 const int * permuteBack = permuteBack_;
00896 if (packed) {
00897 for (i=0;i<number;i++) {
00898 int iRow=regionIndex[i];
00899 double value = region[iRow];
00900 region[iRow]=0.0;
00901 iRow=permuteBack[iRow];
00902 vector[i]=value;
00903 index[i]=iRow;
00904 }
00905 } else {
00906 for (i=0;i<number;i++) {
00907 int iRow=regionIndex[i];
00908 double value = region[iRow];
00909 region[iRow]=0.0;
00910 iRow=permuteBack[iRow];
00911 vector[iRow]=value;
00912 index[i]=iRow;
00913 }
00914 }
00915 regionSparse->setNumElements(0);
00916 regionSparse2->setNumElements(number);
00917 #ifdef COIN_DEBUG
00918 for (i=0;i<numberRowsExtra_;i++) {
00919 assert (!region[i]);
00920 }
00921 #endif
00922 return number;
00923 }
00924
00925
00926
00927 void
00928 CoinFactorization::updateColumnTransposeUDensish
00929 ( CoinIndexedVector * regionSparse,
00930 int smallestIndex) const
00931 {
00932 double *region = regionSparse->denseVector ( );
00933 int numberNonZero = regionSparse->getNumElements ( );
00934 double tolerance = zeroTolerance_;
00935
00936 int *regionIndex = regionSparse->getIndices ( );
00937
00938 int i,j;
00939
00940 const CoinBigIndex *startRow = startRowU_;
00941
00942 const CoinBigIndex *convertRowToColumn = convertRowToColumnU_;
00943 const int *indexColumn = indexColumnU_;
00944
00945 const double * element = elementU_;
00946 int last = numberU_;
00947
00948 double pivotValue;
00949
00950 const int *numberInRow = numberInRow_;
00951
00952 numberNonZero = 0;
00953 for (i=smallestIndex ; i < last; i++ ) {
00954 pivotValue = region[i];
00955 if ( fabs ( pivotValue ) > tolerance ) {
00956 CoinBigIndex start = startRow[i];
00957 int numberIn = numberInRow[i];
00958 CoinBigIndex end = start + numberIn;
00959 for (j = start ; j < end; j ++ ) {
00960 int iRow = indexColumn[j];
00961 CoinBigIndex getElement = convertRowToColumn[j];
00962 double value = element[getElement];
00963
00964 region[iRow] -= value * pivotValue;
00965 }
00966 regionIndex[numberNonZero++] = i;
00967 } else {
00968 region[i] = 0.0;
00969 }
00970 }
00971
00972 regionSparse->setNumElements ( numberNonZero );
00973 }
00974
00975
00976 void
00977 CoinFactorization::updateColumnTransposeUSparsish
00978 ( CoinIndexedVector * regionSparse,
00979 int smallestIndex) const
00980 {
00981 double *region = regionSparse->denseVector ( );
00982 int numberNonZero = regionSparse->getNumElements ( );
00983 double tolerance = zeroTolerance_;
00984
00985 int *regionIndex = regionSparse->getIndices ( );
00986
00987 int i,j;
00988
00989 const CoinBigIndex *startRow = startRowU_;
00990
00991 const CoinBigIndex *convertRowToColumn = convertRowToColumnU_;
00992 const int *indexColumn = indexColumnU_;
00993
00994 const double * element = elementU_;
00995 int last = numberU_;
00996
00997 double pivotValue;
00998
00999 const int *numberInRow = numberInRow_;
01000
01001
01002
01003 int * stack = sparse_;
01004 int * list = stack + maximumRowsExtra_;
01005 int * next = list + maximumRowsExtra_;
01006 CoinCheckZero * mark = (CoinCheckZero *) (next + maximumRowsExtra_);
01007 int nMarked=0;
01008
01009 for (i=0;i<numberNonZero;i++) {
01010 int iPivot=regionIndex[i];
01011 int iWord = iPivot>>CHECK_SHIFT;
01012 int iBit = iPivot-(iWord<<CHECK_SHIFT);
01013 if (mark[iWord]) {
01014 mark[iWord] |= 1<<iBit;
01015 } else {
01016 mark[iWord] = 1<<iBit;
01017 stack[nMarked++]=iWord;
01018 }
01019 }
01020
01021 numberNonZero = 0;
01022
01023 smallestIndex = smallestIndex >> CHECK_SHIFT;
01024 int kLast = last>>CHECK_SHIFT;
01025
01026 int k;
01027
01028 for (k=smallestIndex;k<kLast;k++) {
01029 unsigned int iMark = mark[k];
01030 if (iMark) {
01031
01032 i = k<<CHECK_SHIFT;
01033 int iLast = i+BITS_PER_CHECK;
01034 for ( ; i < iLast; i++ ) {
01035 pivotValue = region[i];
01036 if ( fabs ( pivotValue ) > tolerance ) {
01037 CoinBigIndex start = startRow[i];
01038 int numberIn = numberInRow[i];
01039 CoinBigIndex end = start + numberIn;
01040 for (j = start ; j < end; j ++ ) {
01041 int iRow = indexColumn[j];
01042 CoinBigIndex getElement = convertRowToColumn[j];
01043 double value = element[getElement];
01044 int iWord = iRow>>CHECK_SHIFT;
01045 int iBit = iRow-(iWord<<CHECK_SHIFT);
01046 if (mark[iWord]) {
01047 mark[iWord] |= 1<<iBit;
01048 } else {
01049 mark[iWord] = 1<<iBit;
01050 stack[nMarked++]=iWord;
01051 }
01052 region[iRow] -= value * pivotValue;
01053 }
01054 regionIndex[numberNonZero++] = i;
01055 } else {
01056 region[i] = 0.0;
01057 }
01058 }
01059 mark[k]=0;
01060 }
01061 }
01062 i = kLast<<CHECK_SHIFT;
01063 mark[kLast]=0;
01064 for (; i < last; i++ ) {
01065 pivotValue = region[i];
01066 if ( fabs ( pivotValue ) > tolerance ) {
01067 CoinBigIndex start = startRow[i];
01068 int numberIn = numberInRow[i];
01069 CoinBigIndex end = start + numberIn;
01070 for (j = start ; j < end; j ++ ) {
01071 int iRow = indexColumn[j];
01072 CoinBigIndex getElement = convertRowToColumn[j];
01073 double value = element[getElement];
01074
01075 region[iRow] -= value * pivotValue;
01076 }
01077 regionIndex[numberNonZero++] = i;
01078 } else {
01079 region[i] = 0.0;
01080 }
01081 }
01082 #ifdef COIN_DEBUG
01083 for (i=0;i<maximumRowsExtra_;i++) {
01084 assert (!mark[i]);
01085 }
01086 #endif
01087
01088 regionSparse->setNumElements ( numberNonZero );
01089 }
01090
01091
01092 void
01093 CoinFactorization::updateColumnTransposeUSparse (
01094 CoinIndexedVector * regionSparse) const
01095 {
01096 double *region = regionSparse->denseVector ( );
01097 int numberNonZero = regionSparse->getNumElements ( );
01098 double tolerance = zeroTolerance_;
01099
01100 int *regionIndex = regionSparse->getIndices ( );
01101
01102 int i;
01103
01104 const CoinBigIndex *startRow = startRowU_;
01105
01106 const CoinBigIndex *convertRowToColumn = convertRowToColumnU_;
01107 const int *indexColumn = indexColumnU_;
01108
01109 const double * element = elementU_;
01110
01111 double pivotValue;
01112
01113 int *numberInRow = numberInRow_;
01114
01115
01116
01117 int * stack = sparse_;
01118 int * list = stack + maximumRowsExtra_;
01119 int * next = list + maximumRowsExtra_;
01120 char * mark = (char *) (next + maximumRowsExtra_);
01121 int nList;
01122 int iPivot;
01123 #ifdef COIN_DEBUG
01124 for (i=0;i<maximumRowsExtra_;i++) {
01125 assert (!mark[i]);
01126 }
01127 #endif
01128 #if 0
01129 {
01130 int i;
01131 for (i=0;i<numberRowsExtra_;i++) {
01132 int krs = startRow[i];
01133 int kre = krs + numberInRow[i];
01134 int k;
01135 for (k=krs;k<kre;k++)
01136 assert (indexColumn[k]>i);
01137 }
01138 }
01139 #endif
01140 int k,nStack;
01141 nList=0;
01142 for (k=0;k<numberNonZero;k++) {
01143 iPivot=regionIndex[k];
01144 if(!mark[iPivot]) {
01145 stack[0]=iPivot;
01146 mark[iPivot]=1;
01147 int j=startRow[iPivot]+numberInRow[iPivot];
01148 if (j>startRow[iPivot]) {
01149 int kPivot=indexColumn[--j];
01150
01151 next[0] =j;
01152
01153 if (!mark[kPivot]) {
01154 stack[1]=kPivot;
01155 mark[kPivot]=1;
01156 next[1]=startRow[kPivot]+numberInRow[kPivot];
01157 nStack=1;
01158 } else {
01159 nStack=0;
01160 }
01161 while (nStack>=0) {
01162 int kPivot,j;
01163
01164 kPivot=stack[nStack];
01165 j=next[nStack];
01166
01167 #if 0
01168 if (j==startRow[kPivot]) {
01169
01170 --nStack;
01171 list[nList++]=kPivot;
01172 } else {
01173 kPivot=indexColumn[--j];
01174
01175 next[nStack] =j;
01176 if (!mark[kPivot]) {
01177
01178 stack[++nStack]=kPivot;
01179 mark[kPivot]=1;
01180 next[nStack]=startRow[kPivot]+numberInRow[kPivot];
01181 }
01182 }
01183 #else
01184 bool finished=true;
01185 while (j!=startRow[kPivot]) {
01186 int jPivot=indexColumn[--j];
01187 if (!mark[jPivot]) {
01188
01189 next[nStack] =j;
01190
01191 stack[++nStack]=jPivot;
01192 mark[jPivot]=1;
01193 next[nStack]=startRow[jPivot]+numberInRow[jPivot];
01194 finished=false;
01195 break;
01196 }
01197 }
01198 if (finished) {
01199 --nStack;
01200 list[nList++]=kPivot;
01201 }
01202 #endif
01203 }
01204 } else {
01205
01206 list[nList++]=iPivot;
01207 mark[iPivot]=1;
01208 }
01209 }
01210 }
01211 numberNonZero=0;
01212 for (i=nList-1;i>=0;i--) {
01213 iPivot = list[i];
01214 mark[iPivot]=0;
01215 pivotValue = region[iPivot];
01216 if ( fabs ( pivotValue ) > tolerance ) {
01217 CoinBigIndex start = startRow[iPivot];
01218 int numberIn = numberInRow[iPivot];
01219 CoinBigIndex end = start + numberIn;
01220 CoinBigIndex j;
01221 for (j=start ; j < end; j ++ ) {
01222 int iRow = indexColumn[j];
01223 CoinBigIndex getElement = convertRowToColumn[j];
01224 double value = element[getElement];
01225
01226 region[iRow] = region[iRow]
01227 - value * pivotValue;
01228 }
01229 regionIndex[numberNonZero++] = iPivot;
01230 } else {
01231 region[iPivot] = 0.0;
01232 }
01233 }
01234
01235 regionSparse->setNumElements ( numberNonZero );
01236 }
01237
01238
01239
01240 void
01241 CoinFactorization::updateColumnTransposeU ( CoinIndexedVector * regionSparse,
01242 int smallestIndex) const
01243 {
01244 int number = regionSparse->getNumElements ( );
01245 int goSparse;
01246
01247 if (sparseThreshold_>0) {
01248 if (btranAverageAfterU_) {
01249 int newNumber = (int) (number*btranAverageAfterU_);
01250 if (newNumber< sparseThreshold_)
01251 goSparse = 2;
01252 else if (newNumber< sparseThreshold2_)
01253 goSparse = 1;
01254 else
01255 goSparse = 0;
01256 } else {
01257 if (number<sparseThreshold_)
01258 goSparse = 2;
01259 else
01260 goSparse = 0;
01261 }
01262 } else {
01263 goSparse=0;
01264 }
01265 switch (goSparse) {
01266 case 0:
01267 updateColumnTransposeUDensish(regionSparse,smallestIndex);
01268 break;
01269 case 1:
01270 updateColumnTransposeUSparsish(regionSparse,smallestIndex);
01271 break;
01272 case 2:
01273 updateColumnTransposeUSparse(regionSparse);
01274 break;
01275 }
01276 }
01277
01278
01279
01280 void
01281 CoinFactorization::updateColumnTransposeLDensish
01282 ( CoinIndexedVector * regionSparse ) const
01283 {
01284 double *region = regionSparse->denseVector ( );
01285 int *regionIndex = regionSparse->getIndices ( );
01286 int numberNonZero;
01287 double tolerance = zeroTolerance_;
01288 int base;
01289 int first = -1;
01290
01291 numberNonZero=0;
01292
01293 for (first=numberRows_-1;first>=0;first--) {
01294 if (region[first])
01295 break;
01296 }
01297 if ( first >= 0 ) {
01298 base = baseL_;
01299 const CoinBigIndex *startColumn = startColumnL_;
01300 const int *indexRow = indexRowL_;
01301 const double *element = elementL_;
01302 int last = baseL_ + numberL_;
01303
01304 if ( first >= last ) {
01305 first = last - 1;
01306 }
01307 int i;
01308 double pivotValue;
01309 for (i = first ; i >= base; i-- ) {
01310 CoinBigIndex j;
01311 pivotValue = region[i];
01312 for ( j= startColumn[i] ; j < startColumn[i+1]; j++ ) {
01313 int iRow = indexRow[j];
01314 double value = element[j];
01315 pivotValue -= value * region[iRow];
01316 }
01317 if ( fabs ( pivotValue ) > tolerance ) {
01318 region[i] = pivotValue;
01319 regionIndex[numberNonZero++] = i;
01320 } else {
01321 region[i] = 0.0;
01322 }
01323 }
01324
01325 if ( first < base ) {
01326 base = first + 1;
01327 }
01328 for (i = base -1 ; i >= 0; i-- ) {
01329 pivotValue = region[i];
01330 if ( fabs ( pivotValue ) > tolerance ) {
01331 region[i] = pivotValue;
01332 regionIndex[numberNonZero++] = i;
01333 } else {
01334 region[i] = 0.0;
01335 }
01336 }
01337 }
01338
01339 regionSparse->setNumElements ( numberNonZero );
01340 }
01341
01342
01343 void
01344 CoinFactorization::updateColumnTransposeLByRow
01345 ( CoinIndexedVector * regionSparse ) const
01346 {
01347 double *region = regionSparse->denseVector ( );
01348 int *regionIndex = regionSparse->getIndices ( );
01349 int numberNonZero;
01350 double tolerance = zeroTolerance_;
01351 int first = -1;
01352
01353
01354 const double * element = elementByRowL_;
01355 const CoinBigIndex * startRow = startRowL_;
01356 const int * column = indexColumnL_;
01357 int i;
01358 CoinBigIndex j;
01359 for (first=numberRows_-1;first>=0;first--) {
01360 if (region[first])
01361 break;
01362 }
01363 numberNonZero=0;
01364 for (i=first;i>=0;i--) {
01365 double pivotValue = region[i];
01366 if ( fabs ( pivotValue ) > tolerance ) {
01367 regionIndex[numberNonZero++] = i;
01368 for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
01369 int iRow = column[j];
01370 double value = element[j];
01371 region[iRow] -= pivotValue*value;
01372 }
01373 } else {
01374 region[i] = 0.0;
01375 }
01376 }
01377
01378 regionSparse->setNumElements ( numberNonZero );
01379 }
01380
01381 void
01382 CoinFactorization::updateColumnTransposeLSparsish
01383 ( CoinIndexedVector * regionSparse ) const
01384 {
01385 double *region = regionSparse->denseVector ( );
01386 int *regionIndex = regionSparse->getIndices ( );
01387 int numberNonZero = regionSparse->getNumElements();
01388 double tolerance = zeroTolerance_;
01389
01390
01391 const double * element = elementByRowL_;
01392 const CoinBigIndex * startRow = startRowL_;
01393 const int * column = indexColumnL_;
01394 int i;
01395 CoinBigIndex j;
01396
01397
01398 int * stack = sparse_;
01399 int * list = stack + maximumRowsExtra_;
01400 int * next = list + maximumRowsExtra_;
01401 CoinCheckZero * mark = (CoinCheckZero *) (next + maximumRowsExtra_);
01402 int nMarked=0;
01403 #if 1
01404 for (i=0;i<numberNonZero;i++) {
01405 int iPivot=regionIndex[i];
01406 int iWord = iPivot>>CHECK_SHIFT;
01407 int iBit = iPivot-(iWord<<CHECK_SHIFT);
01408 if (mark[iWord]) {
01409 mark[iWord] |= 1<<iBit;
01410 } else {
01411 mark[iWord] = 1<<iBit;
01412 stack[nMarked++]=iWord;
01413 }
01414 }
01415 numberNonZero = 0;
01416
01417 int jLast = (numberRows_-1)>>CHECK_SHIFT;
01418 jLast = (jLast<<CHECK_SHIFT);
01419 for (i=numberRows_-1;i>=jLast;i--) {
01420 double pivotValue = region[i];
01421 if ( fabs ( pivotValue ) > tolerance ) {
01422 regionIndex[numberNonZero++] = i;
01423 for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
01424 int iRow = column[j];
01425 double value = element[j];
01426 int iWord = iRow>>CHECK_SHIFT;
01427 int iBit = iRow-(iWord<<CHECK_SHIFT);
01428 if (mark[iWord]) {
01429 mark[iWord] |= 1<<iBit;
01430 } else {
01431 mark[iWord] = 1<<iBit;
01432 stack[nMarked++]=iWord;
01433 }
01434 region[iRow] -= pivotValue*value;
01435 }
01436 } else {
01437 region[i] = 0.0;
01438 }
01439 }
01440
01441 jLast = jLast>>CHECK_SHIFT;
01442 int k ;
01443 for (k=jLast-1;k>=0;k--) {
01444 unsigned int iMark = mark[k];
01445 if (iMark) {
01446
01447 int iLast = k<<CHECK_SHIFT;
01448 i = iLast+BITS_PER_CHECK-1;
01449 for ( ; i >= iLast; i-- ) {
01450 double pivotValue = region[i];
01451 if ( fabs ( pivotValue ) > tolerance ) {
01452 regionIndex[numberNonZero++] = i;
01453 for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
01454 int iRow = column[j];
01455 double value = element[j];
01456 int iWord = iRow>>CHECK_SHIFT;
01457 int iBit = iRow-(iWord<<CHECK_SHIFT);
01458 if (mark[iWord]) {
01459 mark[iWord] |= 1<<iBit;
01460 } else {
01461 mark[iWord] = 1<<iBit;
01462 stack[nMarked++]=iWord;
01463 }
01464 region[iRow] -= pivotValue*value;
01465 }
01466 } else {
01467 region[i] = 0.0;
01468 }
01469 }
01470 mark[k]=0;
01471 }
01472 }
01473 mark[jLast]=0;
01474 #ifdef COIN_DEBUG
01475 for (i=0;i<maximumRowsExtra_;i++) {
01476 assert (!mark[i]);
01477 }
01478 #endif
01479 #else
01480 for (first=numberRows_-1;first>=0;first--) {
01481 if (region[first])
01482 break;
01483 }
01484 numberNonZero=0;
01485 for (i=first;i>=0;i--) {
01486 double pivotValue = region[i];
01487 if ( fabs ( pivotValue ) > tolerance ) {
01488 regionIndex[numberNonZero++] = i;
01489 for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
01490 int iRow = column[j];
01491 double value = element[j];
01492 region[iRow] -= pivotValue*value;
01493 }
01494 } else {
01495 region[i] = 0.0;
01496 }
01497 }
01498 #endif
01499
01500 regionSparse->setNumElements ( numberNonZero );
01501 }
01502
01503
01504 void
01505 CoinFactorization::updateColumnTransposeLSparse
01506 ( CoinIndexedVector * regionSparse ) const
01507 {
01508 double *region = regionSparse->denseVector ( );
01509 int *regionIndex = regionSparse->getIndices ( );
01510 int numberNonZero = regionSparse->getNumElements ( );
01511 double tolerance = zeroTolerance_;
01512
01513
01514 const double * element = elementByRowL_;
01515 const CoinBigIndex * startRow = startRowL_;
01516 const int * column = indexColumnL_;
01517 int i;
01518 CoinBigIndex j;
01519
01520
01521 int * stack = sparse_;
01522 int * list = stack + maximumRowsExtra_;
01523 int * next = list + maximumRowsExtra_;
01524 char * mark = (char *) (next + maximumRowsExtra_);
01525 int nList;
01526 int number = numberNonZero;
01527 int k, iPivot;
01528 #ifdef COIN_DEBUG
01529 for (i=0;i<maximumRowsExtra_;i++) {
01530 assert (!mark[i]);
01531 }
01532 #endif
01533 int nStack;
01534 nList=0;
01535 for (k=0;k<number;k++) {
01536 iPivot=regionIndex[k];
01537 if(!mark[iPivot]) {
01538 stack[0]=iPivot;
01539 int j=startRow[iPivot+1]-1;
01540 if (j>=startRow[iPivot]) {
01541 int kPivot=column[j];
01542
01543 next[0] =j-1;
01544
01545 if (!mark[kPivot]) {
01546 stack[1]=kPivot;
01547 mark[kPivot]=2;
01548 next[1]=startRow[kPivot+1]-1;
01549 nStack=2;
01550 } else {
01551 nStack=1;
01552 }
01553 while (nStack) {
01554 int kPivot,j;
01555
01556 kPivot=stack[--nStack];
01557 j=next[nStack];
01558 if (j<startRow[kPivot]) {
01559
01560 list[nList++]=kPivot;
01561 mark[kPivot]=1;
01562 } else {
01563 kPivot=column[j];
01564
01565 next[nStack++] --;
01566 if (!mark[kPivot]) {
01567
01568 stack[nStack]=kPivot;
01569 mark[kPivot]=2;
01570 next[nStack++]=startRow[kPivot+1]-1;
01571 }
01572 }
01573 }
01574 } else {
01575
01576 list[nList++]=iPivot;
01577 mark[iPivot]=1;
01578 }
01579 }
01580 }
01581 numberNonZero=0;
01582 for (i=nList-1;i>=0;i--) {
01583 iPivot = list[i];
01584 mark[iPivot]=0;
01585 double pivotValue = region[iPivot];
01586 if ( fabs ( pivotValue ) > tolerance ) {
01587 regionIndex[numberNonZero++] = iPivot;
01588 for ( j = startRow[iPivot]; j < startRow[iPivot+1]; j ++ ) {
01589 int iRow = column[j];
01590 double value = element[j];
01591 region[iRow] -= value * pivotValue;
01592 }
01593 } else {
01594 region[iPivot]=0.0;
01595 }
01596 }
01597
01598 regionSparse->setNumElements ( numberNonZero );
01599 }
01600
01601 void
01602 CoinFactorization::updateColumnTransposeL ( CoinIndexedVector * regionSparse ) const
01603 {
01604 int number = regionSparse->getNumElements ( );
01605 if (!numberL_&&!numberDense_) {
01606 if (sparse_||number<numberRows_)
01607 return;
01608 }
01609 int goSparse;
01610
01611
01612 if (sparseThreshold_>0) {
01613 if (btranAverageAfterL_) {
01614 int newNumber = (int) (number*btranAverageAfterL_);
01615 if (newNumber< sparseThreshold_)
01616 goSparse = 2;
01617 else if (newNumber< sparseThreshold2_)
01618 goSparse = 1;
01619 else
01620 goSparse = 0;
01621 } else {
01622 if (number<sparseThreshold_)
01623 goSparse = 2;
01624 else
01625 goSparse = 0;
01626 }
01627 } else {
01628 goSparse=-1;
01629 }
01630 #ifdef DENSE_CODE
01631 if (numberDense_) {
01632
01633 int lastSparse = numberRows_-numberDense_;
01634 int number = regionSparse->getNumElements();
01635 double *region = regionSparse->denseVector ( );
01636 int *regionIndex = regionSparse->getIndices ( );
01637 int i=0;
01638 bool doDense=false;
01639 if (number<=numberRows_) {
01640 while (i<number) {
01641 int iRow = regionIndex[i];
01642 if (iRow>=lastSparse) {
01643 doDense=true;
01644 regionIndex[i] = regionIndex[--number];
01645 } else {
01646 i++;
01647 }
01648 }
01649 } else {
01650 for (i=numberRows_-1;i>=lastSparse;i--) {
01651 if (region[i]) {
01652 doDense=true;
01653 break;
01654 }
01655 }
01656 if (sparseThreshold_)
01657 goSparse=0;
01658 else
01659 goSparse=-1;
01660 }
01661 if (doDense) {
01662 regionSparse->setNumElements(number);
01663 char trans = 'T';
01664 int ione=1;
01665 int info;
01666 dgetrs_(&trans,&numberDense_,&ione,denseArea_,&numberDense_,
01667 densePermute_,region+lastSparse,&numberDense_,&info);
01668
01669 if (goSparse>0||!numberL_)
01670 regionSparse->scan(lastSparse,numberRows_,zeroTolerance_);
01671 }
01672 if (!numberL_)
01673 return;
01674 }
01675 #endif
01676 switch (goSparse) {
01677 case -1:
01678 updateColumnTransposeLDensish(regionSparse);
01679 break;
01680 case 0:
01681 updateColumnTransposeLByRow(regionSparse);
01682 break;
01683 case 1:
01684 updateColumnTransposeLSparsish(regionSparse);
01685 break;
01686 case 2:
01687 updateColumnTransposeLSparse(regionSparse);
01688 break;
01689 }
01690 }
01691
01692
01693
01694
01695 bool
01696 CoinFactorization::getRowSpaceIterate ( int iRow,
01697 int extraNeeded )
01698 {
01699 int number = numberInRow_[iRow];
01700 CoinBigIndex *startRow = startRowU_;
01701 int *indexColumn = indexColumnU_;
01702 CoinBigIndex *convertRowToColumn = convertRowToColumnU_;
01703 CoinBigIndex space = lengthAreaU_ - startRow[maximumRowsExtra_];
01704 if ( space < extraNeeded + number + 2 ) {
01705
01706 int iRow = nextRow_[maximumRowsExtra_];
01707 CoinBigIndex put = 0;
01708 while ( iRow != maximumRowsExtra_ ) {
01709
01710 CoinBigIndex get = startRow[iRow];
01711 CoinBigIndex getEnd = startRow[iRow] + numberInRow_[iRow];
01712
01713 startRow[iRow] = put;
01714 CoinBigIndex i;
01715 for ( i = get; i < getEnd; i++ ) {
01716 indexColumn[put] = indexColumn[i];
01717 convertRowToColumn[put] = convertRowToColumn[i];
01718 put++;
01719 }
01720 iRow = nextRow_[iRow];
01721 }
01722 numberCompressions_++;
01723 startRow[maximumRowsExtra_] = put;
01724 space = lengthAreaU_ - put;
01725 if ( space < extraNeeded + number + 2 ) {
01726
01727
01728
01729 status_ = -99;
01730 return false; }
01731 }
01732 CoinBigIndex put = startRow[maximumRowsExtra_];
01733 int next = nextRow_[iRow];
01734 int last = lastRow_[iRow];
01735
01736
01737 nextRow_[last] = next;
01738 lastRow_[next] = last;
01739
01740 last = lastRow_[maximumRowsExtra_];
01741 nextRow_[last] = iRow;
01742 lastRow_[maximumRowsExtra_] = iRow;
01743 lastRow_[iRow] = last;
01744 nextRow_[iRow] = maximumRowsExtra_;
01745
01746 CoinBigIndex get = startRow[iRow];
01747
01748 startRow[iRow] = put;
01749 while ( number ) {
01750 number--;
01751 indexColumnU_[put] = indexColumnU_[get];
01752 convertRowToColumn[put] = convertRowToColumn[get];
01753 put++;
01754 get++;
01755 }
01756
01757 startRow[maximumRowsExtra_] = put + extraNeeded + 4;
01758 return true;
01759 }
01760
01761
01762
01763
01764 bool
01765 CoinFactorization::getColumnSpaceIterateR ( int iColumn, double value,
01766 int iRow)
01767 {
01768 double * elementR = elementR_ + lengthAreaR_;
01769 int * indexRowR = indexRowR_ + lengthAreaR_;
01770 int * startR = startColumnR_+maximumPivots_+1;
01771 int number = numberInColumnPlus_[iColumn];
01772
01773
01774 if (lengthAreaR_-startR[maximumColumnsExtra_]<number+1) {
01775
01776 int jColumn = nextColumn_[maximumColumnsExtra_];
01777 CoinBigIndex put = 0;
01778 while ( jColumn != maximumColumnsExtra_ ) {
01779
01780 CoinBigIndex get;
01781 CoinBigIndex getEnd;
01782
01783 get = startR[jColumn];
01784 getEnd = get + numberInColumnPlus_[jColumn];
01785 startR[jColumn] = put;
01786 CoinBigIndex i;
01787 for ( i = get; i < getEnd; i++ ) {
01788 indexRowR[put] = indexRowR[i];
01789 elementR[put] = elementR[i];
01790 put++;
01791 }
01792 jColumn = nextColumn_[jColumn];
01793 }
01794 numberCompressions_++;
01795 startR[maximumColumnsExtra_]=put;
01796 }
01797
01798 if (lengthAreaR_-startR[maximumColumnsExtra_]<number+1)
01799 return false;
01800
01801 int next = nextColumn_[iColumn];
01802 int last = lastColumn_[iColumn];
01803
01804
01805 nextColumn_[last] = next;
01806 lastColumn_[next] = last;
01807
01808 CoinBigIndex put = startR[maximumColumnsExtra_];
01809
01810 last = lastColumn_[maximumColumnsExtra_];
01811 nextColumn_[last] = iColumn;
01812 lastColumn_[maximumColumnsExtra_] = iColumn;
01813 lastColumn_[iColumn] = last;
01814 nextColumn_[iColumn] = maximumColumnsExtra_;
01815
01816
01817 CoinBigIndex get = startR[iColumn];
01818 startR[iColumn] = put;
01819 int i = 0;
01820 for (i=0 ; i < number; i ++ ) {
01821 elementR[put]= elementR[get];
01822 indexRowR[put++] = indexRowR[get++];
01823 }
01824
01825 elementR[put]=value;
01826 indexRowR[put++]=iRow;
01827 numberInColumnPlus_[iColumn]++;
01828
01829 startR[maximumColumnsExtra_] = min(put + 4,lengthAreaR_);
01830 return true;
01831 }
01832
01833
01834
01835
01836
01837 CoinBigIndex
01838 CoinFactorization::getColumnSpaceIterate ( int iColumn, double value,
01839 int iRow)
01840 {
01841 if (numberInColumnPlus_) {
01842 delete [] numberInColumnPlus_;
01843 numberInColumnPlus_=NULL;
01844 }
01845 int number = numberInColumn_[iColumn];
01846 int iNext=nextColumn_[iColumn];
01847 CoinBigIndex space = startColumnU_[iNext]-startColumnU_[iColumn];
01848 CoinBigIndex put;
01849 if ( space < number + 1 ) {
01850
01851 if (lengthAreaU_-startColumnU_[maximumColumnsExtra_]<number+1) {
01852
01853 int jColumn = nextColumn_[maximumColumnsExtra_];
01854 CoinBigIndex put = 0;
01855 while ( jColumn != maximumColumnsExtra_ ) {
01856
01857 CoinBigIndex get;
01858 CoinBigIndex getEnd;
01859
01860 get = startColumnU_[jColumn];
01861 getEnd = get + numberInColumn_[jColumn];
01862 startColumnU_[jColumn] = put;
01863 CoinBigIndex i;
01864 for ( i = get; i < getEnd; i++ ) {
01865 double value = elementU_[i];
01866 if (value) {
01867 indexRowU_[put] = indexRowU_[i];
01868 elementU_[put] = value;
01869 put++;
01870 } else {
01871 numberInColumn_[jColumn]--;
01872 }
01873 }
01874 jColumn = nextColumn_[jColumn];
01875 }
01876 numberCompressions_++;
01877 startColumnU_[maximumColumnsExtra_]=put;
01878
01879 CoinBigIndex *convertRowToColumn = convertRowToColumnU_;
01880 CoinBigIndex j = 0;
01881 CoinBigIndex *startRow = startRowU_;
01882
01883 int iRow;
01884 for ( iRow = 0; iRow < numberRowsExtra_; iRow++ ) {
01885 startRow[iRow] = j;
01886 j += numberInRow_[iRow];
01887 }
01888 factorElements_=j;
01889
01890 CoinZeroN ( numberInRow_, numberRowsExtra_ );
01891 int i;
01892 for ( i = 0; i < numberRowsExtra_; i++ ) {
01893 CoinBigIndex start = startColumnU_[i];
01894 CoinBigIndex end = start + numberInColumn_[i];
01895
01896 CoinBigIndex j;
01897 for ( j = start; j < end; j++ ) {
01898 int iRow = indexRowU_[j];
01899 int iLook = numberInRow_[iRow];
01900
01901 numberInRow_[iRow] = iLook + 1;
01902 CoinBigIndex k = startRow[iRow] + iLook;
01903
01904 indexColumnU_[k] = i;
01905 convertRowToColumn[k] = j;
01906 }
01907 }
01908 }
01909
01910 if (lengthAreaU_-startColumnU_[maximumColumnsExtra_]>=number+1) {
01911 int next = nextColumn_[iColumn];
01912 int last = lastColumn_[iColumn];
01913
01914
01915 nextColumn_[last] = next;
01916 lastColumn_[next] = last;
01917
01918 put = startColumnU_[maximumColumnsExtra_];
01919
01920 last = lastColumn_[maximumColumnsExtra_];
01921 nextColumn_[last] = iColumn;
01922 lastColumn_[maximumColumnsExtra_] = iColumn;
01923 lastColumn_[iColumn] = last;
01924 nextColumn_[iColumn] = maximumColumnsExtra_;
01925
01926
01927 CoinBigIndex get = startColumnU_[iColumn];
01928 startColumnU_[iColumn] = put;
01929 int i = 0;
01930 for (i=0 ; i < number; i ++ ) {
01931 double value = elementU_[get];
01932 int iRow=indexRowU_[get++];
01933 if (value) {
01934 elementU_[put]= value;
01935 int n=numberInRow_[iRow];
01936 int start = startRowU_[iRow];
01937 int j;
01938 for (j=start;j<start+n;j++) {
01939 if (indexColumnU_[j]==iColumn) {
01940 convertRowToColumnU_[j]=put;
01941 break;
01942 }
01943 }
01944 assert (j<start+n);
01945 indexRowU_[put++] = iRow;
01946 } else {
01947 assert (!numberInRow_[iRow]);
01948 numberInColumn_[iColumn]--;
01949 }
01950 }
01951
01952 int n=numberInRow_[iRow];
01953 int start = startRowU_[iRow];
01954 int j;
01955 for (j=start;j<start+n;j++) {
01956 if (indexColumnU_[j]==iColumn) {
01957 convertRowToColumnU_[j]=put;
01958 break;
01959 }
01960 }
01961 assert (j<start+n);
01962 elementU_[put]=value;
01963 indexRowU_[put]=iRow;
01964 numberInColumn_[iColumn]++;
01965
01966 startColumnU_[maximumColumnsExtra_] = min(put + 4,lengthAreaU_);
01967 } else {
01968
01969 put=-1;
01970 }
01971 } else {
01972
01973 put=startColumnU_[iColumn]+numberInColumn_[iColumn];
01974 int n=numberInRow_[iRow];
01975 int start = startRowU_[iRow];
01976 int j;
01977 for (j=start;j<start+n;j++) {
01978 if (indexColumnU_[j]==iColumn) {
01979 convertRowToColumnU_[j]=put;
01980 break;
01981 }
01982 }
01983 assert (j<start+n);
01984 elementU_[put]=value;
01985 indexRowU_[put]=iRow;
01986 numberInColumn_[iColumn]++;
01987 }
01988 return put;
01989 }