00001
00002
00003
00004 #if defined(_MSC_VER)
00005
00006 # pragma warning(disable:4786)
00007 #endif
00008
00009 #include <cassert>
00010 #include <stdio.h>
00011 #include "CoinFactorization.hpp"
00012 #include "CoinIndexedVector.hpp"
00013 #include "CoinHelperFunctions.hpp"
00014 #if DENSE_CODE==1
00015
00016 extern "C" int dgetrf_(int *m, int *n, double *a, int * lda,
00017 int *ipiv, int *info);
00018 #endif
00019
00020
00021 int
00022 CoinFactorization::factorSparse ( )
00023 {
00024 int *indexRow = indexRowU_;
00025 int *indexColumn = indexColumnU_;
00026 double *element = elementU_;
00027 int count = 1;
00028 double *workArea = new double [ numberRows_ ];
00029
00030
00031 int denseThreshold=denseThreshold_;
00032
00033 CoinZeroN ( workArea, numberRows_ );
00034
00035 CoinBigIndex workSize = 1000;
00036 unsigned int *workArea2 = ( unsigned int * ) new int [ workSize ];
00037
00038 int larger;
00039
00040 if ( numberRows_ < numberColumns_ ) {
00041 larger = numberColumns_;
00042 } else {
00043 larger = numberRows_;
00044 }
00045 int status = 0;
00046
00047 if (biasLU_<3) {
00048 int pivotColumn;
00049 for ( pivotColumn = 0; pivotColumn < numberColumns_;
00050 pivotColumn++ ) {
00051 if ( numberInColumn_[pivotColumn] == 1 ) {
00052 CoinBigIndex start = startColumnU_[pivotColumn];
00053 double value = element[start];
00054 if ( value == slackValue_ && numberInColumnPlus_[pivotColumn] == 0 ) {
00055
00056 int iRow = indexRowU_[start];
00057
00058
00059 totalElements_ -= numberInRow_[iRow];
00060 if ( !pivotColumnSingleton ( iRow, pivotColumn ) ) {
00061 status = -99;
00062 count=biggerDimension_+1;
00063 break;
00064 }
00065 pivotColumn_[numberGoodU_] = pivotColumn;
00066 numberGoodU_++;
00067 }
00068 }
00069 }
00070 }
00071 numberSlacks_ = numberGoodU_;
00072 int *nextCount = nextCount_;
00073 int *numberInRow = numberInRow_;
00074 int *numberInColumn = numberInColumn_;
00075 CoinBigIndex *startRow = startRowU_;
00076 CoinBigIndex *startColumn = startColumnU_;
00077 double pivotTolerance = pivotTolerance_;
00078 int numberTrials = numberTrials_;
00079 int numberRows = numberRows_;
00080
00081 separateLinks(1,(biasLU_>1));
00082 while ( count <= biggerDimension_ ) {
00083 CoinBigIndex minimumCount = INT_MAX;
00084 CoinBigIndex minimumCost = INT_MAX;
00085
00086 count = 1;
00087 bool stopping = false;
00088 int pivotRow = -1;
00089 int pivotColumn = -1;
00090 int pivotRowPosition = -1;
00091 int pivotColumnPosition = -1;
00092 int look = firstCount_[count];
00093 int trials = 0;
00094
00095 while ( !stopping ) {
00096 if ( count == 1 && firstCount_[1] >= 0 &&!biasLU_) {
00097
00098 while ( look >= 0 ) {
00099 if ( look < numberRows_ ) {
00100 look = nextCount_[look];
00101 } else {
00102 int iColumn = look - numberRows_;
00103
00104 #if COIN_DEBUG
00105 if ( numberInColumn_[iColumn] != count ) {
00106 abort ( );
00107 }
00108 #endif
00109 CoinBigIndex start = startColumnU_[iColumn];
00110 int iRow = indexRow[start];
00111
00112 pivotRow = iRow;
00113 pivotRowPosition = start;
00114 pivotColumn = iColumn;
00115 pivotColumnPosition = -1;
00116 stopping = true;
00117 look = -1;
00118 break;
00119 }
00120 }
00121 if ( !stopping ) {
00122
00123 look = firstCount_[1];
00124 }
00125 }
00126 while ( look >= 0 ) {
00127 if ( look < numberRows_ ) {
00128 int iRow = look;
00129
00130 #if COIN_DEBUG
00131 if ( numberInRow[iRow] != count ) {
00132 abort ( );
00133 }
00134 #endif
00135 look = nextCount[look];
00136 bool rejected = false;
00137 CoinBigIndex start = startRow[iRow];
00138 CoinBigIndex end = start + count;
00139
00140 CoinBigIndex i;
00141 for ( i = start; i < end; i++ ) {
00142 int iColumn = indexColumn[i];
00143 CoinBigIndex cost = ( count - 1 ) * numberInColumn[iColumn];
00144
00145 if ( cost < minimumCost ) {
00146 CoinBigIndex where = startColumn[iColumn];
00147 double minimumValue = element[where];
00148
00149 minimumValue = fabs ( minimumValue ) * pivotTolerance;
00150 while ( indexRow[where] != iRow ) {
00151 where++;
00152 }
00153 #if COIN_DEBUG
00154 {
00155 CoinBigIndex end_debug = startColumn[iColumn] +
00156 numberInColumn[iColumn];
00157
00158 if ( where >= end_debug ) {
00159 abort ( );
00160 }
00161 }
00162 #endif
00163 double value = element[where];
00164
00165 value = fabs ( value );
00166 if ( value >= minimumValue ) {
00167 minimumCost = cost;
00168 minimumCount = numberInColumn[iColumn];
00169 pivotRow = iRow;
00170 pivotRowPosition = -1;
00171 pivotColumn = iColumn;
00172 pivotColumnPosition = i;
00173 if ( minimumCount < count ) {
00174 stopping = true;
00175 look = -1;
00176 break;
00177 }
00178 } else if ( pivotRow == -1 ) {
00179 rejected = true;
00180 }
00181 }
00182 }
00183 trials++;
00184 if ( trials >= numberTrials && pivotRow >= 0 ) {
00185 stopping = true;
00186 look = -1;
00187 break;
00188 }
00189 if ( rejected ) {
00190
00191
00192 deleteLink ( iRow );
00193 addLink ( iRow, biggerDimension_ + 1 );
00194 }
00195 } else {
00196 int iColumn = look - numberRows;
00197
00198 #if COIN_DEBUG
00199 if ( numberInColumn[iColumn] != count ) {
00200 abort ( );
00201 }
00202 #endif
00203 look = nextCount[look];
00204 CoinBigIndex start = startColumn[iColumn];
00205 CoinBigIndex end = start + numberInColumn[iColumn];
00206 double minimumValue = element[start];
00207
00208 minimumValue = fabs ( minimumValue ) * pivotTolerance;
00209 CoinBigIndex i;
00210 for ( i = start; i < end; i++ ) {
00211 double value = element[i];
00212
00213 value = fabs ( value );
00214 if ( value >= minimumValue ) {
00215 int iRow = indexRow[i];
00216 CoinBigIndex cost = ( count - 1 ) * numberInRow[iRow];
00217
00218 if ( cost < minimumCost ) {
00219 minimumCost = cost;
00220 minimumCount = numberInRow[iRow];
00221 pivotRow = iRow;
00222 pivotRowPosition = i;
00223 pivotColumn = iColumn;
00224 pivotColumnPosition = -1;
00225 if ( minimumCount <= count + 1 ) {
00226 stopping = true;
00227 look = -1;
00228 break;
00229 }
00230 }
00231 }
00232 }
00233 trials++;
00234 if ( trials >= numberTrials && pivotRow >= 0 ) {
00235 stopping = true;
00236 look = -1;
00237 break;
00238 }
00239 }
00240 }
00241
00242 if ( !stopping ) {
00243 count++;
00244 if ( count <= biggerDimension_ ) {
00245 look = firstCount_[count];
00246 } else {
00247 stopping = true;
00248 }
00249 } else {
00250 if ( pivotRow >= 0 ) {
00251 int numberDoRow = numberInRow_[pivotRow] - 1;
00252 int numberDoColumn = numberInColumn_[pivotColumn] - 1;
00253
00254 totalElements_ -= ( numberDoRow + numberDoColumn + 1 );
00255 if ( numberDoColumn > 0 ) {
00256 if ( numberDoRow > 0 ) {
00257 if ( numberDoColumn > 1 ) {
00258
00259
00260
00261 int increment = numberDoColumn + 1 + 4;
00262
00263 if ( increment & 15 ) {
00264 increment = increment & ( ~15 );
00265 increment += 16;
00266 }
00267 int increment2 =
00268
00269 ( increment + COINFACTORIZATION_BITS_PER_INT - 1 ) >> COINFACTORIZATION_SHIFT_PER_INT;
00270 CoinBigIndex size = increment2 * numberDoRow;
00271
00272 if ( size > workSize ) {
00273 workSize = size;
00274 delete [] workArea2 ;
00275 workArea2 = ( unsigned int * ) new int [ workSize ];
00276 }
00277 bool goodPivot;
00278
00279 if ( larger < 32766 ) {
00280
00281 goodPivot = pivot ( pivotRow, pivotColumn,
00282 pivotRowPosition, pivotColumnPosition,
00283 workArea, workArea2, increment,
00284 increment2, ( short * ) markRow_ ,
00285 32767);
00286 } else {
00287
00288 goodPivot = pivot ( pivotRow, pivotColumn,
00289 pivotRowPosition, pivotColumnPosition,
00290 workArea, workArea2, increment,
00291 increment2, ( int * ) markRow_ ,
00292 INT_MAX);
00293 }
00294 if ( !goodPivot ) {
00295 status = -99;
00296 count=biggerDimension_+1;
00297 break;
00298 }
00299 } else {
00300 if ( !pivotOneOtherRow ( pivotRow, pivotColumn ) ) {
00301 status = -99;
00302 count=biggerDimension_+1;
00303 break;
00304 }
00305 }
00306 } else {
00307 if ( !pivotRowSingleton ( pivotRow, pivotColumn ) ) {
00308 status = -99;
00309 count=biggerDimension_+1;
00310 break;
00311 }
00312 }
00313 } else {
00314 if ( !pivotColumnSingleton ( pivotRow, pivotColumn ) ) {
00315 status = -99;
00316 count=biggerDimension_+1;
00317 break;
00318 }
00319 }
00320 pivotColumn_[numberGoodU_] = pivotColumn;
00321 numberGoodU_++;
00322 }
00323 }
00324 }
00325 #if COIN_DEBUG==2
00326 checkConsistency ( );
00327 #endif
00328 if (denseThreshold) {
00329
00330 int leftRows=numberRows_-numberGoodU_;
00331 double full = leftRows;
00332 full *= full;
00333 assert (full>=0.0);
00334 double leftElements = totalElements_;
00335
00336
00337 if ((1.5*leftElements>full&&leftRows>denseThreshold_)) {
00338
00339 if (status!=0)
00340 break;
00341 #ifdef DENSE_CODE
00342 int check=0;
00343 for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
00344 if (numberInColumn_[iColumn])
00345 check++;
00346 }
00347 if (check!=leftRows) {
00348 printf("** mismatch %d columns left, %d rows\n",check,leftRows);
00349 denseThreshold=0;
00350 } else {
00351 status=2;
00352 if ((messageLevel_&4)!=0)
00353 std::cout<<" Went dense at "<<leftRows<<" rows "<<
00354 totalElements_<<" "<<full<<" "<<leftElements<<std::endl;
00355 break;
00356 }
00357 #endif
00358 }
00359 }
00360 }
00361 delete [] workArea ;
00362 delete [] workArea2 ;
00363 return status;
00364 }
00365
00366
00367 int CoinFactorization::factorDense()
00368 {
00369 int status=0;
00370 #ifdef DENSE_CODE
00371 numberDense_=numberRows_-numberGoodU_;
00372 if (sizeof(CoinBigIndex)==4&&numberDense_>=2<<15) {
00373 abort();
00374 }
00375 CoinBigIndex full = numberDense_*numberDense_;
00376 totalElements_=full;
00377 denseArea_= new double [full];
00378 memset(denseArea_,0,full*sizeof(double));
00379 densePermute_= new int [numberDense_];
00380
00381 int i;
00382 for (i=0;i<numberRows_;i++)
00383 lastRow_[i]=0;
00384 int * indexRow = indexRowU_;
00385 double * element = elementU_;
00386 for (int i=0;i<numberGoodU_;i++) {
00387 int iRow=pivotRowL_[i];
00388 lastRow_[iRow]=-1;
00389 }
00390 int which=0;
00391 for (i=0;i<numberRows_;i++) {
00392 if (!lastRow_[i]) {
00393 lastRow_[i]=which;
00394 nextRow_[i]=numberGoodU_+which;
00395 densePermute_[which]=i;
00396 which++;
00397 }
00398 }
00399
00400 CoinBigIndex endL=startColumnL_[numberGoodL_];
00401
00402 double * column = denseArea_;
00403 int rowsDone=0;
00404
00405 for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
00406 if (numberInColumn_[iColumn]) {
00407
00408 CoinBigIndex start = startColumnU_[iColumn];
00409 int number = numberInColumn_[iColumn];
00410 CoinBigIndex end = start+number;
00411 for (CoinBigIndex i=start;i<end;i++) {
00412 int iRow=indexRow[i];
00413 iRow=lastRow_[iRow];
00414 column[iRow]=element[i];
00415 }
00416 column+=numberDense_;
00417 while (lastRow_[rowsDone]<0) {
00418 rowsDone++;
00419 }
00420 pivotRowL_[numberGoodU_]=rowsDone;
00421 startColumnL_[numberGoodU_+1]=endL;
00422 numberInColumn_[iColumn]=0;
00423 pivotColumn_[numberGoodU_]=iColumn;
00424 pivotRegion_[numberGoodU_]=1.0;
00425 numberGoodU_++;
00426 }
00427 }
00428 assert(numberGoodU_==numberRows_);
00429 #ifndef NEWSTYLE
00430 numberGoodL_=numberRows_;
00431
00432
00433 int info;
00434 dgetrf_(&numberDense_,&numberDense_,denseArea_,&numberDense_,densePermute_,
00435 &info);
00436
00437 if(info)
00438 status = -1;
00439 #else
00440 numberGoodU_ = numberRows_-numberDense_;
00441 int base = numberGoodU_;
00442 int iDense;
00443 double tolerance = zeroTolerance_;
00444 tolerance = 1.0e-30;
00445
00446 for (iDense=0;iDense<numberDense_;iDense++) {
00447
00448 int iColumn = pivotColumn_[base+iDense];
00449 int next = nextColumn_[iColumn];
00450 int numberInPivotColumn = iDense;
00451 CoinBigIndex space = startColumnU_[next]
00452 - startColumnU_[iColumn]
00453 - numberInColumnPlus_[next];
00454
00455 if ( numberInPivotColumn > space ) {
00456
00457 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
00458 return -99;
00459 }
00460 }
00461
00462 numberInColumn_[iColumn]=numberInPivotColumn;
00463 }
00464 if ( lengthL_ + full*0.5 > lengthAreaL_ ) {
00465
00466 std::cout << "more memory needed in middle of invert" << std::endl;
00467 return -99;
00468 }
00469 for (iDense=0;iDense<numberDense_;iDense++) {
00470 int iRow;
00471 int jDense;
00472 int pivotRow=-1;
00473 double * element = denseArea_+iDense*numberDense_;
00474 double largest = 1.0e-12;
00475 for (iRow=iDense;iRow<numberDense_;iRow++) {
00476 if (fabs(element[iRow])>largest) {
00477 largest = fabs(element[iRow]);
00478 pivotRow = iRow;
00479 }
00480 }
00481 if ( pivotRow >= 0 ) {
00482 int iColumn = pivotColumn_[base+iDense];
00483 double pivotElement=element[pivotRow];
00484
00485 int originalRow = densePermute_[pivotRow];
00486
00487 nextRow_[originalRow] = numberGoodU_;
00488
00489 densePermute_[pivotRow]=densePermute_[iDense];
00490 densePermute_[iDense] = originalRow;
00491 for (jDense=iDense;jDense<numberDense_;jDense++) {
00492 double value = element[iDense];
00493 element[iDense] = element[pivotRow];
00494 element[pivotRow] = value;
00495 element += numberDense_;
00496 }
00497 double pivotMultiplier = 1.0 / pivotElement;
00498
00499 pivotRegion_[numberGoodU_] = pivotMultiplier;
00500
00501 element = denseArea_+iDense*numberDense_;
00502 CoinBigIndex l = lengthL_;
00503 startColumnL_[numberGoodL_] = l;
00504 for (iRow=iDense+1;iRow<numberDense_;iRow++) {
00505 double value = element[iRow]*pivotMultiplier;
00506 element[iRow] = value;
00507 if (fabs(value)>tolerance) {
00508 indexRowL_[l] = densePermute_[iRow];
00509 elementL_[l++] = value;
00510 }
00511 }
00512 pivotRowL_[numberGoodL_++] = originalRow;
00513 lengthL_ = l;
00514 startColumnL_[numberGoodL_] = l;
00515
00516 CoinBigIndex start = startColumnU_[iColumn];
00517 for (iRow=0;iRow<iDense;iRow++) {
00518 if (fabs(element[iRow])>tolerance) {
00519 indexRowU_[start] = densePermute_[iRow];
00520 elementU_[start++] = element[iRow];
00521 }
00522 }
00523 numberInColumn_[iColumn]=0;
00524 numberInColumnPlus_[iColumn] += start-startColumnU_[iColumn];
00525 startColumnU_[iColumn]=start;
00526
00527 double * element2 = element+numberDense_;
00528 for (jDense=iDense+1;jDense<numberDense_;jDense++) {
00529 double value = element2[iDense];
00530 for (iRow=iDense+1;iRow<numberDense_;iRow++) {
00531 double oldValue=element2[iRow];
00532 element2[iRow] -= value*element[iRow];
00533 if (oldValue&&!element2[iRow]) {
00534 printf("Updated element for column %d, row %d old %g",
00535 pivotColumn_[base+jDense],densePermute_[iRow],oldValue);
00536 printf(" new %g\n",element2[iRow]);
00537 }
00538 }
00539 element2 += numberDense_;
00540 }
00541 numberGoodU_++;
00542 } else {
00543 return -1;
00544 }
00545 }
00546
00547 delete [] denseArea_;
00548 denseArea_ = NULL;
00549
00550 delete [] densePermute_;
00551 densePermute_ = NULL;
00552 numberDense_=0;
00553 #endif
00554 #endif
00555 return status;
00556 }
00557
00558 void
00559 CoinFactorization::separateLinks(int count,bool rowsFirst)
00560 {
00561 int next = firstCount_[count];
00562 int firstRow=-1;
00563 int firstColumn=-1;
00564 int lastRow=-1;
00565 int lastColumn=-1;
00566 while(next>=0) {
00567 int next2=nextCount_[next];
00568 if (next>=numberRows_) {
00569 nextCount_[next]=-1;
00570
00571 if (firstColumn>=0) {
00572 lastCount_[next]=lastColumn;
00573 nextCount_[lastColumn]=next;
00574 } else {
00575 lastCount_[next]= -2 - count;
00576 firstColumn=next;
00577 }
00578 lastColumn=next;
00579 } else {
00580
00581 if (firstRow>=0) {
00582 lastCount_[next]=lastRow;
00583 nextCount_[lastRow]=next;
00584 } else {
00585 lastCount_[next]= -2 - count;
00586 firstRow=next;
00587 }
00588 lastRow=next;
00589 }
00590 next=next2;
00591 }
00592 if (rowsFirst&&firstRow>=0) {
00593 firstCount_[count]=firstRow;
00594 nextCount_[lastRow]=firstColumn;
00595 if (firstColumn>=0)
00596 lastCount_[firstColumn]=lastRow;
00597 } else if (firstRow<0) {
00598 firstCount_[count]=firstColumn;
00599 } else if (firstColumn>=0) {
00600 firstCount_[count]=firstColumn;
00601 nextCount_[lastColumn]=firstRow;
00602 if (firstRow>=0)
00603 lastCount_[firstRow]=lastColumn;
00604 }
00605 }