Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/CoinAbcBaseFactorization2.cpp

    r1910 r2385  
    1717//return code is <0 error, 0= finished
    1818CoinSimplexInt
    19 CoinAbcTypeFactorization::factorSparse (  )
     19CoinAbcTypeFactorization::factorSparse()
    2020{
    21   CoinSimplexInt * COIN_RESTRICT indexRow = indexRowUAddress_;
    22   CoinSimplexInt * COIN_RESTRICT indexColumn = indexColumnUAddress_;
    23   CoinFactorizationDouble * COIN_RESTRICT element = elementUAddress_;
     21  CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
     22  CoinSimplexInt *COIN_RESTRICT indexColumn = indexColumnUAddress_;
     23  CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
    2424  workArea_.conditionalNew(numberRows_);
    25   workAreaAddress_=workArea_.array();
    26   CoinFactorizationDouble *  COIN_RESTRICT workArea = workAreaAddress_;
     25  workAreaAddress_ = workArea_.array();
     26  CoinFactorizationDouble *COIN_RESTRICT workArea = workAreaAddress_;
    2727  double initialLargest = preProcess3();
    2828  // 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;
     29  initialLargest = 1.0e-10 / initialLargest;
     30#if ABC_NORMAL_DEBUG > 0
     31  double largestPivot = 1.0;
     32  double smallestPivot = 1.0;
    3333#endif
    3434  CoinSimplexInt status = 0;
    3535  //do slacks first
    36   CoinSimplexInt *  COIN_RESTRICT numberInRow = numberInRowAddress_;
    37   CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
     36  CoinSimplexInt *COIN_RESTRICT numberInRow = numberInRowAddress_;
     37  CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
    3838  //CoinSimplexInt * numberInColumnPlus = numberInColumnPlusAddress_;
    39   CoinBigIndex *  COIN_RESTRICT startColumnU = startColumnUAddress_;
     39  CoinBigIndex *COIN_RESTRICT startColumnU = startColumnUAddress_;
    4040  //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;
     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;
    4646#if 0
    4747        {
     
    5858  //separateLinks();
    5959  // while just singletons - do singleton rows first
    60   CoinSimplexInt *  COIN_RESTRICT pivotColumn = pivotColumnAddress_;
     60  CoinSimplexInt *COIN_RESTRICT pivotColumn = pivotColumnAddress_;
    6161  // 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_;
     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_;
    7070  startColumnL[numberGoodL_] = lengthL_; //for luck and first time
    7171#ifdef ABC_USE_FUNCTION_POINTERS
    72   CoinSimplexInt *  COIN_RESTRICT nextColumn = nextColumnAddress_;
    73   CoinSimplexInt *  COIN_RESTRICT lastColumn = lastColumnAddress_;
    74   scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
     72  CoinSimplexInt *COIN_RESTRICT nextColumn = nextColumnAddress_;
     73  CoinSimplexInt *COIN_RESTRICT lastColumn = lastColumnAddress_;
     74  scatterStruct *COIN_RESTRICT scatter = scatterUColumn();
    7575#if ABC_USE_FUNCTION_POINTERS
    7676  extern scatterUpdate AbcScatterLowSubtract[9];
    7777  extern scatterUpdate AbcScatterHighSubtract[4];
    7878#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_;
     79  CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
     80#endif
     81#if CONVERTROW > 1
     82  CoinBigIndex *COIN_RESTRICT convertRowToColumn = convertRowToColumnUAddress_;
     83  CoinBigIndex *COIN_RESTRICT convertColumnToRow = convertColumnToRowUAddress_;
    8484#endif
    8585  //#define PRINT_INVERT_STATS
    8686#ifdef PRINT_INVERT_STATS
    87   int numberGoodUX=numberGoodU_;
    88   int numberTakenOutOfUSpace=0;
     87  int numberGoodUX = numberGoodU_;
     88  int numberTakenOutOfUSpace = 0;
    8989#endif
    9090  CoinSimplexInt look = firstCount[1];
    91   while (look>=0) {
     91  while (look >= 0) {
    9292    checkLinks(1);
    9393    CoinSimplexInt iPivotRow;
    94     CoinSimplexInt iPivotColumn=-1;
     94    CoinSimplexInt iPivotColumn = -1;
    9595#ifdef SMALL_PERMUTE
    9696    int realPivotRow;
    97     int realPivotColumn=-1;
    98 #endif
    99     CoinFactorizationDouble pivotMultiplier=0.0;
    100     if ( look < numberRows_ ) {
    101       assert (numberInRow[look]==1);
     97    int realPivotColumn = -1;
     98#endif
     99    CoinFactorizationDouble pivotMultiplier = 0.0;
     100    if (look < numberRows_) {
     101      assert(numberInRow[look] == 1);
    102102      iPivotRow = look;
    103103      CoinBigIndex start = startRow[iPivotRow];
    104104      iPivotColumn = indexColumn[start];
    105105#ifdef SMALL_PERMUTE
    106       realPivotRow=fromSmallToBigRow_[iPivotRow];
    107       realPivotColumn=fromSmallToBigColumn_[iPivotColumn];
     106      realPivotRow = fromSmallToBigRow_[iPivotRow];
     107      realPivotColumn = fromSmallToBigColumn_[iPivotColumn];
    108108      //printf("singleton row realPivotRow %d realPivotColumn %d\n",
    109109      //     realPivotRow,realPivotColumn);
    110110#endif
    111       assert (numberInColumn[iPivotColumn]>0);
     111      assert(numberInColumn[iPivotColumn] > 0);
    112112      CoinBigIndex startC = startColumn[iPivotColumn];
    113       double minimumValue = fabs(element[startC]*tolerance );
    114 #if CONVERTROW<2
     113      double minimumValue = fabs(element[startC] * tolerance);
     114#if CONVERTROW < 2
    115115      CoinBigIndex where = startC;
    116       while ( indexRow[where] != iPivotRow ) {
    117         where++;
     116      while (indexRow[where] != iPivotRow) {
     117        where++;
    118118      }
    119119#else
    120120      CoinBigIndex where = convertRowToColumn[start];
    121121#endif
    122       assert ( where < startC+numberInColumn[iPivotColumn]);
     122      assert(where < startC + numberInColumn[iPivotColumn]);
    123123      CoinFactorizationDouble value = element[where];
    124124      // 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];
     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];
    137137#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];
     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];
    160160            convertRowToColumn[whereRow] = whereColumnEntry;
    161             convertColumnToRow[whereColumnEntry] = whereRow;
    162 #endif
    163             iNumberInRow--;
    164             numberInRow[iRow] = iNumberInRow;
    165             modifyLink ( iRow, iNumberInRow );
    166             checkLinks();
    167           }
    168         }
     161            convertColumnToRow[whereColumnEntry] = whereRow;
     162#endif
     163            iNumberInRow--;
     164            numberInRow[iRow] = iNumberInRow;
     165            modifyLink(iRow, iNumberInRow);
     166            checkLinks();
     167          }
     168        }
    169169      } else {
    170         //printf("Rejecting as element %g largest %g\n",value,element[startC]);
     170        //printf("Rejecting as element %g largest %g\n",value,element[startC]);
    171171      }
    172172    } else {
    173173      iPivotColumn = look - numberRows_;
    174       assert ( numberInColumn[iPivotColumn] == 1 );
     174      assert(numberInColumn[iPivotColumn] == 1);
    175175      {
    176         CoinBigIndex startC = startColumn[iPivotColumn];
    177         iPivotRow = indexRow[startC];
     176        CoinBigIndex startC = startColumn[iPivotColumn];
     177        iPivotRow = indexRow[startC];
    178178#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]++;
     179        realPivotRow = fromSmallToBigRow_[iPivotRow];
     180        realPivotColumn = fromSmallToBigColumn_[iPivotColumn];
     181        int realPivotRow;
     182        realPivotRow = fromSmallToBigRow_[iPivotRow];
     183        indexRow[startC] = realPivotRow;
     184        //printf("singleton column realPivotRow %d realPivotColumn %d\n",
     185        //   realPivotRow,realPivotColumn);
     186#endif
     187        //store pivot columns (so can easily compress)
     188        assert(numberInColumn[iPivotColumn] == 1);
     189        CoinSimplexInt numberDoRow = numberInRow[iPivotRow];
     190        int numberDoColumn = numberInColumn[iPivotColumn] - 1;
     191        totalElements_ -= (numberDoRow + numberDoColumn);
     192        CoinBigIndex startR = startRow[iPivotRow];
     193        CoinBigIndex endR = startR + numberDoRow;
     194        //clean up counts
     195        pivotMultiplier = 1.0 / element[startC];
     196        // would I be better off doing other way first??
     197        for (CoinBigIndex i = startR; i < endR; i++) {
     198          CoinSimplexInt iColumn = indexColumn[i];
     199          if (iColumn == iPivotColumn)
     200            continue;
     201          assert(numberInColumn[iColumn]);
     202          CoinSimplexInt number = numberInColumn[iColumn] - 1;
     203          //modify linked list
     204          modifyLink(iColumn + numberRows_, number);
     205          CoinBigIndex start = startColumn[iColumn];
     206          //move pivot row element
     207          if (number) {
     208#if CONVERTROW < 2
     209            CoinBigIndex pivot = start;
     210            CoinSimplexInt iRow = indexRow[pivot];
     211            while (iRow != iPivotRow) {
     212              pivot++;
     213              iRow = indexRow[pivot];
     214            }
     215#else
     216            CoinBigIndex pivot = convertRowToColumn[i];
     217#endif
     218            assert(pivot < startColumn[iColumn] + numberInColumn[iColumn]);
     219            if (pivot != start) {
     220              //move largest one up
     221              CoinFactorizationDouble value = element[start];
     222              int iRow = indexRow[start];
     223              element[start] = element[pivot];
     224              indexRow[start] = indexRow[pivot];
     225              element[pivot] = element[start + 1];
     226              indexRow[pivot] = indexRow[start + 1];
     227#if CONVERTROW > 1
     228              CoinBigIndex whereRowEntry = convertColumnToRow[start + 1];
     229              CoinBigIndex whereRowEntry2 = convertColumnToRow[start];
     230              convertRowToColumn[whereRowEntry] = pivot;
     231              convertColumnToRow[pivot] = whereRowEntry;
     232              convertRowToColumn[whereRowEntry2] = start + 1;
     233              convertColumnToRow[start + 1] = whereRowEntry2;
     234#endif
     235              element[start + 1] = value;
     236              indexRow[start + 1] = iRow;
     237            } else {
     238              //find new largest element
     239              CoinSimplexInt iRowSave = indexRow[start + 1];
     240              CoinFactorizationDouble valueSave = element[start + 1];
     241              CoinFactorizationDouble valueLargest = fabs(valueSave);
     242              CoinBigIndex end = start + numberInColumn[iColumn];
     243              CoinBigIndex largest = start + 1;
     244              for (CoinBigIndex k = start + 2; k < end; k++) {
     245                CoinFactorizationDouble value = element[k];
     246                CoinFactorizationDouble valueAbs = fabs(value);
     247                if (valueAbs > valueLargest) {
     248                  valueLargest = valueAbs;
     249                  largest = k;
     250                }
     251              }
     252              indexRow[start + 1] = indexRow[largest];
     253              element[start + 1] = element[largest];
     254#if CONVERTROW > 1
     255              CoinBigIndex whereRowEntry = convertColumnToRow[largest];
     256              CoinBigIndex whereRowEntry2 = convertColumnToRow[start + 1];
     257              convertRowToColumn[whereRowEntry] = start + 1;
     258              convertColumnToRow[start + 1] = whereRowEntry;
     259              convertRowToColumn[whereRowEntry2] = largest;
     260              convertColumnToRow[largest] = whereRowEntry2;
     261#endif
     262              indexRow[largest] = iRowSave;
     263              element[largest] = valueSave;
     264            }
     265          }
     266          //clean up counts
     267          numberInColumn[iColumn]--;
     268          numberInColumnPlus[iColumn]++;
    270269#ifdef SMALL_PERMUTE
    271           indexRow[start]=realPivotRow;
    272 #endif
    273           startColumn[iColumn]++;
    274         }
     270          indexRow[start] = realPivotRow;
     271#endif
     272          startColumn[iColumn]++;
     273        }
    275274      }
    276275    }
     
    283282      nextRow[last] = next;
    284283      lastRow[next] = last;
    285       lastRow[iPivotRow] =-2; //mark
     284      lastRow[iPivotRow] = -2; //mark
    286285#ifdef SMALL_PERMUTE
    287286      // mark
    288       fromSmallToBigRow_[iPivotRow]=-1;
    289       fromSmallToBigColumn_[iPivotColumn]=-1;
    290       permuteAddress_[realPivotRow]=numberGoodU_;
    291 #else
    292       permuteAddress_[iPivotRow]=numberGoodU_;
     287      fromSmallToBigRow_[iPivotRow] = -1;
     288      fromSmallToBigColumn_[iPivotColumn] = -1;
     289      permuteAddress_[realPivotRow] = numberGoodU_;
     290#else
     291      permuteAddress_[iPivotRow] = numberGoodU_;
    293292#endif
    294293#ifdef ABC_USE_FUNCTION_POINTERS
     
    296295#ifdef SMALL_PERMUTE
    297296#if ABC_USE_FUNCTION_POINTERS
    298       if (number<9) {
    299         scatter[realPivotRow].functionPointer=AbcScatterLowSubtract[number];
     297      if (number < 9) {
     298        scatter[realPivotRow].functionPointer = AbcScatterLowSubtract[number];
    300299      } else {
    301         scatter[realPivotRow].functionPointer=AbcScatterHighSubtract[number&3];
    302       }
    303 #endif
    304       scatter[realPivotRow].offset=lastEntryByColumnUPlus_;
    305       scatter[realPivotRow].number=number;
     300        scatter[realPivotRow].functionPointer = AbcScatterHighSubtract[number & 3];
     301      }
     302#endif
     303      scatter[realPivotRow].offset = lastEntryByColumnUPlus_;
     304      scatter[realPivotRow].number = number;
    306305#else
    307306#if ABC_USE_FUNCTION_POINTERS
    308       if (number<9) {
    309         scatter[iPivotRow].functionPointer=AbcScatterLowSubtract[number];
     307      if (number < 9) {
     308        scatter[iPivotRow].functionPointer = AbcScatterLowSubtract[number];
    310309      } 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;
     310        scatter[iPivotRow].functionPointer = AbcScatterHighSubtract[number & 3];
     311      }
     312#endif
     313      scatter[iPivotRow].offset = lastEntryByColumnUPlus_;
     314      scatter[iPivotRow].number = number;
     315#endif
     316      CoinBigIndex startC = startColumn[iPivotColumn] - number;
     317      for (int i = startC; i < startC + number; i++)
     318        elementUColumnPlus[lastEntryByColumnUPlus_++] = element[i] * pivotMultiplier;
     319      CoinAbcMemcpy(reinterpret_cast< CoinSimplexInt * >(elementUColumnPlus + lastEntryByColumnUPlus_), indexRow + startC, number);
     320      lastEntryByColumnUPlus_ += (number + 1) >> 1;
    322321#endif
    323322      numberInColumn[iPivotColumn] = 0;
    324323      numberInRow[iPivotRow] = 0;
    325324      //modify linked list for pivots
    326       deleteLink ( iPivotRow );
    327       deleteLink ( iPivotColumn + numberRows_ );
    328     checkLinks();
     325      deleteLink(iPivotRow);
     326      deleteLink(iPivotColumn + numberRows_);
     327      checkLinks();
    329328#ifdef SMALL_PERMUTE
    330       assert (permuteAddress_[realPivotRow]==numberGoodU_);
     329      assert(permuteAddress_[realPivotRow] == numberGoodU_);
    331330      pivotColumn[numberGoodU_] = realPivotColumn;
    332331#else
    333       assert (permuteAddress_[iPivotRow]==numberGoodU_);
     332      assert(permuteAddress_[iPivotRow] == numberGoodU_);
    334333      pivotColumn[numberGoodU_] = iPivotColumn;
    335334#endif
     
    338337    } else {
    339338      //take out for moment
    340       modifyLink ( iPivotRow, numberRows_ + 1 );
    341     }
    342     look=-1; //nextCount[look];
    343     if (look<0) {
     339      modifyLink(iPivotRow, numberRows_ + 1);
     340    }
     341    look = -1; //nextCount[look];
     342    if (look < 0) {
    344343      // start again
    345344      look = firstCount[1];
    346345    }
    347     assert (iPivotColumn>=0);
     346    assert(iPivotColumn >= 0);
    348347#ifdef ABC_USE_FUNCTION_POINTERS
    349348    if (!numberInColumn[iPivotColumn]) {
    350       int iNext=nextColumn[iPivotColumn];
    351       int iLast=lastColumn[iPivotColumn];
    352       assert (iLast>=0);
     349      int iNext = nextColumn[iPivotColumn];
     350      int iLast = lastColumn[iPivotColumn];
     351      assert(iLast >= 0);
    353352#ifdef PRINT_INVERT_STATS
    354353      numberTakenOutOfUSpace++;
    355354#endif
    356       lastColumn[iNext]=iLast;
    357       nextColumn[iLast]=iNext;
     355      lastColumn[iNext] = iLast;
     356      nextColumn[iLast] = iNext;
    358357    }
    359358#endif
    360359  }
    361360  // put back all rejected
    362   look = firstCount[numberRows_+1];
    363   while (look>=0) {
     361  look = firstCount[numberRows_ + 1];
     362  while (look >= 0) {
    364363#ifndef NDEBUG
    365     if ( look < numberRows_ ) {
    366       assert (numberInRow[look]==1);
     364    if (look < numberRows_) {
     365      assert(numberInRow[look] == 1);
    367366    } else {
    368367      int iColumn = look - numberRows_;
    369       assert ( numberInColumn[iColumn] == 1 );
    370     }
    371 #endif
    372     int nextLook=nextCount[look];
    373     modifyLink ( look, 1 );
    374     look=nextLook;
     368      assert(numberInColumn[iColumn] == 1);
     369    }
     370#endif
     371    int nextLook = nextCount[look];
     372    modifyLink(look, 1);
     373    look = nextLook;
    375374  }
    376375#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;
     376  if (numberGoodU_ < numberRows_ && numberGoodU_ > numberSlacks_ + (numberRows_ >> 3)) {
     377    CoinSimplexInt *COIN_RESTRICT fromBigToSmallRow = reinterpret_cast< CoinSimplexInt * >(workArea_.array());
     378    CoinSimplexInt *COIN_RESTRICT fromBigToSmallColumn = fromBigToSmallRow + numberRows_;
     379    int nSmall = 0;
     380    for (int i = 0; i < numberRowsSmall_; i++) {
     381      int realRow = fromSmallToBigRow_[i];
     382      if (realRow < 0) {
     383        fromBigToSmallRow[i] = -1;
    385384      } 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;
     385        // in
     386        fromBigToSmallRow[i] = nSmall;
     387        numberInRow[nSmall] = numberInRow[i];
     388        startRow[nSmall] = startRow[i];
     389        fromSmallToBigRow_[nSmall++] = realRow;
     390      }
     391    }
     392    nSmall = 0;
     393    for (int i = 0; i < numberRowsSmall_; i++) {
     394      int realColumn = fromSmallToBigColumn_[i];
     395      if (realColumn < 0) {
     396        fromBigToSmallColumn[i] = -1;
    398397      } 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);
     398        // in
     399        fromBigToSmallColumn[i] = nSmall;
     400        numberInColumn[nSmall] = numberInColumn[i];
     401        numberInColumnPlus[nSmall] = numberInColumnPlus[i];
     402        startColumn[nSmall] = startColumn[i];
     403        fromSmallToBigColumn_[nSmall++] = realColumn;
     404      }
     405    }
     406    CoinAbcMemset0(numberInRow + nSmall, numberRowsSmall_ - nSmall);
     407    CoinAbcMemset0(numberInColumn + nSmall, numberRowsSmall_ - nSmall);
    409408    // 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]];
     409    for (int i = 0; i < numberRowsSmall_; i++) {
     410      CoinBigIndex start = startRow[i];
     411      CoinBigIndex end = start + numberInRow[i];
     412      for (CoinBigIndex j = start; j < end; j++) {
     413        indexColumn[j] = fromBigToSmallColumn[indexColumn[j]];
     414      }
     415    }
     416    for (int i = 0; i < numberRowsSmall_; i++) {
     417      CoinBigIndex start = startColumn[i];
     418      CoinBigIndex end = start + numberInColumn[i];
     419      for (CoinBigIndex j = start; j < end; j++) {
     420        indexRow[j] = fromBigToSmallRow[indexRow[j]];
    422421      }
    423422    }
    424423    // find area somewhere
    425     int * temp = fromSmallToBigColumn_+nSmall;
    426     CoinSimplexInt * nextFake=temp;
     424    int *temp = fromSmallToBigColumn_ + nSmall;
     425    CoinSimplexInt *nextFake = temp;
    427426    //CoinSimplexInt lastFake=temp+numberRows_;
    428     CoinAbcMemcpy(nextFake,nextRow,numberRowsSmall_);
    429     nextFake[numberRows_]=nextRow[numberRows_];
     427    CoinAbcMemcpy(nextFake, nextRow, numberRowsSmall_);
     428    nextFake[numberRows_] = nextRow[numberRows_];
    430429    //CoinAbcMemcpy(lastFake,lastRow,numberRowsSmall_);
    431430    // 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;
     431    CoinSimplexInt lastBigRow = numberRows_;
     432    CoinSimplexInt lastSmallRow = numberRows_;
     433    for (CoinSimplexInt i = 0; i < nSmall; i++) {
     434      CoinSimplexInt bigRow = nextFake[lastBigRow];
     435      assert(bigRow >= 0 && bigRow < numberRowsSmall_);
     436      CoinSimplexInt smallRow = fromBigToSmallRow[bigRow];
     437      nextRow[lastSmallRow] = smallRow;
     438      lastRow[smallRow] = lastSmallRow;
     439      lastBigRow = bigRow;
     440      lastSmallRow = smallRow;
     441    }
     442    assert(nextFake[lastBigRow] == numberRows_);
     443    nextRow[lastSmallRow] = numberRows_;
     444    lastRow[numberRows_] = lastSmallRow;
    446445    // 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;
     446    CoinAbcMemcpy(nextFake, nextColumn, numberRowsSmall_);
     447    nextFake[maximumRowsExtra_] = nextColumn[maximumRowsExtra_];
     448    CoinSimplexInt lastBigColumn = maximumRowsExtra_;
     449    CoinSimplexInt lastSmallColumn = maximumRowsExtra_;
     450    for (CoinSimplexInt i = 0; i < nSmall; i++) {
     451      CoinSimplexInt bigColumn = nextFake[lastBigColumn];
     452      assert(bigColumn >= 0 && bigColumn < numberRowsSmall_);
     453      CoinSimplexInt smallColumn = fromBigToSmallColumn[bigColumn];
     454      nextColumn[lastSmallColumn] = smallColumn;
     455      lastColumn[smallColumn] = lastSmallColumn;
     456      lastBigColumn = bigColumn;
     457      lastSmallColumn = smallColumn;
     458    }
     459    assert(nextFake[lastBigColumn] == maximumRowsExtra_);
     460    nextColumn[lastSmallColumn] = maximumRowsExtra_;
     461    lastColumn[maximumRowsExtra_] = lastSmallColumn;
    463462    // 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;
     463    CoinSimplexInt *COIN_RESTRICT lastCount = this->lastCountAddress_;
     464    CoinAbcMemcpy(temp, nextCount, numberRows_ + numberRowsSmall_);
     465    for (int iCount = 0; iCount <= nSmall; iCount++) {
     466      CoinSimplexInt nextBig = firstCount[iCount];
     467      if (nextBig >= 0) {
     468        //CoinSimplexInt next=firstCount[iCount];
     469        CoinSimplexInt nextSmall;
     470        if (nextBig < numberRows_)
     471          nextSmall = fromBigToSmallRow[nextBig];
     472        else
     473          nextSmall = fromBigToSmallColumn[nextBig - numberRows_] + numberRows_;
     474        firstCount[iCount] = nextSmall;
     475        CoinSimplexInt *where = &lastCount[nextSmall];
     476        while (nextBig >= 0) {
     477          CoinSimplexInt thisSmall = nextSmall;
     478          nextBig = temp[nextBig];
     479          if (nextBig >= 0) {
     480            if (nextBig < numberRows_)
     481              nextSmall = fromBigToSmallRow[nextBig];
     482            else
     483              nextSmall = fromBigToSmallColumn[nextBig - numberRows_] + numberRows_;
     484            lastCount[nextSmall] = thisSmall;
     485            nextCount[thisSmall] = nextSmall;
     486          } else {
     487            nextCount[thisSmall] = nextBig;
     488          }
     489        }
     490        assert(nextSmall >= 0);
     491        // fill in odd one
     492        *where = iCount - numberRows_ - 2;
    494493      }
    495494    }
    496495    //printf("%d rows small1 %d small2 %d\n",numberRows_,numberRowsSmall_,nSmall);
    497     numberRowsSmall_=nSmall;
     496    numberRowsSmall_ = nSmall;
    498497    //exit(0);
    499498  }
     
    501500  // Put .hpp functions in separate file for profiling
    502501#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_;
     502  int numberGoodUY = numberGoodU_;
     503  int numberElsLeftX = 0;
     504  for (int i = 0; i < numberRows_; i++)
     505    numberElsLeftX += numberInColumn[i];
     506  int saveN1X = totalElements_;
    508507#endif
    509508#if ABC_DENSE_CODE
     
    514513  CoinBigIndex workSize = 1000;
    515514  workArea2_.conditionalNew(workSize);
    516   workArea2Address_=workArea2_.array();
    517   CoinSimplexUnsignedInt *  COIN_RESTRICT workArea2 = workArea2Address_;
     515  workArea2Address_ = workArea2_.array();
     516  CoinSimplexUnsignedInt *COIN_RESTRICT workArea2 = workArea2Address_;
    518517
    519518  //set markRow so no rows updated
    520519  //set markRow so no rows updated
    521   CoinSimplexInt *  COIN_RESTRICT markRow = markRow_.array();
    522   CoinFillN ( markRow, numberRowsSmall_, LARGE_UNSET);
    523   CoinZeroN ( workArea, numberRowsSmall_ );
     520  CoinSimplexInt *COIN_RESTRICT markRow = markRow_.array();
     521  CoinFillN(markRow, numberRowsSmall_, LARGE_UNSET);
     522  CoinZeroN(workArea, numberRowsSmall_);
    524523  CoinFactorizationDouble pivotTolerance = pivotTolerance_;
    525524  CoinSimplexInt numberTrials = numberTrials_;
     
    528527  CoinSimplexInt count = 1;
    529528  startColumnL[numberGoodL_] = lengthL_; //for luck and first time
    530   numberRowsLeft_=numberRows_-numberGoodU_;
    531   while ( count <= numberRowsLeft_ ) {
     529  numberRowsLeft_ = numberRows_ - numberGoodU_;
     530  while (count <= numberRowsLeft_) {
    532531    CoinBigIndex minimumCount = COIN_INT_MAX;
    533532    CoinFactorizationDouble minimumCost = COIN_DBL_MAX;
     
    553552    CoinSimplexInt trials = 0;
    554553    //CoinSimplexInt *  COIN_RESTRICT pivotColumn = pivotColumnAddress_;
    555     while ( look >= 0 ) {
     554    while (look >= 0) {
    556555      checkLinks(1);
    557       if ( look < numberRows_ ) {
     556      if (look < numberRows_) {
    558557        CoinSimplexInt iRow = look;
    559558        look = nextCount[look];
     
    561560        CoinBigIndex start = startRow[iRow];
    562561        CoinBigIndex end = start + count;
    563        
     562
    564563        CoinBigIndex i;
    565         for ( i = start; i < end; i++ ) {
     564        for (i = start; i < end; i++) {
    566565          CoinSimplexInt iColumn = indexColumn[i];
    567           assert (numberInColumn[iColumn]>0);
    568           CoinFactorizationDouble cost = ( count - 1 ) * numberInColumn[iColumn] + 0.1;
    569          
    570           if ( cost < minimumCost ) {
     566          assert(numberInColumn[iColumn] > 0);
     567          CoinFactorizationDouble cost = (count - 1) * numberInColumn[iColumn] + 0.1;
     568
     569          if (cost < minimumCost) {
    571570            CoinBigIndex where = startColumn[iColumn];
    572571            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 ) {
     572
     573            minimumValue = fabs(minimumValue) * pivotTolerance;
     574            if (count == 1)
     575              minimumValue = CoinMin(minimumValue, 0.999999);
     576            while (indexRow[where] != iRow) {
    578577              where++;
    579             }                   /* endwhile */
    580             assert ( where < startColumn[iColumn] +
    581                      numberInColumn[iColumn]);
     578            } /* endwhile */
     579            assert(where < startColumn[iColumn] + numberInColumn[iColumn]);
    582580            CoinFactorizationDouble value = element[where];
    583            
    584             value = fabs ( value );
    585             if ( value >= minimumValue ) {
     581
     582            value = fabs(value);
     583            if (value >= minimumValue) {
    586584              minimumCost = cost;
    587585              minimumCount = numberInColumn[iColumn];
     
    589587              pivotRowPosition = -1;
    590588              iPivotColumn = iColumn;
    591               assert (iPivotRow>=0&&iPivotColumn>=0);
     589              assert(iPivotRow >= 0 && iPivotColumn >= 0);
    592590              pivotColumnPosition = i;
    593               rejected=false;
    594             } else if ( iPivotRow == -1 && count > 1) {
     591              rejected = false;
     592            } else if (iPivotRow == -1 && count > 1) {
    595593              rejected = true;
    596594            }
     
    598596        }
    599597        trials++;
    600         if ( iPivotRow >= 0 && (trials >= numberTrials||minimumCount<count) ) {
     598        if (iPivotRow >= 0 && (trials >= numberTrials || minimumCount < count)) {
    601599          look = -1;
    602600          break;
    603601        }
    604         if ( rejected ) {
     602        if (rejected) {
    605603          //take out for moment
    606604          //eligible when row changes
    607           modifyLink ( iRow, numberRows_ + 1 );
     605          modifyLink(iRow, numberRows_ + 1);
    608606        }
    609607      } else {
    610608        CoinSimplexInt iColumn = look - numberRows;
    611        
     609
    612610#if 0
    613611        if ( numberInColumn[iColumn] != count ) {
     
    617615        }
    618616#endif
    619         assert ( numberInColumn[iColumn] == count );
     617        assert(numberInColumn[iColumn] == count);
    620618        look = nextCount[look];
    621619        CoinBigIndex start = startColumn[iColumn];
    622620        CoinBigIndex end = start + numberInColumn[iColumn];
    623621        CoinFactorizationDouble minimumValue = element[start];
    624        
    625         minimumValue = fabs ( minimumValue ) * pivotTolerance;
     622
     623        minimumValue = fabs(minimumValue) * pivotTolerance;
    626624        CoinBigIndex i;
    627         for ( i = start; i < end; i++ ) {
     625        for (i = start; i < end; i++) {
    628626          CoinFactorizationDouble value = element[i];
    629          
    630           value = fabs ( value );
    631           if ( value >= minimumValue ) {
     627
     628          value = fabs(value);
     629          if (value >= minimumValue) {
    632630            CoinSimplexInt iRow = indexRow[i];
    633631            CoinSimplexInt nInRow = numberInRow[iRow];
    634             assert (nInRow>0);
    635             CoinFactorizationDouble cost = ( count - 1 ) * nInRow;
    636            
    637             if ( cost < minimumCost ) {
     632            assert(nInRow > 0);
     633            CoinFactorizationDouble cost = (count - 1) * nInRow;
     634
     635            if (cost < minimumCost) {
    638636              minimumCost = cost;
    639637              minimumCount = nInRow;
     
    641639              pivotRowPosition = i;
    642640              iPivotColumn = iColumn;
    643               assert (iPivotRow>=0&&iPivotColumn>=0);
     641              assert(iPivotRow >= 0 && iPivotColumn >= 0);
    644642              pivotColumnPosition = -1;
    645643            }
     
    647645        }
    648646        trials++;
    649         if ( iPivotRow >= 0 && (trials >= numberTrials||minimumCount<count) ) {
     647        if (iPivotRow >= 0 && (trials >= numberTrials || minimumCount < count)) {
    650648          look = -1;
    651649          break;
    652650        }
    653651      }
    654     }                           /* endwhile */
    655     if (iPivotRow>=0) {
    656       assert (iPivotRow<numberRows_);
     652    } /* endwhile */
     653    if (iPivotRow >= 0) {
     654      assert(iPivotRow < numberRows_);
    657655      CoinSimplexInt numberDoRow = numberInRow[iPivotRow] - 1;
    658656      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);
     657
     658      //printf("realPivotRow %d (%d elements) realPivotColumn %d (%d elements)\n",
     659      //   fromSmallToBigRow[iPivotRow],numberDoRow+1,fromSmallToBigColumn[iPivotColumn],numberDoColumn+1);
    662660      //      printf("nGoodU %d pivRow %d (num %d) pivCol %d (num %d) - %d elements left\n",
    663661      //   numberGoodU_,iPivotRow,numberDoRow,iPivotColumn,numberDoColumn,totalElements_);
     
    681679      assert (nn==totalElements_);
    682680#endif
    683       totalElements_ -= ( numberDoRow + numberDoColumn + 1 );
     681      totalElements_ -= (numberDoRow + numberDoColumn + 1);
    684682#if 0
    685683      if (numberInColumn[5887]==0&&numberRowsSmall_>5887) {
     
    690688      }
    691689#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
     690      if (numberDoColumn > 0) {
     691        if (numberDoRow > 0) {
     692          if (numberDoColumn > 1) {
     693            //  if (1) {
     694            //need to adjust more for cache and SMP
     695            //allow at least 4 extra
     696            CoinSimplexInt increment = numberDoColumn + 1 + 4;
     697
     698            if (increment & 15) {
     699              increment = increment & (~15);
     700              increment += 16;
     701            }
     702            CoinSimplexInt increment2 =
     703
     704              (increment + COINFACTORIZATION_BITS_PER_INT - 1) >> COINFACTORIZATION_SHIFT_PER_INT;
     705            CoinBigIndex size = increment2 * numberDoRow;
     706
     707            if (size > workSize) {
     708              workSize = size;
     709              workArea2_.conditionalNew(workSize);
     710              workArea2Address_ = workArea2_.array();
     711              workArea2 = workArea2Address_;
     712            }
     713            //branch out to best pivot routine
    716714#define ABC_PARALLEL_PIVOT
    717715#ifndef ABC_USE_FUNCTION_POINTERS
     
    719717#endif
    720718#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 );
     719            //if (numberRowsSmall_==7202&&numberGoodU_==15609) {
     720            //printf("NumberGoodU %d\n",numberGoodU_);
     721            //}
     722            if (numberDoRow < 20)
     723              status = pivot(iPivotRow, iPivotColumn,
     724                pivotRowPosition, pivotColumnPosition,
     725                workArea, workArea2,
     726                increment2, markRow);
     727            else
     728              status = pivot(iPivotRow, iPivotColumn,
     729                pivotRowPosition, pivotColumnPosition,
     730                markRow);
     731#else
     732            status = pivot(iPivotRow, iPivotColumn,
     733              pivotRowPosition, pivotColumnPosition,
     734              workArea, workArea2,
     735              increment2, markRow);
    738736#endif
    739737#if 0
     
    741739              assert (markRow[i]==LARGE_UNSET);
    742740#endif
    743             //#define CHECK_SIZE
     741              //#define CHECK_SIZE
    744742#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             }
     743            {
     744              double largestU = 0.0;
     745              int iU = -1;
     746              int jU = -1;
     747              for (int i = 0; i < numberRowsSmall_; i++) {
     748                CoinBigIndex start = startColumn[i];
     749                CoinBigIndex end = start + numberInColumn[i];
     750                for (int j = start; j < end; j++) {
     751                  if (fabs(element[j]) > largestU) {
     752                    iU = i;
     753                    jU = j;
     754                    largestU = fabs(element[j]);
     755                  }
     756                }
     757              }
     758              printf("%d largest el %g at i=%d,j=%d start %d end %d count %d\n", numberGoodU_, element[jU],
     759                iU, jU, startColumn[iU], startColumn[iU] + numberInColumn[iU], count);
     760            }
    763761#endif
    764762#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             }
     763            checkLinks();
     764            if (status < 0) {
     765              count = numberRows_ + 1;
     766              break;
     767            }
     768          } else {
     769            if (!pivotOneOtherRow(iPivotRow, iPivotColumn)) {
     770              status = -99;
     771              count = numberRows_ + 1;
     772              break;
     773            }
    776774#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           }
     775            {
     776              double largestU = 0.0;
     777              int iU = -1;
     778              int jU = -1;
     779              for (int i = 0; i < numberRowsSmall_; i++) {
     780                CoinBigIndex start = startColumn[i];
     781                CoinBigIndex end = start + numberInColumn[i];
     782                for (int j = start; j < end; j++) {
     783                  if (fabs(element[j]) > largestU) {
     784                    iU = i;
     785                    jU = j;
     786                    largestU = fabs(element[j]);
     787                  }
     788                }
     789              }
     790              printf("B%d largest el %g at i=%d,j=%d\n", numberGoodU_, element[jU],
     791                iU, jU);
     792            }
     793#endif
     794            checkLinks();
     795          }
     796        } else {
     797          assert(!numberDoRow);
     798          checkLinks();
     799          if (!pivotRowSingleton(iPivotRow, iPivotColumn)) {
     800            status = -99;
     801            count = numberRows_ + 1;
     802            break;
     803          }
    806804#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         }
     805          {
     806            double largestU = 0.0;
     807            int iU = -1;
     808            int jU = -1;
     809            for (int i = 0; i < numberRowsSmall_; i++) {
     810              CoinBigIndex start = startColumn[i];
     811              CoinBigIndex end = start + numberInColumn[i];
     812              for (int j = start; j < end; j++) {
     813                if (fabs(element[j]) > largestU) {
     814                  iU = i;
     815                  jU = j;
     816                  largestU = fabs(element[j]);
     817                }
     818              }
     819            }
     820            printf("C%d largest el %g at i=%d,j=%d\n", numberGoodU_, element[jU],
     821              iU, jU);
     822          }
     823#endif
     824          checkLinks();
     825        }
    828826      } else {
    829         assert (!numberDoColumn);
    830         pivotColumnSingleton ( iPivotRow, iPivotColumn );
     827        assert(!numberDoColumn);
     828        pivotColumnSingleton(iPivotRow, iPivotColumn);
    831829#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_);
     830        {
     831          double largestU = 0.0;
     832          int iU = -1;
     833          int jU = -1;
     834          for (int i = 0; i < numberRowsSmall_; i++) {
     835            CoinBigIndex start = startColumn[i];
     836            CoinBigIndex end = start + numberInColumn[i];
     837            for (int j = start; j < end; j++) {
     838              if (fabs(element[j]) > largestU) {
     839                iU = i;
     840                jU = j;
     841                largestU = fabs(element[j]);
     842              }
     843            }
     844          }
     845          printf("D%d largest el %g at i=%d,j=%d\n", numberGoodU_, element[jU],
     846            iU, jU);
     847        }
     848#endif
     849        checkLinks();
     850      }
     851      double pivotValue = fabs(pivotRegion[numberGoodU_]);
     852#if ABC_NORMAL_DEBUG > 0
     853      largestPivot = CoinMax(largestPivot, pivotValue);
     854      smallestPivot = CoinMin(smallestPivot, pivotValue);
     855#endif
     856      if (pivotValue < initialLargest) {
     857        if (pivotTolerance_ < 0.95) {
     858          pivotTolerance_ = CoinMin(1.1 * pivotTolerance_, 0.99);
     859#if ABC_NORMAL_DEBUG > 0
     860          printf("PPPivot value of %g increasing pivot tolerance to %g\n",
     861            pivotValue, pivotTolerance_);
     862#endif
     863          status = -97;
     864          break;
     865        }
     866      }
     867      afterPivot(iPivotRow, iPivotColumn);
     868      assert(numberGoodU_ <= numberRows_);
     869      assert(startRow[numberRows_] == lengthAreaU_);
    872870#if 0
    873871      static int ixxxxxx=0;
     
    944942#if ABC_DENSE_CODE
    945943      if (!status) {
    946         status=wantToGoDense();
     944        status = wantToGoDense();
    947945      }
    948946#endif
    949947      if (status)
    950         break;
     948        break;
    951949      // start at 1 again
    952950      count = 1;
     
    954952      //end of this - onto next
    955953      count++;
    956     } 
    957   }                             /* endwhile */
    958 #if ABC_NORMAL_DEBUG>0
    959 #if ABC_SMALL<2
    960   int lenU=2*(lastEntryByColumnUPlus_/3);
     954    }
     955  } /* endwhile */
     956#if ABC_NORMAL_DEBUG > 0
     957#if ABC_SMALL < 2
     958  int lenU = 2 * (lastEntryByColumnUPlus_ / 3);
    961959  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;
     960    1.0e-10 / initialLargest, smallestPivot, largestPivot, numberRows_ - numberGoodU_,
     961    lengthL_, lenU);
     962#endif
     963#endif
     964  workArea2_.conditionalDelete();
     965  workArea2Address_ = NULL;
    968966#ifdef PRINT_INVERT_STATS
    969   int numberDense=numberRows_-numberGoodU_;
     967  int numberDense = numberRows_ - numberGoodU_;
    970968  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);
     969    numberRows_, numberRows_ - numberGoodUY, numberGoodUX, numberGoodUY - numberGoodUX,
     970    numberDense, numberElsLeftX, saveN1X, totalElements_, numberTakenOutOfUSpace);
    973971#endif
    974972  return status;
     
    979977CoinSimplexInt CoinAbcTypeFactorization::factorDense()
    980978{
    981   CoinSimplexInt status=0;
    982   numberDense_=numberRows_-numberGoodU_;
    983   if (sizeof(CoinBigIndex)==4&&numberDense_>=(2<<15)) {
     979  CoinSimplexInt status = 0;
     980  numberDense_ = numberRows_ - numberGoodU_;
     981  if (sizeof(CoinBigIndex) == 4 && numberDense_ >= (2 << 15)) {
    984982    abort();
    985   } 
     983  }
    986984  CoinBigIndex full;
    987   full = numberDense_*numberDense_;
    988   totalElements_=full;
     985  full = numberDense_ * numberDense_;
     986  totalElements_ = full;
    989987  // 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)));
     988#if ABC_DENSE_CODE == 1
     989  leadingDimension_ = ((numberDense_ >> 4) + 1) << 4;
     990#else
     991  leadingDimension_ = numberDense_;
     992#endif
     993  CoinBigIndex newSize = (leadingDimension_ + FACTOR_CPU) * numberDense_;
     994  newSize += (numberDense_ + 1) / (sizeof(CoinFactorizationDouble) / sizeof(CoinSimplexInt));
     995#if 1
     996  newSize += 2 * ((numberDense_ + 3) / (sizeof(CoinFactorizationDouble) / sizeof(short)));
     997  newSize += ((numberRows_ + 3) / (sizeof(CoinFactorizationDouble) / sizeof(short)));
    1000998  // so we can align on 256 byte
    1001   newSize+=32;
     999  newSize += 32;
    10021000  //newSize += (numberDense_+1)/(sizeof(CoinFactorizationDouble)/sizeof(CoinSimplexInt));
    10031001#endif
    10041002#ifdef SMALL_PERMUTE
    10051003  // 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_);
     1004  CoinSimplexInt *COIN_RESTRICT fromSmallToBigRow = reinterpret_cast< CoinSimplexInt * >(workArea_.array());
     1005  CoinSimplexInt *COIN_RESTRICT fromSmallToBigColumn = fromSmallToBigRow + numberRowsSmall_;
     1006  CoinAbcMemcpy(fromSmallToBigRow, fromSmallToBigRow_, numberRowsSmall_);
     1007  CoinAbcMemcpy(fromSmallToBigColumn, fromSmallToBigColumn_, numberRowsSmall_);
    10101008#endif
    10111009  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;
     1010  denseArea_.conditionalNew(newSize);
     1011  denseAreaAddress_ = denseArea_.array();
     1012  CoinInt64 xx = reinterpret_cast< CoinInt64 >(denseAreaAddress_);
     1013  int iBottom = static_cast< int >(xx & 63);
     1014  int offset = (256 - iBottom) >> 3;
    10171015  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_;
     1016  CoinFactorizationDouble *COIN_RESTRICT denseArea = denseAreaAddress_;
     1017  CoinZeroN(denseArea, newSize - 32);
     1018  CoinSimplexInt *COIN_RESTRICT densePermute = reinterpret_cast< CoinSimplexInt * >(denseArea + (leadingDimension_ + FACTOR_CPU) * numberDense_);
     1019#if ABC_DENSE_CODE < 0
     1020  CoinSimplexInt *COIN_RESTRICT indexRowU = indexRowUAddress_;
     1021  CoinSimplexInt *COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
    10251022#endif
    10261023  //mark row lookup using lastRow
    10271024  CoinSimplexInt i;
    1028   CoinSimplexInt *  COIN_RESTRICT lastRow = lastRowAddress_;
    1029   CoinSimplexInt *  COIN_RESTRICT numberInColumn = numberInColumnAddress_;
     1025  CoinSimplexInt *COIN_RESTRICT lastRow = lastRowAddress_;
     1026  CoinSimplexInt *COIN_RESTRICT numberInColumn = numberInColumnAddress_;
    10301027#if 0
    10311028  char xxx[17000];
     
    10601057    abort();
    10611058#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++) {
     1059  for (i = 0; i < numberRowsSmall_; i++) {
     1060    if (lastRow[i] >= 0) {
     1061      lastRow[i] = 0;
     1062    }
     1063  }
     1064#endif
     1065  CoinSimplexInt *COIN_RESTRICT indexRow = indexRowUAddress_;
     1066  CoinFactorizationDouble *COIN_RESTRICT element = elementUAddress_;
     1067  CoinSimplexInt which = 0;
     1068  for (i = 0; i < numberRowsSmall_; i++) {
    10721069    if (!lastRow[i]) {
    1073       lastRow[i]=which;
     1070      lastRow[i] = which;
    10741071#ifdef SMALL_PERMUTE
    1075       int realRow=fromSmallToBigRow[i];
     1072      int realRow = fromSmallToBigRow[i];
    10761073      //if (xxx[realRow]!=0) abort();
    10771074      //xxx[realRow]=-1;
    1078       permuteAddress_[realRow]=numberGoodU_+which;
    1079       densePermute[which]=realRow;
    1080 #else   
    1081       permuteAddress_[i]=numberGoodU_+which;
    1082       densePermute[which]=i;
     1075      permuteAddress_[realRow] = numberGoodU_ + which;
     1076      densePermute[which] = realRow;
     1077#else
     1078      permuteAddress_[i] = numberGoodU_ + which;
     1079      densePermute[which] = i;
    10831080#endif
    10841081      which++;
    10851082    }
    1086   } 
     1083  }
    10871084  //for (i=0;i<numberRowsSmall_;i++) {
    10881085  //int realRow=fromSmallToBigRow[i];
     
    10901087  //}
    10911088  //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_];
     1089  CoinBigIndex *COIN_RESTRICT startColumnL = startColumnLAddress_;
     1090#if ABC_DENSE_CODE < 0
     1091  CoinFactorizationDouble *COIN_RESTRICT elementL = elementLAddress_;
     1092  CoinSimplexInt *COIN_RESTRICT indexRowL = indexRowLAddress_;
     1093#endif
     1094  CoinBigIndex endL = startColumnL[numberGoodL_];
    10981095  //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_;
     1096  CoinFactorizationDouble *COIN_RESTRICT column = denseArea;
     1097  CoinSimplexInt rowsDone = 0;
     1098  CoinSimplexInt iColumn = 0;
     1099  CoinSimplexInt *COIN_RESTRICT pivotColumn = pivotColumnAddress_;
     1100  CoinFactorizationDouble *COIN_RESTRICT pivotRegion = pivotRegionAddress_;
     1101  CoinBigIndex *startColumnU = startColumnUAddress_;
    11051102#ifdef ABC_USE_FUNCTION_POINTERS
    1106   scatterStruct * COIN_RESTRICT scatter = scatterUColumn();
     1103  scatterStruct *COIN_RESTRICT scatter = scatterUColumn();
    11071104#if ABC_USE_FUNCTION_POINTERS
    11081105  extern scatterUpdate AbcScatterLowSubtract[9];
    11091106  extern scatterUpdate AbcScatterHighSubtract[4];
    11101107#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++) {
     1108  CoinFactorizationDouble *COIN_RESTRICT elementUColumnPlus = elementUColumnPlusAddress_;
     1109  CoinSimplexInt *COIN_RESTRICT numberInColumnPlus = numberInColumnPlusAddress_;
     1110#endif
     1111#if ABC_DENSE_CODE == 2
     1112  int nDense = 0;
     1113#endif
     1114  for (iColumn = 0; iColumn < numberRowsSmall_; iColumn++) {
    11181115    if (numberInColumn[iColumn]) {
    1119 #if ABC_DENSE_CODE==2
     1116#if ABC_DENSE_CODE == 2
    11201117      nDense++;
    11211118#endif
     
    11231120      CoinBigIndex start = startColumnU[iColumn];
    11241121      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];
     1122      CoinBigIndex end = start + number;
     1123      for (CoinBigIndex i = start; i < end; i++) {
     1124        CoinSimplexInt iRow = indexRow[i];
     1125        iRow = lastRow[iRow];
     1126        assert(iRow >= 0 && iRow < numberDense_);
     1127#if ABC_DENSE_CODE != 2
     1128        column[iRow] = element[i];
     1129#else
     1130#if BLOCKING8 == 8
     1131        int iBlock = iRow >> 3;
     1132#elif BLOCKING8 == 4
     1133        int iBlock = iRow >> 2;
     1134#elif BLOCKING8 == 2
     1135        int iBlock = iRow >> 1;
     1136#else
     1137        abort();
     1138#endif
     1139        iRow = iRow & (BLOCKING8 - 1);
     1140        column[iRow + BLOCKING8X8 * iBlock] = element[i];
    11441141#endif
    11451142      } /* endfor */
    1146 #if ABC_DENSE_CODE!=2
    1147       column+=leadingDimension_;
    1148 #else
    1149       if ((nDense&(BLOCKING8-1))!=0)
    1150         column += BLOCKING8;
     1143#if ABC_DENSE_CODE != 2
     1144      column += leadingDimension_;
     1145#else
     1146      if ((nDense & (BLOCKING8 - 1)) != 0)
     1147        column += BLOCKING8;
    11511148      else
    1152         column += leadingDimension_*BLOCKING8-(BLOCKING8-1)*BLOCKING8;
    1153 #endif
    1154       while (lastRow[rowsDone]<0) {
     1149        column += leadingDimension_ * BLOCKING8 - (BLOCKING8 - 1) * BLOCKING8;
     1150#endif
     1151      while (lastRow[rowsDone] < 0) {
    11551152        rowsDone++;
    11561153      } /* endwhile */
     
    11581155      int numberIn = numberInColumnPlus[iColumn];
    11591156#ifdef SMALL_PERMUTE
    1160       int realRowsDone=fromSmallToBigRow[rowsDone];
     1157      int realRowsDone = fromSmallToBigRow[rowsDone];
    11611158#if ABC_USE_FUNCTION_POINTERS
    1162       if (numberIn<9) {
    1163         scatter[realRowsDone].functionPointer=AbcScatterLowSubtract[numberIn];
     1159      if (numberIn < 9) {
     1160        scatter[realRowsDone].functionPointer = AbcScatterLowSubtract[numberIn];
    11641161      } else {
    1165         scatter[realRowsDone].functionPointer=AbcScatterHighSubtract[numberIn&3];
    1166       }
    1167 #endif
    1168       scatter[realRowsDone].offset=lastEntryByColumnUPlus_;
    1169       scatter[realRowsDone].number=numberIn;
     1162        scatter[realRowsDone].functionPointer = AbcScatterHighSubtract[numberIn & 3];
     1163      }
     1164#endif
     1165      scatter[realRowsDone].offset = lastEntryByColumnUPlus_;
     1166      scatter[realRowsDone].number = numberIn;
    11701167#else
    11711168#if ABC_USE_FUNCTION_POINTERS
    1172       if (numberIn<9) {
    1173         scatter[rowsDone].functionPointer=AbcScatterLowSubtract[numberIn];
     1169      if (numberIn < 9) {
     1170        scatter[rowsDone].functionPointer = AbcScatterLowSubtract[numberIn];
    11741171      } 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;
     1172        scatter[rowsDone].functionPointer = AbcScatterHighSubtract[numberIn & 3];
     1173      }
     1174#endif
     1175      scatter[rowsDone].offset = lastEntryByColumnUPlus_;
     1176      scatter[rowsDone].number = numberIn;
     1177#endif
     1178      CoinBigIndex startC = start - numberIn;
     1179      for (int i = startC; i < startC + numberIn; i++)
     1180        elementUColumnPlus[lastEntryByColumnUPlus_++] = element[i];
     1181      CoinAbcMemcpy(reinterpret_cast< CoinSimplexInt * >(elementUColumnPlus + lastEntryByColumnUPlus_), indexRow + startC, numberIn);
     1182      lastEntryByColumnUPlus_ += (numberIn + 1) >> 1;
    11861183#endif
    11871184#ifdef SMALL_PERMUTE
    1188       permuteAddress_[realRowsDone]=numberGoodU_;
    1189       pivotColumn[numberGoodU_]=fromSmallToBigColumn[iColumn];
    1190 #else
    1191       permuteAddress_[rowsDone]=numberGoodU_;
    1192       pivotColumn[numberGoodU_]=iColumn;
     1185      permuteAddress_[realRowsDone] = numberGoodU_;
     1186      pivotColumn[numberGoodU_] = fromSmallToBigColumn[iColumn];
     1187#else
     1188      permuteAddress_[rowsDone] = numberGoodU_;
     1189      pivotColumn[numberGoodU_] = iColumn;
    11931190#endif
    11941191      rowsDone++;
    1195       startColumnL[numberGoodU_+1]=endL;
    1196       numberInColumn[iColumn]=0;
    1197       pivotRegion[numberGoodU_]=1.0;
     1192      startColumnL[numberGoodU_ + 1] = endL;
     1193      numberInColumn[iColumn] = 0;
     1194      pivotRegion[numberGoodU_] = 1.0;
    11981195      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_)
     1196    }
     1197  }
     1198#if ABC_DENSE_CODE > 0
     1199  //printf("Going dense at %d\n",numberDense_);
     1200  if (denseThreshold_ > 0) {
     1201    if (numberGoodU_ != numberRows_)
    12051202      return -1; // something odd
    1206     numberGoodL_=numberRows_;
     1203    numberGoodL_ = numberRows_;
    12071204    //now factorize
    12081205    //dgef(denseArea,&numberDense_,&numberDense_,densePermute);
    1209 #if ABC_DENSE_CODE!=2
     1206#if ABC_DENSE_CODE != 2
    12101207#ifndef CLAPACK
    12111208    CoinSimplexInt info;
    1212     F77_FUNC(dgetrf,DGETRF)(&numberDense_,&numberDense_,denseArea,&leadingDimension_,densePermute,
    1213                             &info);
    1214    
     1209    F77_FUNC(dgetrf, DGETRF)
     1210    (&numberDense_, &numberDense_, denseArea, &leadingDimension_, densePermute,
     1211      &info);
     1212
    12151213    // need to check size of pivots
    1216     if(info)
     1214    if (info)
    12171215      status = -1;
    12181216#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 );
     1217    status = clapack_dgetrf(CblasColMajor, numberDense_, numberDense_,
     1218      denseArea, leadingDimension_, densePermute);
     1219    assert(!status);
     1220#endif
     1221#else
     1222    status = CoinAbcDgetrf(numberDense_, numberDense_, denseArea, numberDense_, densePermute
     1223#if ABC_PARALLEL == 2
     1224      ,
     1225      parallelMode_
     1226#endif
     1227    );
    12291228    if (status) {
    1230       status=-1;
     1229      status = -1;
    12311230      printf("singular in dense\n");
    12321231    }
    12331232#endif
    12341233    return status;
    1235   } 
    1236 #else
    1237   numberGoodU_ = numberRows_-numberDense_;
     1234  }
     1235#else
     1236  numberGoodU_ = numberRows_ - numberDense_;
    12381237  CoinSimplexInt base = numberGoodU_;
    12391238  CoinSimplexInt iDense;
    1240   CoinSimplexInt numberToDo=numberDense_;
    1241   denseThreshold_=-abs(denseThreshold_); //0;
    1242   CoinSimplexInt *  COIN_RESTRICT nextColumn = nextColumnAddress_;
    1243   const CoinSimplexInt *  COIN_RESTRICT pivotColumnConst = pivotColumnAddress_;
     1239  CoinSimplexInt numberToDo = numberDense_;
     1240  denseThreshold_ = -abs(denseThreshold_); //0;
     1241  CoinSimplexInt *COIN_RESTRICT nextColumn = nextColumnAddress_;
     1242  const CoinSimplexInt *COIN_RESTRICT pivotColumnConst = pivotColumnAddress_;
    12441243  // make sure we have enough space in L and U
    1245   for (iDense=0;iDense<numberToDo;iDense++) {
     1244  for (iDense = 0; iDense < numberToDo; iDense++) {
    12461245    //how much space have we got
    1247     iColumn = pivotColumnConst[base+iDense];
     1246    iColumn = pivotColumnConst[base + iDense];
    12481247    CoinSimplexInt next = nextColumn[iColumn];
    12491248    CoinSimplexInt numberInPivotColumn = iDense;
    1250     CoinBigIndex space = startColumnU[next] 
     1249    CoinBigIndex space = startColumnU[next]
    12511250      - startColumnU[iColumn]
    12521251      - numberInColumnPlus[next];
    12531252    //assume no zero elements
    1254     if ( numberInPivotColumn > space ) {
     1253    if (numberInPivotColumn > space) {
    12551254      //getColumnSpace also moves fixed part
    1256       if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
    1257         return -99;
     1255      if (!getColumnSpace(iColumn, numberInPivotColumn)) {
     1256        return -99;
    12581257      }
    12591258    }
    12601259    // set so further moves will work
    1261     numberInColumn[iColumn]=numberInPivotColumn;
     1260    numberInColumn[iColumn] = numberInPivotColumn;
    12621261  }
    12631262  // 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_ ) {
     1263  for (iColumn = numberGoodU_ + numberToDo; iColumn < numberRows_; iColumn++) {
     1264    nextRow[iColumn] = iColumn;
     1265    startColumnL[iColumn + 1] = endL;
     1266    pivotRegion[iColumn] = 1.0;
     1267  }
     1268  if (lengthL_ + full * 0.5 > lengthAreaL_) {
    12701269    //need more memory
    1271     if ((messageLevel_&4)!=0)
     1270    if ((messageLevel_ & 4) != 0)
    12721271      std::cout << "more memory needed in middle of invert" << std::endl;
    12731272    return -99;
    12741273  }
    1275   CoinFactorizationDouble * COIN_RESTRICT elementU = elementUAddress_;
    1276   CoinSimplexInt *  COIN_RESTRICT ipiv = densePermute-numberDense_;
    1277   int returnCode=CoinAbcDgetrf(numberDense_,numberDense_,denseArea,numberDense_,ipiv);
     1274  CoinFactorizationDouble *COIN_RESTRICT elementU = elementUAddress_;
     1275  CoinSimplexInt *COIN_RESTRICT ipiv = densePermute - numberDense_;
     1276  int returnCode = CoinAbcDgetrf(numberDense_, numberDense_, denseArea, numberDense_, ipiv);
    12781277  if (!returnCode) {
    1279     CoinSimplexDouble *  COIN_RESTRICT element = denseArea;
    1280     for (int iDense=0;iDense<numberToDo;iDense++) {
    1281       int pivotRow=ipiv[iDense];
     1278    CoinSimplexDouble *COIN_RESTRICT element = denseArea;
     1279    for (int iDense = 0; iDense < numberToDo; iDense++) {
     1280      int pivotRow = ipiv[iDense];
    12821281      // get original row
    12831282      CoinSimplexInt originalRow = densePermute[pivotRow];
    12841283      // swap
    1285       densePermute[pivotRow]=densePermute[iDense];
     1284      densePermute[pivotRow] = densePermute[iDense];
    12861285      densePermute[iDense] = originalRow;
    12871286    }
    1288     for (int iDense=0;iDense<numberToDo;iDense++) {
    1289       int iColumn = pivotColumnConst[base+iDense];
     1287    for (int iDense = 0; iDense < numberToDo; iDense++) {
     1288      int iColumn = pivotColumnConst[base + iDense];
    12901289      // get original row
    12911290      CoinSimplexInt originalRow = densePermute[iDense];
    12921291      // do nextRow
    12931292      lastRow[originalRow] = -2; //mark
    1294       permuteAddress_[originalRow]=numberGoodU_;
     1293      permuteAddress_[originalRow] = numberGoodU_;
    12951294      CoinFactorizationDouble pivotMultiplier = element[iDense];
    12961295      //printf("pivotMultiplier %g\n",pivotMultiplier);
     
    12981297      // Do L
    12991298      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         }
     1299      startColumnL[numberGoodL_] = l; //for luck and first time
     1300      for (int iRow = iDense + 1; iRow < numberDense_; iRow++) {
     1301        CoinFactorizationDouble value = element[iRow];
     1302        if (!TEST_LESS_THAN_TOLERANCE(value)) {
     1303          //printf("L %d row %d value %g\n",l,densePermute[iRow],value);
     1304          indexRowL[l] = densePermute[iRow];
     1305          elementL[l++] = value;
     1306        }
    13081307      }
    13091308      numberGoodL_++;
     
    13121311      // update U column
    13131312      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         }
     1313      for (int iRow = 0; iRow < iDense; iRow++) {
     1314        if (!TEST_LESS_THAN_TOLERANCE(element[iRow])) {
     1315          //printf("U %d row %d value %g\n",start,densePermute[iRow],element[iRow]);
     1316          indexRowU[start] = densePermute[iRow];
     1317          elementU[start++] = element[iRow];
     1318        }
    13201319      }
    13211320      element += numberDense_;
    1322       numberInColumn[iColumn]=0;
    1323       numberInColumnPlus[iColumn] += start-startColumnU[iColumn];
    1324       startColumnU[iColumn]=start;
     1321      numberInColumn[iColumn] = 0;
     1322      numberInColumnPlus[iColumn] += start - startColumnU[iColumn];
     1323      startColumnU[iColumn] = start;
    13251324      numberGoodU_++;
    13261325    }
     
    13281327    return -1;
    13291328  }
    1330   numberDense_=0;
     1329  numberDense_ = 0;
    13311330#endif
    13321331  return status;
     
    13891388#endif
    13901389#endif
     1390
     1391/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     1392*/
Note: See TracChangeset for help on using the changeset viewer.