source: trunk/Clp/src/CoinAbcBaseFactorization2.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 3 months ago

formatting

  • Property svn:keywords set to Id
File size: 47.9 KB
Line 
1/* $Id: CoinAbcBaseFactorization2.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5#ifdef ABC_JUST_ONE_FACTORIZATION
6#include "CoinAbcCommonFactorization.hpp"
7#define CoinAbcTypeFactorization CoinAbcBaseFactorization
8#define ABC_SMALL -2
9#include "CoinAbcBaseFactorization.hpp"
10#endif
11#ifdef CoinAbcTypeFactorization
12
13#include "CoinIndexedVector.hpp"
14#include "CoinHelperFunctions.hpp"
15#include "CoinFinite.hpp"
16//  factorSparse.  Does sparse phase of factorization
17//return code is <0 error, 0= finished
18CoinSimplexInt
19CoinAbcTypeFactorization::factorSparse()
20{
21  CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
22  CoinSimplexInt *COIN_RESTRICT indexColumn = indexColumnUAddress_;
23  CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
24  workArea_.conditionalNew(numberRows_);
25  workAreaAddress_ = workArea_.array();
26  CoinFactorizationDouble *COIN_RESTRICT workArea = workAreaAddress_;
27  double initialLargest = preProcess3();
28  // adjust for test
29  initialLargest = 1.0e-10 / initialLargest;
30#if ABC_NORMAL_DEBUG > 0
31  double largestPivot = 1.0;
32  double smallestPivot = 1.0;
33#endif
34  CoinSimplexInt status = 0;
35  //do slacks first
36  CoinSimplexInt *COIN_RESTRICT numberInRow = numberInRowAddress_;
37  CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
38  //CoinSimplexInt * numberInColumnPlus = numberInColumnPlusAddress_;
39  CoinBigIndex *COIN_RESTRICT startColumnU = startColumnUAddress_;
40  //CoinBigIndex *  COIN_RESTRICT startColumnL = startColumnLAddress_;
41  CoinZeroN(workArea, numberRows_);
42  CoinSimplexInt *COIN_RESTRICT nextCount = this->nextCountAddress_;
43  CoinSimplexInt *COIN_RESTRICT firstCount = this->firstCountAddress_;
44  CoinBigIndex *COIN_RESTRICT startRow = startRowUAddress_;
45  CoinBigIndex *COIN_RESTRICT startColumn = startColumnU;
46#if 0
47        {
48          int checkEls=0;
49          for (CoinSimplexInt iColumn=0;iColumn<numberRows_;iColumn++) {
50            if (numberInColumn[iColumn]) {
51              checkEls+=numberInColumn[iColumn];
52            }
53          }
54          assert (checkEls==totalElements_);
55        }
56#endif
57  // Put column singletons first - (if false)
58  //separateLinks();
59  // while just singletons - do singleton rows first
60  CoinSimplexInt *COIN_RESTRICT pivotColumn = pivotColumnAddress_;
61  // allow a bit more on tolerance
62  double tolerance = 0.5 * pivotTolerance_;
63  CoinSimplexInt *COIN_RESTRICT nextRow = nextRowAddress_;
64  CoinSimplexInt *COIN_RESTRICT lastRow = lastRowAddress_;
65  CoinBigIndex *COIN_RESTRICT startColumnL = startColumnLAddress_;
66  CoinFactorizationDouble *COIN_RESTRICT elementL = elementLAddress_;
67  CoinSimplexInt *COIN_RESTRICT indexRowL = indexRowLAddress_;
68  CoinSimplexInt *COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
69  CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
70  startColumnL[numberGoodL_] = lengthL_; //for luck and first time
71#ifdef ABC_USE_FUNCTION_POINTERS
72  CoinSimplexInt *COIN_RESTRICT nextColumn = nextColumnAddress_;
73  CoinSimplexInt *COIN_RESTRICT lastColumn = lastColumnAddress_;
74  scatterStruct *COIN_RESTRICT scatter = scatterUColumn();
75#if ABC_USE_FUNCTION_POINTERS
76  extern scatterUpdate AbcScatterLowSubtract[9];
77  extern scatterUpdate AbcScatterHighSubtract[4];
78#endif
79  CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
80#endif
81#if CONVERTROW > 1
82  CoinBigIndex *COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
83  CoinBigIndex *COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
84#endif
85  //#define PRINT_INVERT_STATS
86#ifdef PRINT_INVERT_STATS
87  int numberGoodUX = numberGoodU_;
88  int numberTakenOutOfUSpace = 0;
89#endif
90  CoinSimplexInt look = firstCount[1];
91  while (look >= 0) {
92    checkLinks(1);
93    CoinSimplexInt iPivotRow;
94    CoinSimplexInt iPivotColumn = -1;
95#ifdef SMALL_PERMUTE
96    int realPivotRow;
97    int realPivotColumn = -1;
98#endif
99    CoinFactorizationDouble pivotMultiplier = 0.0;
100    if (look < numberRows_) {
101      assert(numberInRow[look] == 1);
102      iPivotRow = look;
103      CoinBigIndex start = startRow[iPivotRow];
104      iPivotColumn = indexColumn[start];
105#ifdef SMALL_PERMUTE
106      realPivotRow = fromSmallToBigRow_[iPivotRow];
107      realPivotColumn = fromSmallToBigColumn_[iPivotColumn];
108      //printf("singleton row realPivotRow %d realPivotColumn %d\n",
109      //     realPivotRow,realPivotColumn);
110#endif
111      assert(numberInColumn[iPivotColumn] > 0);
112      CoinBigIndex startC = startColumn[iPivotColumn];
113      double minimumValue = fabs(element[startC] * tolerance);
114#if CONVERTROW < 2
115      CoinBigIndex where = startC;
116      while (indexRow[where] != iPivotRow) {
117        where++;
118      }
119#else
120      CoinBigIndex where = convertRowToColumn[start];
121#endif
122      assert(where < startC + numberInColumn[iPivotColumn]);
123      CoinFactorizationDouble value = element[where];
124      // size shouldn't matter as would be taken if matrix flipped but ...
125      if (fabs(value) > minimumValue) {
126        CoinSimplexInt numberDoColumn = numberInColumn[iPivotColumn];
127        totalElements_ -= numberDoColumn;
128        // L is always big enough here
129        //store column in L, compress in U and take column out
130        CoinBigIndex l = lengthL_;
131        lengthL_ += numberDoColumn - 1;
132        pivotMultiplier = 1.0 / value;
133        CoinBigIndex endC = startC + numberDoColumn;
134        for (CoinBigIndex i = startC; i < endC; i++) {
135          if (i != where) {
136            CoinSimplexInt iRow = indexRow[i];
137#ifdef SMALL_PERMUTE
138            indexRowL[l] = fromSmallToBigRow_[iRow];
139#else
140            indexRowL[l] = iRow;
141#endif
142            elementL[l] = element[i] * pivotMultiplier;
143            l++;
144            //take out of row list
145            CoinBigIndex startR = startRow[iRow];
146            CoinSimplexInt iNumberInRow = numberInRow[iRow];
147            CoinBigIndex endR = startR + iNumberInRow;
148#if CONVERTROW < 2
149            CoinBigIndex whereRow = startR;
150            while (indexColumn[whereRow] != iPivotColumn)
151              whereRow++;
152            assert(whereRow < endR);
153#else
154            CoinBigIndex whereRow = convertColumnToRow[i];
155#endif
156            int iColumn = indexColumn[endR - 1];
157            indexColumn[whereRow] = iColumn;
158#if CONVERTROW > 1
159            CoinBigIndex whereColumnEntry = convertRowToColumn[endR - 1];
160            convertRowToColumn[whereRow] = whereColumnEntry;
161            convertColumnToRow[whereColumnEntry] = whereRow;
162#endif
163            iNumberInRow--;
164            numberInRow[iRow] = iNumberInRow;
165            modifyLink(iRow, iNumberInRow);
166            checkLinks();
167          }
168        }
169      } else {
170        //printf("Rejecting as element %g largest %g\n",value,element[startC]);
171      }
172    } else {
173      iPivotColumn = look - numberRows_;
174      assert(numberInColumn[iPivotColumn] == 1);
175      {
176        CoinBigIndex startC = startColumn[iPivotColumn];
177        iPivotRow = indexRow[startC];
178#ifdef SMALL_PERMUTE
179        realPivotRow = fromSmallToBigRow_[iPivotRow];
180        realPivotColumn = fromSmallToBigColumn_[iPivotColumn];
181        int realPivotRow;
182        realPivotRow = fromSmallToBigRow_[iPivotRow];
183        indexRow[startC] = realPivotRow;
184        //printf("singleton column realPivotRow %d realPivotColumn %d\n",
185        //   realPivotRow,realPivotColumn);
186#endif
187        //store pivot columns (so can easily compress)
188        assert(numberInColumn[iPivotColumn] == 1);
189        CoinSimplexInt numberDoRow = numberInRow[iPivotRow];
190        int numberDoColumn = numberInColumn[iPivotColumn] - 1;
191        totalElements_ -= (numberDoRow + numberDoColumn);
192        CoinBigIndex startR = startRow[iPivotRow];
193        CoinBigIndex endR = startR + numberDoRow;
194        //clean up counts
195        pivotMultiplier = 1.0 / element[startC];
196        // would I be better off doing other way first??
197        for (CoinBigIndex i = startR; i < endR; i++) {
198          CoinSimplexInt iColumn = indexColumn[i];
199          if (iColumn == iPivotColumn)
200            continue;
201          assert(numberInColumn[iColumn]);
202          CoinSimplexInt number = numberInColumn[iColumn] - 1;
203          //modify linked list
204          modifyLink(iColumn + numberRows_, number);
205          CoinBigIndex start = startColumn[iColumn];
206          //move pivot row element
207          if (number) {
208#if CONVERTROW < 2
209            CoinBigIndex pivot = start;
210            CoinSimplexInt iRow = indexRow[pivot];
211            while (iRow != iPivotRow) {
212              pivot++;
213              iRow = indexRow[pivot];
214            }
215#else
216            CoinBigIndex pivot = convertRowToColumn[i];
217#endif
218            assert(pivot < startColumn[iColumn] + numberInColumn[iColumn]);
219            if (pivot != start) {
220              //move largest one up
221              CoinFactorizationDouble value = element[start];
222              int iRow = indexRow[start];
223              element[start] = element[pivot];
224              indexRow[start] = indexRow[pivot];
225              element[pivot] = element[start + 1];
226              indexRow[pivot] = indexRow[start + 1];
227#if CONVERTROW > 1
228              CoinBigIndex whereRowEntry = convertColumnToRow[start + 1];
229              CoinBigIndex whereRowEntry2 = convertColumnToRow[start];
230              convertRowToColumn[whereRowEntry] = pivot;
231              convertColumnToRow[pivot] = whereRowEntry;
232              convertRowToColumn[whereRowEntry2] = start + 1;
233              convertColumnToRow[start + 1] = whereRowEntry2;
234#endif
235              element[start + 1] = value;
236              indexRow[start + 1] = iRow;
237            } else {
238              //find new largest element
239              CoinSimplexInt iRowSave = indexRow[start + 1];
240              CoinFactorizationDouble valueSave = element[start + 1];
241              CoinFactorizationDouble valueLargest = fabs(valueSave);
242              CoinBigIndex end = start + numberInColumn[iColumn];
243              CoinBigIndex largest = start + 1;
244              for (CoinBigIndex k = start + 2; k < end; k++) {
245                CoinFactorizationDouble value = element[k];
246                CoinFactorizationDouble valueAbs = fabs(value);
247                if (valueAbs > valueLargest) {
248                  valueLargest = valueAbs;
249                  largest = k;
250                }
251              }
252              indexRow[start + 1] = indexRow[largest];
253              element[start + 1] = element[largest];
254#if CONVERTROW > 1
255              CoinBigIndex whereRowEntry = convertColumnToRow[largest];
256              CoinBigIndex whereRowEntry2 = convertColumnToRow[start + 1];
257              convertRowToColumn[whereRowEntry] = start + 1;
258              convertColumnToRow[start + 1] = whereRowEntry;
259              convertRowToColumn[whereRowEntry2] = largest;
260              convertColumnToRow[largest] = whereRowEntry2;
261#endif
262              indexRow[largest] = iRowSave;
263              element[largest] = valueSave;
264            }
265          }
266          //clean up counts
267          numberInColumn[iColumn]--;
268          numberInColumnPlus[iColumn]++;
269#ifdef SMALL_PERMUTE
270          indexRow[start] = realPivotRow;
271#endif
272          startColumn[iColumn]++;
273        }
274      }
275    }
276    if (pivotMultiplier) {
277      numberGoodL_++;
278      startColumnL[numberGoodL_] = lengthL_;
279      //take out this bit of indexColumnU
280      CoinSimplexInt next = nextRow[iPivotRow];
281      CoinSimplexInt last = lastRow[iPivotRow];
282      nextRow[last] = next;
283      lastRow[next] = last;
284      lastRow[iPivotRow] = -2; //mark
285#ifdef SMALL_PERMUTE
286      // mark
287      fromSmallToBigRow_[iPivotRow] = -1;
288      fromSmallToBigColumn_[iPivotColumn] = -1;
289      permuteAddress_[realPivotRow] = numberGoodU_;
290#else
291      permuteAddress_[iPivotRow] = numberGoodU_;
292#endif
293#ifdef ABC_USE_FUNCTION_POINTERS
294      int number = numberInColumnPlus[iPivotColumn];
295#ifdef SMALL_PERMUTE
296#if ABC_USE_FUNCTION_POINTERS
297      if (number < 9) {
298        scatter[realPivotRow].functionPointer = AbcScatterLowSubtract[number];
299      } else {
300        scatter[realPivotRow].functionPointer = AbcScatterHighSubtract[number & 3];
301      }
302#endif
303      scatter[realPivotRow].offset = lastEntryByColumnUPlus_;
304      scatter[realPivotRow].number = number;
305#else
306#if ABC_USE_FUNCTION_POINTERS
307      if (number < 9) {
308        scatter[iPivotRow].functionPointer = AbcScatterLowSubtract[number];
309      } else {
310        scatter[iPivotRow].functionPointer = AbcScatterHighSubtract[number & 3];
311      }
312#endif
313      scatter[iPivotRow].offset = lastEntryByColumnUPlus_;
314      scatter[iPivotRow].number = number;
315#endif
316      CoinBigIndex startC = startColumn[iPivotColumn] - number;
317      for (int i = startC; i < startC + number; i++)
318        elementUColumnPlus[lastEntryByColumnUPlus_++] = element[i] * pivotMultiplier;
319      CoinAbcMemcpy(reinterpret_cast< CoinSimplexInt * >(elementUColumnPlus + lastEntryByColumnUPlus_), indexRow + startC, number);
320      lastEntryByColumnUPlus_ += (number + 1) >> 1;
321#endif
322      numberInColumn[iPivotColumn] = 0;
323      numberInRow[iPivotRow] = 0;
324      //modify linked list for pivots
325      deleteLink(iPivotRow);
326      deleteLink(iPivotColumn + numberRows_);
327      checkLinks();
328#ifdef SMALL_PERMUTE
329      assert(permuteAddress_[realPivotRow] == numberGoodU_);
330      pivotColumn[numberGoodU_] = realPivotColumn;
331#else
332      assert(permuteAddress_[iPivotRow] == numberGoodU_);
333      pivotColumn[numberGoodU_] = iPivotColumn;
334#endif
335      pivotRegion[numberGoodU_] = pivotMultiplier;
336      numberGoodU_++;
337    } else {
338      //take out for moment
339      modifyLink(iPivotRow, numberRows_ + 1);
340    }
341    look = -1; //nextCount[look];
342    if (look < 0) {
343      // start again
344      look = firstCount[1];
345    }
346    assert(iPivotColumn >= 0);
347#ifdef ABC_USE_FUNCTION_POINTERS
348    if (!numberInColumn[iPivotColumn]) {
349      int iNext = nextColumn[iPivotColumn];
350      int iLast = lastColumn[iPivotColumn];
351      assert(iLast >= 0);
352#ifdef PRINT_INVERT_STATS
353      numberTakenOutOfUSpace++;
354#endif
355      lastColumn[iNext] = iLast;
356      nextColumn[iLast] = iNext;
357    }
358#endif
359  }
360  // put back all rejected
361  look = firstCount[numberRows_ + 1];
362  while (look >= 0) {
363#ifndef NDEBUG
364    if (look < numberRows_) {
365      assert(numberInRow[look] == 1);
366    } else {
367      int iColumn = look - numberRows_;
368      assert(numberInColumn[iColumn] == 1);
369    }
370#endif
371    int nextLook = nextCount[look];
372    modifyLink(look, 1);
373    look = nextLook;
374  }
375#ifdef SMALL_PERMUTE
376  if (numberGoodU_ < numberRows_ && numberGoodU_ > numberSlacks_ + (numberRows_ >> 3)) {
377    CoinSimplexInt *COIN_RESTRICT fromBigToSmallRow = reinterpret_cast< CoinSimplexInt * >(workArea_.array());
378    CoinSimplexInt *COIN_RESTRICT fromBigToSmallColumn = fromBigToSmallRow + numberRows_;
379    int nSmall = 0;
380    for (int i = 0; i < numberRowsSmall_; i++) {
381      int realRow = fromSmallToBigRow_[i];
382      if (realRow < 0) {
383        fromBigToSmallRow[i] = -1;
384      } else {
385        // in
386        fromBigToSmallRow[i] = nSmall;
387        numberInRow[nSmall] = numberInRow[i];
388        startRow[nSmall] = startRow[i];
389        fromSmallToBigRow_[nSmall++] = realRow;
390      }
391    }
392    nSmall = 0;
393    for (int i = 0; i < numberRowsSmall_; i++) {
394      int realColumn = fromSmallToBigColumn_[i];
395      if (realColumn < 0) {
396        fromBigToSmallColumn[i] = -1;
397      } else {
398        // in
399        fromBigToSmallColumn[i] = nSmall;
400        numberInColumn[nSmall] = numberInColumn[i];
401        numberInColumnPlus[nSmall] = numberInColumnPlus[i];
402        startColumn[nSmall] = startColumn[i];
403        fromSmallToBigColumn_[nSmall++] = realColumn;
404      }
405    }
406    CoinAbcMemset0(numberInRow + nSmall, numberRowsSmall_ - nSmall);
407    CoinAbcMemset0(numberInColumn + nSmall, numberRowsSmall_ - nSmall);
408    // indices
409    for (int i = 0; i < numberRowsSmall_; i++) {
410      CoinBigIndex start = startRow[i];
411      CoinBigIndex end = start + numberInRow[i];
412      for (CoinBigIndex j = start; j < end; j++) {
413        indexColumn[j] = fromBigToSmallColumn[indexColumn[j]];
414      }
415    }
416    for (int i = 0; i < numberRowsSmall_; i++) {
417      CoinBigIndex start = startColumn[i];
418      CoinBigIndex end = start + numberInColumn[i];
419      for (CoinBigIndex j = start; j < end; j++) {
420        indexRow[j] = fromBigToSmallRow[indexRow[j]];
421      }
422    }
423    // find area somewhere
424    int *temp = fromSmallToBigColumn_ + nSmall;
425    CoinSimplexInt *nextFake = temp;
426    //CoinSimplexInt lastFake=temp+numberRows_;
427    CoinAbcMemcpy(nextFake, nextRow, numberRowsSmall_);
428    nextFake[numberRows_] = nextRow[numberRows_];
429    //CoinAbcMemcpy(lastFake,lastRow,numberRowsSmall_);
430    // nextRow has end at numberRows_
431    CoinSimplexInt lastBigRow = numberRows_;
432    CoinSimplexInt lastSmallRow = numberRows_;
433    for (CoinSimplexInt i = 0; i < nSmall; i++) {
434      CoinSimplexInt bigRow = nextFake[lastBigRow];
435      assert(bigRow >= 0 && bigRow < numberRowsSmall_);
436      CoinSimplexInt smallRow = fromBigToSmallRow[bigRow];
437      nextRow[lastSmallRow] = smallRow;
438      lastRow[smallRow] = lastSmallRow;
439      lastBigRow = bigRow;
440      lastSmallRow = smallRow;
441    }
442    assert(nextFake[lastBigRow] == numberRows_);
443    nextRow[lastSmallRow] = numberRows_;
444    lastRow[numberRows_] = lastSmallRow;
445    // nextColumn has end at maximumRowsExtra_
446    CoinAbcMemcpy(nextFake, nextColumn, numberRowsSmall_);
447    nextFake[maximumRowsExtra_] = nextColumn[maximumRowsExtra_];
448    CoinSimplexInt lastBigColumn = maximumRowsExtra_;
449    CoinSimplexInt lastSmallColumn = maximumRowsExtra_;
450    for (CoinSimplexInt i = 0; i < nSmall; i++) {
451      CoinSimplexInt bigColumn = nextFake[lastBigColumn];
452      assert(bigColumn >= 0 && bigColumn < numberRowsSmall_);
453      CoinSimplexInt smallColumn = fromBigToSmallColumn[bigColumn];
454      nextColumn[lastSmallColumn] = smallColumn;
455      lastColumn[smallColumn] = lastSmallColumn;
456      lastBigColumn = bigColumn;
457      lastSmallColumn = smallColumn;
458    }
459    assert(nextFake[lastBigColumn] == maximumRowsExtra_);
460    nextColumn[lastSmallColumn] = maximumRowsExtra_;
461    lastColumn[maximumRowsExtra_] = lastSmallColumn;
462    // now equal counts (could redo - but for now get exact copy)
463    CoinSimplexInt *COIN_RESTRICT lastCount = this->lastCountAddress_;
464    CoinAbcMemcpy(temp, nextCount, numberRows_ + numberRowsSmall_);
465    for (int iCount = 0; iCount <= nSmall; iCount++) {
466      CoinSimplexInt nextBig = firstCount[iCount];
467      if (nextBig >= 0) {
468        //CoinSimplexInt next=firstCount[iCount];
469        CoinSimplexInt nextSmall;
470        if (nextBig < numberRows_)
471          nextSmall = fromBigToSmallRow[nextBig];
472        else
473          nextSmall = fromBigToSmallColumn[nextBig - numberRows_] + numberRows_;
474        firstCount[iCount] = nextSmall;
475        CoinSimplexInt *where = &lastCount[nextSmall];
476        while (nextBig >= 0) {
477          CoinSimplexInt thisSmall = nextSmall;
478          nextBig = temp[nextBig];
479          if (nextBig >= 0) {
480            if (nextBig < numberRows_)
481              nextSmall = fromBigToSmallRow[nextBig];
482            else
483              nextSmall = fromBigToSmallColumn[nextBig - numberRows_] + numberRows_;
484            lastCount[nextSmall] = thisSmall;
485            nextCount[thisSmall] = nextSmall;
486          } else {
487            nextCount[thisSmall] = nextBig;
488          }
489        }
490        assert(nextSmall >= 0);
491        // fill in odd one
492        *where = iCount - numberRows_ - 2;
493      }
494    }
495    //printf("%d rows small1 %d small2 %d\n",numberRows_,numberRowsSmall_,nSmall);
496    numberRowsSmall_ = nSmall;
497    //exit(0);
498  }
499#endif
500  // Put .hpp functions in separate file for profiling
501#ifdef PRINT_INVERT_STATS
502  int numberGoodUY = numberGoodU_;
503  int numberElsLeftX = 0;
504  for (int i = 0; i < numberRows_; i++)
505    numberElsLeftX += numberInColumn[i];
506  int saveN1X = totalElements_;
507#endif
508#if ABC_DENSE_CODE
509  // when to go dense
510  //CoinSimplexInt denseThreshold=abs(denseThreshold_);
511#endif
512  //get space for bit work area
513  CoinBigIndex workSize = 1000;
514  workArea2_.conditionalNew(workSize);
515  workArea2Address_ = workArea2_.array();
516  CoinSimplexUnsignedInt *COIN_RESTRICT workArea2 = workArea2Address_;
517
518  //set markRow so no rows updated
519  //set markRow so no rows updated
520  CoinSimplexInt *COIN_RESTRICT markRow = markRow_.array();
521  CoinFillN(markRow, numberRowsSmall_, LARGE_UNSET);
522  CoinZeroN(workArea, numberRowsSmall_);
523  CoinFactorizationDouble pivotTolerance = pivotTolerance_;
524  CoinSimplexInt numberTrials = numberTrials_;
525  CoinSimplexInt numberRows = numberRows_;
526  // while just singletons - do singleton rows first
527  CoinSimplexInt count = 1;
528  startColumnL[numberGoodL_] = lengthL_; //for luck and first time
529  numberRowsLeft_ = numberRows_ - numberGoodU_;
530  while (count <= numberRowsLeft_) {
531    CoinBigIndex minimumCount = COIN_INT_MAX;
532    CoinFactorizationDouble minimumCost = COIN_DBL_MAX;
533
534    CoinSimplexInt iPivotRow = -1;
535    CoinSimplexInt iPivotColumn = -1;
536    CoinBigIndex pivotRowPosition = -1;
537    CoinBigIndex pivotColumnPosition = -1;
538    CoinSimplexInt look = firstCount[count];
539#if 0
540    if (numberRowsSmall_==2744&&!numberInColumn[1919]) {
541      int look2=look;
542      while (look2>=0) {
543        if (look2==numberRows_+1919) {
544          printf("*** 1919 on list of count %d\n",count);
545          abort();
546        }
547        look2 = nextCount[look2];
548      }
549    }
550#endif
551    //printf("At count %d look %d\n",count,look);
552    CoinSimplexInt trials = 0;
553    //CoinSimplexInt *  COIN_RESTRICT pivotColumn = pivotColumnAddress_;
554    while (look >= 0) {
555      checkLinks(1);
556      if (look < numberRows_) {
557        CoinSimplexInt iRow = look;
558        look = nextCount[look];
559        bool rejected = false;
560        CoinBigIndex start = startRow[iRow];
561        CoinBigIndex end = start + count;
562
563        CoinBigIndex i;
564        for (i = start; i < end; i++) {
565          CoinSimplexInt iColumn = indexColumn[i];
566          assert(numberInColumn[iColumn] > 0);
567          CoinFactorizationDouble cost = (count - 1) * numberInColumn[iColumn] + 0.1;
568
569          if (cost < minimumCost) {
570            CoinBigIndex where = startColumn[iColumn];
571            double minimumValue = element[where];
572
573            minimumValue = fabs(minimumValue) * pivotTolerance;
574            if (count == 1)
575              minimumValue = CoinMin(minimumValue, 0.999999);
576            while (indexRow[where] != iRow) {
577              where++;
578            } /* endwhile */
579            assert(where < startColumn[iColumn] + numberInColumn[iColumn]);
580            CoinFactorizationDouble value = element[where];
581
582            value = fabs(value);
583            if (value >= minimumValue) {
584              minimumCost = cost;
585              minimumCount = numberInColumn[iColumn];
586              iPivotRow = iRow;
587              pivotRowPosition = -1;
588              iPivotColumn = iColumn;
589              assert(iPivotRow >= 0 && iPivotColumn >= 0);
590              pivotColumnPosition = i;
591              rejected = false;
592            } else if (iPivotRow == -1 && count > 1) {
593              rejected = true;
594            }
595          }
596        }
597        trials++;
598        if (iPivotRow >= 0 && (trials >= numberTrials || minimumCount < count)) {
599          look = -1;
600          break;
601        }
602        if (rejected) {
603          //take out for moment
604          //eligible when row changes
605          modifyLink(iRow, numberRows_ + 1);
606        }
607      } else {
608        CoinSimplexInt iColumn = look - numberRows;
609
610#if 0
611        if ( numberInColumn[iColumn] != count ) {
612          printf("column %d has %d elements but in count list of %d\n",
613                 iColumn,numberInColumn[iColumn],count);
614          abort();
615        }
616#endif
617        assert(numberInColumn[iColumn] == count);
618        look = nextCount[look];
619        CoinBigIndex start = startColumn[iColumn];
620        CoinBigIndex end = start + numberInColumn[iColumn];
621        CoinFactorizationDouble minimumValue = element[start];
622
623        minimumValue = fabs(minimumValue) * pivotTolerance;
624        CoinBigIndex i;
625        for (i = start; i < end; i++) {
626          CoinFactorizationDouble value = element[i];
627
628          value = fabs(value);
629          if (value >= minimumValue) {
630            CoinSimplexInt iRow = indexRow[i];
631            CoinSimplexInt nInRow = numberInRow[iRow];
632            assert(nInRow > 0);
633            CoinFactorizationDouble cost = (count - 1) * nInRow;
634
635            if (cost < minimumCost) {
636              minimumCost = cost;
637              minimumCount = nInRow;
638              iPivotRow = iRow;
639              pivotRowPosition = i;
640              iPivotColumn = iColumn;
641              assert(iPivotRow >= 0 && iPivotColumn >= 0);
642              pivotColumnPosition = -1;
643            }
644          }
645        }
646        trials++;
647        if (iPivotRow >= 0 && (trials >= numberTrials || minimumCount < count)) {
648          look = -1;
649          break;
650        }
651      }
652    } /* endwhile */
653    if (iPivotRow >= 0) {
654      assert(iPivotRow < numberRows_);
655      CoinSimplexInt numberDoRow = numberInRow[iPivotRow] - 1;
656      CoinSimplexInt numberDoColumn = numberInColumn[iPivotColumn] - 1;
657
658      //printf("realPivotRow %d (%d elements) realPivotColumn %d (%d elements)\n",
659      //   fromSmallToBigRow[iPivotRow],numberDoRow+1,fromSmallToBigColumn[iPivotColumn],numberDoColumn+1);
660      //      printf("nGoodU %d pivRow %d (num %d) pivCol %d (num %d) - %d elements left\n",
661      //   numberGoodU_,iPivotRow,numberDoRow,iPivotColumn,numberDoColumn,totalElements_);
662#if 0 //ndef NDEBUG
663      static int testXXXX=0;
664      testXXXX++;
665      if ((testXXXX%100)==0)
666        printf("check consistent totalElements_\n");
667      int nn=0;
668      for (int i=0;i<numberRowsSmall_;i++) {
669        int n=numberInColumn[i];
670        if (n) {
671          int start=startColumn[i];
672          double largest=fabs(element[start]);
673          for (int j=1;j<n;j++) {
674            assert (fabs(element[start+j])<=largest);
675          }
676          nn+=n;
677        }
678      }
679      assert (nn==totalElements_);
680#endif
681      totalElements_ -= (numberDoRow + numberDoColumn + 1);
682#if 0
683      if (numberInColumn[5887]==0&&numberRowsSmall_>5887) {
684        int start=startRow[1181];
685        int end=start+numberInRow[1181];
686        for (int i=start;i<end;i++)
687          assert(indexColumn[i]!=5887);
688      }
689#endif
690      if (numberDoColumn > 0) {
691        if (numberDoRow > 0) {
692          if (numberDoColumn > 1) {
693            //  if (1) {
694            //need to adjust more for cache and SMP
695            //allow at least 4 extra
696            CoinSimplexInt increment = numberDoColumn + 1 + 4;
697
698            if (increment & 15) {
699              increment = increment & (~15);
700              increment += 16;
701            }
702            CoinSimplexInt increment2 =
703
704              (increment + COINFACTORIZATION_BITS_PER_INT - 1) >> COINFACTORIZATION_SHIFT_PER_INT;
705            CoinBigIndex size = increment2 * numberDoRow;
706
707            if (size > workSize) {
708              workSize = size;
709              workArea2_.conditionalNew(workSize);
710              workArea2Address_ = workArea2_.array();
711              workArea2 = workArea2Address_;
712            }
713            //branch out to best pivot routine
714#define ABC_PARALLEL_PIVOT
715#ifndef ABC_USE_FUNCTION_POINTERS
716#undef ABC_PARALLEL_PIVOT
717#endif
718#ifdef ABC_PARALLEL_PIVOT
719            //if (numberRowsSmall_==7202&&numberGoodU_==15609) {
720            //printf("NumberGoodU %d\n",numberGoodU_);
721            //}
722            if (numberDoRow < 20)
723              status = pivot(iPivotRow, iPivotColumn,
724                pivotRowPosition, pivotColumnPosition,
725                workArea, workArea2,
726                increment2, markRow);
727            else
728              status = pivot(iPivotRow, iPivotColumn,
729                pivotRowPosition, pivotColumnPosition,
730                markRow);
731#else
732            status = pivot(iPivotRow, iPivotColumn,
733              pivotRowPosition, pivotColumnPosition,
734              workArea, workArea2,
735              increment2, markRow);
736#endif
737#if 0
738            for (int i=0;i<numberRowsSmall_;i++)
739              assert (markRow[i]==LARGE_UNSET);
740#endif
741              //#define CHECK_SIZE
742#ifdef CHECK_SIZE
743            {
744              double largestU = 0.0;
745              int iU = -1;
746              int jU = -1;
747              for (int i = 0; i < numberRowsSmall_; i++) {
748                CoinBigIndex start = startColumn[i];
749                CoinBigIndex end = start + numberInColumn[i];
750                for (int j = start; j < end; j++) {
751                  if (fabs(element[j]) > largestU) {
752                    iU = i;
753                    jU = j;
754                    largestU = fabs(element[j]);
755                  }
756                }
757              }
758              printf("%d largest el %g at i=%d,j=%d start %d end %d count %d\n", numberGoodU_, element[jU],
759                iU, jU, startColumn[iU], startColumn[iU] + numberInColumn[iU], count);
760            }
761#endif
762#undef CHECK_SIZE
763            checkLinks();
764            if (status < 0) {
765              count = numberRows_ + 1;
766              break;
767            }
768          } else {
769            if (!pivotOneOtherRow(iPivotRow, iPivotColumn)) {
770              status = -99;
771              count = numberRows_ + 1;
772              break;
773            }
774#ifdef CHECK_SIZE
775            {
776              double largestU = 0.0;
777              int iU = -1;
778              int jU = -1;
779              for (int i = 0; i < numberRowsSmall_; i++) {
780                CoinBigIndex start = startColumn[i];
781                CoinBigIndex end = start + numberInColumn[i];
782                for (int j = start; j < end; j++) {
783                  if (fabs(element[j]) > largestU) {
784                    iU = i;
785                    jU = j;
786                    largestU = fabs(element[j]);
787                  }
788                }
789              }
790              printf("B%d largest el %g at i=%d,j=%d\n", numberGoodU_, element[jU],
791                iU, jU);
792            }
793#endif
794            checkLinks();
795          }
796        } else {
797          assert(!numberDoRow);
798          checkLinks();
799          if (!pivotRowSingleton(iPivotRow, iPivotColumn)) {
800            status = -99;
801            count = numberRows_ + 1;
802            break;
803          }
804#ifdef CHECK_SIZE
805          {
806            double largestU = 0.0;
807            int iU = -1;
808            int jU = -1;
809            for (int i = 0; i < numberRowsSmall_; i++) {
810              CoinBigIndex start = startColumn[i];
811              CoinBigIndex end = start + numberInColumn[i];
812              for (int j = start; j < end; j++) {
813                if (fabs(element[j]) > largestU) {
814                  iU = i;
815                  jU = j;
816                  largestU = fabs(element[j]);
817                }
818              }
819            }
820            printf("C%d largest el %g at i=%d,j=%d\n", numberGoodU_, element[jU],
821              iU, jU);
822          }
823#endif
824          checkLinks();
825        }
826      } else {
827        assert(!numberDoColumn);
828        pivotColumnSingleton(iPivotRow, iPivotColumn);
829#ifdef CHECK_SIZE
830        {
831          double largestU = 0.0;
832          int iU = -1;
833          int jU = -1;
834          for (int i = 0; i < numberRowsSmall_; i++) {
835            CoinBigIndex start = startColumn[i];
836            CoinBigIndex end = start + numberInColumn[i];
837            for (int j = start; j < end; j++) {
838              if (fabs(element[j]) > largestU) {
839                iU = i;
840                jU = j;
841                largestU = fabs(element[j]);
842              }
843            }
844          }
845          printf("D%d largest el %g at i=%d,j=%d\n", numberGoodU_, element[jU],
846            iU, jU);
847        }
848#endif
849        checkLinks();
850      }
851      double pivotValue = fabs(pivotRegion[numberGoodU_]);
852#if ABC_NORMAL_DEBUG > 0
853      largestPivot = CoinMax(largestPivot, pivotValue);
854      smallestPivot = CoinMin(smallestPivot, pivotValue);
855#endif
856      if (pivotValue < initialLargest) {
857        if (pivotTolerance_ < 0.95) {
858          pivotTolerance_ = CoinMin(1.1 * pivotTolerance_, 0.99);
859#if ABC_NORMAL_DEBUG > 0
860          printf("PPPivot value of %g increasing pivot tolerance to %g\n",
861            pivotValue, pivotTolerance_);
862#endif
863          status = -97;
864          break;
865        }
866      }
867      afterPivot(iPivotRow, iPivotColumn);
868      assert(numberGoodU_ <= numberRows_);
869      assert(startRow[numberRows_] == lengthAreaU_);
870#if 0
871      static int ixxxxxx=0;
872      {
873        int nLeft=0;
874        for (int i=0;i<numberRows_;i++) {
875          if (permuteAddress_[i]<0) {
876            nLeft++;
877          }
878        }
879        assert (nLeft==numberRowsLeft_);
880        if (nLeft!=numberRowsLeft_) {
881          printf("mismatch nleft %d numberRowsleft_ %d\n",nLeft,numberRowsLeft_);
882          abort();
883        }
884      }
885      if (numberRowsSmall_==-262/*2791*/) {
886        bool bad=false;
887#if 0
888        int cols[2800];
889        for (int i=0;i<2791;i++)
890          cols[i]=-1;
891        int n=numberInRow[1714];
892        int start = startRow[1714];
893        int end=start+n;
894        for (int i=start;i<end;i++) {
895          int iColumn=indexColumn[i];
896          if (iColumn<0||iColumn>=2744) {
897            printf("Bad column %d at %d\n",iColumn,i);
898            bad=true;
899          } else if (cols[iColumn]>=0) {
900            printf("Duplicate column %d at %d was at %d\n",iColumn,i,cols[iColumn]);
901            bad=true;
902          } else {
903            cols[iColumn]=i;
904          }
905        }
906#else
907        {
908          int n=numberInRow[1160];
909          int start = startRow[1160];
910          int end=start+n;
911          for (int i=start;i<end;i++) {
912            int iColumn=indexColumn[i];
913            if (iColumn==2111&&!numberInColumn[2111]) {
914              printf("Bad column %d at %d\n",iColumn,i);
915              abort();
916            }
917          }
918        }
919#endif
920        ixxxxxx++;
921        printf("pivrow %d pivcol %d count in 2111 %d in row 2111 %d xxx %d ngood %d\n",iPivotRow,
922               iPivotColumn,numberInColumn[2111],numberInRow[2111],ixxxxxx,numberGoodU_);
923        //printf("pivrow %d pivcol %d numberGoodU_ %d xxx %d\n",iPivotRow,
924        //     iPivotColumn,numberGoodU_,ixxxxxx);
925        if (bad)
926          abort();
927#if 0
928        if (numberInRow[1714]>numberRows_-numberGoodU_) {
929          printf("numberGoodU %d nrow-ngood %d\n",numberGoodU_,numberRows_-numberGoodU_);
930          abort();
931        }
932        if (ixxxxxx==2198||ixxxxxx==1347) {
933          printf("Trouble ahead\n");
934        }
935#else
936        if (ixxxxxx==1756) {
937          printf("Trouble ahead\n");
938        }
939#endif
940      }
941#endif
942#if ABC_DENSE_CODE
943      if (!status) {
944        status = wantToGoDense();
945      }
946#endif
947      if (status)
948        break;
949      // start at 1 again
950      count = 1;
951    } else {
952      //end of this - onto next
953      count++;
954    }
955  } /* endwhile */
956#if ABC_NORMAL_DEBUG > 0
957#if ABC_SMALL < 2
958  int lenU = 2 * (lastEntryByColumnUPlus_ / 3);
959  printf("largest initial element %g, smallestPivot %g largest %g (%d dense rows) - %d in L, %d in U\n",
960    1.0e-10 / initialLargest, smallestPivot, largestPivot, numberRows_ - numberGoodU_,
961    lengthL_, lenU);
962#endif
963#endif
964  workArea2_.conditionalDelete();
965  workArea2Address_ = NULL;
966#ifdef PRINT_INVERT_STATS
967  int numberDense = numberRows_ - numberGoodU_;
968  printf("XX %d rows, %d in bump (%d slacks,%d triangular), % d dense - startels %d (%d) now %d - taken out of uspace (triang) %d\n",
969    numberRows_, numberRows_ - numberGoodUY, numberGoodUX, numberGoodUY - numberGoodUX,
970    numberDense, numberElsLeftX, saveN1X, totalElements_, numberTakenOutOfUSpace);
971#endif
972  return status;
973}
974#if ABC_DENSE_CODE
975//:method factorDense.  Does dense phase of factorization
976//return code is <0 error, 0= finished
977CoinSimplexInt CoinAbcTypeFactorization::factorDense()
978{
979  CoinSimplexInt status = 0;
980  numberDense_ = numberRows_ - numberGoodU_;
981  if (sizeof(CoinBigIndex) == 4 && numberDense_ >= (2 << 15)) {
982    abort();
983  }
984  CoinBigIndex full;
985  full = numberDense_ * numberDense_;
986  totalElements_ = full;
987  // Add coding to align an object
988#if ABC_DENSE_CODE == 1
989  leadingDimension_ = ((numberDense_ >> 4) + 1) << 4;
990#else
991  leadingDimension_ = numberDense_;
992#endif
993  CoinBigIndex newSize = (leadingDimension_ + FACTOR_CPU) * numberDense_;
994  newSize += (numberDense_ + 1) / (sizeof(CoinFactorizationDouble) / sizeof(CoinSimplexInt));
995#if 1
996  newSize += 2 * ((numberDense_ + 3) / (sizeof(CoinFactorizationDouble) / sizeof(short)));
997  newSize += ((numberRows_ + 3) / (sizeof(CoinFactorizationDouble) / sizeof(short)));
998  // so we can align on 256 byte
999  newSize += 32;
1000  //newSize += (numberDense_+1)/(sizeof(CoinFactorizationDouble)/sizeof(CoinSimplexInt));
1001#endif
1002#ifdef SMALL_PERMUTE
1003  // densePermute has fromSmallToBig!!!
1004  CoinSimplexInt *COIN_RESTRICT fromSmallToBigRow = reinterpret_cast< CoinSimplexInt * >(workArea_.array());
1005  CoinSimplexInt *COIN_RESTRICT fromSmallToBigColumn = fromSmallToBigRow + numberRowsSmall_;
1006  CoinAbcMemcpy(fromSmallToBigRow, fromSmallToBigRow_, numberRowsSmall_);
1007  CoinAbcMemcpy(fromSmallToBigColumn, fromSmallToBigColumn_, numberRowsSmall_);
1008#endif
1009  denseArea_.conditionalDelete();
1010  denseArea_.conditionalNew(newSize);
1011  denseAreaAddress_ = denseArea_.array();
1012  CoinInt64 xx = reinterpret_cast< CoinInt64 >(denseAreaAddress_);
1013  int iBottom = static_cast< int >(xx & 63);
1014  int offset = (256 - iBottom) >> 3;
1015  denseAreaAddress_ += offset;
1016  CoinFactorizationDouble *COIN_RESTRICT denseArea = denseAreaAddress_;
1017  CoinZeroN(denseArea, newSize - 32);
1018  CoinSimplexInt *COIN_RESTRICT densePermute = reinterpret_cast< CoinSimplexInt * >(denseArea + (leadingDimension_ + FACTOR_CPU) * numberDense_);
1019#if ABC_DENSE_CODE < 0
1020  CoinSimplexInt *COIN_RESTRICT indexRowU = indexRowUAddress_;
1021  CoinSimplexInt *COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
1022#endif
1023  //mark row lookup using lastRow
1024  CoinSimplexInt i;
1025  CoinSimplexInt *COIN_RESTRICT lastRow = lastRowAddress_;
1026  CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
1027#if 0
1028  char xxx[17000];
1029  assert (numberRows_<=17000);
1030  memset(xxx,0,numberRows_);
1031  int nLeft=0;
1032  for (i=0;i<numberRows_;i++) {
1033    if (permuteAddress_[i]>=0) {
1034      xxx[i]=1;
1035    } else {
1036      nLeft++;
1037    }
1038  }
1039  printf("%d rows left\n",nLeft);
1040  bool bad=false;
1041  for (i=0;i<numberRowsSmall_;i++) {
1042    int iBigRow=fromSmallToBigRow[i];
1043    if (!xxx[iBigRow])
1044      printf("Row %d (big %d) in dense\n",i,iBigRow);
1045    if(iBigRow<0) abort();
1046    if (lastRow[i]>=0) {
1047      lastRow[i]=0;
1048      printf("dense row %d translates to %d permute %d\n",i,iBigRow,permuteAddress_[iBigRow]);
1049    } else {
1050      if (permuteAddress_[iBigRow]<0)
1051      printf("sparse row %d translates to %d permute %d\n",i,iBigRow,permuteAddress_[iBigRow]);
1052      if (xxx[iBigRow]!=1) bad=true;
1053      xxx[iBigRow]=0;
1054    }
1055  }
1056  if (bad)
1057    abort();
1058#else
1059  for (i = 0; i < numberRowsSmall_; i++) {
1060    if (lastRow[i] >= 0) {
1061      lastRow[i] = 0;
1062    }
1063  }
1064#endif
1065  CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
1066  CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
1067  CoinSimplexInt which = 0;
1068  for (i = 0; i < numberRowsSmall_; i++) {
1069    if (!lastRow[i]) {
1070      lastRow[i] = which;
1071#ifdef SMALL_PERMUTE
1072      int realRow = fromSmallToBigRow[i];
1073      //if (xxx[realRow]!=0) abort();
1074      //xxx[realRow]=-1;
1075      permuteAddress_[realRow] = numberGoodU_ + which;
1076      densePermute[which] = realRow;
1077#else
1078      permuteAddress_[i] = numberGoodU_ + which;
1079      densePermute[which] = i;
1080#endif
1081      which++;
1082    }
1083  }
1084  //for (i=0;i<numberRowsSmall_;i++) {
1085  //int realRow=fromSmallToBigRow[i];
1086  //if(xxx[realRow]>0) abort();
1087  //}
1088  //for L part
1089  CoinBigIndex *COIN_RESTRICT startColumnL = startColumnLAddress_;
1090#if ABC_DENSE_CODE < 0
1091  CoinFactorizationDouble *COIN_RESTRICT elementL = elementLAddress_;
1092  CoinSimplexInt *COIN_RESTRICT indexRowL = indexRowLAddress_;
1093#endif
1094  CoinBigIndex endL = startColumnL[numberGoodL_];
1095  //take out of U
1096  CoinFactorizationDouble *COIN_RESTRICT column = denseArea;
1097  CoinSimplexInt rowsDone = 0;
1098  CoinSimplexInt iColumn = 0;
1099  CoinSimplexInt *COIN_RESTRICT pivotColumn = pivotColumnAddress_;
1100  CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
1101  CoinBigIndex *startColumnU = startColumnUAddress_;
1102#ifdef ABC_USE_FUNCTION_POINTERS
1103  scatterStruct *COIN_RESTRICT scatter = scatterUColumn();
1104#if ABC_USE_FUNCTION_POINTERS
1105  extern scatterUpdate AbcScatterLowSubtract[9];
1106  extern scatterUpdate AbcScatterHighSubtract[4];
1107#endif
1108  CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
1109  CoinSimplexInt *COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
1110#endif
1111#if ABC_DENSE_CODE == 2
1112  int nDense = 0;
1113#endif
1114  for (iColumn = 0; iColumn < numberRowsSmall_; iColumn++) {
1115    if (numberInColumn[iColumn]) {
1116#if ABC_DENSE_CODE == 2
1117      nDense++;
1118#endif
1119      //move
1120      CoinBigIndex start = startColumnU[iColumn];
1121      CoinSimplexInt number = numberInColumn[iColumn];
1122      CoinBigIndex end = start + number;
1123      for (CoinBigIndex i = start; i < end; i++) {
1124        CoinSimplexInt iRow = indexRow[i];
1125        iRow = lastRow[iRow];
1126        assert(iRow >= 0 && iRow < numberDense_);
1127#if ABC_DENSE_CODE != 2
1128        column[iRow] = element[i];
1129#else
1130#if BLOCKING8 == 8
1131        int iBlock = iRow >> 3;
1132#elif BLOCKING8 == 4
1133        int iBlock = iRow >> 2;
1134#elif BLOCKING8 == 2
1135        int iBlock = iRow >> 1;
1136#else
1137        abort();
1138#endif
1139        iRow = iRow & (BLOCKING8 - 1);
1140        column[iRow + BLOCKING8X8 * iBlock] = element[i];
1141#endif
1142      } /* endfor */
1143#if ABC_DENSE_CODE != 2
1144      column += leadingDimension_;
1145#else
1146      if ((nDense & (BLOCKING8 - 1)) != 0)
1147        column += BLOCKING8;
1148      else
1149        column += leadingDimension_ * BLOCKING8 - (BLOCKING8 - 1) * BLOCKING8;
1150#endif
1151      while (lastRow[rowsDone] < 0) {
1152        rowsDone++;
1153      } /* endwhile */
1154#ifdef ABC_USE_FUNCTION_POINTERS
1155      int numberIn = numberInColumnPlus[iColumn];
1156#ifdef SMALL_PERMUTE
1157      int realRowsDone = fromSmallToBigRow[rowsDone];
1158#if ABC_USE_FUNCTION_POINTERS
1159      if (numberIn < 9) {
1160        scatter[realRowsDone].functionPointer = AbcScatterLowSubtract[numberIn];
1161      } else {
1162        scatter[realRowsDone].functionPointer = AbcScatterHighSubtract[numberIn & 3];
1163      }
1164#endif
1165      scatter[realRowsDone].offset = lastEntryByColumnUPlus_;
1166      scatter[realRowsDone].number = numberIn;
1167#else
1168#if ABC_USE_FUNCTION_POINTERS
1169      if (numberIn < 9) {
1170        scatter[rowsDone].functionPointer = AbcScatterLowSubtract[numberIn];
1171      } else {
1172        scatter[rowsDone].functionPointer = AbcScatterHighSubtract[numberIn & 3];
1173      }
1174#endif
1175      scatter[rowsDone].offset = lastEntryByColumnUPlus_;
1176      scatter[rowsDone].number = numberIn;
1177#endif
1178      CoinBigIndex startC = start - numberIn;
1179      for (int i = startC; i < startC + numberIn; i++)
1180        elementUColumnPlus[lastEntryByColumnUPlus_++] = element[i];
1181      CoinAbcMemcpy(reinterpret_cast< CoinSimplexInt * >(elementUColumnPlus + lastEntryByColumnUPlus_), indexRow + startC, numberIn);
1182      lastEntryByColumnUPlus_ += (numberIn + 1) >> 1;
1183#endif
1184#ifdef SMALL_PERMUTE
1185      permuteAddress_[realRowsDone] = numberGoodU_;
1186      pivotColumn[numberGoodU_] = fromSmallToBigColumn[iColumn];
1187#else
1188      permuteAddress_[rowsDone] = numberGoodU_;
1189      pivotColumn[numberGoodU_] = iColumn;
1190#endif
1191      rowsDone++;
1192      startColumnL[numberGoodU_ + 1] = endL;
1193      numberInColumn[iColumn] = 0;
1194      pivotRegion[numberGoodU_] = 1.0;
1195      numberGoodU_++;
1196    }
1197  }
1198#if ABC_DENSE_CODE > 0
1199  //printf("Going dense at %d\n",numberDense_);
1200  if (denseThreshold_ > 0) {
1201    if (numberGoodU_ != numberRows_)
1202      return -1; // something odd
1203    numberGoodL_ = numberRows_;
1204    //now factorize
1205    //dgef(denseArea,&numberDense_,&numberDense_,densePermute);
1206#if ABC_DENSE_CODE != 2
1207#ifndef CLAPACK
1208    CoinSimplexInt info;
1209    F77_FUNC(dgetrf, DGETRF)
1210    (&numberDense_, &numberDense_, denseArea, &leadingDimension_, densePermute,
1211      &info);
1212
1213    // need to check size of pivots
1214    if (info)
1215      status = -1;
1216#else
1217    status = clapack_dgetrf(CblasColMajor, numberDense_, numberDense_,
1218      denseArea, leadingDimension_, densePermute);
1219    assert(!status);
1220#endif
1221#else
1222    status = CoinAbcDgetrf(numberDense_, numberDense_, denseArea, numberDense_, densePermute
1223#if ABC_PARALLEL == 2
1224      ,
1225      parallelMode_
1226#endif
1227    );
1228    if (status) {
1229      status = -1;
1230      printf("singular in dense\n");
1231    }
1232#endif
1233    return status;
1234  }
1235#else
1236  numberGoodU_ = numberRows_ - numberDense_;
1237  CoinSimplexInt base = numberGoodU_;
1238  CoinSimplexInt iDense;
1239  CoinSimplexInt numberToDo = numberDense_;
1240  denseThreshold_ = -abs(denseThreshold_); //0;
1241  CoinSimplexInt *COIN_RESTRICT nextColumn = nextColumnAddress_;
1242  const CoinSimplexInt *COIN_RESTRICT pivotColumnConst = pivotColumnAddress_;
1243  // make sure we have enough space in L and U
1244  for (iDense = 0; iDense < numberToDo; iDense++) {
1245    //how much space have we got
1246    iColumn = pivotColumnConst[base + iDense];
1247    CoinSimplexInt next = nextColumn[iColumn];
1248    CoinSimplexInt numberInPivotColumn = iDense;
1249    CoinBigIndex space = startColumnU[next]
1250      - startColumnU[iColumn]
1251      - numberInColumnPlus[next];
1252    //assume no zero elements
1253    if (numberInPivotColumn > space) {
1254      //getColumnSpace also moves fixed part
1255      if (!getColumnSpace(iColumn, numberInPivotColumn)) {
1256        return -99;
1257      }
1258    }
1259    // set so further moves will work
1260    numberInColumn[iColumn] = numberInPivotColumn;
1261  }
1262  // Fill in ?
1263  for (iColumn = numberGoodU_ + numberToDo; iColumn < numberRows_; iColumn++) {
1264    nextRow[iColumn] = iColumn;
1265    startColumnL[iColumn + 1] = endL;
1266    pivotRegion[iColumn] = 1.0;
1267  }
1268  if (lengthL_ + full * 0.5 > lengthAreaL_) {
1269    //need more memory
1270    if ((messageLevel_ & 4) != 0)
1271      std::cout << "more memory needed in middle of invert" << std::endl;
1272    return -99;
1273  }
1274  CoinFactorizationDouble *COIN_RESTRICT elementU = elementUAddress_;
1275  CoinSimplexInt *COIN_RESTRICT ipiv = densePermute - numberDense_;
1276  int returnCode = CoinAbcDgetrf(numberDense_, numberDense_, denseArea, numberDense_, ipiv);
1277  if (!returnCode) {
1278    CoinSimplexDouble *COIN_RESTRICT element = denseArea;
1279    for (int iDense = 0; iDense < numberToDo; iDense++) {
1280      int pivotRow = ipiv[iDense];
1281      // get original row
1282      CoinSimplexInt originalRow = densePermute[pivotRow];
1283      // swap
1284      densePermute[pivotRow] = densePermute[iDense];
1285      densePermute[iDense] = originalRow;
1286    }
1287    for (int iDense = 0; iDense < numberToDo; iDense++) {
1288      int iColumn = pivotColumnConst[base + iDense];
1289      // get original row
1290      CoinSimplexInt originalRow = densePermute[iDense];
1291      // do nextRow
1292      lastRow[originalRow] = -2; //mark
1293      permuteAddress_[originalRow] = numberGoodU_;
1294      CoinFactorizationDouble pivotMultiplier = element[iDense];
1295      //printf("pivotMultiplier %g\n",pivotMultiplier);
1296      pivotRegion[numberGoodU_] = pivotMultiplier;
1297      // Do L
1298      CoinBigIndex l = lengthL_;
1299      startColumnL[numberGoodL_] = l; //for luck and first time
1300      for (int iRow = iDense + 1; iRow < numberDense_; iRow++) {
1301        CoinFactorizationDouble value = element[iRow];
1302        if (!TEST_LESS_THAN_TOLERANCE(value)) {
1303          //printf("L %d row %d value %g\n",l,densePermute[iRow],value);
1304          indexRowL[l] = densePermute[iRow];
1305          elementL[l++] = value;
1306        }
1307      }
1308      numberGoodL_++;
1309      lengthL_ = l;
1310      startColumnL[numberGoodL_] = l;
1311      // update U column
1312      CoinBigIndex start = startColumnU[iColumn];
1313      for (int iRow = 0; iRow < iDense; iRow++) {
1314        if (!TEST_LESS_THAN_TOLERANCE(element[iRow])) {
1315          //printf("U %d row %d value %g\n",start,densePermute[iRow],element[iRow]);
1316          indexRowU[start] = densePermute[iRow];
1317          elementU[start++] = element[iRow];
1318        }
1319      }
1320      element += numberDense_;
1321      numberInColumn[iColumn] = 0;
1322      numberInColumnPlus[iColumn] += start - startColumnU[iColumn];
1323      startColumnU[iColumn] = start;
1324      numberGoodU_++;
1325    }
1326  } else {
1327    return -1;
1328  }
1329  numberDense_ = 0;
1330#endif
1331  return status;
1332}
1333#endif
1334#if 0
1335// Separate out links with same row/column count
1336void
1337CoinAbcTypeFactorization::separateLinks()
1338{
1339  const CoinSimplexInt count=1;
1340  CoinSimplexInt * COIN_RESTRICT nextCount = this->nextCountAddress_;
1341  CoinSimplexInt * COIN_RESTRICT firstCount = this->firstCountAddress_;
1342  CoinSimplexInt * COIN_RESTRICT lastCount = this->lastCountAddress_;
1343  CoinSimplexInt next = firstCount[count];
1344  CoinSimplexInt firstRow=-1;
1345  CoinSimplexInt firstColumn=-1;
1346  CoinSimplexInt lastRow=-1;
1347  CoinSimplexInt lastColumn=-1;
1348  while(next>=0) {
1349    CoinSimplexInt next2=nextCount[next];
1350    if (next>=numberRows_) {
1351      nextCount[next]=-1;
1352      // Column
1353      if (firstColumn>=0) {
1354        lastCount[next]=lastColumn;
1355        nextCount[lastColumn]=next;
1356      } else {
1357        lastCount[next]= -2 - count;
1358        firstColumn=next;
1359      }
1360      lastColumn=next;
1361    } else {
1362      // Row
1363      if (firstRow>=0) {
1364        lastCount[next]=lastRow;
1365        nextCount[lastRow]=next;
1366      } else {
1367        lastCount[next]= -2 - count;
1368        firstRow=next;
1369      }
1370      lastRow=next;
1371    }
1372    next=next2;
1373  }
1374  if (firstRow>=0) {
1375    firstCount[count]=firstRow;
1376    nextCount[lastRow]=firstColumn;
1377    if (firstColumn>=0)
1378      lastCount[firstColumn]=lastRow;
1379  } else if (firstRow<0) {
1380    firstCount[count]=firstColumn;
1381  } else if (firstColumn>=0) {
1382    firstCount[count]=firstColumn;
1383    nextCount[lastColumn]=firstRow;
1384    if (firstRow>=0)
1385      lastCount[firstRow]=lastColumn;
1386  }
1387}
1388#endif
1389#endif
1390
1391/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1392*/
Note: See TracBrowser for help on using the repository browser.