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