1 /* $Id: CoinFactorization2.cpp 1448 2011-06-19 15:34:41Z stefan $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #if defined(_MSC_VER)
7 // Turn off compiler warning about long names
8 # pragma warning(disable:4786)
9 #endif
10
11 #include "CoinUtilsConfig.h"
12
13 #include <cassert>
14 #include <cfloat>
15 #include <stdio.h>
16 #include "CoinFactorization.hpp"
17 #include "CoinIndexedVector.hpp"
18 #include "CoinHelperFunctions.hpp"
19 #include "CoinFinite.hpp"
20 #if DENSE_CODE==1
21 // using simple lapack interface
22 extern "C"
23 {
24 /** LAPACK Fortran subroutine DGETRF. */
25 void F77_FUNC(dgetrf,DGETRF)(ipfint * m, ipfint *n,
26 double *A, ipfint *ldA,
27 ipfint * ipiv, ipfint *info);
28 }
29 #endif
30 #ifndef NDEBUG
31 static int counter1=0;
32 #endif
33 // factorSparse. Does sparse phase of factorization
34 //return code is <0 error, 0= finished
35 int
factorSparse()36 CoinFactorization::factorSparse ( )
37 {
38 int larger;
39
40 if ( numberRows_ < numberColumns_ ) {
41 larger = numberColumns_;
42 } else {
43 larger = numberRows_;
44 }
45 int returnCode;
46 #define LARGELIMIT 65530
47 #define SMALL_SET 65531
48 #define SMALL_UNSET (SMALL_SET+1)
49 #define LARGE_SET COIN_INT_MAX-10
50 #define LARGE_UNSET (LARGE_SET+1)
51 if ( larger < LARGELIMIT )
52 returnCode = factorSparseSmall();
53 else
54 returnCode = factorSparseLarge();
55 return returnCode;
56 }
57 // factorSparse. Does sparse phase of factorization
58 //return code is <0 error, 0= finished
59 int
factorSparseSmall()60 CoinFactorization::factorSparseSmall ( )
61 {
62 int *indexRow = indexRowU_.array();
63 int *indexColumn = indexColumnU_.array();
64 CoinFactorizationDouble *element = elementU_.array();
65 int count = 1;
66 workArea_.conditionalNew(numberRows_);
67 CoinFactorizationDouble * workArea = workArea_.array();
68 #ifndef NDEBUG
69 counter1++;
70 #endif
71 // when to go dense
72 int denseThreshold=denseThreshold_;
73
74 CoinZeroN ( workArea, numberRows_ );
75 //get space for bit work area
76 CoinBigIndex workSize = 1000;
77 workArea2_.conditionalNew(workSize);
78 unsigned int * workArea2 = workArea2_.array();
79
80 //set markRow so no rows updated
81 unsigned short * markRow = reinterpret_cast<unsigned short *> (markRow_.array());
82 CoinFillN ( markRow, numberRows_, static_cast<unsigned short> (SMALL_UNSET));
83 int status = 0;
84 //do slacks first
85 int * numberInRow = numberInRow_.array();
86 int * numberInColumn = numberInColumn_.array();
87 int * numberInColumnPlus = numberInColumnPlus_.array();
88 CoinBigIndex * startColumnU = startColumnU_.array();
89 CoinBigIndex * startColumnL = startColumnL_.array();
90 if (biasLU_<3&&numberColumns_==numberRows_) {
91 int iPivotColumn;
92 int * pivotColumn = pivotColumn_.array();
93 int * nextRow = nextRow_.array();
94 int * lastRow = lastRow_.array();
95 for ( iPivotColumn = 0; iPivotColumn < numberColumns_;
96 iPivotColumn++ ) {
97 if ( numberInColumn[iPivotColumn] == 1 ) {
98 CoinBigIndex start = startColumnU[iPivotColumn];
99 CoinFactorizationDouble value = element[start];
100 if ( value == slackValue_ && numberInColumnPlus[iPivotColumn] == 0 ) {
101 // treat as slack
102 int iRow = indexRow[start];
103 // but only if row not marked
104 if (numberInRow[iRow]>0) {
105 totalElements_ -= numberInRow[iRow];
106 //take out this bit of indexColumnU
107 int next = nextRow[iRow];
108 int last = lastRow[iRow];
109
110 nextRow[last] = next;
111 lastRow[next] = last;
112 nextRow[iRow] = numberGoodU_; //use for permute
113 lastRow[iRow] = -2; //mark
114 //modify linked list for pivots
115 deleteLink ( iRow );
116 numberInRow[iRow]=-1;
117 numberInColumn[iPivotColumn]=0;
118 numberGoodL_++;
119 startColumnL[numberGoodL_] = 0;
120 pivotColumn[numberGoodU_] = iPivotColumn;
121 numberGoodU_++;
122 }
123 }
124 }
125 }
126 // redo
127 preProcess(4);
128 CoinFillN ( markRow, numberRows_, static_cast<unsigned short> (SMALL_UNSET));
129 }
130 numberSlacks_ = numberGoodU_;
131 int *nextCount = nextCount_.array();
132 int *firstCount = firstCount_.array();
133 CoinBigIndex *startRow = startRowU_.array();
134 CoinBigIndex *startColumn = startColumnU;
135 //#define UGLY_COIN_FACTOR_CODING
136 #ifdef UGLY_COIN_FACTOR_CODING
137 CoinFactorizationDouble *elementL = elementL_.array();
138 int *indexRowL = indexRowL_.array();
139 int *saveColumn = saveColumn_.array();
140 int *nextRow = nextRow_.array();
141 int *lastRow = lastRow_.array() ;
142 #endif
143 double pivotTolerance = pivotTolerance_;
144 int numberTrials = numberTrials_;
145 int numberRows = numberRows_;
146 // Put column singletons first - (if false)
147 separateLinks(1,(biasLU_>1));
148 #ifndef NDEBUG
149 int counter2=0;
150 #endif
151 while ( count <= biggerDimension_ ) {
152 #ifndef NDEBUG
153 counter2++;
154 int badRow=-1;
155 if (counter1==-1&&counter2>=0) {
156 // check counts consistent
157 for (int iCount=1;iCount<numberRows_;iCount++) {
158 int look = firstCount[iCount];
159 while ( look >= 0 ) {
160 if ( look < numberRows_ ) {
161 int iRow = look;
162 if (iRow==badRow)
163 printf("row count for row %d is %d\n",iCount,iRow);
164 if ( numberInRow[iRow] != iCount ) {
165 printf("failed debug on %d entry to factorSparse and %d try\n",
166 counter1,counter2);
167 printf("row %d - count %d number %d\n",iRow,iCount,numberInRow[iRow]);
168 abort();
169 }
170 look = nextCount[look];
171 } else {
172 int iColumn = look - numberRows;
173 if ( numberInColumn[iColumn] != iCount ) {
174 printf("failed debug on %d entry to factorSparse and %d try\n",
175 counter1,counter2);
176 printf("column %d - count %d number %d\n",iColumn,iCount,numberInColumn[iColumn]);
177 abort();
178 }
179 look = nextCount[look];
180 }
181 }
182 }
183 }
184 #endif
185 CoinBigIndex minimumCount = COIN_INT_MAX;
186 double minimumCost = COIN_DBL_MAX;
187
188 int iPivotRow = -1;
189 int iPivotColumn = -1;
190 int pivotRowPosition = -1;
191 int pivotColumnPosition = -1;
192 int look = firstCount[count];
193 int trials = 0;
194 int * pivotColumn = pivotColumn_.array();
195
196 if ( count == 1 && firstCount[1] >= 0 &&!biasLU_) {
197 //do column singletons first to put more in U
198 while ( look >= 0 ) {
199 if ( look < numberRows_ ) {
200 look = nextCount[look];
201 } else {
202 int iColumn = look - numberRows_;
203
204 assert ( numberInColumn[iColumn] == count );
205 CoinBigIndex start = startColumnU[iColumn];
206 int iRow = indexRow[start];
207
208 iPivotRow = iRow;
209 pivotRowPosition = start;
210 iPivotColumn = iColumn;
211 assert (iPivotRow>=0&&iPivotColumn>=0);
212 pivotColumnPosition = -1;
213 look = -1;
214 break;
215 }
216 } /* endwhile */
217 if ( iPivotRow < 0 ) {
218 //back to singletons
219 look = firstCount[1];
220 }
221 }
222 while ( look >= 0 ) {
223 if ( look < numberRows_ ) {
224 int iRow = look;
225 #ifndef NDEBUG
226 if ( numberInRow[iRow] != count ) {
227 printf("failed on %d entry to factorSparse and %d try\n",
228 counter1,counter2);
229 printf("row %d - count %d number %d\n",iRow,count,numberInRow[iRow]);
230 abort();
231 }
232 #endif
233 look = nextCount[look];
234 bool rejected = false;
235 CoinBigIndex start = startRow[iRow];
236 CoinBigIndex end = start + count;
237
238 CoinBigIndex i;
239 for ( i = start; i < end; i++ ) {
240 int iColumn = indexColumn[i];
241 assert (numberInColumn[iColumn]>0);
242 double cost = ( count - 1 ) * numberInColumn[iColumn];
243
244 if ( cost < minimumCost ) {
245 CoinBigIndex where = startColumn[iColumn];
246 double minimumValue = element[where];
247
248 minimumValue = fabs ( minimumValue ) * pivotTolerance;
249 while ( indexRow[where] != iRow ) {
250 where++;
251 } /* endwhile */
252 assert ( where < startColumn[iColumn] +
253 numberInColumn[iColumn]);
254 CoinFactorizationDouble value = element[where];
255
256 value = fabs ( value );
257 if ( value >= minimumValue ) {
258 minimumCost = cost;
259 minimumCount = numberInColumn[iColumn];
260 iPivotRow = iRow;
261 pivotRowPosition = -1;
262 iPivotColumn = iColumn;
263 assert (iPivotRow>=0&&iPivotColumn>=0);
264 pivotColumnPosition = i;
265 rejected=false;
266 if ( minimumCount < count ) {
267 look = -1;
268 break;
269 }
270 } else if ( iPivotRow == -1 ) {
271 rejected = true;
272 }
273 }
274 }
275 trials++;
276 if ( trials >= numberTrials && iPivotRow >= 0 ) {
277 look = -1;
278 break;
279 }
280 if ( rejected ) {
281 //take out for moment
282 //eligible when row changes
283 deleteLink ( iRow );
284 addLink ( iRow, biggerDimension_ + 1 );
285 }
286 } else {
287 int iColumn = look - numberRows;
288
289 assert ( numberInColumn[iColumn] == count );
290 look = nextCount[look];
291 CoinBigIndex start = startColumn[iColumn];
292 CoinBigIndex end = start + numberInColumn[iColumn];
293 CoinFactorizationDouble minimumValue = element[start];
294
295 minimumValue = fabs ( minimumValue ) * pivotTolerance;
296 CoinBigIndex i;
297 for ( i = start; i < end; i++ ) {
298 CoinFactorizationDouble value = element[i];
299
300 value = fabs ( value );
301 if ( value >= minimumValue ) {
302 int iRow = indexRow[i];
303 int nInRow = numberInRow[iRow];
304 assert (nInRow>0);
305 double cost = ( count - 1 ) * nInRow;
306
307 if ( cost < minimumCost ) {
308 minimumCost = cost;
309 minimumCount = nInRow;
310 iPivotRow = iRow;
311 pivotRowPosition = i;
312 iPivotColumn = iColumn;
313 assert (iPivotRow>=0&&iPivotColumn>=0);
314 pivotColumnPosition = -1;
315 if ( minimumCount <= count + 1 ) {
316 look = -1;
317 break;
318 }
319 }
320 }
321 }
322 trials++;
323 if ( trials >= numberTrials && iPivotRow >= 0 ) {
324 look = -1;
325 break;
326 }
327 }
328 } /* endwhile */
329 if (iPivotRow>=0) {
330 assert (iPivotRow<numberRows_);
331 int numberDoRow = numberInRow[iPivotRow] - 1;
332 int numberDoColumn = numberInColumn[iPivotColumn] - 1;
333
334 totalElements_ -= ( numberDoRow + numberDoColumn + 1 );
335 if ( numberDoColumn > 0 ) {
336 if ( numberDoRow > 0 ) {
337 if ( numberDoColumn > 1 ) {
338 // if (1) {
339 //need to adjust more for cache and SMP
340 //allow at least 4 extra
341 int increment = numberDoColumn + 1 + 4;
342
343 if ( increment & 15 ) {
344 increment = increment & ( ~15 );
345 increment += 16;
346 }
347 int increment2 =
348
349 ( increment + COINFACTORIZATION_BITS_PER_INT - 1 ) >> COINFACTORIZATION_SHIFT_PER_INT;
350 CoinBigIndex size = increment2 * numberDoRow;
351
352 if ( size > workSize ) {
353 workSize = size;
354 workArea2_.conditionalNew(workSize);
355 workArea2 = workArea2_.array();
356 }
357 bool goodPivot;
358 #ifndef UGLY_COIN_FACTOR_CODING
359 //branch out to best pivot routine
360 goodPivot = pivot ( iPivotRow, iPivotColumn,
361 pivotRowPosition, pivotColumnPosition,
362 workArea, workArea2,
363 increment2, markRow ,
364 SMALL_SET);
365 #else
366 #define FAC_SET SMALL_SET
367 #include "CoinFactorization.hpp"
368 #undef FAC_SET
369 #undef UGLY_COIN_FACTOR_CODING
370 #endif
371 if ( !goodPivot ) {
372 status = -99;
373 count=biggerDimension_+1;
374 break;
375 }
376 } else {
377 if ( !pivotOneOtherRow ( iPivotRow, iPivotColumn ) ) {
378 status = -99;
379 count=biggerDimension_+1;
380 break;
381 }
382 }
383 } else {
384 assert (!numberDoRow);
385 if ( !pivotRowSingleton ( iPivotRow, iPivotColumn ) ) {
386 status = -99;
387 count=biggerDimension_+1;
388 break;
389 }
390 }
391 } else {
392 assert (!numberDoColumn);
393 if ( !pivotColumnSingleton ( iPivotRow, iPivotColumn ) ) {
394 status = -99;
395 count=biggerDimension_+1;
396 break;
397 }
398 }
399 assert (nextRow_.array()[iPivotRow]==numberGoodU_);
400 pivotColumn[numberGoodU_] = iPivotColumn;
401 numberGoodU_++;
402 // This should not need to be trapped here - but be safe
403 if (numberGoodU_==numberRows_)
404 count=biggerDimension_+1;
405 #if COIN_DEBUG==2
406 checkConsistency ( );
407 #endif
408 #if 0
409 // Even if no dense code may be better to use naive dense
410 if (!denseThreshold_&&totalElements_>1000) {
411 int leftRows=numberRows_-numberGoodU_;
412 double full = leftRows;
413 full *= full;
414 assert (full>=0.0);
415 double leftElements = totalElements_;
416 double ratio;
417 if (leftRows>2000)
418 ratio=4.0;
419 else
420 ratio=3.0;
421 if (ratio*leftElements>full)
422 denseThreshold=1;
423 }
424 #endif
425 if (denseThreshold) {
426 // see whether to go dense
427 int leftRows=numberRows_-numberGoodU_;
428 double full = leftRows;
429 full *= full;
430 assert (full>=0.0);
431 double leftElements = totalElements_;
432 //if (leftRows==100)
433 //printf("at 100 %d elements\n",totalElements_);
434 double ratio;
435 if (leftRows>2000)
436 ratio=4.0;
437 else if (leftRows>800)
438 ratio=3.0;
439 else if (leftRows>256)
440 ratio=2.0;
441 else
442 ratio=1.5;
443 if ((ratio*leftElements>full&&leftRows>denseThreshold_)) {
444 //return to do dense
445 if (status!=0)
446 break;
447 int check=0;
448 for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
449 if (numberInColumn[iColumn])
450 check++;
451 }
452 if (check!=leftRows&&denseThreshold_) {
453 //printf("** mismatch %d columns left, %d rows\n",check,leftRows);
454 denseThreshold=0;
455 } else {
456 status=2;
457 if ((messageLevel_&4)!=0)
458 std::cout<<" Went dense at "<<leftRows<<" rows "<<
459 totalElements_<<" "<<full<<" "<<leftElements<<std::endl;
460 if (!denseThreshold_)
461 denseThreshold_=-check; // say how many
462 break;
463 }
464 }
465 }
466 // start at 1 again
467 count = 1;
468 } else {
469 //end of this - onto next
470 count++;
471 }
472 } /* endwhile */
473 workArea_.conditionalDelete() ;
474 workArea2_.conditionalDelete() ;
475 return status;
476 }
477
478 //:method factorDense. Does dense phase of factorization
479 //return code is <0 error, 0= finished
factorDense()480 int CoinFactorization::factorDense()
481 {
482 int status=0;
483 numberDense_=numberRows_-numberGoodU_;
484 if (sizeof(CoinBigIndex)==4&&numberDense_>=2<<15) {
485 abort();
486 }
487 CoinBigIndex full;
488 if (denseThreshold_>0)
489 full = numberDense_*numberDense_;
490 else
491 full = - denseThreshold_*numberDense_;
492 totalElements_=full;
493 denseArea_= new double [full];
494 CoinZeroN(denseArea_,full);
495 densePermute_= new int [numberDense_];
496 int * indexRowU = indexRowU_.array();
497 //mark row lookup using lastRow
498 int i;
499 int * nextRow = nextRow_.array();
500 int * lastRow = lastRow_.array();
501 int * numberInColumn = numberInColumn_.array();
502 int * numberInColumnPlus = numberInColumnPlus_.array();
503 for (i=0;i<numberRows_;i++) {
504 if (lastRow[i]>=0)
505 lastRow[i]=0;
506 }
507 int * indexRow = indexRowU_.array();
508 CoinFactorizationDouble * element = elementU_.array();
509 int which=0;
510 for (i=0;i<numberRows_;i++) {
511 if (!lastRow[i]) {
512 lastRow[i]=which;
513 nextRow[i]=numberGoodU_+which;
514 densePermute_[which]=i;
515 which++;
516 }
517 }
518 //for L part
519 CoinBigIndex * startColumnL = startColumnL_.array();
520 CoinFactorizationDouble * elementL = elementL_.array();
521 int * indexRowL = indexRowL_.array();
522 CoinBigIndex endL=startColumnL[numberGoodL_];
523 //take out of U
524 double * column = denseArea_;
525 int rowsDone=0;
526 int iColumn=0;
527 int * pivotColumn = pivotColumn_.array();
528 CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
529 CoinBigIndex * startColumnU = startColumnU_.array();
530 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
531 if (numberInColumn[iColumn]) {
532 //move
533 CoinBigIndex start = startColumnU[iColumn];
534 int number = numberInColumn[iColumn];
535 CoinBigIndex end = start+number;
536 for (CoinBigIndex i=start;i<end;i++) {
537 int iRow=indexRow[i];
538 iRow=lastRow[iRow];
539 assert (iRow>=0&&iRow<numberDense_);
540 column[iRow]=element[i];
541 } /* endfor */
542 column+=numberDense_;
543 while (lastRow[rowsDone]<0) {
544 rowsDone++;
545 } /* endwhile */
546 nextRow[rowsDone]=numberGoodU_;
547 rowsDone++;
548 startColumnL[numberGoodU_+1]=endL;
549 numberInColumn[iColumn]=0;
550 pivotColumn[numberGoodU_]=iColumn;
551 pivotRegion[numberGoodU_]=1.0;
552 numberGoodU_++;
553 }
554 }
555 #ifdef DENSE_CODE
556 if (denseThreshold_>0) {
557 assert(numberGoodU_==numberRows_);
558 numberGoodL_=numberRows_;
559 //now factorize
560 //dgef(denseArea_,&numberDense_,&numberDense_,densePermute_);
561 int info;
562 F77_FUNC(dgetrf,DGETRF)(&numberDense_,&numberDense_,denseArea_,&numberDense_,densePermute_,
563 &info);
564 // need to check size of pivots
565 if(info)
566 status = -1;
567 return status;
568 }
569 #endif
570 numberGoodU_ = numberRows_-numberDense_;
571 int base = numberGoodU_;
572 int iDense;
573 int numberToDo=-denseThreshold_;
574 denseThreshold_=0;
575 double tolerance = zeroTolerance_;
576 tolerance = 1.0e-30;
577 int * nextColumn = nextColumn_.array();
578 const int * pivotColumnConst = pivotColumn_.array();
579 // make sure we have enough space in L and U
580 for (iDense=0;iDense<numberToDo;iDense++) {
581 //how much space have we got
582 iColumn = pivotColumnConst[base+iDense];
583 int next = nextColumn[iColumn];
584 int numberInPivotColumn = iDense;
585 CoinBigIndex space = startColumnU[next]
586 - startColumnU[iColumn]
587 - numberInColumnPlus[next];
588 //assume no zero elements
589 if ( numberInPivotColumn > space ) {
590 //getColumnSpace also moves fixed part
591 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
592 return -99;
593 }
594 }
595 // set so further moves will work
596 numberInColumn[iColumn]=numberInPivotColumn;
597 }
598 // Fill in ?
599 for (iColumn=numberGoodU_+numberToDo;iColumn<numberRows_;iColumn++) {
600 nextRow[iColumn]=iColumn;
601 startColumnL[iColumn+1]=endL;
602 pivotRegion[iColumn]=1.0;
603 }
604 if ( lengthL_ + full*0.5 > lengthAreaL_ ) {
605 //need more memory
606 if ((messageLevel_&4)!=0)
607 std::cout << "more memory needed in middle of invert" << std::endl;
608 return -99;
609 }
610 CoinFactorizationDouble *elementU = elementU_.array();
611 for (iDense=0;iDense<numberToDo;iDense++) {
612 int iRow;
613 int jDense;
614 int pivotRow=-1;
615 double * element = denseArea_+iDense*numberDense_;
616 CoinFactorizationDouble largest = 1.0e-12;
617 for (iRow=iDense;iRow<numberDense_;iRow++) {
618 if (fabs(element[iRow])>largest) {
619 largest = fabs(element[iRow]);
620 pivotRow = iRow;
621 }
622 }
623 if ( pivotRow >= 0 ) {
624 iColumn = pivotColumnConst[base+iDense];
625 CoinFactorizationDouble pivotElement=element[pivotRow];
626 // get original row
627 int originalRow = densePermute_[pivotRow];
628 // do nextRow
629 nextRow[originalRow] = numberGoodU_;
630 lastRow[originalRow] = -2; //mark
631 // swap
632 densePermute_[pivotRow]=densePermute_[iDense];
633 densePermute_[iDense] = originalRow;
634 for (jDense=iDense;jDense<numberDense_;jDense++) {
635 CoinFactorizationDouble value = element[iDense];
636 element[iDense] = element[pivotRow];
637 element[pivotRow] = value;
638 element += numberDense_;
639 }
640 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
641 //printf("pivotMultiplier %g\n",pivotMultiplier);
642 pivotRegion[numberGoodU_] = pivotMultiplier;
643 // Do L
644 element = denseArea_+iDense*numberDense_;
645 CoinBigIndex l = lengthL_;
646 startColumnL[numberGoodL_] = l; //for luck and first time
647 for (iRow=iDense+1;iRow<numberDense_;iRow++) {
648 CoinFactorizationDouble value = element[iRow]*pivotMultiplier;
649 element[iRow] = value;
650 if (fabs(value)>tolerance) {
651 indexRowL[l] = densePermute_[iRow];
652 elementL[l++] = value;
653 }
654 }
655 numberGoodL_++;
656 lengthL_ = l;
657 startColumnL[numberGoodL_] = l;
658 // update U column
659 CoinBigIndex start = startColumnU[iColumn];
660 for (iRow=0;iRow<iDense;iRow++) {
661 if (fabs(element[iRow])>tolerance) {
662 indexRowU[start] = densePermute_[iRow];
663 elementU[start++] = element[iRow];
664 }
665 }
666 numberInColumn[iColumn]=0;
667 numberInColumnPlus[iColumn] += start-startColumnU[iColumn];
668 startColumnU[iColumn]=start;
669 // update other columns
670 double * element2 = element+numberDense_;
671 for (jDense=iDense+1;jDense<numberToDo;jDense++) {
672 CoinFactorizationDouble value = element2[iDense];
673 for (iRow=iDense+1;iRow<numberDense_;iRow++) {
674 //double oldValue=element2[iRow];
675 element2[iRow] -= value*element[iRow];
676 //if (oldValue&&!element2[iRow]) {
677 //printf("Updated element for column %d, row %d old %g",
678 // pivotColumnConst[base+jDense],densePermute_[iRow],oldValue);
679 //printf(" new %g\n",element2[iRow]);
680 //}
681 }
682 element2 += numberDense_;
683 }
684 numberGoodU_++;
685 } else {
686 return -1;
687 }
688 }
689 // free area (could use L?)
690 delete [] denseArea_;
691 denseArea_ = NULL;
692 // check if can use another array for densePermute_
693 delete [] densePermute_;
694 densePermute_ = NULL;
695 numberDense_=0;
696 return status;
697 }
698 // Separate out links with same row/column count
699 void
separateLinks(int count,bool rowsFirst)700 CoinFactorization::separateLinks(int count,bool rowsFirst)
701 {
702 int *nextCount = nextCount_.array();
703 int *firstCount = firstCount_.array();
704 int *lastCount = lastCount_.array();
705 int next = firstCount[count];
706 int firstRow=-1;
707 int firstColumn=-1;
708 int lastRow=-1;
709 int lastColumn=-1;
710 while(next>=0) {
711 int next2=nextCount[next];
712 if (next>=numberRows_) {
713 nextCount[next]=-1;
714 // Column
715 if (firstColumn>=0) {
716 lastCount[next]=lastColumn;
717 nextCount[lastColumn]=next;
718 } else {
719 lastCount[next]= -2 - count;
720 firstColumn=next;
721 }
722 lastColumn=next;
723 } else {
724 // Row
725 if (firstRow>=0) {
726 lastCount[next]=lastRow;
727 nextCount[lastRow]=next;
728 } else {
729 lastCount[next]= -2 - count;
730 firstRow=next;
731 }
732 lastRow=next;
733 }
734 next=next2;
735 }
736 if (rowsFirst&&firstRow>=0) {
737 firstCount[count]=firstRow;
738 nextCount[lastRow]=firstColumn;
739 if (firstColumn>=0)
740 lastCount[firstColumn]=lastRow;
741 } else if (firstRow<0) {
742 firstCount[count]=firstColumn;
743 } else if (firstColumn>=0) {
744 firstCount[count]=firstColumn;
745 nextCount[lastColumn]=firstRow;
746 if (firstRow>=0)
747 lastCount[firstRow]=lastColumn;
748 }
749 }
750 // Debug - save on file
751 int
saveFactorization(const char * file) const752 CoinFactorization::saveFactorization (const char * file ) const
753 {
754 FILE * fp = fopen(file,"wb");
755 if (fp) {
756 // Save so we can pick up scalars
757 const char * first = reinterpret_cast<const char *> ( &pivotTolerance_);
758 const char * last = reinterpret_cast<const char *> ( &biasLU_);
759 // increment
760 last += sizeof(int);
761 if (fwrite(first,last-first,1,fp)!=1)
762 return 1;
763 // Now arrays
764 if (CoinToFile(elementU_.array(),lengthAreaU_ , fp ))
765 return 1;
766 if (CoinToFile(indexRowU_.array(),lengthAreaU_ , fp ))
767 return 1;
768 if (CoinToFile(indexColumnU_.array(),lengthAreaU_ , fp ))
769 return 1;
770 if (CoinToFile(convertRowToColumnU_.array(),lengthAreaU_ , fp ))
771 return 1;
772 if (CoinToFile(elementByRowL_.array(),lengthAreaL_ , fp ))
773 return 1;
774 if (CoinToFile(indexColumnL_.array(),lengthAreaL_ , fp ))
775 return 1;
776 if (CoinToFile(startRowL_.array() , numberRows_+1, fp ))
777 return 1;
778 if (CoinToFile(elementL_.array(),lengthAreaL_ , fp ))
779 return 1;
780 if (CoinToFile(indexRowL_.array(),lengthAreaL_ , fp ))
781 return 1;
782 if (CoinToFile(startColumnL_.array(),numberRows_ + 1 , fp ))
783 return 1;
784 if (CoinToFile(markRow_.array(),numberRows_ , fp))
785 return 1;
786 if (CoinToFile(saveColumn_.array(),numberColumns_ , fp))
787 return 1;
788 if (CoinToFile(startColumnR_.array() , maximumPivots_ + 1 , fp ))
789 return 1;
790 if (CoinToFile(startRowU_.array(),maximumRowsExtra_ + 1 , fp ))
791 return 1;
792 if (CoinToFile(numberInRow_.array(),maximumRowsExtra_ + 1 , fp ))
793 return 1;
794 if (CoinToFile(nextRow_.array(),maximumRowsExtra_ + 1 , fp ))
795 return 1;
796 if (CoinToFile(lastRow_.array(),maximumRowsExtra_ + 1 , fp ))
797 return 1;
798 if (CoinToFile(pivotRegion_.array(),maximumRowsExtra_ + 1 , fp ))
799 return 1;
800 if (CoinToFile(permuteBack_.array(),maximumRowsExtra_ + 1 , fp ))
801 return 1;
802 if (CoinToFile(permute_.array(),maximumRowsExtra_ + 1 , fp ))
803 return 1;
804 if (CoinToFile(pivotColumnBack_.array(),maximumRowsExtra_ + 1 , fp ))
805 return 1;
806 if (CoinToFile(startColumnU_.array(),maximumColumnsExtra_ + 1 , fp ))
807 return 1;
808 if (CoinToFile(numberInColumn_.array(),maximumColumnsExtra_ + 1 , fp ))
809 return 1;
810 if (CoinToFile(numberInColumnPlus_.array(),maximumColumnsExtra_ + 1 , fp ))
811 return 1;
812 if (CoinToFile(firstCount_.array(),biggerDimension_ + 2 , fp ))
813 return 1;
814 if (CoinToFile(nextCount_.array(),numberRows_ + numberColumns_ , fp ))
815 return 1;
816 if (CoinToFile(lastCount_.array(),numberRows_ + numberColumns_ , fp ))
817 return 1;
818 if (CoinToFile(pivotRowL_.array(),numberRows_ + 1 , fp ))
819 return 1;
820 if (CoinToFile(pivotColumn_.array(),maximumColumnsExtra_ + 1 , fp ))
821 return 1;
822 if (CoinToFile(nextColumn_.array(),maximumColumnsExtra_ + 1 , fp ))
823 return 1;
824 if (CoinToFile(lastColumn_.array(),maximumColumnsExtra_ + 1 , fp ))
825 return 1;
826 if (CoinToFile(denseArea_ , numberDense_*numberDense_, fp ))
827 return 1;
828 if (CoinToFile(densePermute_ , numberDense_, fp ))
829 return 1;
830 fclose(fp);
831 }
832 return 0;
833 }
834 // Debug - restore from file
835 int
restoreFactorization(const char * file,bool factorIt)836 CoinFactorization::restoreFactorization (const char * file , bool factorIt )
837 {
838 FILE * fp = fopen(file,"rb");
839 if (fp) {
840 // Get rid of current
841 gutsOfDestructor();
842 CoinBigIndex newSize=0; // for checking - should be same
843 // Restore so we can pick up scalars
844 char * first = reinterpret_cast<char *> ( &pivotTolerance_);
845 char * last = reinterpret_cast<char *> ( &biasLU_);
846 // increment
847 last += sizeof(int);
848 if (fread(first,last-first,1,fp)!=1)
849 return 1;
850 CoinBigIndex space = lengthAreaL_ - lengthL_;
851 // Now arrays
852 CoinFactorizationDouble *elementU = elementU_.array();
853 if (CoinFromFile(elementU,lengthAreaU_ , fp, newSize )==1)
854 return 1;
855 assert (newSize==lengthAreaU_);
856 int * indexRowU = indexRowU_.array();
857 if (CoinFromFile(indexRowU,lengthAreaU_ , fp, newSize )==1)
858 return 1;
859 assert (newSize==lengthAreaU_);
860 int * indexColumnU = indexColumnU_.array();
861 if (CoinFromFile(indexColumnU,lengthAreaU_ , fp, newSize )==1)
862 return 1;
863 assert (newSize==lengthAreaU_);
864 CoinBigIndex *convertRowToColumnU = convertRowToColumnU_.array();
865 if (CoinFromFile(convertRowToColumnU,lengthAreaU_ , fp, newSize )==1)
866 return 1;
867 assert (newSize==lengthAreaU_||(newSize==0&&!convertRowToColumnU_.array()));
868 CoinFactorizationDouble * elementByRowL = elementByRowL_.array();
869 if (CoinFromFile(elementByRowL,lengthAreaL_ , fp, newSize )==1)
870 return 1;
871 assert (newSize==lengthAreaL_||(newSize==0&&!elementByRowL_.array()));
872 int * indexColumnL = indexColumnL_.array();
873 if (CoinFromFile(indexColumnL,lengthAreaL_ , fp, newSize )==1)
874 return 1;
875 assert (newSize==lengthAreaL_||(newSize==0&&!indexColumnL_.array()));
876 CoinBigIndex * startRowL = startRowL_.array();
877 if (CoinFromFile(startRowL , numberRows_+1, fp, newSize )==1)
878 return 1;
879 assert (newSize==numberRows_+1||(newSize==0&&!startRowL_.array()));
880 CoinFactorizationDouble * elementL = elementL_.array();
881 if (CoinFromFile(elementL,lengthAreaL_ , fp, newSize )==1)
882 return 1;
883 assert (newSize==lengthAreaL_);
884 int * indexRowL = indexRowL_.array();
885 if (CoinFromFile(indexRowL,lengthAreaL_ , fp, newSize )==1)
886 return 1;
887 assert (newSize==lengthAreaL_);
888 CoinBigIndex * startColumnL = startColumnL_.array();
889 if (CoinFromFile(startColumnL,numberRows_ + 1 , fp, newSize )==1)
890 return 1;
891 assert (newSize==numberRows_+1);
892 int * markRow = markRow_.array();
893 if (CoinFromFile(markRow,numberRows_ , fp, newSize )==1)
894 return 1;
895 assert (newSize==numberRows_);
896 int * saveColumn = saveColumn_.array();
897 if (CoinFromFile(saveColumn,numberColumns_ , fp, newSize )==1)
898 return 1;
899 assert (newSize==numberColumns_);
900 CoinBigIndex * startColumnR = startColumnR_.array();
901 if (CoinFromFile(startColumnR , maximumPivots_ + 1 , fp, newSize )==1)
902 return 1;
903 assert (newSize==maximumPivots_+1||(newSize==0&&!startColumnR_.array()));
904 CoinBigIndex * startRowU = startRowU_.array();
905 if (CoinFromFile(startRowU,maximumRowsExtra_ + 1 , fp, newSize )==1)
906 return 1;
907 assert (newSize==maximumRowsExtra_+1||(newSize==0&&!startRowU_.array()));
908 int * numberInRow = numberInRow_.array();
909 if (CoinFromFile(numberInRow,maximumRowsExtra_ + 1 , fp, newSize )==1)
910 return 1;
911 assert (newSize==maximumRowsExtra_+1);
912 int * nextRow = nextRow_.array();
913 if (CoinFromFile(nextRow,maximumRowsExtra_ + 1 , fp, newSize )==1)
914 return 1;
915 assert (newSize==maximumRowsExtra_+1);
916 int * lastRow = lastRow_.array();
917 if (CoinFromFile(lastRow,maximumRowsExtra_ + 1 , fp, newSize )==1)
918 return 1;
919 assert (newSize==maximumRowsExtra_+1);
920 CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
921 if (CoinFromFile(pivotRegion,maximumRowsExtra_ + 1 , fp, newSize )==1)
922 return 1;
923 assert (newSize==maximumRowsExtra_+1);
924 int * permuteBack = permuteBack_.array();
925 if (CoinFromFile(permuteBack,maximumRowsExtra_ + 1 , fp, newSize )==1)
926 return 1;
927 assert (newSize==maximumRowsExtra_+1||(newSize==0&&!permuteBack_.array()));
928 int * permute = permute_.array();
929 if (CoinFromFile(permute,maximumRowsExtra_ + 1 , fp, newSize )==1)
930 return 1;
931 assert (newSize==maximumRowsExtra_+1||(newSize==0&&!permute_.array()));
932 int * pivotColumnBack = pivotColumnBack_.array();
933 if (CoinFromFile(pivotColumnBack,maximumRowsExtra_ + 1 , fp, newSize )==1)
934 return 1;
935 assert (newSize==maximumRowsExtra_+1||(newSize==0&&!pivotColumnBack_.array()));
936 CoinBigIndex * startColumnU = startColumnU_.array();
937 if (CoinFromFile(startColumnU,maximumColumnsExtra_ + 1 , fp, newSize )==1)
938 return 1;
939 assert (newSize==maximumColumnsExtra_+1);
940 int * numberInColumn = numberInColumn_.array();
941 if (CoinFromFile(numberInColumn,maximumColumnsExtra_ + 1 , fp, newSize )==1)
942 return 1;
943 assert (newSize==maximumColumnsExtra_+1);
944 int * numberInColumnPlus = numberInColumnPlus_.array();
945 if (CoinFromFile(numberInColumnPlus,maximumColumnsExtra_ + 1 , fp, newSize )==1)
946 return 1;
947 assert (newSize==maximumColumnsExtra_+1);
948 int *firstCount = firstCount_.array();
949 if (CoinFromFile(firstCount,biggerDimension_ + 2 , fp, newSize )==1)
950 return 1;
951 assert (newSize==biggerDimension_+2);
952 int *nextCount = nextCount_.array();
953 if (CoinFromFile(nextCount,numberRows_ + numberColumns_ , fp, newSize )==1)
954 return 1;
955 assert (newSize==numberRows_+numberColumns_);
956 int *lastCount = lastCount_.array();
957 if (CoinFromFile(lastCount,numberRows_ + numberColumns_ , fp, newSize )==1)
958 return 1;
959 assert (newSize==numberRows_+numberColumns_);
960 int * pivotRowL = pivotRowL_.array();
961 if (CoinFromFile(pivotRowL,numberRows_ + 1 , fp, newSize )==1)
962 return 1;
963 assert (newSize==numberRows_+1);
964 int * pivotColumn = pivotColumn_.array();
965 if (CoinFromFile(pivotColumn,maximumColumnsExtra_ + 1 , fp, newSize )==1)
966 return 1;
967 assert (newSize==maximumColumnsExtra_+1);
968 int * nextColumn = nextColumn_.array();
969 if (CoinFromFile(nextColumn,maximumColumnsExtra_ + 1 , fp, newSize )==1)
970 return 1;
971 assert (newSize==maximumColumnsExtra_+1);
972 int * lastColumn = lastColumn_.array();
973 if (CoinFromFile(lastColumn,maximumColumnsExtra_ + 1 , fp, newSize )==1)
974 return 1;
975 assert (newSize==maximumColumnsExtra_+1);
976 if (CoinFromFile(denseArea_ , numberDense_*numberDense_, fp, newSize )==1)
977 return 1;
978 assert (newSize==numberDense_*numberDense_);
979 if (CoinFromFile(densePermute_ , numberDense_, fp, newSize )==1)
980 return 1;
981 assert (newSize==numberDense_);
982 lengthAreaR_ = space;
983 elementR_ = elementL_.array() + lengthL_;
984 indexRowR_ = indexRowL_.array() + lengthL_;
985 fclose(fp);
986 if (factorIt) {
987 if (biasLU_>=3||numberRows_!=numberColumns_)
988 preProcess ( 2 );
989 else
990 preProcess ( 3 ); // no row copy
991 factor ( );
992 }
993 }
994 return 0;
995 }
996 // factorSparse. Does sparse phase of factorization
997 //return code is <0 error, 0= finished
998 int
factorSparseLarge()999 CoinFactorization::factorSparseLarge ( )
1000 {
1001 int *indexRow = indexRowU_.array();
1002 int *indexColumn = indexColumnU_.array();
1003 CoinFactorizationDouble *element = elementU_.array();
1004 int count = 1;
1005 workArea_.conditionalNew(numberRows_);
1006 CoinFactorizationDouble * workArea = workArea_.array();
1007 #ifndef NDEBUG
1008 counter1++;
1009 #endif
1010 // when to go dense
1011 int denseThreshold=denseThreshold_;
1012
1013 CoinZeroN ( workArea, numberRows_ );
1014 //get space for bit work area
1015 CoinBigIndex workSize = 1000;
1016 workArea2_.conditionalNew(workSize);
1017 unsigned int * workArea2 = workArea2_.array();
1018
1019 //set markRow so no rows updated
1020 int * markRow = markRow_.array();
1021 CoinFillN ( markRow, numberRows_, COIN_INT_MAX-10+1);
1022 int status = 0;
1023 //do slacks first
1024 int * numberInRow = numberInRow_.array();
1025 int * numberInColumn = numberInColumn_.array();
1026 int * numberInColumnPlus = numberInColumnPlus_.array();
1027 CoinBigIndex * startColumnU = startColumnU_.array();
1028 CoinBigIndex * startColumnL = startColumnL_.array();
1029 if (biasLU_<3&&numberColumns_==numberRows_) {
1030 int iPivotColumn;
1031 int * pivotColumn = pivotColumn_.array();
1032 int * nextRow = nextRow_.array();
1033 int * lastRow = lastRow_.array();
1034 for ( iPivotColumn = 0; iPivotColumn < numberColumns_;
1035 iPivotColumn++ ) {
1036 if ( numberInColumn[iPivotColumn] == 1 ) {
1037 CoinBigIndex start = startColumnU[iPivotColumn];
1038 CoinFactorizationDouble value = element[start];
1039 if ( value == slackValue_ && numberInColumnPlus[iPivotColumn] == 0 ) {
1040 // treat as slack
1041 int iRow = indexRow[start];
1042 // but only if row not marked
1043 if (numberInRow[iRow]>0) {
1044 totalElements_ -= numberInRow[iRow];
1045 //take out this bit of indexColumnU
1046 int next = nextRow[iRow];
1047 int last = lastRow[iRow];
1048
1049 nextRow[last] = next;
1050 lastRow[next] = last;
1051 nextRow[iRow] = numberGoodU_; //use for permute
1052 lastRow[iRow] = -2; //mark
1053 //modify linked list for pivots
1054 deleteLink ( iRow );
1055 numberInRow[iRow]=-1;
1056 numberInColumn[iPivotColumn]=0;
1057 numberGoodL_++;
1058 startColumnL[numberGoodL_] = 0;
1059 pivotColumn[numberGoodU_] = iPivotColumn;
1060 numberGoodU_++;
1061 }
1062 }
1063 }
1064 }
1065 // redo
1066 preProcess(4);
1067 CoinFillN ( markRow, numberRows_, COIN_INT_MAX-10+1);
1068 }
1069 numberSlacks_ = numberGoodU_;
1070 int *nextCount = nextCount_.array();
1071 int *firstCount = firstCount_.array();
1072 CoinBigIndex *startRow = startRowU_.array();
1073 CoinBigIndex *startColumn = startColumnU;
1074 //double *elementL = elementL_.array();
1075 //int *indexRowL = indexRowL_.array();
1076 //int *saveColumn = saveColumn_.array();
1077 //int *nextRow = nextRow_.array();
1078 //int *lastRow = lastRow_.array() ;
1079 double pivotTolerance = pivotTolerance_;
1080 int numberTrials = numberTrials_;
1081 int numberRows = numberRows_;
1082 // Put column singletons first - (if false)
1083 separateLinks(1,(biasLU_>1));
1084 #ifndef NDEBUG
1085 int counter2=0;
1086 #endif
1087 while ( count <= biggerDimension_ ) {
1088 #ifndef NDEBUG
1089 counter2++;
1090 int badRow=-1;
1091 if (counter1==-1&&counter2>=0) {
1092 // check counts consistent
1093 for (int iCount=1;iCount<numberRows_;iCount++) {
1094 int look = firstCount[iCount];
1095 while ( look >= 0 ) {
1096 if ( look < numberRows_ ) {
1097 int iRow = look;
1098 if (iRow==badRow)
1099 printf("row count for row %d is %d\n",iCount,iRow);
1100 if ( numberInRow[iRow] != iCount ) {
1101 printf("failed debug on %d entry to factorSparse and %d try\n",
1102 counter1,counter2);
1103 printf("row %d - count %d number %d\n",iRow,iCount,numberInRow[iRow]);
1104 abort();
1105 }
1106 look = nextCount[look];
1107 } else {
1108 int iColumn = look - numberRows;
1109 if ( numberInColumn[iColumn] != iCount ) {
1110 printf("failed debug on %d entry to factorSparse and %d try\n",
1111 counter1,counter2);
1112 printf("column %d - count %d number %d\n",iColumn,iCount,numberInColumn[iColumn]);
1113 abort();
1114 }
1115 look = nextCount[look];
1116 }
1117 }
1118 }
1119 }
1120 #endif
1121 CoinBigIndex minimumCount = COIN_INT_MAX;
1122 double minimumCost = COIN_DBL_MAX;
1123
1124 int iPivotRow = -1;
1125 int iPivotColumn = -1;
1126 int pivotRowPosition = -1;
1127 int pivotColumnPosition = -1;
1128 int look = firstCount[count];
1129 int trials = 0;
1130 int * pivotColumn = pivotColumn_.array();
1131
1132 if ( count == 1 && firstCount[1] >= 0 &&!biasLU_) {
1133 //do column singletons first to put more in U
1134 while ( look >= 0 ) {
1135 if ( look < numberRows_ ) {
1136 look = nextCount[look];
1137 } else {
1138 int iColumn = look - numberRows_;
1139
1140 assert ( numberInColumn[iColumn] == count );
1141 CoinBigIndex start = startColumnU[iColumn];
1142 int iRow = indexRow[start];
1143
1144 iPivotRow = iRow;
1145 pivotRowPosition = start;
1146 iPivotColumn = iColumn;
1147 assert (iPivotRow>=0&&iPivotColumn>=0);
1148 pivotColumnPosition = -1;
1149 look = -1;
1150 break;
1151 }
1152 } /* endwhile */
1153 if ( iPivotRow < 0 ) {
1154 //back to singletons
1155 look = firstCount[1];
1156 }
1157 }
1158 while ( look >= 0 ) {
1159 if ( look < numberRows_ ) {
1160 int iRow = look;
1161 #ifndef NDEBUG
1162 if ( numberInRow[iRow] != count ) {
1163 printf("failed on %d entry to factorSparse and %d try\n",
1164 counter1,counter2);
1165 printf("row %d - count %d number %d\n",iRow,count,numberInRow[iRow]);
1166 abort();
1167 }
1168 #endif
1169 look = nextCount[look];
1170 bool rejected = false;
1171 CoinBigIndex start = startRow[iRow];
1172 CoinBigIndex end = start + count;
1173
1174 CoinBigIndex i;
1175 for ( i = start; i < end; i++ ) {
1176 int iColumn = indexColumn[i];
1177 assert (numberInColumn[iColumn]>0);
1178 double cost = ( count - 1 ) * numberInColumn[iColumn];
1179
1180 if ( cost < minimumCost ) {
1181 CoinBigIndex where = startColumn[iColumn];
1182 CoinFactorizationDouble minimumValue = element[where];
1183
1184 minimumValue = fabs ( minimumValue ) * pivotTolerance;
1185 while ( indexRow[where] != iRow ) {
1186 where++;
1187 } /* endwhile */
1188 assert ( where < startColumn[iColumn] +
1189 numberInColumn[iColumn]);
1190 CoinFactorizationDouble value = element[where];
1191
1192 value = fabs ( value );
1193 if ( value >= minimumValue ) {
1194 minimumCost = cost;
1195 minimumCount = numberInColumn[iColumn];
1196 iPivotRow = iRow;
1197 pivotRowPosition = -1;
1198 iPivotColumn = iColumn;
1199 assert (iPivotRow>=0&&iPivotColumn>=0);
1200 pivotColumnPosition = i;
1201 rejected=false;
1202 if ( minimumCount < count ) {
1203 look = -1;
1204 break;
1205 }
1206 } else if ( iPivotRow == -1 ) {
1207 rejected = true;
1208 }
1209 }
1210 }
1211 trials++;
1212 if ( trials >= numberTrials && iPivotRow >= 0 ) {
1213 look = -1;
1214 break;
1215 }
1216 if ( rejected ) {
1217 //take out for moment
1218 //eligible when row changes
1219 deleteLink ( iRow );
1220 addLink ( iRow, biggerDimension_ + 1 );
1221 }
1222 } else {
1223 int iColumn = look - numberRows;
1224
1225 assert ( numberInColumn[iColumn] == count );
1226 look = nextCount[look];
1227 CoinBigIndex start = startColumn[iColumn];
1228 CoinBigIndex end = start + numberInColumn[iColumn];
1229 CoinFactorizationDouble minimumValue = element[start];
1230
1231 minimumValue = fabs ( minimumValue ) * pivotTolerance;
1232 CoinBigIndex i;
1233 for ( i = start; i < end; i++ ) {
1234 CoinFactorizationDouble value = element[i];
1235
1236 value = fabs ( value );
1237 if ( value >= minimumValue ) {
1238 int iRow = indexRow[i];
1239 int nInRow = numberInRow[iRow];
1240 assert (nInRow>0);
1241 double cost = ( count - 1 ) * nInRow;
1242
1243 if ( cost < minimumCost ) {
1244 minimumCost = cost;
1245 minimumCount = nInRow;
1246 iPivotRow = iRow;
1247 pivotRowPosition = i;
1248 iPivotColumn = iColumn;
1249 assert (iPivotRow>=0&&iPivotColumn>=0);
1250 pivotColumnPosition = -1;
1251 if ( minimumCount <= count + 1 ) {
1252 look = -1;
1253 break;
1254 }
1255 }
1256 }
1257 }
1258 trials++;
1259 if ( trials >= numberTrials && iPivotRow >= 0 ) {
1260 look = -1;
1261 break;
1262 }
1263 }
1264 } /* endwhile */
1265 if (iPivotRow>=0) {
1266 if ( iPivotRow >= 0 ) {
1267 int numberDoRow = numberInRow[iPivotRow] - 1;
1268 int numberDoColumn = numberInColumn[iPivotColumn] - 1;
1269
1270 totalElements_ -= ( numberDoRow + numberDoColumn + 1 );
1271 if ( numberDoColumn > 0 ) {
1272 if ( numberDoRow > 0 ) {
1273 if ( numberDoColumn > 1 ) {
1274 // if (1) {
1275 //need to adjust more for cache and SMP
1276 //allow at least 4 extra
1277 int increment = numberDoColumn + 1 + 4;
1278
1279 if ( increment & 15 ) {
1280 increment = increment & ( ~15 );
1281 increment += 16;
1282 }
1283 int increment2 =
1284
1285 ( increment + COINFACTORIZATION_BITS_PER_INT - 1 ) >> COINFACTORIZATION_SHIFT_PER_INT;
1286 CoinBigIndex size = increment2 * numberDoRow;
1287
1288 if ( size > workSize ) {
1289 workSize = size;
1290 workArea2_.conditionalNew(workSize);
1291 workArea2 = workArea2_.array();
1292 }
1293 bool goodPivot;
1294
1295 //might be able to do better by permuting
1296 #ifndef UGLY_COIN_FACTOR_CODING
1297 //branch out to best pivot routine
1298 goodPivot = pivot ( iPivotRow, iPivotColumn,
1299 pivotRowPosition, pivotColumnPosition,
1300 workArea, workArea2,
1301 increment2, markRow ,
1302 LARGE_SET);
1303 #else
1304 #define FAC_SET LARGE_SET
1305 #include "CoinFactorization.hpp"
1306 #undef FAC_SET
1307 #endif
1308 if ( !goodPivot ) {
1309 status = -99;
1310 count=biggerDimension_+1;
1311 break;
1312 }
1313 } else {
1314 if ( !pivotOneOtherRow ( iPivotRow, iPivotColumn ) ) {
1315 status = -99;
1316 count=biggerDimension_+1;
1317 break;
1318 }
1319 }
1320 } else {
1321 assert (!numberDoRow);
1322 if ( !pivotRowSingleton ( iPivotRow, iPivotColumn ) ) {
1323 status = -99;
1324 count=biggerDimension_+1;
1325 break;
1326 }
1327 }
1328 } else {
1329 assert (!numberDoColumn);
1330 if ( !pivotColumnSingleton ( iPivotRow, iPivotColumn ) ) {
1331 status = -99;
1332 count=biggerDimension_+1;
1333 break;
1334 }
1335 }
1336 assert (nextRow_.array()[iPivotRow]==numberGoodU_);
1337 pivotColumn[numberGoodU_] = iPivotColumn;
1338 numberGoodU_++;
1339 // This should not need to be trapped here - but be safe
1340 if (numberGoodU_==numberRows_)
1341 count=biggerDimension_+1;
1342 }
1343 #if COIN_DEBUG==2
1344 checkConsistency ( );
1345 #endif
1346 #if 0
1347 // Even if no dense code may be better to use naive dense
1348 if (!denseThreshold_&&totalElements_>1000) {
1349 int leftRows=numberRows_-numberGoodU_;
1350 double full = leftRows;
1351 full *= full;
1352 assert (full>=0.0);
1353 double leftElements = totalElements_;
1354 double ratio;
1355 if (leftRows>2000)
1356 ratio=4.0;
1357 else
1358 ratio=3.0;
1359 if (ratio*leftElements>full)
1360 denseThreshold=1;
1361 }
1362 #endif
1363 if (denseThreshold) {
1364 // see whether to go dense
1365 int leftRows=numberRows_-numberGoodU_;
1366 double full = leftRows;
1367 full *= full;
1368 assert (full>=0.0);
1369 double leftElements = totalElements_;
1370 //if (leftRows==100)
1371 //printf("at 100 %d elements\n",totalElements_);
1372 double ratio;
1373 if (leftRows>2000)
1374 ratio=4.0;
1375 else if (leftRows>800)
1376 ratio=3.0;
1377 else if (leftRows>256)
1378 ratio=2.0;
1379 else
1380 ratio=1.5;
1381 if ((ratio*leftElements>full&&leftRows>denseThreshold_)) {
1382 //return to do dense
1383 if (status!=0)
1384 break;
1385 int check=0;
1386 for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
1387 if (numberInColumn[iColumn])
1388 check++;
1389 }
1390 if (check!=leftRows&&denseThreshold_) {
1391 //printf("** mismatch %d columns left, %d rows\n",check,leftRows);
1392 denseThreshold=0;
1393 } else {
1394 status=2;
1395 if ((messageLevel_&4)!=0)
1396 std::cout<<" Went dense at "<<leftRows<<" rows "<<
1397 totalElements_<<" "<<full<<" "<<leftElements<<std::endl;
1398 if (!denseThreshold_)
1399 denseThreshold_=-check; // say how many
1400 break;
1401 }
1402 }
1403 }
1404 // start at 1 again
1405 count = 1;
1406 } else {
1407 //end of this - onto next
1408 count++;
1409 }
1410 } /* endwhile */
1411 workArea_.conditionalDelete() ;
1412 workArea2_.conditionalDelete() ;
1413 return status;
1414 }
1415