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