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