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

Last change on this file since 1879 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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