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

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